70 views

Skip to first unread message

Dec 12, 2011, 10:11:01 PM12/12/11

to sage-...@googlegroups.com

I'm doing some integrals:

sage: a, b, t = var('a b t')

sage: f(a,b,t) = sin(t)^2/(a + b*cos(t))^2

sage: integrate(f(3/2,1,t), (t,0,2*pi))

-2/5*(sqrt(5) - 3)*pi*sqrt(5)

Okay, that's fine. But

sage: integrate(f(1.5,1,t), (t,0,2*pi))

blows up with:

RuntimeError: ECL says: Error executing code in Maxima: CRECIP:

attempted inverse of zero (mod 3)

This is with 4.7.2 on 64-bit Ubuntu. I get the same error on sage.math.

How on earth is it trying to invert something mod 3?

Dan

--

--- Dan Drake

----- http://mathsci.kaist.ac.kr/~drake

-------

Dec 13, 2011, 9:48:57 AM12/13/11

to sage-devel

On Dec 12, 10:11 pm, Dan Drake <dr...@kaist.edu> wrote:

> I'm doing some integrals:

>

> sage: a, b, t = var('a b t')

> sage: f(a,b,t) = sin(t)^2/(a + b*cos(t))^2

> sage: integrate(f(3/2,1,t), (t,0,2*pi))

> -2/5*(sqrt(5) - 3)*pi*sqrt(5)

>

> Okay, that's fine. But

>

> sage: integrate(f(1.5,1,t), (t,0,2*pi))

>

> blows up with:

>

> RuntimeError: ECL says: Error executing code in Maxima: CRECIP:

> attempted inverse of zero (mod 3)

>

> This is with 4.7.2 on 64-bit Ubuntu. I get the same error on sage.math.

>

> How on earth is it trying to invert something mod 3?

Hmm, not obvious to me. I can get it to happen on Mac, but I can't

get it to happen in the Maxima console.

(%i1) display2d:false;

(%o1) false

(%i2) h:sin(t)^2/(cos(t)+1.5)^2;

(%o2) sin(t)^2/(cos(t)+1.5)^2

(%i6) integrate(h,t,0,2*%pi);

rat: replaced 1.5 by 3/2 = 1.5

rat: replaced 1.5 by 3/2 = 1.5

rat: replaced 1.5 by 3/2 = 1.5

rat: replaced 1.5 by 3/2 = 1.5

rat: replaced 1.5 by 3/2 = 1.5

(%o6) -2*%pi

Can you open a ticket and cc: me?

- kcrisman

Dec 13, 2011, 12:23:06 PM12/13/11

to sage-devel

Maxima 5.25.1

(%i2) display2d:false;

(%o2) false

(%i3) h:sin(t)^2/(cos(t)+1.5)^2;

(%o3) sin(t)^2/(cos(t)+1.5)^2

(%i4) integrate(h,t,0,2*%pi);

rat: replaced 1.5 by 3/2 = 1.5

rat: replaced 1.5 by 3/2 = 1.5

rat: replaced 1.5 by 3/2 = 1.5

rat: replaced 1.5 by 3/2 = 1.5

rat: replaced 1.5 by 3/2 = 1.5

(%o4) -(2*sqrt(5)-6)*%pi/sqrt(5)

