Series expansion around float fails with NotImplementedError

[oscarbenjamin]

This works fine with rational S.Half but fails with float 0.5:

In [5]: series(1/sqrt(1-x**2), x, S.Half)
Out[5]: 
                                      2                   3                   4                    
2⋅√3   4⋅√3⋅(x - 1/2)   8⋅√3⋅(x - 1/2)    112⋅√3⋅(x - 1/2)    608⋅√3⋅(x - 1/2)    1088⋅√3⋅(x - 1/2)
──── + ────────────── + ─────────────── + ───────────────── + ───────────────── + ─────────────────
 3           9                 9                  81                 243                 243       

5                         
     ⎛         6         ⎞
─ + O⎝(x - 1/2) ; x → 1/2⎠
                          

In [6]: series(1/sqrt(1-x**2), x, 0.5)
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-6-ab154a1ba9fb> in <module>
----> 1 series(1/sqrt(1-x**2), x, 0.5)

~/current/sympy/sympy.git/sympy/series/series.py in series(expr, x, x0, n, dir)
     61     """
     62     expr = sympify(expr)
---> 63     return expr.series(x, x0, n, dir)

~/current/sympy/sympy.git/sympy/core/expr.py in series(self, x, x0, n, dir, logx, cdir)
   2956                 rep2 = x
   2957                 rep2b = -x0
-> 2958             s = self.subs(x, rep).series(x, x0=0, n=n, dir='+', logx=logx, cdir=cdir)
   2959             if n is None:  # lseries...
   2960                 return (si.subs(x, rep2 + rep2b) for si in s)

~/current/sympy/sympy.git/sympy/core/expr.py in series(self, x, x0, n, dir, logx, cdir)
   2966             # replace x with an x that has a positive assumption
   2967             xpos = Dummy('x', positive=True)
-> 2968             rv = self.subs(x, xpos).series(xpos, x0, n, dir, logx=logx, cdir=cdir)
   2969             if n is None:
   2970                 return (s.subs(xpos, x) for s in rv)

~/current/sympy/sympy.git/sympy/core/expr.py in series(self, x, x0, n, dir, logx, cdir)
   2974         from sympy.series.order import Order
   2975         if n is not None:  # nseries handling
-> 2976             s1 = self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
   2977             o = s1.getO() or S.Zero
   2978             if o:

~/current/sympy/sympy.git/sympy/core/power.py in _eval_nseries(self, x, n, logx, cdir)
   1735             _, d = g.leadterm(x)
   1736             if not d.is_positive:
-> 1737                 raise NotImplementedError()
   1738 
   1739         from sympy.functions.elementary.integers import ceiling

NotImplementedError: 

[praneethratna]

I think whenever series encounters a Float, it should be converted to Sympy Rational.The following change fixes this issue,

diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index e607234d0e..51f1283dd8 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -2933,6 +2933,10 @@ def series(self, x=None, x0=0, n=6, dir="+", logx=None, cdir=0):
             If "x0" is an infinity object

         """
+        x0 = sympify(x0)
+        if isinstance(x0, Float):
+            x0 = Rational(x0)
+
         if x is None:
             syms = self.free_symbols
             if not syms:

[jksuom]

To avoid issues with inexact numbers all Floats should probably be converted to Rationals. In this case, the problem is that g = b/f - 1 has a tiny nonzero constant part (-5.55111512312578e-17) resulting from the inexact division by f, It seems that it could be made to disappear if g were expanded on this line:

sympy/sympy/core/power.py

Line 1689 in 0322154

g = (b/f - S.One).cancel(expand=False)

Expansion was prevented for some reason. (I forget why but maybe @0sidharth would remember.) Expansion could be reintroduced on this "last resort" line:

sympy/sympy/core/power.py

Line 1734 in 0322154

g = g.simplify()

simplify does not help. (I does not make the rounding error disappear.) Instead, g could be recomputed on this line with expansion like in

g = ((b - f)/f).cancel()  # division by f is the last operation to avoid error

or possibly only

g = ((b - f)/f).expand()

[praneethratna]

g = ((b - f)/f).cancel()  # division by f is the last operation to avoid error

or possibly only

g = ((b - f)/f).expand()

With this change,

diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index e607234d0e..51f1283dd8 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -2933,6 +2933,10 @@ def series(self, x=None, x0=0, n=6, dir="+", logx=None, cdir=0):
             If "x0" is an infinity object

         """
+        x0 = sympify(x0)
+        if isinstance(x0, Float):
+            x0 = Rational(x0)
+
         if x is None:
             syms = self.free_symbols
             if not syms:
diff --git a/sympy/core/power.py b/sympy/core/power.py
index b62d839811..5f2ff67321 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1734,7 +1734,7 @@ def mul(d1, d2):
             else:
                 raise NotImplementedError()
         if not d.is_positive:
-            g = g.simplify()
+            g = ((b - f)/f).cancel()
             _, d = g.leadterm(x)
             if not d.is_positive:
                 raise NotImplementedError()

The output is now,

>>> series(1/sqrt(1-x**2), x, 0.5)
2*sqrt(3)/3 + 4*sqrt(3)*(x - 1/2)/9 + 8*sqrt(3)*(x - 1/2)**2/9 + 112*sqrt(3)*(x - 1/2)**3/81 + 608*sqrt(3)*(x - 1/2)**4/243 + 1088*sqrt(3)*(x - 1/2)**5/243 + O((x - 1/2)**6, (x, 1/2))

[praneethratna]

But with the above change there are some failing tests in sympy/series/tests/test_series.py,

________________________________________________________________________________________________________________________
__________________________________ sympy/series/tests/test_series.py:test_issue_12791 __________________________________
Traceback (most recent call last):
  File "/home/sympy/sympy/series/tests/test_series.py", line 237, in test_issue_12791
    assert expr.series(beta, 0.5, 2).trigsimp() == sol
AssertionError
________________________________________________________________________________________________________________________
__________________________________ sympy/series/tests/test_series.py:test_issue_18842 __________________________________
Traceback (most recent call last):
  File "/home/sympy/sympy/series/tests/test_series.py", line 280, in test_issue_18842
    assert f.series(x, 0.491, n=1).removeO().nsimplify() ==  \
AssertionError

[jksuom]

Maybe the simplify line should be left as it is because of these errors. Then the recomputed cancel line could added here

sympy/sympy/core/power.py

Lines 1736 to 1737 in 0322154

	if not d.is_positive:
	raise NotImplementedError()

together with a yet another if not d.is_positive: test before raising NotImplementedError.

[praneethratna]

Do you mean a change something like this? But this change still raises the same errors.

diff --git a/sympy/core/power.py b/sympy/core/power.py
index b62d839811..e4505b274d 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1736,6 +1736,8 @@ def mul(d1, d2):
         if not d.is_positive:
             g = g.simplify()
             _, d = g.leadterm(x)
+            if not d.is_positive:
+                g = ((b - f)/f).expand()
             if not d.is_positive:
                 raise NotImplementedError()

[jksuom]

Also d has to be recomputed.

[praneethratna]

Yes, this change is working, But should'nt all floats be converted to Rationals?

diff --git a/sympy/core/power.py b/sympy/core/power.py
index b62d839811..4b54b19980 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1736,6 +1736,9 @@ def mul(d1, d2):
         if not d.is_positive:
             g = g.simplify()
             _, d = g.leadterm(x)
+            if not d.is_positive:
+                g = ((b - f)/f).cancel()
+                _, d = g.leadterm(x)
             if not d.is_positive:
                 raise NotImplementedError()

Adding code for converting floats to rationals is raising errors,

diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index e607234d0e..51f1283dd8 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -2933,6 +2933,10 @@ def series(self, x=None, x0=0, n=6, dir="+", logx=None, cdir=0):
             If "x0" is an infinity object

         """
+        x0 = sympify(x0)
+        if isinstance(x0, Float):
+            x0 = Rational(x0)
+
         if x is None:
             syms = self.free_symbols
             if not syms:
diff --git a/sympy/core/power.py b/sympy/core/power.py
index b62d839811..4b54b19980 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1736,6 +1736,9 @@ def mul(d1, d2):
         if not d.is_positive:
             g = g.simplify()
             _, d = g.leadterm(x)
+            if not d.is_positive:
+                g = ((b - f)/f).cancel()
+                _, d = g.leadterm(x)
             if not d.is_positive:
                 raise NotImplementedError()

[jksuom]

There might be Floats also elsewhere in the expression. Changing x0 looks like a partial solution. I would recommend the change to Rationals only if nothing else would help.

[praneethratna]

The change that you have suggested:

diff --git a/sympy/core/power.py b/sympy/core/power.py
index b62d839811..12f8f723a6 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1736,6 +1736,9 @@ def mul(d1, d2):
         if not d.is_positive:
             g = g.simplify()
             _, d = g.leadterm(x)
+            if not d.is_positive:
+                g = ((b - f)/f).expand()
+                _, d = g.leadterm(x)
             if not d.is_positive:
                 raise NotImplementedError()

This is causing errors in test_limitseq.py where test_limit_seq() is freezed, this is caused by the _, d = g.leadterm(x) after the recomputed cancel line.

[jksuom]

This is causing errors

This is unfortunately typical when series code is modified. It looks like the expanded g will be different from the earlier version in a bad way. The difference should be investigated.

[praneethratna]

Also on a side note, this output is different from what we get initially, after the change output is

