division by zero from integrate to giac. Where does long giac output come from?

asked 2021-05-30 04:14:38 +0100

Nasser gravatar image

updated 2021-05-30 04:45:35 +0100

I found number of integrals where sagemath says ValueError: power::eval(): division by zero: when using giac algorithm.

But when I try the same integrate command in giac directly, it actually gives no output. Not able to integrate it. giac just says Done.

So I wonder why then sagemath gives division by zero? It also gives very large output below that, which I do not know where it came from, since giac does not show any output.

At first I thought giac may be does not display on the screen an output which is too large? I could not find googling around anything about it. But sagemath claims that giac has returned this vey very long result from integrate. How is that possible?

Using giac 1.7 and sagemath 9.3 on Linux

Here is the command inside giac

10>> integrate((e*x+d)^(9/2)/(a*d*e+(a*e^2+c*d^2)*x+c*d*e*x^2),x)
Warning, need to choose a branch for the root of a polynomial with parameters. This might be wrong.
The choice was done assuming [a,c,d,exp(1),exp(2)]=[60,-94,37,-33,49]
Warning, need to choose a branch for the root of a polynomial with parameters. This might be wrong.
The choice was done assuming [a,c,d,exp(1),exp(2)]=[-3,-98,-46,47,-64]
Warning, need to choose a branch for the root of a polynomial with parameters. This might be wrong.
The choice was done assuming [a,c,d,exp(1),exp(2)]=[-22,22,-65,22,-81]
Warning, need to choose a branch for the root of a polynomial with parameters. This might be wrong.
The choice was done assuming [a,c,d,exp(1),exp(2)]=[1,6,75,-85,-12]

Evaluation time: 10.44
Done
// Time 10.44

And inside sagemath

var('e x d a c')
sage: integrate((e*x+d)^(9/2)/(a*d*e+(a*e^2+c*d^2)*x+c*d*e*x^2),x, algorithm="giac")
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/usr/lib/python3.9/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._div_ (build/cythonized/sage/symbolic/expression.cpp:25839)()
   3894             else:
-> 3895                 x = left._gobj / _right._gobj
   3896             return new_Expression_from_GEx(left._parent, x)

ValueError: power::eval(): division by zero

During handling of the above exception, another exception occurred:

ZeroDivisionError                         Traceback (most recent call last)
/usr/lib/python3.9/site-packages/sage/interfaces/giac.py in _sage_(self, locals)
   1131             try:
-> 1132                 return symbolic_expression_from_string(result, lsymbols,
   1133                     accept_sequence=True)

/usr/lib/python3.9/site-packages/sage/calculus/calculus.py in symbolic_expression_from_string(s, syms, accept_sequence)
   2407                                     if isinstance(v,Function)})
-> 2408     return parse_func(s)
   2409 

/usr/lib/python3.9/site-packages/sage/misc/parser.pyx in sage.misc.parser.Parser.parse_sequence (build/cythonized/sage/misc/parser.c:5837)()
    549 
--> 550     cpdef parse_sequence(self, s):
    551         """

/usr/lib/python3.9/site-packages/sage/misc/parser.pyx in sage.misc.parser.Parser.parse_sequence (build/cythonized/sage/misc/parser.c:5702)()
    565         cdef Tokenizer tokens = Tokenizer(s)
--> 566         all = self.p_sequence(tokens)
    567         if tokens.next() != EOS:

/usr/lib/python3.9/site-packages/sage/misc/parser.pyx in sage.misc.parser.Parser.p_sequence (build/cythonized/sage/misc/parser.c:6450)()
    631             elif token == '(':
--> 632                 obj = self.p_tuple(tokens)
    633             elif token == EOS:

/usr/lib/python3.9/site-packages/sage/misc/parser.pyx in sage.misc.parser.Parser.p_tuple (build/cythonized/sage/misc/parser.c:7237)()
    699                 tokens.reset(start)
--> 700                 return self.p_eqn(tokens)
    701 

/usr/lib/python3.9/site-packages/sage/misc/parser.pyx in sage.misc.parser.Parser.p_eqn (build/cythonized/sage/misc/parser.c:7394)()
    728         """
