Prototype risch_integrate() function ready for testing

[Aaron S. Meurer]

(copied from issue 2010)

I have ready in my integration3 branch a prototype risch_integrate() function, that is a user-level function for the full Risch Algorithm I have been implementing this summer. Pull from http://github.com/asmeurer/sympy/tree/integration3
.

This is NOT ready to go in. It is a prototype function that I am making available so people can try out the new algorithm and hopefully help me to find the bugs in it. Please pass it your favorite non-elementary integrals and see if it can determine that they are not elementary. If you try to pass it a very crazy function at random, the chances are pretty high that it will not be elementary. So a better way to test it is to come up with a crazy function, then differentiate it. Then pass the derivative and see if it can give you your original function back. Note that it will probably not look exactly the same as your original function, and may differ by a constant. You should verify by differentiating the result you get and calling cancel() (or simplify(), but usually cancel() is enough) on the difference.

So you can review the code too, if you like, but just know that things are not stable yet, and this isn't strictly a branch for review.

So far, this function only supports exponentials and logarithms.
Support for trigonometric functions is planned. Algebraic functions are
not supported. If the function returns an unevaluated Integral, it means
that it has proven the integral to be non-elementary. Note that several
cases are still not implemented, so you may get NotImplementedError
instead. Eventually, these will all be eliminated, and the only
NotImplementedError you should see from this function is
NotImplementedError("Algebraic extensions are not supported.")

This function has not been integrated in any way with the already
existing integrate() yet, and you can use it to compare.

Examples:

In [1]: risch_integrate(exp(x**2), x)
Out[1]:
⌠
⎮ ⎛ 2⎞
⎮ ⎝x ⎠
⎮ ℯ dx
⌡

In [2]: risch_integrate(x**100*exp(x), x).diff(x)
Out[2]:
100 x
x ⋅ℯ

In [3]: %timeit risch_integrate(x**100*exp(x), x).diff(x)
1 loops, best of 3: 270 ms per loop

In [4]: integrate(x**100*exp(x), x)
... hangs ...

In [5]: risch_integrate(x/log(x), x)
Out[5]:
⌠
⎮ x
⎮ ────── dx
⎮ log(x)
⌡

In [6]: risch_integrate(log(x)**10, x).diff(x)
Out[6]:
10
log (x)

In [7]: integrate(log(x)**10, x).diff(x)
Out[7]:
10
log (x)

In [8]: %timeit risch_integrate(log(x)**10, x).diff(x)
10 loops, best of 3: 159 ms per loop

In [9]: %timeit integrate(log(x)**10, x).diff(x)
1 loops, best of 3: 2.35 s per loop

Be warned that things are still very buggy and you should always verify
results by differentiating. Usually, cancel(diff(result, x) - result)
should be enough. This should go to 0.

So please, please, PLEASE, try out this function and report any bugs that you find. It is not necessary to report NotImplementedError bugs, because I already know about those (I put them in there), and as I mentioned above, they are all planned to disappear. Also, I am continually updating my branch with fixes, so you should do a "git pull" and try again before you report anything.

Also, I am aware that there are test failures. This is because I had to hack exp._eval_subs() to only do exact substitution (no algebraic substitution). It's just a quick hack workaround, and I should eventually get a real fix.

Finally, I'm thinking there needs to be a way to differentiate between an unevaluated Integral because the integrator failed and an unevaluated Integral because it has proven the integral to be non-elementary. Any ideas?

Aaron Meurer

[Ondrej Certik]

Hi Aaron!

On Thu, Aug 5, 2010 at 2:01 PM, Aaron S. Meurer <asme...@gmail.com> wrote:
> (copied from issue 2010)
>
> I have ready in my integration3 branch a prototype risch_integrate() function, that is a user-level function for the full Risch Algorithm I have been implementing this summer.  Pull from http://github.com/asmeurer/sympy/tree/integration3

This is just excellent!

I would like to invite everyone to try this. (Read Aaron's email above
for things to try and not to try yet.) So here are some tougher cases:

In [1]: risch_integrate(1/(x**8+1), x)
[hangs]

In [4]: cancel(risch_integrate(1/(x**9+1), x).diff(x))
Out[4]:
d ⎛ ⎛ 6 3 ⎞⎞ 3 d ⎛
1 + 3⋅──⎝RootSum⎝531441⋅t + 729⋅t + 1, Λ(t, t⋅log(x + 9⋅t))⎠⎠ + 3⋅x ⋅──⎝Root
dx dx
──────────────────────────────────────────────────────────────────────────────
3
3 + 3⋅x

⎛ 6 3 ⎞⎞
Sum⎝531441⋅t + 729⋅t + 1, Λ(t, t⋅log(x + 9⋅t))⎠⎠

──────────────────────────────────────────────────

It should be equivalent but it's a bit messy. Maybe we need to
implement some symbolic manipulation of RootSum().