>>> expr = 1/sqrt(1 - x**2)
>>> series(expr, x, 0.5)
0.769800358919501 + 1.539600717839*(x - 0.5)**2 + 2.39493444997178*(x - 0.5)**3 + 4.33369090947275*(x - 0.5)**4 + 7.75502583800386*(x - 0.5)**5 + 0.769800358919501*x + O((x - 1/2)**6, (x, 1/2))

But initially the output was,

>>> series(1/sqrt(1-x**2), x, 0.5)
2*sqrt(3)/3 + 4*sqrt(3)*(x - 1/2)/9 + 8*sqrt(3)*(x - 1/2)**2/9 + 112*sqrt(3)*(x - 1/2)**3/81 + 608*sqrt(3)*(x - 1/2)**4/243 + 1088*sqrt(3)*(x - 1/2)**5/243 + O((x - 1/2)**6, (x, 1/2))

[jksuom]

It seems that the only test in test_limit_seq() that will run this code is this one:

sympy/sympy/series/tests/test_limitseq.py

Lines 90 to 91 in 0322154

	e = Sum(harmonic(k)**2/k, (k, 1, 2*n)) / harmonic(n)**3
	assert limit_seq(e, n) == S.One / 3

where g.simplify() is zero. Currently, NotImplementedError is returned, but it would be more correct to return f**e. It looks like g.is_zero should be tested before recomputing g.

[praneethratna]

The error is resolved now, but now the output is in exact form,

>>> series(1/sqrt(1-x**2), x, 0.5)
0.769800358919501 + 1.539600717839*(x - 0.5)**2 + 2.39493444997178*(x - 0.5)**3 + 4.33369090947275*(x - 0.5)**4 + 7.75502583800386*(x - 0.5)**5 + 0.769800358919501*x + O((x - 1/2)**6, (x, 1/2))

instead of something like this,

>>> series(1/sqrt(1-x**2), x, S.Half)
2*sqrt(3)/3 + 4*sqrt(3)*(x - 1/2)/9 + 8*sqrt(3)*(x - 1/2)**2/9 + 112*sqrt(3)*(x - 1/2)**3/81 + 608*sqrt(3)*(x - 1/2)**4/243 + 1088*sqrt(3)*(x - 1/2)**5/243 + O((x - 1/2)**6, (x, 1/2))

Should'nt the outputs for both series(1/sqrt(1-x**2), x, S.Half) and series(1/sqrt(1-x**2), x, 0.5) be same?

[jksuom]

Should'nt the outputs for both series(1/sqrt(1-x**2), x, S.Half) and series(1/sqrt(1-x**2), x, 0.5) be same?

Not necessarily, I think. Rational(0,5) is very simple (Rational(1,2)) but, in general, a Float converted to a Rational has very big numerator and denominator. The output is inexact though it may look formally exact with complicated rational coefficients.

[faze-geek]

The error is resolved now, but now the output is in exact form,

I guess same as #13061 ( another floats issue ) which earlier gave a canonical result ( exp, sqrt terms preserved) and now giving expanded numerical results once the targeted bug is solved .

Maybe someone familiar could help to handle it -

Earlier -

 In [4]: sympy.inverse_laplace_transform( 1.0/(s**4 + 10*s**3 + 34*s**2 + 46*s + 21), s, t )
Out[4]: 0.5*(2*exp(2*t) - 2*exp(sqrt(2)*t) - sqrt(2)*exp(sqrt(2)*t) + 2 - 2*exp(-sqrt(2)*t) + sqrt(2)*exp(-sqrt(2)*t))*exp(-3*t)*Heaviside(t)/((-2 + sqrt(2))**2*(-1 + sqrt(2))*(1 + sqrt(2))*(sqrt(2) + 2)**2)

Current -

>>> inverse_laplace_transform( 1.0/(s**4 + 10*s**3 + 34*s**2 + 46*s + 21), s, t )
(-(0.0732233047033631*exp(1.58578643762691*t) + 0.426776695296637*exp(4.41421356237309*t))*exp(4.0*t)*Heaviside(-(1 - exp(t))*exp_polar(2*I*pi)*exp(-t)) + 0.25*exp(7.0*t)*Heaviside(t) + 0.25*exp(9.0*t)*Heaviside(t))*exp(-10.0*t)

[ThePauliPrinciple]

How about replacing the float with a dummy variable and then substituting back the Float once the series has been expanded? It might supress this behaviour. Additionally, not every user will be happy if their Float is converted into a Rational, since even for 0.5, it isn't the same as 1/2 in every context.

[0sidharth]

(I forget why but maybe @0sidharth would remember.)

I don't recall exactly either - most probably, it was to prevent a test case from failing, or, although unlikely, added to prevent the unnecessary expansion (at that time) and speed up slightly.