Improve the way that sage evaluates symbolic functions on basic types

[sagetrac-bober]

Currently when a symbolic function (e.g. gamma(), exp(), sin()) gets a python float, it passes it off to ginac, and let's ginac evaluate it. This can be very slow, e.g.:

sage: timeit('exp(6.0r)')
625 loops, best of 3: 476 µs per loop

This also leads to possibly undesirable inconsistencies across different platforms: for example, ginac ends up calling the local libc tgammal to evaluate the gamma function, and even when this tgammal is evaluated by the same version of eglibc on two different platforms the answer varies depending on whether 80 bit doubles are available.

To fix this, we make Sage check its types which correspond most closely to python type for the function that should be called, and use the functions for those types if available.

Now we have

sage: timeit('exp(6.0r)')
625 loops, best of 3: 1.02 µs per loop

For the specific case of the gamma function, the is also a change here to make RDF use python's math.gamma(). It is more accurate than gsl. It is perhaps not as accurate as eglibc's gamma(), but it should give reliable results on different platforms.

I think this will fix 3 of the 4 failing numerical accuracy doctests on ARM, though the test for binomial(.5r, 5) had to be weakened slightly. (see http://groups.google.com/group/sage-devel/browse_thread/thread/3c8c61ea113ea60c ) binomial() is currently not computed in a very good way, though, so it is reasonable to weaken this test temporarily. (see #12448 )

Finally, I'll remark that these improvements could probably be much better. The above timings should be compared to

sage: timeit('math.exp(6.0r)')
625 loops, best of 3: 112 ns per loop

Apply:

• attachment: trac12449.patch

Depends on #4498
Depends on #12507
Depends on #9130

CC: @SnarkBoojum @eviatarbach

Component: symbolics

Keywords: gamma function

Author: Jonathan Bober, Burcin Erocal

Reviewer: Burcin Erocal, Jonathan Bober

Issue created by migration from https://trac.sagemath.org/ticket/12449

[sagetrac-bober]

Author: bober

[burcin]

comment:2

Thanks for the patch. I agree that considering float and complex types separately in the fast evaluation path makes sense for speed. However, this causes inconsistencies with the results you get when a value is substituted into the expression.

We still need to modify the py_tgamma() function in sage/symbolic/pynac.pyx to fix this:

sage: gamma(x).subs(x=6.r)          
120.0

AFAICT (without applying the patch and testing), attachment: 12449.patch removes the try/except block around return getattr(args[0], self._name)(), disabling the call to BuiltinFunction.__call__ later on in the function body.

This would be a good occasion to merge my patch from attachment: trac_1173-move_call.patch:ticket:1173I'll try to rebase your patch on top of that.

[burcin]

Changed author from bober to Jonathan Bober

[burcin]

Reviewer: Burcin Erocal

[burcin]

comment:3

Replying to @burcin:

AFAICT (without applying the patch and testing), attachment: 12449.patch removes the try/except block around return getattr(args[0], self._name)(), disabling the call to BuiltinFunction.__call__ later on in the function body.

Sorry. I was looking at the code after applying attachment: trac_1173-move_call.patch:ticket:1173There is no try/except block in the original code.

Regarding the changes to sage/rings/real_double.pyx: IIRC, python's math module just calls the underlying libc functions. On ARM won't this introduce new errors?

[sagetrac-bober]

comment:4

Replying to @burcin:

Replying to @burcin:

AFAICT (without applying the patch and testing), attachment: 12449.patch removes the try/except block around return getattr(args[0], self._name)(), disabling the call to BuiltinFunction.__call__ later on in the function body.

Sorry. I was looking at the code after applying attachment: trac_1173-move_call.patch:ticket:1173There is no try/except block in the original code.

Regarding the changes to sage/rings/real_double.pyx: IIRC, python's math module just calls the underlying libc functions. On ARM won't this introduce new errors?

For the gamma function, at least, python has its own implementation. Of course, that function itself might call some underlying libc functions, so maybe it won't be entirely consistent, now that I think about it.

In the bit of testing that I did, math.gamma was less accurate than tgammal on my computer, but much more accurate than gsl's gamma function.

[burcin]

Changed author from Jonathan Bober to Jonathan Bober, Burcin Erocal

[burcin]

comment:5

It turns out that Python implements the gamma function itself:

http://hg.python.org/cpython/file/58bd6a58365d/Modules/mathmodule.c#l226

I uploaded a new patch that changes py_tgamma() in sage/symbolic/pynac.pyx to use Python's gamma implementation.

I give a positive review to attachment: 12449.patchIf someone can review my patch we can switch this to positive review.

BTW, this will fix #9162. We should close that as a duplicate.

[burcin]

comment:6

I uploaded a new version of my patch since the previous function had a typo that make doctests fail.

[sagetrac-bober]

comment:7

Well, I suppose that I will just go ahead and give it a positive review, since your changes are simple and orthogonal to mine (and I checked that all tests still pass on my machine). I would like to get some confirmation that this fixes things on ARM, but I suppose that we'll find that out eventually; that's a secondary point as far as this patch goes, though, so it isn't that important.

[sagetrac-bober]

Changed reviewer from Burcin Erocal to Burcin Erocal, Jonathan Bober

[jdemeyer]

comment:8

This should be rebased to #4498 + #12507 + #9130.

[jdemeyer]

Dependencies: #4498, #12507, #9130

[sagetrac-bober]

comment:25

Attachment: trac12449.patch.gz

Ok, after a brief discussion at SD 40.5 (after Jeroen went to sleep and wasn't there to complain) I'm restoring this to positive review. A few people seemed to agree that if a function is called on a machine floating point type then it is reasonable for it to give machine dependent answers.

I hope that all the tests actually pass now. Also, I've twice updated the patch to remove trailing whitespace, so I hope that I actually got it all this time.

(Of course, Jeroen will still have an opportunity to object...)

[jdemeyer]

comment:26

Replying to @sagetrac-bober:

A few people seemed to agree that if a function is called on a machine floating point type then it is reasonable for it to give machine dependent answers.

Obviously the following comment doesn't apply then:

# It is more accurate than gsl's gamma function and 
# should provide consistant results across platforms. 

Also: you should add tests for the infinity stuff that you added.

[jdemeyer]

comment:27

*** ping ***

The issues I mentioned in the last comment are fairly trivial.

[rwst]

comment:33

The patch actually contains code that switches FP sage_tgamma1 from GSL to python_gamma which might be a good thing, but requires a separate ticket. The timings are identical, so I suspect that the implementation might be too.

[rwst]

comment:34

@jdemeyer: Is the patch still relevant after the changes you made to BuiltinFunction.__call__?

[jdemeyer]

comment:36

Replying to @rwst:

@jdemeyer: Is the patch still relevant after the changes you made to BuiltinFunction.__call__?

No, this overlaps a lot with #17130.

[rwst]

comment:37

The original issue seems resolved:

sage: %timeit('exp(6.0r)')
100000000 loops, best of 3: 7.31 ns per loop
sage: %timeit('math.exp(6.0r)')
100000000 loops, best of 3: 7.28 ns per loop