Bug:Different results when using preparse

[M M]

I get different results from Sage when I try to get a numerical approximation for an expression and if I use evaluate a preparse of the string. I get different results on different versions of sage as well. Here are samples:

----------------------------------------------------------------------

| Sage Version 5.3, Release Date: 2012-09-08                         |

| Type "notebook()" for the browser-based notebook interface.        |

| Type "help()" for help.                                            |

----------------------------------------------------------------------

sage: eval(preparse("numerical_approx(integral(x/(x^3-x+1), x, 1, 2))"))

0.132008722884722

sage: numerical_approx(integral(x/(x^3-x+1), x, 1, 2))

0.132008722884722

sage: A = integral(x/(x^3-x+1), x, 1, 2)

sage: A_str = str(A)

sage: eval(preparse("numerical_approx("+A_str +")"))

-0.393296585552547

sage: 

┌────────────────────────────────────────────────────────────────────┐
│ Sage Version 6.5, Release Date: 2015-02-17                         │
│ Type "notebook()" for the browser-based notebook interface.        │
│ Type "help()" for help.                                            │
└────────────────────────────────────────────────────────────────────┘
sage: eval(preparse("numerical_approx(integral(x/(x^3-x+1), x, 1, 2))"))
1.45943044687563
sage: numerical_approx(integral(x/(x^3-x+1), x, 1, 2))
1.45943044687563
sage: A = integral(x/(x^3-x+1), x, 1, 2)
sage:
sage: A_str = str(A)
sage:
sage: eval(preparse("numerical_approx("+A_str +")"))
0.159046901967485
sage:

┌────────────────────────────────────────────────────────────────────┐
│ Sage Version 6.2, Release Date: 2014-05-06                         │
│ Type "notebook()" for the browser-based notebook interface.        │
│ Type "help()" for help.                                            │
└────────────────────────────────────────────────────────────────────┘
sage: eval(preparse("numerical_approx(integral(x/(x^3-x+1), x, 1, 2))"))
1.64714767119638
sage: numerical_approx(integral(x/(x^3-x+1), x, 1, 2))
1.64714767119638
sage: A = integral(x/(x^3-x+1), x, 1, 2)
sage: A_str = str(A)
sage: eval(preparse("numerical_approx("+A_str +")"))
0.162416510011260
sage:

On the Notebook on the cloud it gives me the following error although the version is the same as one of the versions I tried locally

version()

'Sage Version 6.5, Release Date: 2015-02-17'

numerical_approx(integral(x/(x^3-x+1), x, 1, 2))

Error in lines 1-1 Traceback (most recent call last): File "/projects/f700a2f3-7f30-4b47-9f18-e0309eb8c48c/.sagemathcloud/sage_server.py", line 875, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "<string>", line 1, in <module> File "/usr/local/sage/sage-6.5/local/lib/python2.7/site-packages/sage/misc/functional.py", line 1298, in numerical_approx return x._numerical_approx(prec, algorithm=algorithm) File "sage/symbolic/expression.pyx", line 4861, in sage.symbolic.expression.Expression._numerical_approx (build/cythonized/sage/symbolic/expression.cpp:27184) x = self._convert(kwds) File "sage/symbolic/expression.pyx", line 1034, in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:7790) cdef GEx res = self._gobj.evalf(0, kwds) ValueError: power::eval(): division by zero

[Dima Pasechnik]

On 2015-03-10, M M <mirette...@gmail.com> wrote:
> I get different results from Sage when I try to get a numerical
> approximation for an expression and if I use evaluate a preparse of the
> string. I get different results on different versions of sage as well. Here
> are samples:

Integration is done by Maxima, and it is a bloody mess; e.g.
(with Sage 6.6.beta3)

sage: numerical_integral(x/(x^3-x+1.), 1., 2.)
(0.565799916456428, 6.281640945809164e-15)
sage: integral(x/(x^3-x+1), x, 1, 2).n()
1.11396858782562
sage: integral(x/(x^3-x+1), x, 1, 2).full_simplify().n()
0.168522351678247
sage:

