What is the status of sympy in sagemath? In particular the integrator.

[Nasser M. Abbasi]

I use sagemath to run the independent CAS integrations tests for Fricas, Giac and Maxima, since it is much easier to use the same script to all CAS systems instead of learning how to use each separate CAS. The result is put on this page.

I found that sympy now can be used from sagemath. 

So I said, great. Instead of having separate script for sympy in python will use the same sagemath script and just change the name of the algorithm to 'sympy'. Makes life easier.

But when I tried this on one test file, I found many integrals now fail, where they work using sympy directly in Python.

I am not sure if this is because sympy is not yet fully yet supported in sagemath or if this is just a bug and overlooked support.  

For example, on this one file,  sympy used to score 84.66% passing score when used directly, but now in sagemath it scores 65.64%.  

This translates to about 30 more integrals failing in this file of 163 integrals.

Below will give one example. All seem to give the same exception

NotImplementedError('conversion to SageMath is not implemented')

Here is one example in sagemath 9.8

var('A B a alpha b beta m n x ')

integrate(x/((b*x^2+a)^m),x, algorithm='sympy')

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In [2], line 1
----> 1 integrate(x/(b*x**Integer(3)+a)**Integer(2),x, algorithm='sympy')

File ~/TMP/sage-9.8/src/sage/misc/functional.py:773, in integral(x, *args, **kwds)
    648 """
    649 Return an indefinite or definite integral of an object ``x``.
    650
   (...)
    770
    771 """
    772 if hasattr(x, 'integral'):
--> 773     return x.integral(*args, **kwds)
    774 else:
    775     from sage.symbolic.ring import SR

File ~/TMP/sage-9.8/src/sage/symbolic/expression.pyx:13211, in sage.symbolic.expression.Expression.integral()
  13209                 R = SR
  13210         return R(integral(f, v, a, b, **kwds))
> 13211     return integral(self, *args, **kwds)
  13212
  13213 integrate = integral

File ~/TMP/sage-9.8/src/sage/symbolic/integration/integral.py:1063, in integrate(expression, v, a, b, algorithm, hold)
   1061     if not integrator:
   1062         raise ValueError("Unknown algorithm: %s" % algorithm)
-> 1063     return integrator(expression, v, a, b)
   1064 if a is None:
   1065     return indefinite_integral(expression, v, hold=hold)

File ~/TMP/sage-9.8/src/sage/symbolic/integration/external.py:69, in sympy_integrator(expression, v, a, b)
     67 else:
     68     result = sympy.integrate(ex, (v, a._sympy_(), b._sympy_()))
---> 69 return result._sage_()

File ~/TMP/sage-9.8/src/sage/interfaces/sympy.py:216, in _sympysage_add(self)
    214 s = 0
    215 for x in self.args:
--> 216     s += x._sage_()
    217 return s

File ~/TMP/sage-9.8/local/var/lib/sage/venv-python3.11.1/lib/python3.11/site-packages/sympy/core/basic.py:1959, in Basic._sage_(self)
   1957 sympy_init()  # may monkey-patch _sage_ method into self's class or superclasses
   1958 if old_method == self._sage_:
-> 1959     raise NotImplementedError('conversion to SageMath is not implemented')
   1960 else:
   1961     # call the freshly monkey-patched method
   1962     return self._sage_()

Here is same integral in sympy itself. You see it works.

>python
Python 3.10.9 (main, Dec 19 2022, 17:35:49) [GCC 12.2.0] on linux
>>> from sympy import *
>>> A,B,a,alpha,b,beta,m,n,x= symbols('A B a alpha b beta m n x ')
>>> integrate(x/(b*x**3+a)**2,x)

x**2/(3*a**2 + 3*a*b*x**3) + RootSum(729*_t**3*a**4*b**2 + 1, Lambda(_t, _t*log(81*_t**2*a**3*b + x)))

The sympy version is 1.11.1 in both cases, all on Linux.

age: ver = installed_packages()
sage: ver['sympy']
'1.11.1'

Will give the list of failed integrals in this one file in a follow up post.

--Nasser

[Dima Pasechnik]

Thanks for showing this. As far as I know, the problem is that Sage does not support RootSum expressions - although they are very useful for integration in particular.

--
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/f756dced-6c0b-41cd-b510-6df90906629an%40googlegroups.com.

[Nasser M. Abbasi]

Not all failed one is due to RootSum. Here is one that fails in sagemage but works in sympy and does not generate RootSum but Piecewise

var('A B a alpha b beta m n x ')

integrate(x^(1/2)/(b*x+a),x, algorithm="sympy")

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)

Cell In [7], line 1
----> 1 integrate(x**(Integer(1)/Integer(2))/(b*x+a),x, algorithm="sympy")

