sage icon indicating copy to clipboard operation
sage copied to clipboard

Use the faster libgiac interface for "giac" integration

Open orlitzky opened this issue 1 year ago • 6 comments

While working on https://github.com/sagemath/sage/issues/38668, I noticed that all of our examples that use giac for integration do so with algorithm='giac'. This uses the pexpect interface that is much slower (and simply sloppier) than the library interface that would be used with algorithm='libgiac'.

To fix this, I could just change those examples to algorithm='libgiac', but since the libgiac interface is better in every way, that seems kind of rude to me. Most users aren't going to know what lib-anything is, and if they just want to use giac, they're going to try algorithm='giac'. So instead this PR makes algorithm='giac' do the right thing, and use the better interface. I've left the algorithm='libgiac' alone for now to avoid breaking any code.

Afterwards I have removed the pexpect giac integrator entirely, because there's no reason to use it, and no other code in sagelib uses it.

orlitzky avatar Sep 20 '24 13:09 orlitzky

Documentation preview for this PR (built with commit 09572b108b0216703d7a77055ad629383bd8d37b; changes) is ready! :tada: This preview will update shortly after each push to this PR.

github-actions[bot] avatar Sep 20 '24 14:09 github-actions[bot]

AFAICS the CI failure is an unrelated timeout:

2024-09-20T15:19:09.1450225Z ----------------------------------------------------------------------
2024-09-20T15:19:09.1451853Z sage -t --long --random-seed=286735480429121101562228604801325644303 src/sage/data_structures/bounded_integer_sequences.pyx  # Timed out
2024-09-20T15:19:09.1453353Z ----------------------------------------------------------------------

orlitzky avatar Sep 20 '24 17:09 orlitzky

This looks like a good idea, thanks.

mkoeppe avatar Sep 29 '24 01:09 mkoeppe

I didn't see the label change, thanks for the review.

orlitzky avatar Oct 02 '24 20:10 orlitzky

I'm getting

sage -t --long --warn-long 45.2 --random-seed=123 src/sage/symbolic/integration/integral.py
**********************************************************************
File "src/sage/symbolic/integration/integral.py", line 64, in sage.symbolic.integration.integral.IndefiniteIntegral.__init__
Failed example:
    integrate(Ex, x, algorithm='giac')  # long time
Expected:
    4*(-2*x^(1/3) + 1)^(3/4) + 6*arctan((-2*x^(1/3) + 1)^(1/4)) - 3*log((-2*x^(1/3) + 1)^(1/4) + 1) + 3*log(abs((-2*x^(1/3) + 1)^(1/4) - 1))
Got:
    4*(-2*x^(1/3) + 1)^(3/4) + 6*arctan((-2*x^(1/3) + 1)^(1/4)) - 3*log(abs((-2*x^(1/3) + 1)^(1/4) + 1)) + 3*log(abs((-2*x^(1/3) + 1)^(1/4) - 1))
**********************************************************************
File "src/sage/symbolic/integration/integral.py", line 209, in sage.symbolic.integration.integral.DefiniteIntegral.__init__
Failed example:
    integral(1/max_symbolic(x, 1)**2, x, 0, oo, algorithm='giac')
Exception raised:
    Traceback (most recent call last):
      File "/home/release/Sage/src/sage/doctest/forker.py", line 714, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/release/Sage/src/sage/doctest/forker.py", line 1144, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.symbolic.integration.integral.DefiniteIntegral.__init__[4]>", line 1, in <module>
        integral(Integer(1)/max_symbolic(x, Integer(1))**Integer(2), x, Integer(0), oo, algorithm='giac')
      File "/home/release/Sage/src/sage/misc/functional.py", line 785, in integral
        return x.integral(*args, **kwds)
               ^^^^^^^^^^^^^^^^^^^^^^^^^
      File "sage/symbolic/expression.pyx", line 13222, in sage.symbolic.expression.Expression.integral
        return integral(self, *args, **kwds)
      File "/home/release/Sage/src/sage/symbolic/integration/integral.py", line 1073, in integrate
        return integrator(expression, v, a, b)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      File "/home/release/Sage/src/sage/symbolic/integration/external.py", line 259, in libgiac_integrator
        result = libgiac.integrate(Pygen(expression), v, a, b)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      File "sage/libs/giac/giac.pyx", line 1934, in sage.libs.giac.giac.GiacFunction.__call__
        return Pygen.__call__(self, *args)
      File "sage/libs/giac/giac.pyx", line 1109, in sage.libs.giac.giac.Pygen.__call__
        raise
      File "sage/libs/giac/giac.pyx", line 1099, in sage.libs.giac.giac.Pygen.__call__
        result = self.gptr[0](right.gptr[0], context_ptr)
    RuntimeError: Unable to check sign: 1>Infinity
