rewrite fast_float to support more datatypes

[sagetrac-cwitty]

I've mentioned several times over the past months that I'm rewriting fast_float. Here's a preliminary patch, that shows the direction I'm taking (automatically generate interpreters for each type, by having Python code write C code).

This version replaces fast_float with the new code, and also adds a new entry point, fast_callable.

fast_callable(EXPR, domain=R)

(where EXPR is a symbolic expression or a polynomial) is essentially equivalent to evaluating EXPR with calls to R() at every node. There's special fast support for domain=RDF or domain=RealField(n), and there's slowish generic code that should work for arbitrary R.

The code is not ready for submission... there are very few doctests, and a lot of the documentation is simply wrong. But if anybody has any comments, I'd be happy to hear them.

CC: @robertwb @jasongrout @sagetrac-bober

Component: basic arithmetic

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

[robertwb]

comment:1

This looks very interesting, and much further reaching than the original fast_float code I mostly wrote on a plane ride home from DC. Wish I had time to review it.

[sagetrac-cwitty]

Attachment: trac5093-fast-callable.patch.gz

[sagetrac-cwitty]

comment:3

OK, I've put up a new patch which I think is ready for review. The big feature here is support for multiple types -- we get accelerated evaluation over RealField(k) as well as RDF, and even the unaccelerated Python-object evaluator (which calls PyNumber_Add and friends) is still faster than the evaluators for multivariate and univariate polynomials over QQ.

We're also typically somewhat faster than the old fast_float, but often not by a lot. (In the three benchmarks below, the speed gain ranges from 2% faster to more than 2x as fast.) There's still a lot of optimization I can do.

Here are the benchmark details:

sage: set_random_seed(0)
sage: K.<x,y,z> = QQ[]
sage: p = K.random_element(10, 20)
sage: p_rr = p.change_ring(RR)
sage: p_rdf = p.change_ring(RDF)
sage: fp_old = fast_float(p, x, y, z, old=True)
sage: fp_rdf = fast_callable(p, domain=RDF)
sage: fp_rr = fast_callable(p, domain=RR)
sage: fp_any = fast_callable(p)
sage: 
sage: timeit('p(1,2,3)')
625 loops, best of 3: 369 µs per loop
sage: timeit('fp_any(1,2,3)')
625 loops, best of 3: 115 µs per loop
sage: 
sage: timeit('p(1.0,2.0,3.0)')
625 loops, best of 3: 460 µs per loop
sage: timeit('p_rr(1.0,2.0,3.0)')
625 loops, best of 3: 341 µs per loop
sage: timeit('fp_rr(1.0,2.0,3.0)')
625 loops, best of 3: 121 µs per loop
sage: 
sage: timeit('p(1.0r,2.0r,3.0r)')
625 loops, best of 3: 334 µs per loop
sage: timeit('p_rdf(1.0r,2.0r,3.0r)')
625 loops, best of 3: 101 µs per loop
sage: timeit('fp_old(1.0r,2.0r,3.0r)')
625 loops, best of 3: 4.27 µs per loop
sage: timeit('fp_rdf(1.0r,2.0r,3.0r)')
625 loops, best of 3: 4.18 µs per loop
sage: 
sage: set_random_seed(0)
sage: K.<x> = QQ[]
sage: p = K.random_element(20)
sage: p_rr = p.change_ring(RR)
sage: p_rdf = p.change_ring(RDF)
sage: fp_old = fast_float(p, x, old=True)
sage: fp_rdf = fast_callable(p, domain=RDF)
sage: fp_rr = fast_callable(p, domain=RR)
sage: fp_any = fast_callable(p)
sage: 
sage: timeit('p(2)')
625 loops, best of 3: 53.3 µs per loop
sage: timeit('fp_any(2)')
625 loops, best of 3: 37 µs per loop
sage: 
sage: timeit('p(2.0)')
625 loops, best of 3: 54.9 µs per loop
sage: timeit('p_rr(2.0)')
625 loops, best of 3: 5.76 µs per loop
sage: timeit('fp_rr(2.0)')
625 loops, best of 3: 7.8 µs per loop
sage: 
sage: timeit('p(2.0r)')
625 loops, best of 3: 153 µs per loop
sage: timeit('p_rdf(2.0r)')
625 loops, best of 3: 20.6 µs per loop
sage: timeit('fp_old(2.0r)')
625 loops, best of 3: 973 ns per loop
sage: timeit('fp_rdf(2.0r)')
625 loops, best of 3: 403 ns per loop
sage: 
sage: var('x,y,z')
(x, y, z)
sage: v = sin(sqrt(x*x + y^3) + z/4) * (x+y+z^5)
sage: fv_old = fast_float(v, x, y, z, old=True)
sage: fv_rdf = fast_callable(v, domain=RDF)
sage: fv_rr = fast_callable(v, domain=RR)
sage: fv_any = fast_callable(v)
sage: 
sage: timeit('v(x=1.0,y=2.0,z=3.0)')
625 loops, best of 3: 438 µs per loop
sage: timeit('fv_rr(1.0,2.0,3.0)')
625 loops, best of 3: 35.5 µs per loop
sage: 
sage: timeit('v(x=2.0r,y=2.0r,z=3.0r)')
625 loops, best of 3: 400 µs per loop
sage: timeit('fv_old(1.0r,2.0r,3.0r)')
625 loops, best of 3: 680 ns per loop
sage: timeit('fv_rdf(1.0r,2.0r,3.0r)')
625 loops, best of 3: 544 ns per loop

[jasongrout]

comment:8

Is anyone looking at reviewing this ticket? Lots of people are getting emailed about updates, so apparently there is interest :).

I think it could probably benefit from more than one pair of eyes since it is a major rewrite of a very key piece of functionality. I might be able to look at at least some of it next week. Who's with me? :)