strange ECL RunTimeError when integrating 
Dan Drake 
12/12/11 7:11 PM 
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 64bit 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 

Re: strange ECL RunTimeError when integrating 
kcrisman 
12/13/11 6:48 AM 
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 64bit 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

Re: strange ECL RunTimeError when integrating 
achrzesz 
12/13/11 9:23 AM 
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

Re: strange ECL RunTimeError when integrating 
kcrisman 
12/13/11 9:42 AM 
> 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

Re: strange ECL RunTimeError when integrating 
Nils Bruin 
12/13/11 10:28 AM 
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: 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.

Re: strange ECL RunTimeError when integrating 
kcrisman 
12/13/11 11:07 AM 
> 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/sagemain/file/9e29a3d84c48/sage/interfaces/maxima_lib.py#l142 and the keepfloat:true. What an unhelpful error message, then!

Re: strange ECL RunTimeError when integrating 
Nils Bruin 
12/13/11 11:27 AM 
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/sagemain/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

Re: strange ECL RunTimeError when integrating 
achrzesz 
12/13/11 11:51 AM 
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

Re: [sagedevel] Re: strange ECL RunTimeError when integrating 
Dan Drake 
12/13/11 4:09 PM 
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. ...especially when the float is 1.5, which can be represented exactly as a binary float! Dan   Dan Drake  http://mathsci.kaist.ac.kr/~drake 

Re: [sagedevel] Re: strange ECL RunTimeError when integrating 
Dan Drake 
12/13/11 4:58 PM 

Re: strange ECL RunTimeError when integrating 
Nils Bruin 
12/13/11 5:42 PM 
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! Actually, I think it is by design that "integrate" does not work right when "keepfloat: true;". The routine "crecip" (reciprocal modulo "modulus") relies one dynamicallyscoped 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/sbbsdsockets.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 
Re: strange ECL RunTimeError when integrating 
kcrisman 
12/13/11 6:38 PM 
> 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... > (%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.

Re: [sagedevel] Re: strange ECL RunTimeError when integrating 
Dan Drake 
12/13/11 7:36 PM 
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... 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. Dan
  Dan Drake  http://mathsci.kaist.ac.kr/~drake 

Re: strange ECL RunTimeError when integrating 
achrzesz 
12/14/11 1:20 AM 
> 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 nonsimplified) answer) Andrzej Chrzeszczyk

Re: strange ECL RunTimeError when integrating 
Nils Bruin 
12/14/11 8:16 AM 
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... 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 a1 zero or nonzero] Of course one could get subtle and wonder what integrate( 1/(x1.0) * 1/(x1.0), x) should be (i.e., do we treat this as 1/(x1.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)
