Simple integral raises AttributeError

157 views
Skip to first unread message

Eric Gourgoulhon

unread,
Dec 4, 2019, 10:38:34 AM12/4/19
to sage-devel
Hi,

In Sage 9.0.beta8 we have

sage: a = var('a')
sage: integrate(1/(x^4 + x^2 + a), x)
...
AttributeError: 'RootSum' object has no attribute '_sage_'

The same error occurs in Sage 8.9, but not in Sage 8.8 (and below). In Sage 8.8, we have instead:

sage: a = var('a')
sage: integrate(1/(x^4 + x^2 + a), x)
integrate(1/(x^4 + x^2 + a), x)

Note that RootSum is a SymPy object. Actually, in Sage 9.0.beta8, forcing the algorithm to 'maxima' yields the same result as in Sage 8.8:

sage: integrate(1/(x^4 + x^2 + a), x, algorithm='maxima')
integrate(1/(x^4 + x^2 + a), x)

So it seems that since Sage 8.9, when integrate() is not capable to find an answer via Maxima, it tries SymPy but is not capable to translate the result back to Sage. I could not find a ticket about this. Shall I open one?

Eric.

PS: for the record, a primitive of 1/(x^4 + x^2 + a) is

sage: b = sqrt(1 - 4*a)
sage: f = sqrt(2)/b*(arctan(sqrt(2)*x/sqrt(1 - b))/sqrt(1 - b)  - arctan(sqrt(2)*x/sqrt(1 + b))/sqrt(1 + b))

as we can check:

sage: diff(f, x).simplify_full()
1/(x^4 + x^2 + a)


Dima Pasechnik

unread,
Dec 4, 2019, 12:37:15 PM12/4/19
to sage-devel
sure, please open a ticket.

--
You received this message because you are subscribed to the Google Groups "sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/5e7fbd5a-f422-4268-b2fb-f77d58c4a2e2%40googlegroups.com.

Emmanuel Charpentier

unread,
Dec 4, 2019, 2:14:48 PM12/4/19
to sage-devel
The old (> 5 years...) Trac#16816 ticket is germane...

BTW, I'm not sure that this would be very helpful:

sage: foo=sympy.integrate(*[sympy.sympify(u) for u in (1/(x^4+x^2+a), x)]); foo
RootSum(_t**4*(256*a**3 - 128*a**2 + 16*a) + _t**2*(4 - 16*a) + 1, Lambda(_t, _t*log(32*_t**3*a**2 - 8*_t**3*a + 4*_t*a - 2*_t + x)))
sage: foo.doit()._sage_().factor().simplify_full()
1/4*sqrt(2)*(sqrt((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))*log(1/4*sqrt(2)*((4*a^2 - a)*((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))^(3/2) + 2*sqrt(2)*x + 2*(2*a - 1)*sqrt((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a)))) - sqrt((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))*log(-1/4*sqrt(2)*((4*a^2 - a)*((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))^(3/2) - 2*sqrt(2)*x + 2*(2*a - 1)*sqrt((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a)))) + sqrt((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))*log(1/4*sqrt(2)*((4*a^2 - a)*((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))^(3/2) + 2*sqrt(2)*x + 2*(2*a - 1)*sqrt((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a)))) - sqrt((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))*log(-1/4*sqrt(2)*((4*a^2 - a)*((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))^(3/2) - 2*sqrt(2)*x + 2*(2*a - 1)*sqrt((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a)))))

Cheers ! ...

However, Sage can come to the rescue. We can compute the root sum *in Sage*, thus getting a "more interesting" result:

sage: S1=sum([u[0]^u[1] for u in foo.args[0]._sage_().roots()]).expand().simplify(); S1
1/2

And (manually) use this result in the second expression. Here, I'm stuck, because I do not (yet) understand how the lambda expression, second argument of RooSum, is supposed to be used.

BTW, fricas can integrate this expression:

sage: IF=f(x).integrate(x,algorithm="fricas").expand().factor().simplify_full(); 
....: IF
1/4*sqrt(2)*(sqrt(((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) + 1)/(4*a^2 - a))*log(1/2*sqrt(2)*(2*sqrt(2)*x + ((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) + 4*a - 1)*sqrt(((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) + 1)/(4*a^2 - a)))) - sqrt(((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) + 1)/(4*a^2 - a))*log(1/2*sqrt(2)*(2*sqrt(2)*x - ((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) + 4*a - 1)*sqrt(((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) + 1)/(4*a^2 - a)))) - sqrt(-((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) - 1)/(4*a^2 - a))*log(1/2*sqrt(2)*(2*sqrt(2)*x + ((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) - 4*a + 1)*sqrt(-((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) - 1)/(4*a^2 - a)))) + sqrt(-((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) - 1)/(4*a^2 - a))*log(1/2*sqrt(2)*(2*sqrt(2)*x - ((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) - 4*a + 1)*sqrt(-((4*a^2 - a)*sqrt(-1/(4*a^3 - a^2)) - 1)/(4*a^2 - a)))))

