Repository: sympy/sympy
Pull request: [integrals] fix substitution formula in meijerint
<!-- 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 #25786 

#### Brief description of what is fixed or changed

fix substitution formula in meijerint


#### 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:


* integrals
  * Fixed a bug with integrating over negative real numbers using meijerint algorithm. Formerly, `integrate(exp(-x**2), (x,-5,0), meijerg=True)` 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 -->
* integrals
  * Fixed a bug with integrating over negative real numbers using meijerint algorithm. Formerly, `integrate(exp(-x**2), (x,-5,0), meijerg=True)` incorrectly gave `-sqrt(pi)/2 * erf(5)` instead of `sqrt(pi)/2 * erf(5)`.
<!-- END RELEASE NOTES -->


Problem statement:
Wrong result for a simple integral
`sympy.integrate` returns a wrong result for a simple integral of a Gaussian function. The following code
```python
import sympy as sp
x = sp.Symbol('x')
f = sp.exp(-0.8*x**2)
sp.integrate(f, (x, -3, -0.2))
```
returns $−0.447288590994492\sqrt{\pi}$, which has the correct magnitude, but the wrong sign. This was observed with version 1.12 of `sympy`.

Interestingly, if the constant in the exponential is defined as a `Rational`:
```python
f = sp.exp(-sp.Rational(8,10)*x**2)
```
then the result is correct. If we use a different constant, for instance
```python
f = sp.exp(-0.4*x**2)
```
the result is also correct.

Hints from issue discussion:
In fact if we calculate the indefinite integral first and then substitute limits, sympy gives a correct result :

``` 
indefinite_integral = sp.integrate(f, x)
result = indefinite_integral.subs(x, -0.2) - indefinite_integral.subs(x, -3)
```

One of the reason for this disparity I believe is because sympy evaluates definite integrals symbolically, while indefinite integrals are more of a step by step approach, _**Although Ideally both should have same answer**_
The correct output is seen if both risch and meijerg are disabled:
```python
In [21]: integrate(f, (x, -3, -0.2), meijerg=False, risch=False)
Out[21]: 0.447288590994492⋅√π
```
It may be related to
    erf.eval  @ sympy/sympy/functions/special/error_functions.py
```
       # Try to pull out factors of -1
        if arg.could_extract_minus_sign():
            return -cls(-arg)
```

Consider the following snippet:
```
from sympy import Symbol, exp
from sympy.integrals.meijerint import meijerint_indefinite

x = Symbol('x', negative=True)
f = exp(-x**2)

print(meijerint_indefinite(f, x))
```
This will print "-sqrt(pi)*erf(x)/2" instead of sqrt(pi)*erf(x)/2. Is this intended?

