strange ECL RunTimeError when integrating
Dan Drake
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
-------
kcrisman
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
achrzesz
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
kcrisman
> 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
Nils Bruin
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.
kcrisman
> 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!
Nils Bruin
> 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
achrzesz
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
Dan Drake
> 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!
Dan Drake
> Can you open a ticket and cc: me?
This is now http://trac.sagemath.org/sage_trac/ticket/12152
Nils Bruin
> ...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
kcrisman
> 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...
Dan Drake
> > 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.
achrzesz
> 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
Nils Bruin
> 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)