--> 729         lhs = self.p_expr(tokens)
    730         cdef int op = tokens.next()

/usr/lib/python3.9/site-packages/sage/misc/parser.pyx in sage.misc.parser.Parser.p_expr (build/cythonized/sage/misc/parser.c:7746)()
    768         cdef int op
--> 769         operand1 = self.p_term(tokens)
    770         op = tokens.next()

/usr/lib/python3.9/site-packages/sage/misc/parser.pyx in sage.misc.parser.Parser.p_term (build/cythonized/sage/misc/parser.c:8150)()
    812             else:
--> 813                 operand1 = operand1 / operand2
    814             op = tokens.next()

/usr/lib/python3.9/site-packages/sage/structure/element.pyx in sage.structure.element.Element.__truediv__ (build/cythonized/sage/structure/element.c:13162)()
   1734         if HAVE_SAME_PARENT(cl):
-> 1735             return (<Element>left)._div_(right)
   1736         if BOTH_ARE_ELEMENT(cl):

/usr/lib/python3.9/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._div_ (build/cythonized/sage/symbolic/expression.cpp:25919)()
   3900             if 'division by zero' in str(msg):
-> 3901                 raise ZeroDivisionError("symbolic division by zero")
   3902             else:

ZeroDivisionError: symbolic division by zero

During handling of the above exception, another exception occurred:

NotImplementedError                       Traceback (most recent call last)
<ipython-input-6-754567723580> in <module>
----> 1 integrate((e*x+d)**(Integer(9)/Integer(2))/(a*d*e+(a*e**Integer(2)+c*d**Integer(2))*x+c*d*e*x**Integer(2)),x, algorithm="giac")

