Repository: sympy/sympy
Pull request: Fix for bug in pow class
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #25165 

#### Brief description of what is fixed or changed
Added fixes for better execution of thee operation to remove the constant term from an expression in _eval_nseries method in Pow class


#### Other comments
Changed some old tests, as the new implementation gives more simplified forms of the expected results.

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* functions
  * Fixed a bug in _eval_nseries method of Pow class to better handle the operation of negating the constant term from the expression.
<!-- END RELEASE NOTES -->


Problem statement:
Series expansion not working
Hello,

In our project we use very large functions and a lot of series expansions throw a "not implemented error". I managed to extract a relatively simple case:

`(1/sqrt(( - y + 1)**2 + (y - 0.23)**4)).series(y,x0=0,n=5)`

When using `apart()` the series expansion works: 

`apart((1/sqrt(( - y + 1)**2 + (y - 0.23)**4))).series(y,x0=0,n=5)`

sadly its not that easy with our larger functions. Any ideas or am I missing something fundamental? 

The error is: 

`Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/home/user/.local/lib/python3.10/site-packages/sympy/core/expr.py", line 3026, in series
    rv = self.subs(x, xpos).series(xpos, x0, n, dir, logx=logx, cdir=cdir)
  File "/home/user/.local/lib/python3.10/site-packages/sympy/core/expr.py", line 3034, in series
    s1 = self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
  File "/home/user/.local/lib/python3.10/site-packages/sympy/core/power.py", line 1758, in _eval_nseries
    raise NotImplementedError()
NotImplementedError`

Hints from issue discussion:
Hey @RagonEbker thanks for reporting this issue. We actually have the necessary tools to address the above case. 

The problem lies here https://github.com/sympy/sympy/blob/aa6db7c688855c181e020929dd34b9639dec64ad/sympy/core/power.py#L1581-L1586

Let's say we want to calculate the series for `(1/sqrt(( - y + 1)**2 + (y - 0.6)**4))`

Hence we encounter a terms `g` which is the following 
```
>>> g
0.885269121813031*(y - 1)**2 + 0.885269121813031*(y - 0.6)**4 - 1.0
```

If we look at this term the constant term here is S.Zero or rather there is no constant term
```
>>> g = 0.885269121813031*(y - 1)**2 + 0.885269121813031*(y - 0.6)**4 - 1.0
>>> (0.885269121813031*(y - 1)**2).expand()
0.885269121813031*y**2 - 1.77053824362606*y + 0.885269121813031
>>> (0.885269121813031*(y - 0.6)**4).expand()
0.885269121813031*y**4 - 2.12464589235127*y**3 + 1.91218130311615*y**2 - 0.764872521246459*y + 0.114730878186969
>>> 0.885269121813031 + 0.114730878186969 - 1.0
0.0
```

But if we calculate the leading term for `g`, we get a constant due to a rounding off error.
```
>>> g.as_leading_term(y)
-2.22044604925031e-16
```
This leads to the wrong answer. @oscarbenjamin what can be done here ? We have all the tools necessary. We just need to expand/simplify `g` such that the constant terms goes to 0 rather than returning a float

I think that `nsimplify` is being called in the wrong place. It is already too late because rounding errors have already happened and so after `nsimplify` the constant term is still not found to be zero.

If this is going to use `nsimplify` to handle floats then it needs to be done right at the beginning of the calculation. Otherwise everything needs to be designed to handle floats throughout.
cc @arnabnandikgp 
As you've been contributing to the series module these days, you might be interested in picking up this issue. We have all the tools required to solve this though we are making some mistakes with simplyfing things.
I was toggling through the git history of the associated lines and realised that the code block 
https://github.com/sympy/sympy/blob/aa6db7c688855c181e020929dd34b9639dec64ad/sympy/core/power.py#L1581-L1586 
was added through #22213  in order to handle the very same negation problems that we are trying to handle in this case but it worked then as the failing test cases weren't beyond the limits of nsimplify(which is 10**-15) but in this example we presumably exceeded that causing us problems.

Commit message: changes added

Signed-off-by: arnabnandikgp <arnabnandi2002@gmail.com>

Patch to review:
diff --git a/sympy/core/power.py b/sympy/core/power.py
index a8dc47bdfe67..d695e7d08320 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1533,7 +1533,8 @@ def _eval_nseries(self, x, n, logx, cdir=0):
                 return res
 
         f = b.as_leading_term(x, logx=logx)
-        g = (b/f - S.One).cancel(expand=False)
+        g = (b.expand() - f).cancel()
+        g = g/f
         if not m.is_number:
             raise NotImplementedError()
         maxpow = n - m*e