**********************************************************************
File "src/sage/symbolic/integration/integral.py", line 697, in sage.symbolic.integration.integral.integrate
Failed example:
    integrate(abs(cos(x)), x, 0, 2*pi, algorithm='giac')
Expected:
    4
Got:
    Warning, integration of abs or sign assumes constant sign by intervals (correct if the argument is real):
    Check [abs(cos(sageVARx))]
    Discontinuities at zeroes of cos(sageVARx) were not checked
    4
**********************************************************************

vbraun avatar Oct 04 '24 22:10 vbraun

Expected: 4*(-2x^(1/3) + 1)^(3/4) + 6arctan((-2x^(1/3) + 1)^(1/4)) - 3log((-2x^(1/3) + 1)^(1/4) + 1) + 3log(abs((-2x^(1/3) + 1)^(1/4) - 1)) Got: 4(-2x^(1/3) + 1)^(3/4) + 6arctan((-2x^(1/3) + 1)^(1/4)) - 3log(abs((-2x^(1/3) + 1)^(1/4) + 1)) + 3log(abs((-2*x^(1/3) + 1)^(1/4) - 1))

These are the same anwer. I updated the expected output.

Expected: 4 Got: Warning, integration of abs or sign assumes constant sign by intervals (correct if the argument is real): Check [abs(cos(sageVARx))] Discontinuities at zeroes of cos(sageVARx) were not checked 4

Nothing too scary here. This doctest already looks for the warning output on the previous line, where integrate is used without an algorithm argument. I tweaked the description a bit and added some ....

File "src/sage/symbolic/integration/integral.py", line 209, in sage.symbolic.integration.integral.DefiniteIntegral.init Failed example: integral(1/max_symbolic(x, 1)**2, x, 0, oo, algorithm='giac') ... RuntimeError: Unable to check sign: 1>Infinity

This one on the other hand is a treat. It's reproducible using libgiac_integrator, but not through integrate(), which (if maxima fails) calls libgiac_integrator with exactly the same arguments. My guess is hidden state within sage.libs.giac.

In any case, I was able to work around it. Don't ask me what the difference is, but giac likes +infinity better than Infinity. You get the latter passing Sage's infinity directly to libgiac, and the former by calling _giac_() on Sage's infinity. So I added a few lines to the integrator to do that when the endpoints are plus/minus infinity.

orlitzky avatar Oct 06 '24 11:10 orlitzky

In any case, I was able to work around it. Don't ask me what the difference is, but giac likes +infinity better than Infinity. You get the latter passing Sage's infinity directly to libgiac, and the former by calling _giac_() on Sage's infinity. So I added a few lines to the integrator to do that when the endpoints are plus/minus infinity.

I went back and reworked this. Infinity in giac is unsigned infinity, and the conversion was being done incorrectly (from signed in sage to unsigned in giac). After adding another case to the Pygen constructor, everything is fine. No hacks needed in the integration routine.

orlitzky avatar Nov 13 '24 02:11 orlitzky

If someone has a moment, I lost my reviewer.

This adds the final set of # needs tags to make giac optional (using meson, for now), but at the same time it switches integrate(..., algorithm="giac") to use the faster libgiac interface. That should have been trivial but there was a bug in sage.libs.giac's handling of infinity that complicated things. That is now fixed too.

orlitzky avatar Dec 04 '24 21:12 orlitzky

Thanks, I might finally be able to pull sage.libs.giac out into a separate package now (if the GCC 15 bugs ever end).

orlitzky avatar Dec 05 '24 15:12 orlitzky