Repository: sympy/sympy
Pull request: Test case added for numerical evaluation in test_simplify
<!-- 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 #11004

#### Brief description of what is fixed or changed


#### Other comments


#### 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 -->
NO ENTRY 
<!-- END RELEASE NOTES -->


Problem statement:
problems with numerical evaluation
I working with a simple statistical mechanics example,
using the Stirling approximation of the factorial in combinatorial functions:

```
from __future__ import division
from sympy.abc import *
from sympy import *
from sympy.physics.units import avogadro_number as N


# Stirling's approximation for factorial
def f(n):
    return sqrt(2*pi*n)*(n/E)**n

# number of k-equipartitions of n particles
def m(n, k):
    return  f(n)/ (f(n/k)**k)

# probability of a k-partition being an equipartition
def p(n,k):
    return m(n, k)/(k**n)

q = (p(n, 2)/p(n, 3)).simplify()

z = log(p(n, k)/p(n, k+1)).expand(force=True)
```

Working tests:

```
>>> print q.n(4)
0.9648*n**0.5
>>> print q.n(4,subs={n: N})
7.487e+11
>>> print z   
k*n*log(n)/(k + 1) - k*n*log(k + 1)/(k + 1) - k*n/(k + 1) + k*log(k)/2 - k*log(k + 1)/2 - n*log(n) + n*log(k + 1) + n + n*log(n)/(k + 1) - n*log(k + 1)/(k + 1) - n/(k + 1) + log(n)/2 - log(k + 1)/2 + log(2)/2 + log(pi)/2
>>> print z.simplify()
k*log(k)/2 - k*log(k + 1)/2 + log(n)/2 - log(k + 1)/2 + log(2)/2 + log(pi)/2
>>> print z.simplify().n(4)
0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189
>>> print z.n(4).simplify()
0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189
>>> print z.subs(n,N).simplify()
k*log(k)/2 - k*log(k + 1)/2 - log(k + 1)/2 + log(2)/2 + log(pi)/2 + log(602214179000000000000000)/2
>>> print z.subs(n,N).simplify().n(4)
0.5*k*log(k) - 0.5*k*log(k + 1) - 0.5*log(k + 1) + 28.3
```

The problem cases:
- Substitution not working `n` is not replaced by `N` - note that `print q.n(4,subs={n: N})` above works

```
>>> print z.n(4, subs={n: N})
k*n*log(n)/(k + 1.0) - k*n*log(k + 1)/(k + 1.0) - k*n/(k + 1.0) + 0.5*k*log(k) - 0.5*k*log(k + 1) - n*log(n) + n*log(k + 1) + n + n*log(n)/(k + 1.0) - n*log(k + 1)/(k + 1.0) - n/(k + 1.0) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189

>>> print z.simplify().n(4, subs={n: N}) 
0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189
```
- This looses the `log(n)` term completely and performs some other invalid transformations

```
>>> print z.subs(n,N).n(4).simplify()
0.5*k*(k*log(k) - k*log(k + 1) + log(k))/(k + 1)
```

Hints from issue discussion:
Nothing stands out to me as being wrong:

```
>>> z.subs(n,2).n(4).simplify()
0.5*k*log(k) - 0.5*k*log(k + 1) - 0.5*log(k + 1) + 1.266
>>> z.simplify().subs(n,2).n(4)
0.5*k*log(k) - 0.5*k*log(k + 1) - 0.5*log(k + 1) + 1.266
>>> z.simplify().n(4, subs={n:2})
0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189
```

