Fixes invert on a symbol returns 0

[1e9abhi1e10]

References to other Issues or PRs

Fixes #24394

Brief description of what is fixed or changed

Other comments

Release Notes

NO ENTRY

[sympy-bot]

✅

Hi, I am the SymPy bot (v169). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

• No release notes entry will be added for this pull request.

Click here to see the pull request description that was parsed.

<!-- 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 #24394 


#### 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.

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 -->

[oscarbenjamin]

The slow tests failed showing that something somewhere depends on this. That needs further investigation.

[asmeurer]

This is the code in question:

            if invert:
                h_lc = Poly(h.as_poly(DE.t).LC(), DE.t, field=True, expand=False)
                inv, coeffs = h_lc.as_poly(z, field=True).invert(s), [S.One]

                for coeff in h.coeffs()[1:]:
                    L = reduced(inv*coeff.as_poly(inv.gens), [s])[1]
                    coeffs.append(L.as_expr())

                h = Poly(dict(list(zip(h.monoms(), coeffs))), DE.t)

In this example, h.coeffs() only has one term, so the inv value is not actually used. I'll check the book to see if the code here looks correct from the pseudocode.

[asmeurer]

The code comes from 8079e44 and isn't strictly part of the pseudocode in the book. I don't remember exactly why this was added, but it seems likely enough to me that the logic there is simply incorrect.

[asmeurer]

I guess the book does state, "we can make $pp_i(R_m)(\alpha, t)$ monic in order to simplify the answer", but it doesn't include that in the pseudocode (this is all on page 153). I wouldn't be surprised if there is a more direct way of achieving this than those 6 lines of code.

[asmeurer]

This is in residue_reduce in sympy/integrals/risch.py btw

[asmeurer]

The example from the book, risch_integrate((2*log(x)**2 - log(x) - x**2)/(log(x)**3 - x**2*log(x)), x), shows why it is "needed". Right now it gives

>>> pprint(risch_integrate((2*log(x)**2 - log(x) - x**2)/(log(x)**3 - x**2*log(x)), x))
                                       ⌠
  log(-x + log(x))   log(x + log(x))   ⎮   1
- ──────────────── + ─────────────── + ⎮ ────── dx
         2                  2          ⎮ log(x)
                                       ⌡

But if you change the default to invert=False, you instead get

>>> pprint(risch_integrate((2*log(x)**2 - log(x) - x**2)/(log(x)**3 - x**2*log(x)), x))
   ⎛          2                              ⎞      ⎛        2                              ⎞
   ⎜   3   3⋅x    x   ⎛   2   3⋅x   1⎞       ⎟      ⎜ 3   3⋅x    x   ⎛   2   3⋅x   1⎞       ⎟   ⌠
log⎜- x  + ──── - ─ + ⎜- x  + ─── - ─⎟⋅log(x)⎟   log⎜x  + ──── + ─ + ⎜- x  - ─── - ─⎟⋅log(x)⎟   ⎮    4      2             2
   ⎝        2     2   ⎝        2    2⎠       ⎠      ⎝      2     2   ⎝        2    2⎠       ⎠   ⎮ 4⋅x  - 6⋅x ⋅log(x) - 5⋅x  + 3⋅log(x) + 1
────────────────────────────────────────────── - ──────────────────────────────────────────── + ⎮ ──────────────────────────────────────── dx
                      2                                               2                         ⎮       4             2
                                                                                                ⎮    4⋅x ⋅log(x) - 5⋅x ⋅log(x) + log(x)
                                                                                                ⌡

which doesn't even split out the non-elementary integral from the rational integral correctly (differentiating confirms both answers are correct).

Apparently this integral isn't tested in test_risch.py. When I changed the default, the only non-slow test in sympy/integrals that failed was a direct test for the residue_reduce function. So it certainly seems to be under tested, suggesting the logic there could indeed be wrong.

[1e9abhi1e10]

The example from the book, risch_integrate((2*log(x)**2 - log(x) - x**2)/(log(x)**3 - x**2*log(x)), x), shows why it is "needed". Right now it gives

>>> pprint(risch_integrate((2*log(x)**2 - log(x) - x**2)/(log(x)**3 - x**2*log(x)), x))
                                       ⌠
  log(-x + log(x))   log(x + log(x))   ⎮   1
- ──────────────── + ─────────────── + ⎮ ────── dx
         2                  2          ⎮ log(x)
                                       ⌡

But if you change the default to invert=False, you instead get

>>> pprint(risch_integrate((2*log(x)**2 - log(x) - x**2)/(log(x)**3 - x**2*log(x)), x))
   ⎛          2                              ⎞      ⎛        2                              ⎞
   ⎜   3   3⋅x    x   ⎛   2   3⋅x   1⎞       ⎟      ⎜ 3   3⋅x    x   ⎛   2   3⋅x   1⎞       ⎟   ⌠
log⎜- x  + ──── - ─ + ⎜- x  + ─── - ─⎟⋅log(x)⎟   log⎜x  + ──── + ─ + ⎜- x  - ─── - ─⎟⋅log(x)⎟   ⎮    4      2             2
   ⎝        2     2   ⎝        2    2⎠       ⎠      ⎝      2     2   ⎝        2    2⎠       ⎠   ⎮ 4⋅x  - 6⋅x ⋅log(x) - 5⋅x  + 3⋅log(x) + 1
────────────────────────────────────────────── - ──────────────────────────────────────────── + ⎮ ──────────────────────────────────────── dx
                      2                                               2                         ⎮       4             2
                                                                                                ⎮    4⋅x ⋅log(x) - 5⋅x ⋅log(x) + log(x)
                                                                                                ⌡

which doesn't even split out the non-elementary integral from the rational integral correctly (differentiating confirms both answers are correct).

Apparently this integral isn't tested in test_risch.py. When I changed the default, the only non-slow test in sympy/integrals that failed was a direct test for the residue_reduce function. So it certainly seems to be under tested, suggesting the logic there could indeed be wrong.

So can we add tests here of this example

[oscarbenjamin]

@asmeurer I think that unless you say precisely how the risch code could be fixed this PR will stall.

[jksuom]

This looks like one of those issues caused by the gcd of polynomials being different over field and non-field coefficients. This is a possible fix:

--- a/sympy/integrals/risch.py
+++ b/sympy/integrals/risch.py
@@ -1312,7 +1312,7 @@ def residue_reduce(a, d, DE, z=None, invert=True):
     for i in R:
         R_map[i.degree()] = i
 
-    r = Poly(r, z)
+    r = Poly(r, z, field=True)
     Np, Sp = splitfactor_sqf(r, DE, coefficientD=True, z=z)
 
     for s, i in Sp:

[jksuom]

There are assertion failures in test_residue_reduce, but that is to be expected. When a gcd of polynomials is computed with field coefficients, (non-zero) coefficients are treated as units and not taken into account. Hence the polynomials output by splitfactor_sqf differ by constant factors in general. Those polynomials are also correct (gcd's are not unique), and the tests can be fixed by changing their results.