File ~/TMP/sage-9.8/src/sage/interfaces/sympy.py:603, in _sympysage_piecewise(self)
    588 """
    589 EXAMPLES::
    590
   (...)
    600     -y*z + cases(((log(x) != 0, x^y/log(x)), (1, y)))
    601 """
    602 from sage.functions.other import cases
--> 603 return cases([(p.cond._sage_(),p.expr._sage_()) for p in self.args])

File ~/TMP/sage-9.8/src/sage/interfaces/sympy.py:603, in <listcomp>(.0)
    588 """
    589 EXAMPLES::
    590
   (...)
    600     -y*z + cases(((log(x) != 0, x^y/log(x)), (1, y)))
    601 """
    602 from sage.functions.other import cases
--> 603 return cases([(p.cond._sage_(),p.expr._sage_()) for p in self.args])

File ~/TMP/sage-9.8/local/var/lib/sage/venv-python3.11.1/lib/python3.11/site-packages/sympy/core/basic.py:1959, in Basic._sage_(self)
   1957 sympy_init()  # may monkey-patch _sage_ method into self's class or superclasses
   1958 if old_method == self._sage_:
-> 1959     raise NotImplementedError('conversion to SageMath is not implemented')
   1960 else:
   1961     # call the freshly monkey-patched method
   1962     return self._sage_()

NotImplementedError: conversion to SageMath is not implemented
sage:

Here is same integral in sympy

>>> A,B,a,alpha,b,beta,m,n,x= symbols('A B a alpha b beta m n x ')

>>> integrate(S("x**(1/2)/(b*x+a)"),x)
Piecewise((zoo*sqrt(x), Eq(a, 0) & Eq(b, 0)), (2*x**(3/2)/(3*a), Eq(b, 0)), (2*sqrt(x)/b, Eq(a, 0)), (-a*log(sqrt(x) - sqrt(-a/b))/(b**2*sqrt(-a/b)) + a*log(sqrt(x) + sqrt(-a/b))/(b**2*sqrt(-a/b)) + 2*sqrt(x)/b, True))
>>>

Or without using S

>>> integrate(x**Rational(1/2)/(b*x+a),x)
Piecewise((zoo*sqrt(x), Eq(a, 0) & Eq(b, 0)), (2*x**(3/2)/(3*a), Eq(b, 0)), (2*sqrt(x)/b, Eq(a, 0)), (-a*log(sqrt(x) - sqrt(-a/b))/(b**2*sqrt(-a/b)) + a*log(sqrt(x) + sqrt(-a/b))/(b**2*sqrt(-a/b)) + 2*sqrt(x)/b, True))

I did not have the time to check each integral to find why it fails but will post list of all failed ones from this one file next.

--Nasser

[Dima Pasechnik]

looks like Piecewise () is another non-implemented conversion - although it's probably easy to fix.

To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/38d80cba-9e43-4f8b-afde-5d971f159121n%40googlegroups.com.

[Oscar Benjamin]

One thing Sage could do with SymPy's RootSum is to call doit which
will expand using radical formulae if possible:

x**2/(3*a**2 + 3*a*b*x**3) + RootSum(729*_t**3*a**4*b**2 + 1,
Lambda(_t, _t*log(81*_t**2*a**3*b + x)))

In [37]: x, a, b, _t = symbols('x, a, b, _t')

In [38]: expr = x**2/(3*a**2 + 3*a*b*x**3) +
RootSum(Poly(729*_t**3*a**4*b**2 + 1, _t), Lambda(_t,
_t*log(81*_t**2*a**3*b + x)))

In [39]: print(expr)

x**2/(3*a**2 + 3*a*b*x**3) + RootSum(729*_t**3*a**4*b**2 + 1,
Lambda(_t, _t*log(81*_t**2*a**3*b + x)))

In [40]: print(expr.doit())
x**2/(3*a**2 + 3*a*b*x**3) +
(-1/(a**4*b**2))**(1/3)*log(a**3*b*(-1/(a**4*b**2))**(2/3) + x)/9 +
(-(-1/(a**4*b**2))**(1/3)/18 -
sqrt(3)*I*(-1/(a**4*b**2))**(1/3)/18)*log(81*a**3*b*(-(-1/(a**4*b**2))**(1/3)/18
- sqrt(3)*I*(-1/(a**4*b**2))**(1/3)/18)**2 + x) +
(-(-1/(a**4*b**2))**(1/3)/18 +
sqrt(3)*I*(-1/(a**4*b**2))**(1/3)/18)*log(81*a**3*b*(-(-1/(a**4*b**2))**(1/3)/18
+ sqrt(3)*I*(-1/(a**4*b**2))**(1/3)/18)**2 + x)