@@ -1583,8 +1584,8 @@ def mul(d1, d2):
             # Convert floats like 0.5 to exact SymPy numbers like S.Half, to
             # prevent rounding errors which can induce wrong values of d leading
             # to a NotImplementedError being returned from the block below.
-            from sympy.simplify.simplify import nsimplify
-            _, d = nsimplify(g).leadterm(x, logx=logx)
+            g = g.replace(lambda x: x.is_Float, lambda x: Rational(x))
+            _, d = g.leadterm(x, logx=logx)
         if not d.is_positive:
             g = g.simplify()
             if g.is_zero:
diff --git a/sympy/core/tests/test_expr.py b/sympy/core/tests/test_expr.py
index af83b22ab361..a95144d99cc9 100644
--- a/sympy/core/tests/test_expr.py
+++ b/sympy/core/tests/test_expr.py
@@ -479,14 +479,13 @@ def test_as_leading_term():
     # https://github.com/sympy/sympy/issues/21177
     e = -3*x + (x + Rational(3, 2) - sqrt(3)*S.ImaginaryUnit/2)**2\
         - Rational(3, 2) + 3*sqrt(3)*S.ImaginaryUnit/2
-    assert e.as_leading_term(x) == \
-        (12*sqrt(3)*x - 12*S.ImaginaryUnit*x)/(4*sqrt(3) + 12*S.ImaginaryUnit)
+    assert e.as_leading_term(x) == -sqrt(3)*I*x
 
     # https://github.com/sympy/sympy/issues/21245
     e = 1 - x - x**2
     d = (1 + sqrt(5))/2
     assert e.subs(x, y + 1/d).as_leading_term(y) == \
-        (-576*sqrt(5)*y - 1280*y)/(256*sqrt(5) + 576)
+        (-40*y - 16*sqrt(5)*y)/(16 + 8*sqrt(5))
 
 
 def test_leadterm2():
diff --git a/sympy/core/tests/test_power.py b/sympy/core/tests/test_power.py
index 640834ead0a9..f1080e338de7 100644
--- a/sympy/core/tests/test_power.py
+++ b/sympy/core/tests/test_power.py
@@ -18,6 +18,7 @@
 from sympy.testing.pytest import warns, _both_exp_pow
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.abc import a, b, c, x, y
+from sympy.core.numbers import all_close
 
 def test_rational():
     a = Rational(1, 5)
@@ -661,3 +662,9 @@ def test_issue_26546():
     assert Pow(x+I, Rational(1,2)).is_extended_real is False
     assert Pow(x+I, Rational(1,13)).is_extended_real is False
     assert Pow(x+I, Rational(2,3)).is_extended_real is None
+
+
+def test_issue_25165():
+    e1 = (1/sqrt(( - x + 1)**2 + (x - 0.23)**4)).series(x, 0, 2)
+    e2 = 0.998603724830355 + 1.02004923189934*x + O(x**2)
+    assert all_close(e1, e2)
diff --git a/sympy/series/tests/test_series.py b/sympy/series/tests/test_series.py
index 6ae7ed2b848b..3313b8e6ee48 100644
--- a/sympy/series/tests/test_series.py
+++ b/sympy/series/tests/test_series.py
@@ -367,9 +367,9 @@ def test_issue_20697():
 def test_issue_21245():
     fi = (1 + sqrt(5))/2
     assert (1/(1 - x - x**2)).series(x, 1/fi, 1).factor() == \
-        (-4812 - 2152*sqrt(5) + 1686*x + 754*sqrt(5)*x\
-        + O((x - 2/(1 + sqrt(5)))**2, (x, 2/(1 + sqrt(5)))))/((1 + sqrt(5))\
-        *(20 + 9*sqrt(5))**2*(x + sqrt(5)*x - 2))
+        (-37*sqrt(5) - 83 + 13*sqrt(5)*x + 29*x + O((x - 2/(1 + sqrt(5)))**2, (x\
+            , 2/(1 + sqrt(5)))))/((2*sqrt(5) + 5)**2*(x + sqrt(5)*x - 2))
+
 
 
 def test_issue_21938():


Code review comments:
diff_hunk: @@ -1533,7 +1533,8 @@ def _eval_nseries(self, x, n, logx, cdir=0):
                 return res
 
         f = b.as_leading_term(x, logx=logx)
-        g = (b/f - S.One).cancel(expand=False)
+        g = (b.expand() - f).cancel()
path: sympy/core/power.py
text: @user1:
```suggestion
        g = (_mexpand(b) - f).cancel()
```