And fiddling with sympy.cse suggests that the reference given might be correct:

sage: sympy.cse(sympy.sympify(IF), optimizations="basic")
([(x0, sqrt(2)),
  (x1, 4*a),
  (x2, x1 - 1),
  (x3, 1/x2),
  (x4, a*x2*sqrt(-x3/a**2)),
  (x5, x4 + 1),
  (x6, x3/a),
  (x7, sqrt(x5*x6)),
  (x8, 2*x*x0),
  (x9, x7*(x2 + x4)),
  (x10, x0/2),
  (x11, sqrt(-x6*(x4 - 1))),
  (x12, x11*(-x1 + x5))],
 [-x0*(-x11*log(x10*(-x12 + x8)) + x11*log(x10*(x12 + x8)) + x7*log(x10*(x8 - x9)) - x7*log(x10*(x8 + x9)))/4])

HTH,

Emmanuel Charpentier

unread,
Dec 4, 2019, 2:17:51 PM12/4/19
to sage-devel
And, BTW, another data point:

sage: f(x).integrate(x, algorithm="mathematica_free")
-1/2*sqrt(2)*(sqrt(-sqrt(-4*a + 1) + 1)*arctan(sqrt(2)*x/sqrt(sqrt(-4*a + 1) + 1)) - sqrt(sqrt(-4*a + 1) + 1)*arctan(sqrt(2)*x/sqrt(-sqrt(-4*a + 1) + 1)))/sqrt(-(4*a - 1)*a)


HTH,

Eric Gourgoulhon

unread,
Dec 4, 2019, 5:11:56 PM12/4/19
to sage-devel
Le mercredi 4 décembre 2019 18:37:15 UTC+1, Dima Pasechnik a écrit :
sure, please open a ticket.

This is now
Please review.

Eric.

Eric Gourgoulhon

unread,
Dec 4, 2019, 5:17:36 PM12/4/19
to sage-devel
Le mercredi 4 décembre 2019 20:14:48 UTC+1, Emmanuel Charpentier a écrit :
The old (> 5 years...) Trac#16816 ticket is germane...

Thanks for pointing out this ticket!
In #28842, I propose a quick fix until #16816 is solved, i.e. until sum of roots exist in Sage.
This fix restores the behavior prior to Sage 8.9, namely it returns the unevaluated integral.

Eric.

Dima Pasechnik

unread,
Dec 4, 2019, 5:21:39 PM12/4/19
to sage-devel
On Wed, Dec 4, 2019 at 7:14 PM Emmanuel Charpentier
<emanuel.c...@gmail.com> wrote:
>
> The old (> 5 years...) Trac#16816 ticket is germane...
>
> BTW, I'm not sure that this would be very helpful:
>
> sage: foo=sympy.integrate(*[sympy.sympify(u) for u in (1/(x^4+x^2+a), x)]); foo
> RootSum(_t**4*(256*a**3 - 128*a**2 + 16*a) + _t**2*(4 - 16*a) + 1, Lambda(_t, _t*log(32*_t**3*a**2 - 8*_t**3*a + 4*_t*a - 2*_t + x)))
> sage: foo.doit()._sage_().factor().simplify_full()
> 1/4*sqrt(2)*(sqrt((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))*log(1/4*sqrt(2)*((4*a^2 - a)*((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))^(3/2) + 2*sqrt(2)*x + 2*(2*a - 1)*sqrt((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a)))) - sqrt((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))*log(-1/4*sqrt(2)*((4*a^2 - a)*((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))^(3/2) - 2*sqrt(2)*x + 2*(2*a - 1)*sqrt((4*a + sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a)))) + sqrt((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))*log(1/4*sqrt(2)*((4*a^2 - a)*((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))^(3/2) + 2*sqrt(2)*x + 2*(2*a - 1)*sqrt((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a)))) - sqrt((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))*log(-1/4*sqrt(2)*((4*a^2 - a)*((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a))^(3/2) - 2*sqrt(2)*x + 2*(2*a - 1)*sqrt((4*a - sqrt(-64*a^3 + 48*a^2 - 12*a + 1) - 1)/(16*a^3 - 8*a^2 + a)))))
>
> Cheers ! ...
>
> However, Sage can come to the rescue. We can compute the root sum *in Sage*, thus getting a "more interesting" result:
>
> sage: S1=sum([u[0]^u[1] for u in foo.args[0]._sage_().roots()]).expand().simplify(); S1
> 1/2
>
> And (manually) use this result in the second expression. Here, I'm stuck, because I do not (yet) understand how the lambda expression, second argument of RooSum, is supposed to be used.