Do you know which one is actually correct?
It looks like the 2nd one might be right (by looking at
sage: plot(x/(x^3-x+1), x, 1, 2))


PS. .n() is the same as numerical_approx()

[Dima Pasechnik]

On 2015-03-10, Dima Pasechnik <dim...@gmail.com> wrote:
> On 2015-03-10, M M <mirette...@gmail.com> wrote:
>> I get different results from Sage when I try to get a numerical
>> approximation for an expression and if I use evaluate a preparse of the
>> string. I get different results on different versions of sage as well. Here
>> are samples:
>
> Integration is done by Maxima, and it is a bloody mess; e.g.
> (with Sage 6.6.beta3)
>
> sage: numerical_integral(x/(x^3-x+1.), 1., 2.)
> (0.565799916456428, 6.281640945809164e-15)
> sage: integral(x/(x^3-x+1), x, 1, 2).n()
> 1.11396858782562
> sage: integral(x/(x^3-x+1), x, 1, 2).full_simplify().n()
> 0.168522351678247
> sage:
>
> Do you know which one is actually correct?
> It looks like the 2nd one might be right (by looking at
> sage: plot(x/(x^3-x+1), x, 1, 2))

sorry, the 1st looks right (the plot is not from 0 to 2, but from 1 to
2, so the value looks a bit above 1/2).
(and Wolfram Alpha outputs 0.5658)

[Nils Bruin]

On Tuesday, March 10, 2015 at 9:04:14 AM UTC-7, Dima Pasechnik wrote:

On 2015-03-10, M M <mirette...@gmail.com> wrote:
> I get different results from Sage when I try to get a numerical
> approximation for an expression and if I use evaluate a preparse of the
> string. I get different results on different versions of sage as well. Here
> are samples:

Integration is done by Maxima, and it is a bloody mess; e.g.

It is, but I suspect that's not the cause here. I think it's just numerical instability.

sage: sage: integral(x/(x^3-x+1), x, 1, 2).n(100)
0.56579991645642210974671290281
sage: sage: integral(x/(x^3-x+1), x, 1, 2).simplify_full().n(100)
0.56579991645643344322713389324

i.e., approximating the exact result just needs some more digits to be accurate. The default "n" only specifies the precision *used* using evaluation, so the error in the numerical approximation isn't a priori controlled and depends on the particular evaluation scheme chosen for the expression (expect that to be just "left to right evaluate" as you would do on a calculator). Specifying more digits to work with should give a better result for well-behaved expressions.

In principle, using the RealIntervalField you should get guaranteed digits:

sage: I=integral(x/(x^3-x+1), x, 1, 2)
sage: RealIntervalField(53)(I)
[-infinity .. +infinity]
sage: RealIntervalField(60)(I)
1.?
sage: RealIntervalField(70)(I)
0.566?
sage: RealIntervalField(90)(I)
0.565799917?

Needless to say, probably numerical integration is the better way of getting a numerical approximation.

[ssin...@coe.edu]

On Tuesday, March 10, 2015 at 3:18:45 PM UTC-5, Nils Bruin wrote:

On Tuesday, March 10, 2015 at 9:04:14 AM UTC-7, Dima Pasechnik wrote:

On 2015-03-10, M M <mirette...@gmail.com> wrote:

Integration is done by Maxima, and it is a bloody mess; e.g.

It is, but I suspect that's not the cause here. I think it's just numerical instability.

sage: sage: integral(x/(x^3-x+1), x, 1, 2).n(100)
0.56579991645642210974671290281
sage: sage: integral(x/(x^3-x+1), x, 1, 2).simplify_full().n(100)
0.56579991645643344322713389324

For what it's worth...
┌────────────────────────────────────────────────────────────────────┐

│ Sage Version 6.5, Release Date: 2015-02-17                         │