A more precise instruction would be:

expr.replace(lambda e: isinstance(e, RootSum), lambda e: e.doit(deep=False))

--
Oscar

> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/CAAWYfq35yV9jK1mH5%3Dq_AHwr5DgGaYe5JQuwhnwGanjoYK4ajw%40mail.gmail.com.

[Nasser M. Abbasi]

Here is list of integrals that fail in sagemath from this one file.  It seems Piecewise and RootSum are the cause of these failures.

sagemath version

var('A B a alpha b beta m n x ')
integrate(x/((b*x^2+a)^m),x, algorithm="sympy")

integrate(1/(b*x^3+a),x, algorithm="sympy")
integrate(x/(b*x^3+a),x, algorithm="sympy")
integrate(x^3/(b*x^3+a),x, algorithm="sympy")
integrate(x^4/(b*x^3+a),x, algorithm="sympy")
integrate(1/(b*x^3+a)^2,x, algorithm="sympy")
integrate(x/(b*x^3+a)^2,x, algorithm="sympy")
integrate(x^3/(b*x^3+a)^2,x, algorithm="sympy")
integrate(1/x^2/(b*x^3+a),x, algorithm="sympy")
integrate(1/x^3/(b*x^3+a),x, algorithm="sympy")
integrate(1/x^2/(b*x^3+a)^2,x, algorithm="sympy")
integrate(1/x^3/(b*x^3+a)^2,x, algorithm="sympy")
integrate(1/(b*x^4+a),x, algorithm="sympy")
integrate(x^2/(b*x^4+a),x, algorithm="sympy")
integrate(1/(b*x^4+a)^2,x, algorithm="sympy")
integrate(x^2/(b*x^4+a)^2,x, algorithm="sympy")
integrate(1/x^2/(b*x^4+a),x, algorithm="sympy")
integrate(1/(b*x+a)/x^(1/2),x, algorithm="sympy")
integrate(x^(1/2)/(b*x+a),x, algorithm="sympy")
integrate(x^(3/2)/(b*x+a),x, algorithm="sympy")
integrate(x^(5/2)/(b*x+a),x, algorithm="sympy")
integrate(1/(b*x+a)^2/x^(1/2),x, algorithm="sympy")
integrate(x^(1/2)/(b*x+a)^2,x, algorithm="sympy")
integrate(x^(3/2)/(b*x+a)^2,x, algorithm="sympy")
integrate(x^(5/2)/(b*x+a)^2,x, algorithm="sympy")
integrate(1/(b*x+a)^3/x^(1/2),x, algorithm="sympy")
integrate(x^(1/2)/(b*x+a)^3,x, algorithm="sympy")
integrate(x^(3/2)/(b*x+a)^3,x, algorithm="sympy")
integrate(x^(5/2)/(b*x+a)^3,x, algorithm="sympy")
integrate(1/x^(1/2)/(b*x^2+a),x, algorithm="sympy")
integrate(x^(1/2)/(b*x^2+a),x, algorithm="sympy")
integrate(x^(3/2)/(b*x^2+a),x, algorithm="sympy")
integrate(x^(5/2)/(b*x^2+a),x, algorithm="sympy")
integrate(1/x^(1/2)/(b*x^2+a)^2,x, algorithm="sympy")
integrate(x^(1/2)/(b*x^2+a)^2,x, algorithm="sympy")
integrate(x^(3/2)/(b*x^2+a)^2,x, algorithm="sympy")
integrate(x^(5/2)/(b*x^2+a)^2,x, algorithm="sympy")
integrate(1/x^(1/2)/(b*x^2+a)^3,x, algorithm="sympy")
integrate(x^(1/2)/(b*x^2+a)^3,x, algorithm="sympy")

sympy version

>python

from sympy import *

A,B,a,alpha,b,beta,m,n,x= symbols('A B a alpha b beta m n x ')

