polarify meijerg z param prior to hyperexpand #26551
Conversation
|
✅ Hi, I am the SymPy bot. I'm here to help you write a release notes entry. Please read the guide on how to write release notes. Your release notes are in good order. Here is what the release notes will look like:
This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13. Click here to see the pull request description that was parsed. |
|
I have never really understood exactly what Can you explain what the purpose of polarify is especially in relation to meijerg and the integration routines and how the change here improves things (without making other things worse)? |
|
This doesn't address this case: >>> meijerg(((1, 2), ()), ((S(1)/2,), (0,)), 0)
meijerg(((1, 2), ()), ((1/2,), (0,)), 0)
>>> hyperexpand(_) # prior to this pull this results in -4*sqrt(pi)
0
>>> meijerg(((1, 2), ()), ((S(1)/2,), (0,)), x)
meijerg(((1, 2), ()), ((1/2,), (0,)), x)
>>> hyperexpand(_)
-sqrt(pi)*sqrt(x)*(4*(x/2 + 1/2)/sqrt(x + 1) + 2*asinh(sqrt(x))/sqrt(x))
>>> _.subs(x,0)
nanbut this is the same behavior as mathematica: (If it was worthwhile to address this could be handled similar to the "heurisch wrapper" reintegration which mathematica also doesn't do) |
I'm not quite sure what you think should be addressed. The PR gives: In [1]: e = meijerg(((1, 2), ()), ((S(1)/2,), (0,)), x)
In [2]: hyperexpand(e.subs(x, 0))
Out[2]: 0
In [3]: hyperexpand(e).subs(x, 0)
Out[3]: nan
In [4]: hyperexpand(e).limit(x, 0)
Out[4]: 0
In [5]: e.limit(x, 0)
Out[5]: 0
In [6]: hyperexpand(e)
Out[6]:
⎛ ⎛x 1⎞ ⎞
⎜4⋅⎜─ + ─⎟ ⎟
⎜ ⎝2 2⎠ 2⋅asinh(√x)⎟
-√π⋅√x⋅⎜───────── + ───────────⎟
⎜ _______ √x ⎟
⎝╲╱ x + 1 ⎠
In [7]: hyperexpand(e).expand()
Out[7]:
3/2
2⋅√π⋅x 2⋅√π⋅√x
- ───────── - ───────── - 2⋅√π⋅asinh(√x)
_______ _______
╲╱ x + 1 ╲╱ x + 1
In [8]: hyperexpand(e).expand().subs(x, 0)
Out[8]: 0The This is unfortunate (not caused by this PR): In [9]: e.series(x, n=3)
Out[9]:
╭─╮1, 2 ⎛0, 2 │ ⎞ ╭─╮1, 2 ⎛-1, 2 │ ⎞ ╭─╮1, 2 ⎛1, 2 │ ⎞ ⎛ 3⎞
zoo⋅│╶┐ ⎜ │ 0⎟ + zoo⋅│╶┐ ⎜ │ 0⎟ + │╶┐ ⎜ │ 0⎟ + O⎝x ⎠
╰─╯2, 2 ⎝1/2 0 │ ⎠ ╰─╯2, 2 ⎝ 1/2 0 │ ⎠ ╰─╯2, 2 ⎝1/2 0 │ ⎠
In [10]: hyperexpand(e).series(x, n=3)
Out[10]:
3/2 5/2
2⋅√π⋅x √π⋅x ⎛ 3⎞
-4⋅√π⋅√x - ───────── + ─────── + O⎝x ⎠
3 10 |
I don't have much more understanding beyond what's in the script + comments added in #26551 (never taken even the "complex variables" course). I think I'd have a chance of understanding why this works for the z=0 case (and |
I think if we wanted to be fully pedantically correct the result should be a Piecewise |
I haven't followed through the code but I wonder if it should be valid to multiple the In [6]: e1, e2, e3 = sqrt(x), x, 1/sqrt(x)
In [7]: e1*(e2 + e3)
Out[7]:
⎛ 1 ⎞
√x⋅⎜x + ──⎟
⎝ √x⎠
In [8]: e1*e2 + e1*e3
Out[8]:
3/2
x + 1 |
|
Benchmark results from GitHub Actions Lower numbers are good, higher numbers are bad. A ratio less than 1 Significantly changed benchmark results (PR vs master) Significantly changed benchmark results (master vs previous release) | Change | Before [2487dbb5] | After [b37bfdd3] | Ratio | Benchmark (Parameter) |
|----------|----------------------|---------------------|---------|----------------------------------------------------------------------|
| - | 69.0±1ms | 43.6±0.2ms | 0.63 | integrate.TimeIntegrationRisch02.time_doit(10) |
| - | 68.1±0.6ms | 43.0±0.2ms | 0.63 | integrate.TimeIntegrationRisch02.time_doit_risch(10) |
| + | 18.3±0.5μs | 30.9±0.2μs | 1.69 | integrate.TimeIntegrationRisch03.time_doit(1) |
| - | 5.38±0.03ms | 2.86±0.01ms | 0.53 | logic.LogicSuite.time_load_file |
| - | 74.4±0.3ms | 28.6±0.2ms | 0.38 | polys.TimeGCD_GaussInt.time_op(1, 'dense') |
| - | 26.0±0.1ms | 16.6±0.04ms | 0.64 | polys.TimeGCD_GaussInt.time_op(1, 'expr') |
| - | 74.9±0.3ms | 28.7±0.2ms | 0.38 | polys.TimeGCD_GaussInt.time_op(1, 'sparse') |
| - | 260±2ms | 124±0.6ms | 0.48 | polys.TimeGCD_GaussInt.time_op(2, 'dense') |
| - | 264±3ms | 124±0.7ms | 0.47 | polys.TimeGCD_GaussInt.time_op(2, 'sparse') |
| - | 658±6ms | 371±2ms | 0.56 | polys.TimeGCD_GaussInt.time_op(3, 'dense') |
| - | 659±5ms | 371±1ms | 0.56 | polys.TimeGCD_GaussInt.time_op(3, 'sparse') |
| - | 499±2μs | 286±3μs | 0.57 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(1, 'dense') |
| - | 1.79±0.01ms | 1.04±0.01ms | 0.59 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(2, 'dense') |
| - | 5.85±0.04ms | 3.11±0.04ms | 0.53 | polys.TimeGCD_LinearDenseQuadraticGCD.time_op(3, 'dense') |
| - | 452±2μs | 229±2μs | 0.51 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(1, 'dense') |
| - | 1.49±0.01ms | 685±3μs | 0.46 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(2, 'dense') |
| - | 5.01±0.04ms | 1.64±0.02ms | 0.33 | polys.TimeGCD_QuadraticNonMonicGCD.time_op(3, 'dense') |
| - | 378±1μs | 207±1μs | 0.55 | polys.TimeGCD_SparseGCDHighDegree.time_op(1, 'dense') |
| - | 2.51±0ms | 1.24±0.01ms | 0.5 | polys.TimeGCD_SparseGCDHighDegree.time_op(3, 'dense') |
| - | 9.94±0.1ms | 4.31±0.02ms | 0.43 | polys.TimeGCD_SparseGCDHighDegree.time_op(5, 'dense') |
| - | 365±2μs | 169±0.5μs | 0.46 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(1, 'dense') |
| - | 2.53±0.01ms | 902±3μs | 0.36 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(3, 'dense') |
| - | 9.76±0.03ms | 2.63±0.02ms | 0.27 | polys.TimeGCD_SparseNonMonicQuadratic.time_op(5, 'dense') |
| - | 1.04±0.01ms | 421±2μs | 0.41 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'dense') |
| - | 1.76±0.01ms | 507±2μs | 0.29 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| - | 6.01±0.04ms | 1.77±0.01ms | 0.3 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'dense') |
| - | 8.62±0.07ms | 1.51±0ms | 0.17 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse') |
| - | 287±1μs | 64.6±0.09μs | 0.23 | polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse') |
| - | 3.57±0.04ms | 395±2μs | 0.11 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'dense') |
| - | 4.06±0.01ms | 280±1μs | 0.07 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse') |
| - | 7.24±0.08ms | 1.25±0ms | 0.17 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'dense') |
| - | 8.91±0.03ms | 834±3μs | 0.09 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse') |
| - | 5.03±0.02ms | 3.02±0.04ms | 0.6 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse') |
| - | 12.3±0.07ms | 6.51±0.06ms | 0.53 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'dense') |
| - | 22.5±0.1ms | 9.08±0.01ms | 0.4 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| - | 5.21±0.04ms | 875±10μs | 0.17 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse') |
| - | 12.5±0.05ms | 7.05±0.03ms | 0.57 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse') |
| - | 104±0.7ms | 25.7±0.02ms | 0.25 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'dense') |
| - | 169±1ms | 54.2±0.2ms | 0.32 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse') |
| - | 173±1μs | 111±0.5μs | 0.64 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'dense') |
| - | 362±0.8μs | 219±2μs | 0.6 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(1, 'sparse') |
| - | 4.40±0.01ms | 843±8μs | 0.19 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'dense') |
| - | 5.46±0.07ms | 385±4μs | 0.07 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse') |
| - | 20.4±0.2ms | 2.81±0.03ms | 0.14 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'dense') |
| - | 23.1±0.3ms | 639±3μs | 0.03 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse') |
| - | 481±0.8μs | 137±1μs | 0.29 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse') |
| - | 4.85±0.01ms | 619±3μs | 0.13 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'dense') |
| - | 5.41±0.04ms | 140±1μs | 0.03 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse') |
| - | 13.6±0.2ms | 1.28±0ms | 0.09 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'dense') |
| - | 14.4±0.07ms | 143±0.8μs | 0.01 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse') |
| - | 133±2μs | 75.2±0.5μs | 0.56 | solve.TimeMatrixOperations.time_rref(3, 0) |
| - | 250±1μs | 88.8±0.8μs | 0.36 | solve.TimeMatrixOperations.time_rref(4, 0) |
| - | 24.4±0.1ms | 10.4±0.5ms | 0.43 | solve.TimeSolveLinSys189x49.time_solve_lin_sys |
| - | 28.3±0.2ms | 15.5±0.1ms | 0.55 | solve.TimeSparseSystem.time_linsolve_Aaug(20) |
| - | 55.0±0.2ms | 24.8±0.1ms | 0.45 | solve.TimeSparseSystem.time_linsolve_Aaug(30) |
| - | 28.3±0.2ms | 15.2±0.09ms | 0.54 | solve.TimeSparseSystem.time_linsolve_Ab(20) |
| - | 54.8±0.08ms | 24.7±0.2ms | 0.45 | solve.TimeSparseSystem.time_linsolve_Ab(30) |
Full benchmark results can be found as artifacts in GitHub Actions |
ahh, missed this part. |
|
The blog of the original meijerint author @ness01 discusses the design of polarify/unpolarify in multiple posts e.g. https://nessgrh.wordpress.com/2011/08/06/branching-once-more/ Regarding
I would suspect that this change never affects the working of meijerint because for example in the still failing test case (the only one that actually comes from meijerint!): z has already been "polarified". I suspect meijerint uses G functions with already "polarified" arguments because of the hardcoded use of polar_lift in create_lookup_table although on this last point I'm not too sure. One quote from the blog: "we let S be the riemann surface of the logarithm, we understand G to be defined thereon " - this doesn't quite explain why mpmath can correctly evaluate G functions with rectangular complex z. Also the Riemann surface of the logarithm is not defined at 0. Also note that meijerint in the indefinite integral case will never pass a G function with a purely numeric z. However, the above comments don't address why this change seems to fix the first test case added in the commit: and also why simply adding One suspicious thing I did find is the mpmath documentation: https://mpmath.org/doc/current/functions/hypergeometric.html#meijer-g-function "The default series is chosen based on the degree and in order to be consistent with Mathematica’s. This definition of the Meijer G-function has a discontinuity at |z| = 1 for some orders, which can be avoided by explicitly specifying a series." Most of the added testcases have |z| = 1 |
|
Note also the way that mpmath's meijerg function is called by SymPy: sympy/sympy/functions/special/hyper.py Lines 694 to 724 in 37bb755 |
|
I've just never put it all together to understand how the branching is supposed to work but you can see that it goes wrong sometimes e.g. taking the integral mentioned in this blog post: The behaviour today is: In [25]: i = 1/(x*sqrt(1-x**2))
In [26]: i
Out[26]:
1
─────────────
________
╱ 2
x⋅╲╱ 1 - x
In [27]: i.integrate(x)
Out[27]:
⎧ ⎛1⎞ 1
⎪-acosh⎜─⎟ for ──── > 1
⎪ ⎝x⎠ │ 2│
⎪ │x │
⎨
⎪ ⎛1⎞
⎪ⅈ⋅asin⎜─⎟ otherwise
⎪ ⎝x⎠
⎩
In [28]: i.integrate(x).diff(x)
Out[28]:
⎧ 1 1
⎪─────────────────────────── for ──── > 1
⎪ ________ _______ │ 2│
⎪ 2 ╱ 1 ╱ 1 │x │
⎪x ⋅ ╱ -1 + ─ ⋅ ╱ 1 + ─
⎪ ╲╱ x ╲╱ x
⎪
⎨ -ⅈ
⎪ ──────────────── otherwise
⎪ ________
⎪ 2 ╱ 1
⎪ x ⋅ ╱ 1 - ──
⎪ ╱ 2
⎪ ╲╱ x
⎩
In [29]: i.subs(x, 2*I).n()
Out[29]: -0.223606797749979⋅ⅈ
In [31]: i.integrate(x).diff(x).subs(x, 2*I).n()
Out[31]: 0.223606797749979⋅ⅈ |
|
I don't understand enough about this yet to be confident that this is the right fix. I read about half the blog posts which helped (they are well written) but I am not there yet... |
|
No rush, thanks for helping |
References to other Issues or PRs
partially addresses #26525 (but not the first formula mentioned in the issue title)
Brief description of what is fixed or changed
Other comments
Release Notes