(Dan's result for a=3/2)

Andrzej Chrzeszczyk

Dec 13, 2011, 12:42:04 PM12/13/11

to sage-devel

> Maxima 5.25.1

>

> (%i2) display2d:false;

> (%o2) false

> (%i3) h:sin(t)^2/(cos(t)+1.5)^2;

> (%o3) sin(t)^2/(cos(t)+1.5)^2

> (%i4) integrate(h,t,0,2*%pi);

>

> rat: replaced 1.5 by 3/2 = 1.5

>

> rat: replaced 1.5 by 3/2 = 1.5

>

> rat: replaced 1.5 by 3/2 = 1.5

>

> rat: replaced 1.5 by 3/2 = 1.5

>

> rat: replaced 1.5 by 3/2 = 1.5

> (%o4) -(2*sqrt(5)-6)*%pi/sqrt(5)

>

> (Dan's result for a=3/2)

Right, but I'm wondering why I can't reproduce the " Error executing

code in Maxima: CRECIP:

attempted inverse of zero (mod 3) " bit in the Maxima currently in

Sage? The problem is that then upgrading Maxima might not fix this.

- kcrisman

Dec 13, 2011, 1:28:54 PM12/13/11

to sage-devel

The bug can be triggered via the pexpect interface to maxima, so it

probably has something to do with the things that sage preloads in

maxima:

probably has something to do with the things that sage preloads in

maxima:

sage: maxima("integrate(sin(t)^2/(1*cos(t) + 1.5)^2,t,0,2*%pi)")

TypeError: Error executing code in Maxima

CODE:

sage1 : integrate(sin(t)^2/(1*cos(t) + 1.5)^2,t,0,2*%pi)$

Maxima ERROR:

CRECIP: attempted inverse of zero (mod 3)

-- an error. To debug this try: debugmode(true);

versus:

sage: maxima("integrate(sin(t)^2/(1*cos(t) + 3/2)^2,t,0,2*%pi)")

-2*%pi

This shows that the problem does not lie in differences between

library mode maxima and maxima proper (this maxima interface runs the

same executable as the maxima console)

The problem arises from the rational approximation that maxima tries

to do resulting from 1.5. Given that the maxima console seems to do

this "right", my guess would be that we turn off this rational

approximation trick and that it gets triggered later, resulting in

some "mod 3" computations. CRECIP has to do with finding rational

approximations to floats.

Dec 13, 2011, 2:07:37 PM12/13/11

to sage-devel

> The problem arises from the rational approximation that maxima tries

> to do resulting from 1.5. Given that the maxima console seems to do

> this "right", my guess would be that we turn off this rational

> approximation trick and that it gets triggered later, resulting in

Yes, we try to turn it off as much as possible - see

http://hg.sagemath.org/sage-main/file/9e29a3d84c48/sage/interfaces/maxima_lib.py#l142

and the keepfloat:true. What an unhelpful error message, then!

Dec 13, 2011, 2:27:49 PM12/13/11

to sage-devel

On Dec 13, 11:07 am, kcrisman <kcris...@gmail.com> wrote:

> Yes, we try to turn it off as much as possible - seehttp://hg.sagemath.org/sage-main/file/9e29a3d84c48/sage/interfaces/ma...

> and the keepfloat:true. What an unhelpful error message, then!

> Yes, we try to turn it off as much as possible - seehttp://hg.sagemath.org/sage-main/file/9e29a3d84c48/sage/interfaces/ma...

> and the keepfloat:true. What an unhelpful error message, then!

I think the error is triggered by an actual bug in "rationalize"

somewhere, so it might be worth reporting the error upstream. There is

no reason why "replace a float by a nearby rational number" should

ever fail.

Workaround:

sage: integrate(f(1.5,1,t), (t,0,2*pi))

sage: M=sage.calculus.calculus.maxima

sage: M.eval("keepfloat: false;")

'false'

sage: integrate(f(1.5,1,t), (t,0,2*pi))

-2*pi

sage: M.eval("keepfloat: true;") # this is the default

'true'

sage: integrate(f(1.5,1,t), (t,0,2*pi))

RuntimeError

Dec 13, 2011, 2:51:10 PM12/13/11

to sage-devel

The workaround should give a correct answer

with newer maxima. (As I have noticed) with the present maxima

version

the result -2*pi is wrong

Dec 13, 2011, 7:09:55 PM12/13/11

to sage-...@googlegroups.com

On Tue, 13 Dec 2011 at 11:27AM -0800, Nils Bruin wrote:

> I think the error is triggered by an actual bug in "rationalize"

> somewhere, so it might be worth reporting the error upstream. There is

> no reason why "replace a float by a nearby rational number" should

> ever fail.

> I think the error is triggered by an actual bug in "rationalize"

> somewhere, so it might be worth reporting the error upstream. There is

> no reason why "replace a float by a nearby rational number" should

> ever fail.

...especially when the float is 1.5, which can be represented exactly as

a binary float!

Dec 13, 2011, 7:58:10 PM12/13/11

to sage-...@googlegroups.com

On Tue, 13 Dec 2011 at 06:48AM -0800, kcrisman wrote:

> Can you open a ticket and cc: me?

> Can you open a ticket and cc: me?

This is now http://trac.sagemath.org/sage_trac/ticket/12152

Dec 13, 2011, 8:42:44 PM12/13/11

to sage-devel

On Dec 13, 4:09 pm, Dan Drake <dr...@kaist.edu> wrote:

> ...especially when the float is 1.5, which can be represented exactly as

> a binary float!

> ...especially when the float is 1.5, which can be represented exactly as

> a binary float!

Actually, I think it is by design that "integrate" does not work right

when "keepfloat: true;".

The routine "crecip" (reciprocal modulo "modulus") relies one

dynamically-scoped variable modulus. It looks like polynomial

arithmetic (which is conceivably necessary when computing integrals)

can be thrown off when "modulus" is set in some outer scope. Evidence

(transcript below). Verbal description

- we run the example with debugger "on" and "keepfloat: true". This

brings us in the debugger shell, with the scope of the error. When we

try to compute the gcd(x^2+1.5,x), we get an error, clearly showing

that the "modulus" is the problem.

- when we run gcd(x^2+1.5,x) at toplevel, no such error occurs;

presumably because modulus does not have a value in that scope.

Conjecture: The integration error occurs due to a similar problem, but

in a less shielded way, so we get a less intelligible error message

about essentially the same problem.

From a maxima perspective, I don't think anymore this really is an

error. I don't think you can expect the advanced integration routines

to work in the presence of floats, because polynomial arithmetic is

necessary and things like gcds are horrible to compute for polynomials

with float coefficients. Therefore, it is written with the assumption

that keepfloat:false holds, or at least that integration does not

encounter floats. With that assumption, they are free to set

"modulus", because the codepaths that are affected by it should not be

triggered.

Proposed bugfix: symbolic integration should reject floats.

sage: maxima.console()

;;; Loading #P"/usr/local/sage/local/lib/ecl/sb-bsd-sockets.fas"

;;; Loading #P"/usr/local/sage/local/lib/ecl/sockets.fas"

;;; Loading #P"/usr/local/sage/local/lib/ecl/defsystem.fas"

;;; Loading #P"/usr/local/sage/local/lib/ecl/cmp.fas"

Maxima 5.23.2 http://maxima.sourceforge.net

using Lisp ECL 11.1.1

Distributed under the GNU Public License. See the file COPYING.

Dedicated to the memory of William Schelter.

The function bug_report() provides bug reporting information.

(%i1) keepfloat: true;

(%o1) true

(%i2) debugmode(true);

(%o2) true

(%i3) integrate(sin(t)^2/(1*cos(t) + 1.5)^2,t,0,2*%pi);

CRECIP: attempted inverse of zero (mod 3)

-- an error. Entering the Maxima debugger.

Enter ':h' for help.

(dbm:1) gcd(x^2+1.5,x);

rat: can't rationalize 1.5 when modulus = 3

[...]

(%i4) gcd(x^2+1.5,x);

rat: replaced 1.5 by 3/2 = 1.5

1

(%o4) -

2

Dec 13, 2011, 9:38:48 PM12/13/11

to sage-devel

> Conjecture: The integration error occurs due to a similar problem, but

> in a less shielded way, so we get a less intelligible error message

> about essentially the same problem.

>

> From a maxima perspective, I don't think anymore this really is an

> error. I don't think you can expect the advanced integration routines

> to work in the presence of floats, because polynomial arithmetic is

> necessary and things like gcds are horrible to compute for polynomials

> with float coefficients. Therefore, it is written with the assumption

> that keepfloat:false holds, or at least that integration does not

> encounter floats. With that assumption, they are free to set

> "modulus", because the codepaths that are affected by it should not be

> triggered.

Nice analysis.

> Proposed bugfix: symbolic integration should reject floats.

Yikes! So we wouldn't want

sage: integrate(1.0*x,x,0,2)

2

to work any more? (Pace the error messages below.)

Could we not instead try to catch this kind of error and recommend to

the user that they use rationals? Or something...

Dec 13, 2011, 10:36:32 PM12/13/11

to sage-...@googlegroups.com

On Tue, 13 Dec 2011 at 06:38PM -0800, kcrisman wrote:

> > Proposed bugfix: symbolic integration should reject floats.

>

> Yikes! So we wouldn't want

>

> sage: integrate(1.0*x,x,0,2)

> 2

>

> to work any more? (Pace the error messages below.)

>

> Could we not instead try to catch this kind of error and recommend to

> the user that they use rationals? Or something...

> > Proposed bugfix: symbolic integration should reject floats.

>

> Yikes! So we wouldn't want

>

> sage: integrate(1.0*x,x,0,2)

> 2

>

> to work any more? (Pace the error messages below.)

>

> Could we not instead try to catch this kind of error and recommend to

> the user that they use rationals? Or something...

As long as there are informative error messages, I don't think I'd mind

too much. In that case, I would also want some kind of rational

approximation method for symbolic expressions: something like

(1.0*x^2).float_coeffs_to_rational()

becoming

1*x^2.

This could just run exact_rational() or simplest_rational() on the float

coefficients, perhaps with a keyword to select what you want. I would consider

defaulting to simplest_rational() because of things like this:

sage: (1/3).n().exact_rational()

6004799503160661/18014398509481984

I should mention that this entire thing got started when I took my f()

function and tried to use list_plot to see what was happening as you

vary the coefficients: I did something like

[integrate(f(a, 1, t), (t, 0, 2*pi)) for a in srange(1, 2, .1)]

and it failed right away.

Dec 14, 2011, 4:20:27 AM12/14/11

to sage-devel

> signature.asc

> < 1KWyświetlPobierz

Dan's problems can be avoided:

sage: a, b, t = var('a b t')

sage: f(a,b,t) = sin(t)^2/(a + b*cos(t))^2

sage: integrate(f(QQ(1.5),1,t), (t,0,2*pi))

-2/5*(sqrt(5) - 3)*pi*sqrt(5)

#for list_plot:

sage: [numerical_integral(f(a, 1, t), 0, 2*pi)[0] for a in srange(1.1,

2, .1)]

[8.7989525514600722, 5.0835245940142055, 3.5501007957306188,

2.6946635039534597, 2.1465923700692833, 1.7657335988517402,

1.4864008624226634, 1.2734482804526248, 1.1062834223104587]

but how to avoid the problem (with an exact integral)

sage: integrate(f(sqrt(2),1,t), (t,0,2*pi))

RuntimeError: ECL says: Error executing code in Maxima: Division by 0

(hint: maxima 5.25.1 gives the correct (but non-simplified) answer)

Andrzej Chrzeszczyk

Dec 14, 2011, 11:16:00 AM12/14/11

to sage-devel

On Dec 13, 6:38 pm, kcrisman <kcris...@gmail.com> wrote:

> Yikes! So we wouldn't want

>

> sage: integrate(1.0*x,x,0,2)

> 2

>

> to work any more? (Pace the error messages below.)

>

> Could we not instead try to catch this kind of error and recommend to

> the user that they use rationals? Or something...

> Yikes! So we wouldn't want

>

> sage: integrate(1.0*x,x,0,2)

> 2

>

> to work any more? (Pace the error messages below.)

>

> Could we not instead try to catch this kind of error and recommend to

> the user that they use rationals? Or something...

well, we *could* replace floats with dummy variables, do the

integration and evaluate back again. That would at least allow you to

do the integrations where the value of the float is irrelevant (as

above) and would fail in a reasonable way on things like

integrate(x^-1.0,x) [NOW: division by 0 error]

integrate(x^-a,x) [error "Is a-1 zero or nonzero]

Of course one could get subtle and wonder what

integrate( 1/(x-1.0) * 1/(x-1.0), x)

should be (i.e., do we treat this as 1/(x-1.0)^2 [this already happens

for this particular example], or do we handle 1.0 and 1.0 as

independent?)

Symbolic integration and floats really don't mix. It's just a matter

of deciding how far you want to go before you give up.

Maxima's approach of just taking the rational value of the digits of

the float (or something close to it) is perhaps not that bad if you

absolutely want to return an answer. However, things like

sage: integrate(1/(x^2+1.0000001),x)

3/10000*sqrt(11111110)*arctan(3/10000*sqrt(11111110)*x)

suggest to me we should have given up. Incidentally,

sage: integrate(x^-1.00000000000000000000001,x)

RuntimeError: ECL says: Error executing code in Maxima: Division by 0

that's because currently any float gets translated to a "double" in

maxima. Maxima has the bfloat type for arbitrary precision floats.

Especially with the "binary" interface that gets used for integration,

it would be straightforward to translate RealNumber of a precision

that does not fit into a double, into a bfloat instead. I don't know

how well supported bfloats are in maxima, though (plus, we should

probably not be using maxima for anything involving floats anyway)

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu