Improve the way that sage evaluates symbolic functions on basic types

[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 ​http://trac.sagemath.org/sage_trac/ticket/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​

[burcin]

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:1173​. I'll try to rebase your patch on top of that.

[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:1173​. There 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?

[bober]

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:1173​. There 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]

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.patch. If 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]

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

[bober]

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.

[jdemeyer]

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

[bober]

I've attached a patch rebased to 5.0-beta5. I had no problem applying the patches, so I just made sure all tests still pass.

Maybe I should just go ahead and put it back to positive review, but I'll give Burcin a day or so to take a look at it.

[burcin]

Looks good to me.

[jdemeyer]

Rebasing to sage-5.0.beta5 isn't sufficient, since #9130 was only merged in sage-5.0.beta6. So, this patch still needs an extra rebase.

[bober]

Oops. I guess that explains why I didn't actually need to do anything to rebase. I'll try to remember to update this later today.

[bober]

Ok, now it is rebased to 5.0beta7.

[jdemeyer]

There is numerical noise. On hawk (OpenSolaris? 06.2009-32):

sage -t  --long -force_lib devel/sage/sage/symbolic/function.pyx
Warning: divide by zero encountered in divide
**********************************************************************
File "/export/home/buildbot/build/sage/hawk-1/hawk_full/build/sage-5.0.beta8/devel/sage-main/sage/symbolic/function.pyx", line 185:
    sage: cot(complex(1,2))
Expected:
    (0.03279775553375259-0.9843292264581908j)
Got:
    (0.03279775553375261-0.9843292264581911j)
**********************************************************************
sage -t  --long -force_lib devel/sage/sage/functions/hyperbolic.py
Warning: divide by zero encountered in divide
Warning: divide by zero encountered in divide
**********************************************************************
File "/export/home/buildbot/build/sage/hawk-1/hawk_full/build/sage-5.0.beta8/devel/sage-main/sage/functions/hyperbolic.py", line 543:
    sage: float(arccoth(2))
Expected:
    0.5493061443340548
Got:
    0.5493061443340549
**********************************************************************

[jdemeyer]

On bsd (OS X 10.6 64-bit):

sage -t  --long -force_lib devel/sage/sage/functions/hyperbolic.py
Warning: divide by zero encountered in divide
Warning: divide by zero encountered in divide
**********************************************************************
File "/Users/buildbot/build/sage/bsd-1/bsd_64_full/build/sage-5.0.beta8/devel/sage-main/sage/functions/hyperbolic.py", line 543:
    sage: float(arccoth(2))
Expected:
    0.5493061443340548
Got:
    0.5493061443340549
**********************************************************************

[jdemeyer]

On iras (ia64):

sage -t  --long -force_lib devel/sage/sage/rings/arith.py
**********************************************************************
File "/home/buildbot/build/sage/cleo-1/cleo_full/build/sage-5.0.beta8/devel/sage-main/sage/rings/arith.py", line 3066:
    sage: abs(binomial(0.5r, 5) - 0.02734375) < 1e-17r
Expected:
    True
Got:
    False
**********************************************************************

[Snark]

Replying to jdemeyer:

On iras (ia64):

sage -t  --long -force_lib devel/sage/sage/rings/arith.py
**********************************************************************
File "/home/buildbot/build/sage/cleo-1/cleo_full/build/sage-5.0.beta8/devel/sage-main/sage/rings/arith.py", line 3066:
    sage: abs(binomial(0.5r, 5) - 0.02734375) < 1e-17r
Expected:
    True
Got:
    False
**********************************************************************

Three remarks:

1. that test is also one that is failing on ARM ; and there's a trac ticket about binomial's implementation being pretty wrong (#12448).
2. should that test use "#tol abs" instead of doing the computation itself?
3. shouldn't it be even a "#tol rel" ?

[jdemeyer]

On cicero (Linux i386):

sage -t  --long -force_lib devel/sage/sage/symbolic/function.pyx
Warning: divide by zero encountered in divide
**********************************************************************
File "/home/buildbot/build/sage/cicero-1/cicero_full/build/sage-5.0.beta8/devel/sage-main/sage/symbolic/function.pyx", line 185:
    sage: cot(complex(1,2))
Expected:
    (0.03279775553375259-0.9843292264581908j)
Got:
    (0.032797755533752596-0.9843292264581911j)
**********************************************************************

[bober]

Ok, here is a version that should quiet those numerical errors. I'm leaving this as need review for the moment instead of putting it back to positive review in case someone wants to revisit the issue of whether or not we might want to make sure that we don't return different results on different machines.

[jdemeyer]

I think it's better to use RR (= RealField(53)) or CC (=ComplexField(53)) instead of RDF or CDF. That should give more consistent results.

[dsm]

Missing "-":

            sage: abs(float(arccoth(2)) - 0.5493061443340548) < 1e16 

[bober]

Replying to dsm:

Missing "-":

            sage: abs(float(arccoth(2)) - 0.5493061443340548) < 1e16 

Oops. That is fixed, and the 16 is changed to a 15, since 16 is actually too large, I think. I also took the opportunity to write a better commit message.

Of course, I'm still not sure if this is the right way to do things.

[bober]

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]

Replying to 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]

* ping *

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

[rws]

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.

[rws]

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

[jdemeyer]

Replying to rws:

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

No, this overlaps a lot with #17130.

[rws]

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