In [18]: risch_integrate(sqrt(1+exp(x)), x)
NotImplementedError: Couldn't find an elementary transcendental
extension for (1 + exp(x))**(1/2). Try using a manual extension with
the extension flag.

[18] is probably not supported yet.

Otherwise it works great. I wasn't able to make it break.

Ondrej

[Aaron S. Meurer]

So a few words. First, if you just pass it a strictly rational function, it is nothing new. It will return the exact same result as integrate() because it uses the exact same function, which is the already existing ratint() (however, there is a fix in my branch for rational functions with symbolic coefficients that fail in master).

To really test this, you need to pass it a function that has exp and/or log in it. For example, here are some functions that you might try that better test the algorithm:

In [2]: risch_integrate(1/(exp(x)**9 + 1), x)
Out[2]:
⎛ 9⋅x⎞
log⎝1 + ℯ ⎠
x - ─────────────
9

(your example with x replaced with exp(x))

In [18]: risch_integrate(diff(exp(x)*log(x)/(x + 1), x), x)
Out[18]:
x
ℯ ⋅log(x)
─────────
1 + x

This seems to be a simple example, but observe that our current integrate() cannot handle it:

In [19]: integrate(diff(exp(x)*log(x)/(x + 1), x), x)
Out[19]:
⌠
⎮ ⎛ x x x ⎞
⎮ ⎜ℯ ⋅log(x) ℯ ⋅log(x) ℯ ⎟
⎮ ⎜───────── - ───────── + ─────────⎟ dx
⎮ ⎜ 1 + x 2 x⋅(1 + x)⎟
⎮ ⎝ (1 + x) ⎠
⌡

Also, it's fun to pass it functions that you know do not have elementary anti-derivatives, to see if it can verify that fact:

In [20]: risch_integrate(exp(x**2), x)
Out[20]:
⌠
⎮ ⎛ 2⎞
⎮ ⎝x ⎠
⎮ ℯ dx
⌡

In [21]: risch_integrate(1/log(x), x)
Out[21]:
⌠
⎮ 1
⎮ ────── dx
⎮ log(x)
⌡

Also, remember what I said about integrating random functions. If you try to integrate a random function, the chances are pretty good that it will not be elementary, and even if it is, the result could be quite complicated and it could take a long time to compute. Much better is to come up with an expression that is as complex or not complex as you like, then differentiate it and see if risch_integrate() can give you the original thing back again.

On Aug 18, 2010, at 9:05 PM, Ondrej Certik wrote:

> Hi Aaron!
>
> On Thu, Aug 5, 2010 at 2:01 PM, Aaron S. Meurer <asme...@gmail.com> wrote:
>> (copied from issue 2010)
>>
>> I have ready in my integration3 branch a prototype risch_integrate() function, that is a user-level function for the full Risch Algorithm I have been implementing this summer. Pull from http://github.com/asmeurer/sympy/tree/integration3
>
> This is just excellent!
>
> I would like to invite everyone to try this. (Read Aaron's email above
> for things to try and not to try yet.) So here are some tougher cases:
>
> In [1]: risch_integrate(1/(x**8+1), x)
> [hangs]

This is our old friend expand() again (if you break anything that hangs in sympy these days, it usually ends up being in expand in the traceback). I am hoping that Mateusz's continual Poly improvements will make this go away eventually.

>
>
> In [4]: cancel(risch_integrate(1/(x**9+1), x).diff(x))
> Out[4]:
> d ⎛ ⎛ 6 3 ⎞⎞ 3 d ⎛
> 1 + 3⋅──⎝RootSum⎝531441⋅t + 729⋅t + 1, Λ(t, t⋅log(x + 9⋅t))⎠⎠ + 3⋅x ⋅──⎝Root
> dx dx
> ──────────────────────────────────────────────────────────────────────────────
> 3
> 3 + 3⋅x
>
> ⎛ 6 3 ⎞⎞
> Sum⎝531441⋅t + 729⋅t + 1, Λ(t, t⋅log(x + 9⋅t))⎠⎠
>
> ──────────────────────────────────────────────────
>
>
> It should be equivalent but it's a bit messy. Maybe we need to
> implement some symbolic manipulation of RootSum().

Mateusz should look at this.

>
> In [18]: risch_integrate(sqrt(1+exp(x)), x)
> NotImplementedError: Couldn't find an elementary transcendental
> extension for (1 + exp(x))**(1/2). Try using a manual extension with
> the extension flag.

No, and it won't be any time soon. The algorithm I am implementing is only the transcendental case (no algebraic functions like sqrt). Now, the algebraic part does exist, but it is much more difficult to implement, and I won't even start to do it until the transcendental case is finished. If you get the above error, and you believe that the function is really not algebraic, please post it here because it could be a bug in my preparser algorithm, (be aware that you could be wrong, though).

>
>
> [18] is probably not supported yet.
>
>
> Otherwise it works great. I wasn't able to make it break.