│ Enhanced for SageMathCloud.                                        │
└────────────────────────────────────────────────────────────────────┘
sage: integral(x/(x^3-x+1), x, 1, 2, algorithm='sympy').n()
0.565799916456428
sage: integral(x/(x^3-x+1), x, 1, 2, algorithm='sympy').simplify_full().n()
0.565799916456428
sage: integral(x/(x^3-x+1), x, 1, 2, algorithm='sympy').n(100)
0.56579991645642798527545946385
sage: integral(x/(x^3-x+1), x, 1, 2, algorithm='sympy').simplify_full().n(100)
0.56579991645642798527545946385

--Steve

[Dima Pasechnik]

On 2015-03-10, Nils Bruin <nbr...@sfu.ca> wrote:
> On Tuesday, March 10, 2015 at 9:04:14 AM UTC-7, Dima Pasechnik wrote:
>>

>> On 2015-03-10, M M <mirette...@gmail.com <javascript:>> wrote:
>> > I get different results from Sage when I try to get a numerical
>> > approximation for an expression and if I use evaluate a preparse of the
>> > string. I get different results on different versions of sage as well.
>> Here
>> > are samples:
>>
>> Integration is done by Maxima, and it is a bloody mess; e.g.
>>
>
> It is, but I suspect that's not the cause here. I think it's just numerical
> instability.
>
> sage: sage: integral(x/(x^3-x+1), x, 1, 2).n(100)
> 0.56579991645642210974671290281
> sage: sage: integral(x/(x^3-x+1), x, 1, 2).simplify_full().n(100)
> 0.56579991645643344322713389324

Indeed, the integers produced by Maxima in the exact answer are rather long.

I tried this integral directly in Maxima, and taking bfloat of it
outputs nonsense.

It seems that Sage's .n(100) works better here.

Interestingly, Axiom/FriCAS is able to compute the corresponding indefinite
integral:

http://axiom-wiki.newsynthesis.org/ExampleIntegration

I wish there was a more accessible full implementation of Risch algorithm...

[Nils Bruin]

On Wednesday, March 11, 2015 at 2:46:25 AM UTC-7, Dima Pasechnik wrote:

I tried this integral directly in Maxima, and taking bfloat of it
outputs nonsense.

I have noticed before that bfloats aren't infectious enough: operations on bfloats can easily result in a normal "double". I think there are ways to convince maxima to use bfloats more pervasively. Perhaps a global precision setting somewhere?
 

I wish there was a more accessible full implementation of Risch algorithm...

This is a rational function, so a first calculus course would already teach you the relevant part of the Risch algorithm. It's a little more tricky to get an ostensibly real-valued function as an antiderivative. Anyway, sympy produces a reasonable-looking antiderivative.
 
Interestingly, we have:

sage: I=integral(x/(x^3-x+1), x, 1, 2, algorithm='sympy')
sage: RIF(I)
TypeError: unable to simplify to a real interval approximation

The offending subexpression seems to be:

sage: A=(299838966359964800*69^(5/6)*2^(2/3) - 11515081166050000*69^(2/3)*2^(1/3)*(25*sqrt(69) + 207)^(1/3) - 99785894223312000*sqrt(69)*(25*sqrt(69) + 207)^(2/3) + 2271318237097115625*69^(1/3)*2^(2/3) - 99785894223312000*69^(1/6)*2^(1/3)*(25*sqrt(69) + 207)^(1/3) - 497728835949744*9522^(1/3)*(25*sqrt(69) + 207)^(1/3) - 828883890137982336*(25*sqrt(69) + 207)^(2/3) + 219331275901257879*276^(1/3))^(QQ(-1))
sage: RIF(A)
TypeError: unable to simplify to a real interval approximation

Note the *rational* exponent -1. If that's an integer there's no problem. Using RealIntervalField(200) has the same problem. Using RealField(...) seems to work fine.

[Dima Pasechnik]