integrate(S("x/((b*x**2+a)**m)"),x)
integrate(S("1/(b*x**3+a)"),x)
integrate(S("x/(b*x**3+a)"),x)
integrate(S("x**3/(b*x**3+a)"),x)
integrate(S("x**4/(b*x**3+a)"),x)
integrate(S("1/(b*x**3+a)**2"),x)
integrate(S("x/(b*x**3+a)**2"),x)
integrate(S("x**3/(b*x**3+a)**2"),x)
integrate(S("1/x**2/(b*x**3+a)"),x)
integrate(S("1/x**3/(b*x**3+a)"),x)
integrate(S("1/x**2/(b*x**3+a)**2"),x)
integrate(S("1/x**3/(b*x**3+a)**2"),x)
integrate(S("1/(b*x**4+a)"),x)
integrate(S("x**2/(b*x**4+a)"),x)
integrate(S("1/(b*x**4+a)**2"),x)
integrate(S("x**2/(b*x**4+a)**2"),x)
integrate(S("1/x**2/(b*x**4+a)"),x)
integrate(S("1/(b*x+a)/x**(1/2)"),x)
integrate(S("x**(1/2)/(b*x+a)"),x)
integrate(S("x**(3/2)/(b*x+a)"),x)
integrate(S("x**(5/2)/(b*x+a)"),x)
integrate(S("1/(b*x+a)**2/x**(1/2)"),x)
integrate(S("x**(1/2)/(b*x+a)**2"),x)
integrate(S("x**(3/2)/(b*x+a)**2"),x)
integrate(S("x**(5/2)/(b*x+a)**2"),x)
integrate(S("1/(b*x+a)**3/x**(1/2)"),x)
integrate(S("x**(1/2)/(b*x+a)**3"),x)
integrate(S("x**(3/2)/(b*x+a)**3"),x)
integrate(S("x**(5/2)/(b*x+a)**3"),x)
integrate(S("1/x**(1/2)/(b*x**2+a)"),x)
integrate(S("x**(1/2)/(b*x**2+a)"),x)
integrate(S("x**(3/2)/(b*x**2+a)"),x)
integrate(S("x**(5/2)/(b*x**2+a)"),x)
integrate(S("1/x**(1/2)/(b*x**2+a)**2"),x)
integrate(S("x**(1/2)/(b*x**2+a)**2"),x)
integrate(S("x**(3/2)/(b*x**2+a)**2"),x)
integrate(S("x**(5/2)/(b*x**2+a)**2"),x)
integrate(S("1/x**(1/2)/(b*x**2+a)**3"),x)
integrate(S("x**(1/2)/(b*x**2+a)**3"),x)

All using sagemath 9.8 and sympy 1.11.1 on Linux

--Nasser

[Martin R]

If I recall correctly, this is what I did for the FriCAS interface.  It would be nice to factor out any common functionality, if possible.

[Oscar Benjamin]

On Thu, 27 Apr 2023 at 06:25, 'Martin R' via sage-devel
<sage-...@googlegroups.com> wrote:
>
> On Wednesday, 26 April 2023 at 21:06:30 UTC+2 Oscar Benjamin wrote:
>>
>> One thing Sage could do with SymPy's RootSum is to call doit which
>> will expand using radical formulae if possible:
>>
>> x**2/(3*a**2 + 3*a*b*x**3) + RootSum(729*_t**3*a**4*b**2 + 1,
>> Lambda(_t, _t*log(81*_t**2*a**3*b + x)))
>>
>> In [37]: x, a, b, _t = symbols('x, a, b, _t')
>>
>> In [38]: expr = x**2/(3*a**2 + 3*a*b*x**3) +
>> RootSum(Poly(729*_t**3*a**4*b**2 + 1, _t), Lambda(_t,
>> _t*log(81*_t**2*a**3*b + x)))
>>
>> In [39]: print(expr)
>> x**2/(3*a**2 + 3*a*b*x**3) + RootSum(729*_t**3*a**4*b**2 + 1,
>> Lambda(_t, _t*log(81*_t**2*a**3*b + x)))
>>
>> In [40]: print(expr.doit())
>> x**2/(3*a**2 + 3*a*b*x**3) +
>> (-1/(a**4*b**2))**(1/3)*log(a**3*b*(-1/(a**4*b**2))**(2/3) + x)/9 +
>> (-(-1/(a**4*b**2))**(1/3)/18 -
>> sqrt(3)*I*(-1/(a**4*b**2))**(1/3)/18)*log(81*a**3*b*(-(-1/(a**4*b**2))**(1/3)/18
>> - sqrt(3)*I*(-1/(a**4*b**2))**(1/3)/18)**2 + x) +
>> (-(-1/(a**4*b**2))**(1/3)/18 +
>> sqrt(3)*I*(-1/(a**4*b**2))**(1/3)/18)*log(81*a**3*b*(-(-1/(a**4*b**2))**(1/3)/18
>> + sqrt(3)*I*(-1/(a**4*b**2))**(1/3)/18)**2 + x)
>
> If I recall correctly, this is what I did for the FriCAS interface. It would be nice to factor out any common functionality, if possible.

Obviously though the RootSum is better than the radicals which is why
it is used in the first place so the ideal solution would be to
preserve the RootSum. The simple case above already shows that it is
better but in others the radical formulae can explode and in more
complicated cases formulae won't even exist.

--
Oscar

[Martin R]

I fully agree.

https://github.com/sagemath/sage/issues/34420

https://github.com/sagemath/sage/issues/32133

https://github.com/sagemath/sage/issues/32143