Consider the following substitution t=a*x**b:
```
integral G(params, a*x**b) dx 
= 1/(a*b) integral x**(1-b) G(params, a*x**b) * ab x**(b-1) dx 
= 1/(a*b) integral (t/a)**(1/b-1) G(params, t) * dt 
= 1/(a**(1/b) b) integral G(params + 1/b-1, t) * dt 
= 1/(a**(1/b) b) G(1, params + 1/b, 0, t) 
= 1/(a**(1/b) b) G(1, params + 1/b, 0, a*x**b)
```
I think the code in [meijerint.py#L1729](https://github.com/sympy/sympy/blob/ffa18fad4a63d4a3a999b5bdb933bda2af1dce2f/sympy/integrals/meijerint.py#L1729) is based on this idea. But this is not quite right, e.g.
```
integral G(((),()),((0),()),x**2) dx
= integral exp(-x**2) dx
= sqrt(pi)*erf(x)/2
```
but `1/2 G(((1),()),((1/2),(0)),x^2)=sqrt(pi)*erf(x)/2 * sqrt(x**2)/x`
> Consider the following snippet:
> 
> ```
> from sympy import Symbol, exp
> from sympy.integrals.meijerint import meijerint_indefinite
> 
> x = Symbol('x', negative=True)
> f = exp(-x**2)
> 
> print(meijerint_indefinite(f, x))
> ```
> 
> This will print "-sqrt(pi)*erf(x)/2" instead of sqrt(pi)*erf(x)/2. Is this intended?
> 
> Consider the following substitution t=a*x**b:
> 
> ```
> integral G(params, a*x**b) dx 
> = 1/(a*b) integral x**(1-b) G(params, a*x**b) * ab x**(b-1) dx 
> = 1/(a*b) integral (t/a)**(1/b-1) G(params, t) * dt 
> = 1/(a**(1/b) b) integral G(params + 1/b-1, t) * dt 
> = 1/(a**(1/b) b) G(1, params + 1/b, 0, t) 
> = 1/(a**(1/b) b) G(1, params + 1/b, 0, a*x**b)
> ```
> 
> I think the code in [meijerint.py#L1729](https://github.com/sympy/sympy/blob/ffa18fad4a63d4a3a999b5bdb933bda2af1dce2f/sympy/integrals/meijerint.py#L1729) is based on this idea. But this is not quite right, e.g.
> 
> ```
> integral G(((),()),((0),()),x**2) dx
> = integral exp(-x**2) dx
> = sqrt(pi)*erf(x)/2
> ```
> 
> but `1/2 G(((1),()),((1/2),(0)),x^2)=sqrt(pi)*erf(x)/2 * sqrt(x**2)/x`


Regarding the problematic substitution I asked on [math.stackexchange](https://math.stackexchange.com/questions/4831298/antiderivative-of-meijer-g-function) and [mathoverflow](https://mathoverflow.net/questions/461606/antiderivative-of-meijer-g-function)
> I think the code in [meijerint.py#L1729](https://github.com/sympy/sympy/blob/ffa18fad4a63d4a3a999b5bdb933bda2af1dce2f/sympy/integrals/meijerint.py#L1729) is based on this idea. But this is not quite right, e.g.
> 
> ```
> integral G(((),()),((0),()),x**2) dx
> = integral exp(-x**2) dx
> = sqrt(pi)*erf(x)/2
> ```
> 
> but `1/2 G(((1),()),((1/2),(0)),x^2)=sqrt(pi)*erf(x)/2 * sqrt(x**2)/x`

Note the all to `powdenest` and the comment here:
https://github.com/sympy/sympy/blob/ffa18fad4a63d4a3a999b5bdb933bda2af1dce2f/sympy/integrals/meijerint.py#L1738-L1740
```python
In [5]: e = sqrt(x**2)

In [6]: e
Out[6]: 
   ____
  ╱  2 
╲╱  x  

In [7]: powdenest(e)
Out[7]: 
   ____
  ╱  2 
╲╱  x  

In [8]: powdenest(e,  polar=True)
Out[8]: x

In [9]: x = symbols('x', negative=True)

In [10]: e2 = sqrt(x**2)

In [11]: e2
Out[11]: -x

In [13]: powdenest(e2,  polar=True)
Out[13]: -x
```
Is that the problem here?

The premise of the meijerg code seems to be that `sqrt(x**2)` would remain unevaluated so that `powdenest` can simplify it later. When there is a need to take care of branches I think that meijerg uses `exp_polar` which is what the `polar=True` references.

I don't fully understand how meijerg uses `exp_polar` and how to do it correctly but the impression I have is that the code depends on not allowing the sort of assumptions based "simplification" that we see when the symbol has `negative=True`. It is possibly important to prevent the integration variable from having any assumptions set that would cause this.

Commit message: 25786 fix test

Patch to review:
diff --git a/.mailmap b/.mailmap
index 0a73b94bf2b5..7ceb7715d2e9 100644
--- a/.mailmap
+++ b/.mailmap
@@ -945,6 +945,7 @@ Matthew Wardrop <matthew.wardrop@airbnb.com>
 Matthias Bussonnier <bussonniermatthias@gmail.com>
 Matthias Geier <Matthias.Geier@gmail.com>
 Matthias Köppe <mkoeppe@math.ucdavis.edu> Matthias Koeppe <mkoeppe@math.ucdavis.edu>
+Matthias Liesenfeld <116307294+maliesen@users.noreply.github.com>
 Matthias Rettl <matthias.rettl@stud.unileoben.ac.at>
 Matthias Toews <mat.toews@googlemail.com>
 Mauro Garavello <mauro.garavello@unimib.it>
diff --git a/sympy/integrals/meijerint.py b/sympy/integrals/meijerint.py
index aa9dac6db166..5b518fdf2a16 100644
--- a/sympy/integrals/meijerint.py
+++ b/sympy/integrals/meijerint.py
@@ -1706,8 +1706,8 @@ def _meijerint_indefinite_1(f, x):
         c += s
 
         # we do a substitution t=a*x**b, get integrand fac*t**rho*g
-        fac_ = fac * C / (b*a**((1 + c)/b))
-        rho = (c + 1)/b - 1
+        fac_ = fac * C * x**(1 + c) / b
+        rho = (c + 1)/b
 
         # we now use t**rho*G(params, t) = G(params + rho, t)
         # [L, page 150, equation (4)]
@@ -1720,13 +1720,13 @@ def _meijerint_indefinite_1(f, x):
         t = _dummy('t', 'meijerint-indefinite', S.One)
 
         def tr(p):
-            return [a + rho + 1 for a in p]
+            return [a + rho for a in p]
         if any(b.is_integer and (b <= 0) == True for b in tr(g.bm)):
             r = -meijerg(
-                tr(g.an), tr(g.aother) + [1], tr(g.bm) + [0], tr(g.bother), t)
+                list(g.an), list(g.aother) + [1-rho], list(g.bm) + [-rho], list(g.bother), t)
         else:
             r = meijerg(
-                tr(g.an) + [1], tr(g.aother), tr(g.bm), tr(g.bother) + [0], t)
+                list(g.an) + [1-rho], list(g.aother), list(g.bm), list(g.bother) + [-rho], t)
         # The antiderivative is most often expected to be defined
         # in the neighborhood of  x = 0.
         if b.is_extended_nonnegative and not f.subs(x, 0).has(S.NaN, S.ComplexInfinity):
diff --git a/sympy/integrals/tests/test_meijerint.py b/sympy/integrals/tests/test_meijerint.py
index f23975e65c03..8cb7121f1959 100644
--- a/sympy/integrals/tests/test_meijerint.py
+++ b/sympy/integrals/tests/test_meijerint.py
@@ -755,10 +755,15 @@ def test_issue_6462():
 
 def test_indefinite_1_bug():
     assert integrate((b + t)**(-a), t, meijerg=True
-        ) == -b**(1 - a)*(1 + t/b)**(1 - a)/(a - 1)
+        ).equals(-b**(1 - a)*(1 + t/b)**(1 - a)/(a - 1))
 
 
 def test_pr_23583():
     # This result is wrong. Check whether new result is correct when this test fail.
     assert integrate(1/sqrt((x - I)**2-1), meijerg=True) == \
            Piecewise((acosh(x - I), Abs((x - I)**2) > 1), (-I*asin(x - I), True))
+
+
+# 25786
+def test_integrate_function_of_square_over_negatives():
+    assert integrate(exp(-x**2), (x,-5,0), meijerg=True) == sqrt(pi)/2 * erf(5)
diff --git a/sympy/physics/continuum_mechanics/tests/test_beam.py b/sympy/physics/continuum_mechanics/tests/test_beam.py
index 2c33fca5f9ae..3ee6b044b448 100644
--- a/sympy/physics/continuum_mechanics/tests/test_beam.py
+++ b/sympy/physics/continuum_mechanics/tests/test_beam.py
@@ -336,7 +336,7 @@ def test_variable_moment():
     assert b.slope().expand() == ((10*x*SingularityFunction(x, 0, 0)
         - 10*(x - 4)*SingularityFunction(x, 4, 0))/E).expand()
     assert b.deflection().expand() == ((5*x**2*SingularityFunction(x, 0, 0)
-        - 10*Piecewise((0, Abs(x)/4 < 1), (16*meijerg(((3, 1), ()), ((), (2, 0)), x/4), True))
+        - 10*Piecewise((0, Abs(x)/4 < 1), (x**2*meijerg(((-1, 1), ()), ((), (-2, 0)), x/4), True))
         + 40*SingularityFunction(x, 4, 1))/E).expand()
 
     b = Beam(4, E - x, I)


Code review comments:
diff_hunk: @@ -755,10 +755,15 @@ def test_issue_6462():
 
 def test_indefinite_1_bug():
     assert integrate((b + t)**(-a), t, meijerg=True
-        ) == -b**(1 - a)*(1 + t/b)**(1 - a)/(a - 1)
+        ).equals(-b**(1 - a)*(1 + t/b)**(1 - a)/(a - 1))
path: sympy/integrals/tests/test_meijerint.py
text: @user1:
It would be better to change this test answer to show what is now returned.

The old result was:
```python
In [4]: a, b, c = symbols("a b c")

In [5]: integrate((b + t)**(-a), t, meijerg=True)
Out[5]: 
               1 - a 
  1 - a ⎛    t⎞      
-b     ⋅⎜1 + ─⎟      
        ⎝    b⎠      
─────────────────────
        a - 1 
```
The new result is
```python
In [2]: integrate((b + t)**(-a), t, meijerg=True)
Out[2]: 
          1 - a 
   ⎛    t⎞      
-b⋅⎜1 + ─⎟      
   ⎝    b⎠      
────────────────
      a    a    
   a⋅b  - b 
```