On 2015-03-11, Nils Bruin <nbr...@sfu.ca> wrote:
> On Wednesday, March 11, 2015 at 2:46:25 AM UTC-7, Dima Pasechnik wrote:
>>
>> I tried this integral directly in Maxima, and taking bfloat of it
>> outputs nonsense.
>>
>
> I have noticed before that bfloats aren't infectious enough: operations on
> bfloats can easily result in a normal "double". I think there are ways to
> convince maxima to use bfloats more pervasively. Perhaps a global precision
> setting somewhere?
>
>
>> I wish there was a more accessible full implementation of Risch
>> algorithm...
>>
>
> This is a rational function, so a first calculus course would already teach
> you the relevant part of the Risch algorithm. It's a little more tricky to

Risch, as implemented in Axiom, does not do factorisation (i.e. no
partial fractions).
In this example at least it produces much nicer looking antiderivative,
no huge integers.
http://axiom-wiki.newsynthesis.org/ExampleIntegration

Dima

[William Stein]

On Wed, Mar 11, 2015 at 9:35 AM, Dima Pasechnik <dim...@gmail.com> wrote:
> On 2015-03-11, Nils Bruin <nbr...@sfu.ca> wrote:
>> On Wednesday, March 11, 2015 at 2:46:25 AM UTC-7, Dima Pasechnik wrote:
>>>
>>> I tried this integral directly in Maxima, and taking bfloat of it
>>> outputs nonsense.
>>>
>>
>> I have noticed before that bfloats aren't infectious enough: operations on
>> bfloats can easily result in a normal "double". I think there are ways to
>> convince maxima to use bfloats more pervasively. Perhaps a global precision
>> setting somewhere?
>>
>>
>>> I wish there was a more accessible full implementation of Risch
>>> algorithm...
>>>
>>
>> This is a rational function, so a first calculus course would already teach
>> you the relevant part of the Risch algorithm. It's a little more tricky to
>
> Risch, as implemented in Axiom, does not do factorisation (i.e. no
> partial fractions).
> In this example at least it produces much nicer looking antiderivative,
> no huge integers.
> http://axiom-wiki.newsynthesis.org/ExampleIntegration
>
> Dima

For what it's worth, here's how to mostly do that Axiom session, but
in a SageMathCloud worksheet...

https://cloud.sagemath.com/projects/4a5f0542-5873-4eed-a85c-a18c706e8bcd/files/support/2015-03-11-093745-axiom-integral.sagews

William

>
>> get an ostensibly real-valued function as an antiderivative. Anyway, sympy
>> produces a reasonable-looking antiderivative.
>>
>> Interestingly, we have:
>>
>> sage: I=integral(x/(x^3-x+1), x, 1, 2, algorithm='sympy')
>> sage: RIF(I)
>> TypeError: unable to simplify to a real interval approximation
>>
>> The offending subexpression seems to be:
>>
>> sage: A=(299838966359964800*69^(5/6)*2^(2/3) -
>> 11515081166050000*69^(2/3)*2^(1/3)*(25*sqrt(69) + 207)^(1/3) -
>> 99785894223312000*sqrt(69)*(25*sqrt(69) + 207)^(2/3) +
>> 2271318237097115625*69^(1/3)*2^(2/3) -
>> 99785894223312000*69^(1/6)*2^(1/3)*(25*sqrt(69) + 207)^(1/3) -
>> 497728835949744*9522^(1/3)*(25*sqrt(69) + 207)^(1/3) -
>> 828883890137982336*(25*sqrt(69) + 207)^(2/3) +
>> 219331275901257879*276^(1/3))^(QQ(-1))
>> sage: RIF(A)
>> TypeError: unable to simplify to a real interval approximation
>>
>> Note the *rational* exponent -1. If that's an integer there's no problem.
>> Using RealIntervalField(200) has the same problem. Using RealField(...)
>> seems to work fine.
>>
>

> --
> You received this message because you are subscribed to the Google Groups "sage-support" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
> To post to this group, send email to sage-s...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-support.
> For more options, visit https://groups.google.com/d/optout.