The only thing different in case 3 is that the n has not been replaced (because doing so doesn't change the expression to a number that can be evaluated) so it isn't combined in with the numerical term -- thus 0.9189 instead of 1.266.

Is there something that I am missing?

In my last example

```
>>> print z.subs(n,N).n(4).simplify()
0.5*k*(k*log(k) - k*log(k + 1) + log(k))/(k + 1)
```

There is no numerical term that comes from  `N`. It does not depend on `N` at all. Where did the `N` go?
Also when you divide through by `(k + 1)` coefficient of `log(k + 1)` is `0.5*k*k/(k + 1)` which
not equivalent to that in the other evaluations.

In the `z.simplify().n(4, subs={n:2})` case I do not understand the explanation:

> the n has not been replaced (because doing so doesn't change the expression to a number that can be evaluated)

What is the difference between `X.n(4, subs={n:2})` and `X.subs(n,2).n(4)`?
I thought they should be equvalent.

I think the problem is just that evalf doesn't work very well with expressions involving symbols unless all of them are fully substituted for actual numbers.
> What is the difference between `X.n(4, subs={n:2})` and `X.subs(n,2).n(4)`?



> There is no numerical term that comes from `N`. It does not depend on `N` at all. Where did the `N` go?

It is now there in master:
```
>>> simplify(z.subs(n,N).n(4))
0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(avogadro_number) - 0.5*log(k + 1) + 0.9189
```
This issue could probably be closed with a test (or confirmation that a test addressing the root issue has already been added by the commit that fixed this).
I can't find avogadro_number defined in the physics module.

Tried adding test cases to `test_simplify.py`
```python
N=Symbol('N')
def f(n):
    return sqrt(2*pi*n)*(n/E)**n
def m(n, k):
    return  f(n)/ (f(n/k)**k)
def p(n,k):
    return m(n, k)/(k**n)
q = (p(n, 2)/p(n, 3)).simplify()
z = log(p(n, k)/p(n, k+1)).expand(force=True)
assert simplify(z.subs(n,N).n(4)) == 0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(N) - 0.5*log(k + 1) + 0.9189
```
Still, it gives assertion error, but checking for simplify(z.subs(n,N).n(4)), it returns the same expression as above.
It returns `0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(N) - 0.5*log(k + 1) + 0.9189`
What could be the possible reason @smichr @oscarbenjamin @mhlr?

You need to create floats with the correct precision:
```python
In [23]: r = simplify(z.subs(n,N).n(4))

In [24]: half = Float('0.5', 4)

In [25]: r == half*k*log(k) - half*k*log(k + 1) + half*log(N) - half*log(k + 1) + Float(0.9189224, 4)
Out[25]: True
```
Thank you for your assistance.
Kindly review the PR.
> I can't find avogadro_number defined in the physics module.
> 
> Tried adding test cases to `test_simplify.py`
> 
> ```python
> N=Symbol('N')
> def f(n):
>     return sqrt(2*pi*n)*(n/E)**n
> def m(n, k):
>     return  f(n)/ (f(n/k)**k)
> def p(n,k):
>     return m(n, k)/(k**n)
> q = (p(n, 2)/p(n, 3)).simplify()
> z = log(p(n, k)/p(n, k+1)).expand(force=True)
> assert simplify(z.subs(n,N).n(4)) == 0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(N) - 0.5*log(k + 1) + 0.9189
> ```
> 
> Still, it gives assertion error, but checking for simplify(z.subs(n,N).n(4)), it returns the same expression as above. It returns `0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(N) - 0.5*log(k + 1) + 0.9189` What could be the possible reason @smichr @oscarbenjamin @mhlr?

You can find avogadro number using following code
![image](https://user-images.githubusercontent.com/78342519/137574753-a191bc72-708a-45e1-a6ea-dd0e59cb9b7d.png)

Commit message: fixed formatting

Patch to review:
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/tests/test_simplify.py
index 4ed4fc47f6f2..4cf06956408d 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -1036,6 +1036,22 @@ def test_issue_23543():
     x, y, z = symbols("x y z", commutative=False)
     assert (x*(y + z/2)).simplify() == x*(2*y + z)/2
 
+def test_issue_11004():
+    def f(n):
+        return sqrt(2*pi*n) * (n/E)**n
+
+    def m(n, k):
+        return  f(n) / (f(n/k)**k)
+
+    def p(n,k):
+        return m(n, k) / (k**n)
+
+    N, k = symbols('N k')
+    half = Float('0.5', 4)
+    z = log(p(n, k) / p(n, k + 1)).expand(force=True)
+    r = simplify(z.subs(n, N).n(4))
+    assert r == (half*k*log(k) - half*k*log(k + 1) + half*log(N) - half*log(k + 1) + Float(0.9189224, 4))
+
 
 def test_issue_19161():
     polynomial = Poly('x**2').simplify()


Code review comments:
diff_hunk: @@ -1036,6 +1036,22 @@ def test_issue_23543():
     x, y, z = symbols("x y z", commutative=False)
     assert (x*(y + z/2)).simplify() == x*(2*y + z)/2
 
+def test_issue_11004():
+    def f(n):
+        return sqrt(2*pi*n) * (n/E)**n
+
+    def m(n, k):
+        return  f(n) / (f(n/k)**k)
+
+    def p(n,k):
+        return m(n, k) / (k**n)
+
+    N, k = symbols('N k')
+    half = Float('0.5', 4)
+    z = log(p(n, k) / p(n, k + 1)).expand(force=True)
+    r = simplify(z.subs(n, N).n(4))
+    assert r == (half*k*log(k) - half*k*log(k + 1) + half*log(N) - half*log(k + 1) + Float(0.9189224, 4))
path: sympy/simplify/tests/test_simplify.py
text: @user1:
Personally I think this line is too long at over 100 characters. Black (with more SymPy-esque whitespace around operators) would suggest to reformat it like this:
```suggestion
    assert r == (
        half*k*log(k)
        - half*k*log(k + 1)
        + half*log(N)
        - half*log(k + 1)
        + Float(0.9189224, 4)
    )
```
Others may disagree that line length here isn't an issue, but it is considerably longer than anything surrounding it in this tests module.