Well, I know it's possible, because I do it all the time (and then I go and fix the bug).

>
>
> Ondrej

Aaron Meurer

[Christian Muise]

I'm having issues with the branch:

>>> from sympy import *

Traceback (most recent call last):

  File "<input>", line 1, in <module>

  File "/usr/local/lib/python2.6/dist-packages/sympy/__init__.py", line 24, in <module>

    from polys import *

  File "/usr/local/lib/python2.6/dist-packages/sympy/polys/__init__.py", line 3, in <module>

    from polytools import (

  File "/usr/local/lib/python2.6/dist-packages/sympy/polys/polytools.py", line 66, in <module>

    from sympy.polys.domains import FF, QQ

ImportError: No module named domains

>>> 

  I guess it's the version I'm using? Python 2.6.5 is installed.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To post to this group, send email to sy...@googlegroups.com.
To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/sympy?hl=en.

[Christian Muise]

Apparently the setup.py is dated:

cjmuise@haz:~/Projects/sympy$ git diff

diff --git a/setup.py b/setup.py

index bdec104..538904d 100755

--- a/setup.py

+++ b/setup.py

@@ -71,6 +71,7 @@

     'sympy.mpmath.matrices',

     'sympy.mpmath.calculus',

     'sympy.polys',

+    'sympy.polys.domains',

     'sympy.printing',

     'sympy.printing.pretty',

     'sympy.series',

cjmuise@haz:~/Projects/sympy$ 

[Christian Muise]

  So it's been about 6 years since I've done calculus, and in fact the only math I see from day to day is pretty much 100% symbolic. That being said, I've attached a script that may help, or at least satisfy some curiosity. It's crude, but simple enough to toy around with and checks that the integral of the diff is the same as the original (not sure what to do to check this doesn't require constant factors). Just running it a few times finds NotImplemented errors, odd looking results, etc. This is about all I can offer as far as testing the thing goes ;).

  Cheers

[Aaron S. Meurer]

This is an artifact from Mateusz's polys11, which I am merged into.  There are some other bugs from that branch that you might run into, such as if you try using a ground type other than gmpy or a python other than 2.7. Basically, he made some changes to the domains that aren't quite ready yet, but I had to merge into his branch anyway because of a serious bug that was fixed there.  

By the way, why do you install the package?  There is no need to do that.  Just run ./bin/isympy from the sympy directory.

Aaron Meurer

[Christian Muise]

By the way, why do you install the package?  There is no need to do that.  Just run ./bin/isympy from the sympy directory.

  Habit ;). I like running my python files from the command line, and `bin/isympy ri.py` executes the file before the preface for isympy (setting up the environment, etc).

[Aaron S. Meurer]

Good idea.  The ability to generate random expressions is very useful.  We should make this more general and add it to sympy/utilities/iterables.py.

I had to make a few modifications to your script. First off, you were testing if integral(function.diff(x)) == function, which has two problems: first, you don't even pass it to simplify().  As I said, the integral will usually not look like the original function.  

But more importantly, it need not even be true.  The function and the integral of the derivative of the function can differer by a constant.  

Second, it's better to catch NotImplementedError and continue.

Third, you should never use from import * outside of the interpreter (i.e., in a script).  I know it's lazy, but it's a bad habit to get into.

Fourth, I also regular print the original expression, which makes it easier to copy and paste it if it is wrong and I need to debug it.

I attached an improved script.  Thanks for taking the time to look into this.

Aaron Meurer

[Christian Muise]

Good idea.  The ability to generate random expressions is very useful.  We should make this more general and add it to sympy/utilities/iterables.py.

  It's mesmerizing...may make the thing into a screen saver ;). 

I had to make a few modifications to your script. First off, you were testing if integral(function.diff(x)) == function, which has two problems: first, you don't even pass it to simplify().  As I said, the integral will usually not look like the original function.  

  Ahh, right. Does it make sense to simplify the original formula before returning it from build_func? Attached is a copy that simplifies that equation, as well as the diff.

 

Fourth, I also regular print the original expression, which makes it easier to copy and paste it if it is wrong and I need to debug it.

I attached an improved script.  Thanks for taking the time to look into this.

  So the attached alteration (which is just a few edits to yours...we should invent a system to control versions of files or something ;)) also makes the seed explicit. You can copy the printed seed into the source to repeat the same tests.

  Cheers

[Aaron S. Meurer]

It already exists.  It's called GitHub gist.

  Cheers

Aaron Meurer

[Ondrej Certik]

On Thu, Aug 19, 2010 at 12:18 PM, Christian Muise
<christi...@gmail.com> wrote:
>> Good idea.  The ability to generate random expressions is very useful.  We
>> should make this more general and add it to sympy/utilities/iterables.py.

Yes, definitely. We need such things for all kinds of things in sympy
(limits, series expansion) and so on.

So that it's easy to write a simple script that randomly tries all
kinds of formulas and spits all failures.

Ondrej