RootSum(f(t), Lambda(y,g(y)) means

sum_{y: f(y)=0} g(y)

so we sum the function g() in Lambda over all the roots of the 1st
argument, f().
> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/8d0afbd4-0431-451d-8dde-466e19550490%40googlegroups.com.

Emmanuel Charpentier

unread,
Dec 6, 2019, 10:28:40 PM12/6/19
to sage-devel
Okay. That means that we need
  • a wrapper for RootSum, AND
  • a wrapper for Lambda.
ISTR that Mathematica has a similar setup.A rare occasion to kill two birds wit the same (two) stones. However, Mathematica's "pure functions" are, as far as I can tell, totally unknown to Sage...

Thank you, Dima !
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-...@googlegroups.com.

rjf

unread,
Dec 8, 2019, 9:16:10 PM12/8/19
to sage-devel
FYI: 
In Maxima,  integrate_use_rootsof: true
changes the result from an unevaluated integral to
lsum(log(x-%r1)/(4*%r1^3+2*%r1),%r1,rootsof(a+%r1^4+%r1^2,%r1))

which is fairly compact, and to the point, but
which may or may not be what you want.  The inside of  rootsof
can be picked out and the 4 roots solved in terms of radicals
by Maxima's solve().  
If the expression rootsof() is replaced by a list of the roots
returns by solve,  { something like   substpart( map(rhs, solve( Z,3,1,x)), Z 3), lsum   }
you get a large explicit mess of nested radicals and logs.

That Mathematical result with arctan is nicer that the explicit
radical stuff.

RJF

Emmanuel Charpentier

unread,
Dec 9, 2019, 6:32:42 AM12/9/19
to sage-devel

Le lundi 9 décembre 2019 03:16:10 UTC+1, rjf a écrit :

FYI: 
In Maxima,  integrate_use_rootsof: true
changes the result from an unevaluated integral to
lsum(log(x-%r1)/(4*%r1^3+2*%r1),%r1,rootsof(a+%r1^4+%r1^2,%r1))

which is fairly compact, and to the point, but
which may or may not be what you want.  The inside of  rootsof
can be picked out and the 4 roots solved in terms of radicals
by Maxima's solve().  *

Aciddentally ! For example, x^5+x^3+1 doesn’t have a solution in terms of radicals. Trying to solve it in Maxima fails ; and %solve gives a list of numerical roots.

If the expression rootsof() is replaced by a list of the roots
returns by solve,  { something like   substpart( map(rhs, solve( Z,3,1,x)), Z 3), lsum   }
you get a large explicit mess of nested radicals and logs.

Indeed. And one may add that this solution has some drawbacks problems: computing only one value of the resultt primitive takes several seconds, trying to plot it is pointless…

That Mathematical result with arctan is nicer that the explicit
radical stuff.

“Nicer” depends of what you want to do with the result. In some cases, the approximate answer obtained via %solve is preferable (e. g. plotting). In other cases, the explicit log/roots mess is preferable (e. g. trying to prove that the value of the primitive is real for all real values of the aurument).

The “exact” answers given by Sage via (x^5+x^3+1).roots(ring=QQbar) may be preferable: a radical expression, if it exists, may be obtained via the radical_expression() method, and the precision is not limited to any specific floating-point approximation ; some equalities may be provable using them.

OTOH, such values won’t pass to Maxima, and the resulting answers cannot be further simplified…

The more general way to do this is :

  • remark that the polynoimial x^5+x^3+1 has five (possibly confounded) roots, and name them r_1,\dots,r_5 ;
  • rewrite it as \displaystyle\prod_{r\in r_1,\dots,\r_5}(x-r),
  • integrate using this rewriting,
  • substitute r_1\dots,r_5 as needed by your various uses of the result.

HTH,

Thierry

unread,
Dec 9, 2019, 8:55:25 AM12/9/19
to sage-...@googlegroups.com
Hi,

On Wed, Dec 04, 2019 at 07:38:34AM -0800, Eric Gourgoulhon wrote:
> Hi,
>
> In Sage 9.0.beta8 we have
>
> sage: a = var('a')
> sage: integrate(1/(x^4 + x^2 + a), x)
> ...
> AttributeError: 'RootSum' object has no attribute '_sage_'
>
> The same error occurs in Sage 8.9, but not in Sage 8.8 (and below). In Sage
> 8.8, we have instead:
>
> sage: a = var('a')
> sage: integrate(1/(x^4 + x^2 + a), x)
> integrate(1/(x^4 + x^2 + a), x)
>
> Note that RootSum is a SymPy object. Actually, in Sage 9.0.beta8, forcing
> the algorithm to 'maxima' yields the same result as in Sage 8.8:
>
> sage: integrate(1/(x^4 + x^2 + a), x, algorithm='maxima')
> integrate(1/(x^4 + x^2 + a), x)
>
> So it seems that since Sage 8.9, when integrate() is not capable to find an
> answer via Maxima, it tries SymPy but is not capable to translate the
> result back to Sage. I could not find a ticket about this. Shall I open one?

For what it worth, the change was done at
https://trac.sagemath.org/ticket/27958

Ciao,
Thierry


> Eric.
>
> PS: for the record, a primitive of 1/(x^4 + x^2 + a) is
>
> sage: b = sqrt(1 - 4*a)
> sage: f = sqrt(2)/b*(arctan(sqrt(2)*x/sqrt(1 - b))/sqrt(1 - b) -
> arctan(sqrt(2)*x/sqrt(1 + b))/sqrt(1 + b))
>
> as we can check:
>
> sage: diff(f, x).simplify_full()
> 1/(x^4 + x^2 + a)
>
>
> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/5e7fbd5a-f422-4268-b2fb-f77d58c4a2e2%40googlegroups.com.

Dima Pasechnik

unread,
Dec 9, 2019, 12:14:53 PM12/9/19
to sage-devel
On Mon, Dec 9, 2019 at 1:55 PM Thierry <sage-goo...@lma.metelu.net> wrote:
>
> Hi,
>
> On Wed, Dec 04, 2019 at 07:38:34AM -0800, Eric Gourgoulhon wrote:
> > Hi,
> >
> > In Sage 9.0.beta8 we have
> >
> > sage: a = var('a')
> > sage: integrate(1/(x^4 + x^2 + a), x)
> > ...
> > AttributeError: 'RootSum' object has no attribute '_sage_'
> >
> > The same error occurs in Sage 8.9, but not in Sage 8.8 (and below). In Sage
> > 8.8, we have instead:
> >
> > sage: a = var('a')
> > sage: integrate(1/(x^4 + x^2 + a), x)
> > integrate(1/(x^4 + x^2 + a), x)

Fricas just returns the answer in radicals, if you call
sage: integrate(1/(x^4 + x^2 + a), x, algorithm='fricas')

with
sage: integrate(1/(x^5 + x^3 + a), x, algorithm='fricas')
one gets an error, for a reason similar to the RootSum thing for sympy
(i.e. it can express the integral as a sum over roots of a polynomial, it's
just Sage doesn't know how to parse it)
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/20191209135522.n7ymn5yfpcurdppq%40metelu.net.

Emmanuel Charpentier

unread,
Dec 9, 2019, 5:35:05 PM12/9/19
to sage-devel
Well... My solution has problems,at list with 8.9.beta9 and also on Sagecell, which is currently at 9.0) . Consider :
 
L=flatten([[u[0]]*u[1] for u in (x^5+x^3+1).roots(ring=QQbar)]); print("L = ",L)
P=prod([x-u for u in L]);print("P = ",P)
Rac=[var("r_{}".format(u)) for u in range(len(L))];print("R&c = ",Rac)
D=dict(zip(Rac,L));print("D = ",D)
Pr=prod([x-u for u in Rac]).expand();print("Pr = ",Pr)
%time print(Pr.coefficient(x,2).subs(D)==0)

But %time print(bool(Pr.coefficient(x,2).subs(D)==0)) never returns. Other trials show that trying to use SR methods on expressions including QQbar coeficients fails.

Is this known ? It seems to me that it's either a bug or a design flaw.

Advice ?
Reply all
Reply to author
Forward
0 new messages