--
William (http://wstein.org)

[Dima Pasechnik]

On Wednesday, 11 March 2015 16:44:32 UTC, William wrote:

On Wed, Mar 11, 2015 at 9:35 AM, Dima Pasechnik <dim...@gmail.com> wrote:
> On 2015-03-11, Nils Bruin <nbr...@sfu.ca> wrote:
>> On Wednesday, March 11, 2015 at 2:46:25 AM UTC-7, Dima Pasechnik wrote:
>>>
>>> I tried this integral directly in Maxima, and taking bfloat of it
>>> outputs nonsense.
>>>
>>
>> I have noticed before that bfloats aren't infectious enough: operations on
>> bfloats can easily result in a normal "double". I think there are ways to
>> convince maxima to use bfloats more pervasively. Perhaps a global precision
>> setting somewhere?
>>
>>
>>> I wish there was a more accessible full implementation of Risch
>>> algorithm...
>>>
>>
>> This is a rational function, so a first calculus course would already teach
>> you the relevant part of the Risch algorithm. It's a little more tricky to
>
> Risch, as implemented in Axiom, does not do factorisation (i.e. no
> partial fractions).
> In this example at least it produces much nicer looking antiderivative,
> no huge integers.
> http://axiom-wiki.newsynthesis.org/ExampleIntegration
>
> Dima

For what it's worth, here's how to mostly do that Axiom session, but
in a SageMathCloud worksheet...

   https://cloud.sagemath.com/projects/4a5f0542-5873-4eed-a85c-a18c706e8bcd/files/support/2015-03-11-093745-axiom-integral.sagews

here is a better version (all the stuff works):
https://cloud.sagemath.com/projects/bb6fd6ca-6304-4dda-be31-bd2dd5eb3d98/files/support/2015-03-11-093745-axiom-integral.sagews

Dima
 

[William Stein]

Thanks! I replaced mine by yours.

William

[M M]

Thanks so much for all the efforts for making sage output consistent results for the numerical approximation. However, the main problem I had was the fact that Sage returns different answers when I preparse the string of the same sage expression as in the examples of the original post.

[Nils Bruin]

On Thursday, March 12, 2015 at 1:00:28 PM UTC-7, M M wrote:

Thanks so much for all the efforts for making sage output consistent results for the numerical approximation. However, the main problem I had was the fact that Sage returns different answers when I preparse the string of the same sage expression as in the examples of the original post.

I was not able to replicate that behaviour on a recent build.

[M M]

I tried on a different machine with Sage Version 6.5, Release Date: 2015-02-17   and I get the same behaviour described in the original post. Where parsing a string and evaluating returns different results than using the resulted sage expression.

[Nils Bruin]

On Friday, March 13, 2015 at 11:32:11 AM UTC-7, M M wrote:

I tried on a different machine with Sage Version 6.5, Release Date: 2015-02-17   and I get the same behaviour described in the original post. Where parsing a string and evaluating returns different results than using the resulted sage expression.

Ah, right. Yes, evaluating a string representation of an object does necessarily result in the same object. We have already established that the expression in question is ill-conditioned for numerical purposes, so you probably end up with a subtly different internal representation of the object, which has huge numerical consequences. Indeed:

sage: A = integral(x/(x^3-x+1), x, 1, 2)

sage: B = eval(preparse(str(A)))
sage: (A-B).n()
0.894231297231034
sage: (A-B).simplify()
0

This is entirely consistent with float operations not being commutative/associative/etc. It's unfortunate, but if this is a bug then floating point is a bug.

[Nils Bruin]

On Friday, March 13, 2015 at 1:33:57 PM UTC-7, Nils Bruin wrote:

Ah, right. Yes, evaluating a string representation of an object does NOT necessarily result in the same object.

Apologies for the typo.

[Vincent Delecroix]

>> Ah, right. Yes, evaluating a string representation of an object does

>> *NOT* necessarily result in the same object.

Note that we have sage_input for that

sage: m1 = matrix(ZZ, 2, range(4))
sage: m2 = matrix(GF(7), 2, range(4))

sage: str(m1) == str(m2)
True

sage: sage_input(m1)
matrix(ZZ, [[0, 1], [2, 3]])
sage: sage_input(m2)
matrix(GF(7), [[0, 1], [2, 3]])

Vincent