/usr/lib/python3.9/site-packages/sage/misc/functional.py in integral(x, *args, **kwds)
    757     """
    758     if hasattr(x, 'integral'):
--> 759         return x.integral(*args, **kwds)
    760     else:
    761         from sage.symbolic.ring import SR

/usr/lib/python3.9/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.integral (build/cythonized/sage/symbolic/expression.cpp:66867)()
  12645                     R = SR
  12646             return R(integral(f, v, a, b, **kwds))
> 12647         return integral(self, *args, **kwds)
  12648 
  12649     integrate = integral

/usr/lib/python3.9/site-packages/sage/symbolic/integration/integral.py in integrate(expression, v, a, b, algorithm, hold)
    988         if not integrator:
    989             raise ValueError("Unknown algorithm: %s" % algorithm)
--> 990         return integrator(expression, v, a, b)
    991     if a is None:
    992         return indefinite_integral(expression, v, hold=hold)

/usr/lib/python3.9/site-packages/sage/symbolic/integration/external.py in giac_integrator(expression, v, a, b)
    446         return expression.integrate(v, a, b, hold=True)
    447     else:
--> 448         return result._sage_()

/usr/lib/python3.9/site-packages/sage/interfaces/giac.py in _sage_(self, locals)
   1134 
   1135             except Exception:
-> 1136                 raise NotImplementedError("Unable to parse Giac output: %s" % result)
   1137         else:
   1138             return [entry.sage() for entry in self]

NotImplementedError: Unable to parse Giac output: ((-4*a^8*c^2*d^2*exp(2)^8+2*a^8*sqrt(-c^2*d^3-c*d*sqrt(c^2*d^4-4*a*c*d^2*exp(1)^2+a^2*exp(2)^2+2*a*c*d^2*exp(2))+a*c*d*exp(2))*sqrt(2)*sqrt(c^2*d^4-4*a*c*d^2*exp(1)^2+a^2*exp(2)^2+2*a*c*d^2*exp(2))*exp(2)^8+44*a^7*c^3*d^4*exp(1)^2*exp(2)^6-12*a^7*c^3*d^4*exp(2)^7-22*a^7*c*d^2*sqrt(-c^2*d^3-c*d*sqrt(c^2*d^4-4*a*c*d^2*exp(1)^2+a^2*exp(2)^2+2*a*c*d^2*exp(2))+a*c*d*exp(2))*sqrt(2)*sqrt(c^2*d^4-

... pages and pages of output like the above follows....

p(1))*d^9*c^6-4*sqrt(d+x*exp(1))*d^7*exp(1)^2*c^5*a-2*sqrt(d+x*exp(1))*d^7*c^5*a*exp(2)+4*sqrt(d+x*exp(1))*d^5*exp(1)^2*c^4*a^2*exp(2)+2*sqrt(d+x*exp(1))*d^5*c^4*a^2*exp(2)^2-2*sqrt(d+x*exp(1))*d^3*c^3*a^3*exp(2)^3)/d^7/c^7
sage:

My question is: Where did all this output come from? The ones starting after Unable to parse Giac output:.

I know because giac has warning messages, this could have confused sagemath in parsing the return value, and this might explain the division by zero. There is a question on this already, when giac returns warning messages.

But my question here is different. Since giac did not give any output from the integrate command, I was wondering then where did this very long output came from?

edit retag flag offensive close merge delete

Comments

I cannot answer your main question, but have a few remarks. First, this works fine with libgiac. Secondly, there is a problem that the Giac interface cannot distinguish between the number $e$ and the variable $e$, so you may want to avoid using $e$ as a variable; see https://trac.sagemath.org/ticket/30133. Finally, this integral also works fine with Giac 1.6.0.47 (which comes with the Sage 9.3 distribution), apparently because $e$ seems to be treated differently there than with Giac 1.7.0.

mwageringel gravatar imagemwageringel ( 2021-05-31 21:55:34 +0100 )edit

this works fine with libgiac hummm. I am using 9.3 sagemath with giac 1.7 and it does not work for me. It throws an exception

sage: from sage.libs.giac import libgiac
// Giac share root-directory:/usr/share/giac/
// Giac share root-directory:/usr/share/giac/
Added 0 synonyms
sage: var('e x d a c')
(e, x, d, a, c)
sage: integrand=(e*x+d)^(9/2)/(a*d*e+(a*e^2+c*d^2)*x+c*d*e*x^2)
sage: anti=libgiac.integrate(integrand,x).sage()
   NotImplementedError: Unable to parse Giac output: <built-in method type of sage.libs.giac.giac.Pygen object at 0x7f8f47dd19d0>
Result is too big for Display. If you really want to see it use print

And when doing print(anti) nothing is displayed.

Nasser gravatar imageNasser ( 2021-06-01 11:59:54 +0100 )edit

As for changing the letter e, I am running the Rubi integration test suite as is. This is the same input used for all CAS systems and not practical for me to edit these files and change them differently for each CAS and change letters. Also, inside giac itself, it works using e. But as I mentioned, I see no output in giac.

Welcome to giac readline interface, version 1.7.0
1>> integrate((e*x+d)^(9/2)/(a*d*e+(a*e^2+c*d^2)*x+c*d*e*x^2),x)
Warning, need to choose a branch for the root of a polynomial with parameters. This might be wrong.
The choice was done assuming [a,c,d,exp(1),exp(2)]=[92,8,-70,-12,80]
Evaluation time: 9.75
Done
// Time 9.75
2>>
Nasser gravatar imageNasser ( 2021-06-01 12:05:15 +0100 )edit

Indeed, I get the same error with libgiac now. I am not sure what I have done differently before. Though, if I replace e by something else, then it works.

Regarding the letter e, I agree this is a problem. Hopefully it will be fixed by the ticket I linked to above, if someone has enough interest in fixing this bug. Feel free to take a look at the source code.

mwageringel gravatar imagemwageringel ( 2021-06-04 19:12:08 +0100 )edit