Fast Expression Evaluation

For many applications such as numerical integration, differential equation approximation, plotting a 3d surface, optimization problems, Monte-Carlo simulations, etc., one wishes to pass around and evaluate a single algebraic expression many, many times at various floating point values. Other applications may need to evaluate an expression many times in interval arithmetic, or in a finite field. Doing this via recursive calls over a python representation of the object (even if Maxima or other outside packages are not involved) is extremely inefficient.

This module provides a function, fast_callable(), to transform such expressions into a form where they can be evaluated quickly:

sage: f = sin(x) + 3*x^2
sage: ff = fast_callable(f, vars=[x])
sage: ff(3.5)
36.3992167723104
sage: ff(RIF(3.5))
36.39921677231038?

By default, fast_callable() only removes some interpretive overhead from the evaluation, but all of the individual arithmetic operations are done using standard Sage arithmetic. This is still a huge win over sage.calculus, which evidently has a lot of overhead. Compare the cost of evaluating Wilkinson’s polynomial (in unexpanded form) at x=30:

sage: wilk = prod((x-i) for i in [1 .. 20]); wilk
(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10)*(x - 11)*(x - 12)*(x - 13)*(x - 14)*(x - 15)*(x - 16)*(x - 17)*(x - 18)*(x - 19)*(x - 20)
sage: timeit('wilk.subs(x=30)') # random, long time
625 loops, best of 3: 1.43 ms per loop
sage: fc_wilk = fast_callable(wilk, vars=[x])
sage: timeit('fc_wilk(30)') # random, long time
625 loops, best of 3: 9.72 us per loop

You can specify a particular domain for the evaluation using domain=:

sage: fc_wilk_zz = fast_callable(wilk, vars=[x], domain=ZZ)

The meaning of domain=D is that each intermediate and final result is converted to type D. For instance, the previous example of sin(x) + 3*x^2 with domain=D would be equivalent to D(D(sin(D(x))) + D(D(3)*D(D(x)^2))). (This example also demonstrates the one exception to the general rule: if an exponent is an integral constant, then it is not wrapped with D().)

At first glance, this seems like a very bad idea if you want to compute quickly. And it is a bad idea, for types where we don’t have a special interpreter. It’s not too bad of a slowdown, though. To mitigate the costs, we check whether the value already has the correct parent before we call D.

We don’t yet have a special interpreter with domain ZZ, so we can see how that compares to the generic fc_wilk example above:

sage: timeit('fc_wilk_zz(30)') # random, long time
625 loops, best of 3: 15.4 us per loop

However, for other types, using domain=D will get a large speedup, because we have special-purpose interpreters for those types. One example is RDF. Since with domain=RDF we know that every single operation will be floating-point, we can just execute the floating-point operations directly and skip all the Python object creations that you would get from actually using RDF objects:

sage: fc_wilk_rdf = fast_callable(wilk, vars=[x], domain=RDF)
sage: timeit('fc_wilk_rdf(30.0)') # random, long time
625 loops, best of 3: 7 us per loop

The domain does not need to be a Sage type; for instance, domain=float also works. (We actually use the same fast interpreter for domain=float and domain=RDF; the only difference is that when domain=RDF is used, the return value is an RDF element, and when domain=float is used, the return value is a Python float.)

sage: fc_wilk_float = fast_callable(wilk, vars=[x], domain=float)
sage: timeit('fc_wilk_float(30.0)') # random, long time
625 loops, best of 3: 5.04 us per loop

We also have support for RR:

sage: fc_wilk_rr = fast_callable(wilk, vars=[x], domain=RR)
sage: timeit('fc_wilk_rr(30.0)') # random, long time
625 loops, best of 3: 13 us per loop

For CC:

sage: fc_wilk_cc = fast_callable(wilk, vars=[x], domain=CC)
sage: timeit('fc_wilk_cc(30.0)') # random, long time
625 loops, best of 3: 23 us per loop

And support for CDF:

sage: fc_wilk_cdf = fast_callable(wilk, vars=[x], domain=CDF)
sage: timeit('fc_wilk_cdf(30.0)') # random, long time
625 loops, best of 3: 10.2 us per loop

Currently, fast_callable() can accept two kinds of objects: polynomials (univariate and multivariate) and symbolic expressions (elements of the Symbolic Ring). (This list is likely to grow significantly in the near future.) For polynomials, you can omit the ‘vars’ argument; the variables will default to the ring generators (in the order used when creating the ring).

sage: K.<x,y,z> = QQ[]
sage: p = 10*y + 100*z + x
sage: fp = fast_callable(p)
sage: fp(1,2,3)
321

But you can also specify the variable names to override the default ordering (you can include extra variable names here, too).

sage: fp = fast_callable(p, vars=('x','w','z','y'))

For symbolic expressions, you need to specify the variable names, so that fast_callable() knows what order to use.

sage: var('y,z,x')
(y, z, x)
sage: f = 10*y + 100*z + x
sage: ff = fast_callable(f, vars=(x,y,z))
sage: ff(1,2,3)
321

You can also specify extra variable names:

sage: ff = fast_callable(f, vars=('x','w','z','y'))
sage: ff(1,2,3,4)
341

This should be enough for normal use of fast_callable(); let’s discuss some more advanced topics.

Sometimes it may be useful to create a fast version of an expression without going through symbolic expressions or polynomials; perhaps because you want to describe to fast_callable() an expression with common subexpressions.

Internally, fast_callable() works in two stages: it constructs an expression tree from its argument, and then it builds a fast evaluator from that expression tree. You can bypass the first phase by building your own expression tree and passing that directly to fast_callable(), using an ExpressionTreeBuilder.

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=('x','y','z'))

An ExpressionTreeBuilder has three interesting methods: constant(), var(), and call(). All of these methods return ExpressionTree objects.

The var() method takes a string, and returns an expression tree for the corresponding variable.

sage: x = etb.var('x')
sage: y = etb.var('y')
sage: z = etb.var('y')

Expression trees support Python’s numeric operators, so you can easily build expression trees representing arithmetic expressions.

sage: v1 = (x+y)*(y+z) + (y//z)

The constant() method takes a Sage value, and returns an expression tree representing that value.

sage: v2 = etb.constant(3.14159) * x + etb.constant(1729) * y

The call() method takes a sage/Python function and zero or more expression trees, and returns an expression tree representing the function call.

sage: v3 = etb.call(sin, v1+v2)
sage: v3
sin(add(add(mul(add(v_0, v_1), add(v_1, v_1)), floordiv(v_1, v_1)), add(mul(3.14159000000000, v_0), mul(1729, v_1))))

Many sage/Python built-in functions are specially handled; for instance, when evaluating an expression involving sin() over RDF, the C math library function sin() is called. Arbitrary functions are allowed, but will be much slower since they will call back to Python code on every call; for example, the following will work.

sage: def my_sqrt(x): return pow(x, 0.5)
sage: e = etb.call(my_sqrt, v1); e
{my_sqrt}(add(mul(add(v_0, v_1), add(v_1, v_1)), floordiv(v_1, v_1)))
sage: fast_callable(e)(1, 2, 3)
3.60555127546399

To provide fast_callable() for your own class (so that fast_callable(x) works when x is an instance of your class), implement a method _fast_callable_(self, etb) for your class. This method takes an ExpressionTreeBuilder, and returns an expression tree built up using the methods described above.

EXAMPLES:

sage: var('x')
x
sage: f = fast_callable(sqrt(x^7+1), vars=[x], domain=float)

sage: f(1)
1.4142135623730951
sage: f.op_list()
[('load_arg', 0), ('ipow', 7), ('load_const', 1.0), 'add', 'sqrt', 'return']

To interpret that last line, we load argument 0 (‘x’ in this case) onto the stack, push the constant 7.0 onto the stack, call the pow function (which takes 2 arguments from the stack), push the constant 1.0, add the top two arguments of the stack, and then call sqrt.

Here we take sin of the first argument and add it to f:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder('x')
sage: x = etb.var('x')
sage: f = etb.call(sqrt, x^7 + 1)
sage: g = etb.call(sin, x)
sage: fast_callable(f+g).op_list()
[('load_arg', 0), ('ipow', 7), ('load_const', 1), 'add', ('py_call', <function sqrt at ...>, 1), ('load_arg', 0), ('py_call', sin, 1), 'add', 'return']

AUTHOR:

• Carl Witty (2009-02): initial version (heavily inspired by Robert Bradshaw’s fast_eval.pyx)

Todo

The following bits of text were written for the module docstring. They are not true yet, but I hope they will be true someday, at which point I will move them into the main text.

The final interesting method of ExpressionTreeBuilder is choice(). This produces conditional expressions, like the C COND ? T : F expression or the Python T if COND else F. This lets you define piecewise functions using fast_callable().

sage: v4 = etb.choice(v3 >= etb.constant(0), v1, v2)  # not tested

The arguments are (COND, T, F) (the same order as in C), so the above means that if v3 evaluates to a nonnegative number, then v4 will evaluate to the result of v1; otherwise, v4 will evaluate to the result of v2.

Let’s see an example where we see that fast_callable() does not evaluate common subexpressions more than once. We’ll make a fast_callable() expression that gives the result of 16 iterations of the Mandelbrot function.

sage: etb = ExpressionTreeBuilder('c')
sage: z = etb.constant(0)
sage: c = etb.var('c')
sage: for i in range(16):
....:     z = z*z + c
sage: mand = fast_callable(z, domain=CDF)

Now ff does 32 complex arithmetic operations on each call (16 additions and 16 multiplications). However, if z*z produced code that evaluated z twice, then this would do many thousands of arithmetic operations instead.

Note that the handling for common subexpressions only checks whether expression trees are the same Python object; for instance, the following code will evaluate x+1 twice:

sage: etb = ExpressionTreeBuilder('x')
sage: x = etb.var('x')
sage: (x+1)*(x+1)
mul(add(v_0, 1), add(v_0, 1))

but this code will only evaluate x+1 once:

sage: v = x+1; v*v
mul(add(v_0, 1), add(v_0, 1))

class sage.ext.fast_callable.CompilerInstrSpec(n_inputs, n_outputs, parameters)

Bases: object

Describe a single instruction to the fast_callable code generator.

An instruction has a number of stack inputs, a number of stack outputs, and a parameter list describing extra arguments that must be passed to the InstructionStream.instr method (that end up as extra words in the code).

The parameter list is a list of strings. Each string is one of the following:

• ‘args’ - The instruction argument refers to an input argument of the wrapper class; it is just appended to the code.

• ‘constants’, ‘py_constants’ - The instruction argument is a value; the value is added to the corresponding list (if it’s not already there) and the index is appended to the code.

• ‘n_inputs’, ‘n_outputs’ - The instruction actually takes a variable number of inputs or outputs (the n_inputs and n_outputs attributes of this instruction are ignored). The instruction argument specifies the number of inputs or outputs (respectively); it is just appended to the code.

class sage.ext.fast_callable.Expression

Bases: object

Represent an expression for fast_callable.

Supports the standard Python arithmetic operators; if arithmetic is attempted between an Expression and a non-Expression, the non-Expression is converted to an expression (using the __call__ method of the Expression’s ExpressionTreeBuilder).

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: x = etb.var(x)
sage: etb(x)
v_0
sage: etb(3)
3
sage: etb.call(sin, x)
sin(v_0)
sage: (x+1)/(x-1)
div(add(v_0, 1), sub(v_0, 1))
sage: x//5
floordiv(v_0, 5)
sage: -abs(~x)
neg(abs(inv(v_0)))

abs()

Compute the absolute value of an Expression.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: x = etb(x)
sage: abs(x)
abs(v_0)
sage: x.abs()
abs(v_0)
sage: x.__abs__()
abs(v_0)

class sage.ext.fast_callable.ExpressionCall

Bases: Expression

An Expression that represents a function call.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: type(etb.call(sin, x))
<class 'sage.ext.fast_callable.ExpressionCall'>

arguments()

Return the arguments from this ExpressionCall.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: etb.call(sin, x).arguments()
[v_0]

function()

Return the function from this ExpressionCall.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: etb.call(sin, x).function()
sin

class sage.ext.fast_callable.ExpressionChoice

Bases: Expression

A conditional expression.

(It’s possible to create choice nodes, but they don’t work yet.)

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: etb.choice(etb.call(operator.eq, x, 0), 0, 1/x)
(0 if {eq}(v_0, 0) else div(1, v_0))

condition()

Return the condition of an ExpressionChoice.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: x = etb(x)
sage: etb.choice(x, ~x, 0).condition()
v_0

if_false()

Return the false branch of an ExpressionChoice.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: x = etb(x)
sage: etb.choice(x, ~x, 0).if_false()
0

if_true()

Return the true branch of an ExpressionChoice.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: x = etb(x)
sage: etb.choice(x, ~x, 0).if_true()
inv(v_0)

class sage.ext.fast_callable.ExpressionConstant

Bases: Expression

An Expression that represents an arbitrary constant.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: type(etb(3))
<class 'sage.ext.fast_callable.ExpressionConstant'>

value()

Return the constant value of an ExpressionConstant.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: etb(3).value()
3

class sage.ext.fast_callable.ExpressionIPow

Bases: Expression

A power Expression with an integer exponent.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: type(etb.var('x')^17)
<class 'sage.ext.fast_callable.ExpressionIPow'>

base()

Return the base from this ExpressionIPow.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: (etb(33)^42).base()
33

exponent()

Return the exponent from this ExpressionIPow.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: (etb(x)^(-1)).exponent()
-1

class sage.ext.fast_callable.ExpressionTreeBuilder

Bases: object

A class with helper methods for building Expressions.

An instance of this class is passed to _fast_callable_ methods; you can also instantiate it yourself to create your own expressions for fast_callable, bypassing _fast_callable_.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder('x')
sage: x = etb.var('x')
sage: (x+3)*5
mul(add(v_0, 3), 5)

call(fn, *args)

Construct a call node, given a function and a list of arguments.

The arguments will be converted to Expressions using ExpressionTreeBuilder.__call__.

As a special case, notices if the function is operator.pow and the second argument is integral, and constructs an ExpressionIPow instead.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: etb.call(cos, x)
cos(v_0)
sage: etb.call(sin, 1)
sin(1)
sage: etb.call(sin, etb(1))
sin(1)
sage: etb.call(factorial, x+57)
{factorial}(add(v_0, 57))
sage: etb.call(operator.pow, x, 543)
ipow(v_0, 543)

choice(cond, iftrue, iffalse)

Construct a choice node (a conditional expression), given the condition, and the values for the true and false cases.

(It’s possible to create choice nodes, but they don’t work yet.)

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: etb.choice(etb.call(operator.eq, x, 0), 0, 1/x)
(0 if {eq}(v_0, 0) else div(1, v_0))

constant(c)

Turn the argument into an ExpressionConstant, converting it to our domain if we have one.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder('x')
sage: etb.constant(pi)
pi
sage: etb = ExpressionTreeBuilder('x', domain=RealField(200))
sage: etb.constant(pi)
3.1415926535897932384626433832795028841971693993751058209749

var(v)

Turn the argument into an ExpressionVariable. Look it up in the list of variables. (Variables are matched by name.)

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: var('a,b,some_really_long_name')
(a, b, some_really_long_name)
sage: x = polygen(QQ)
sage: etb = ExpressionTreeBuilder(vars=('a','b',some_really_long_name, x))
sage: etb.var(some_really_long_name)
v_2
sage: etb.var('some_really_long_name')
v_2
sage: etb.var(x)
v_3
sage: etb.var('y')
Traceback (most recent call last):
...
ValueError: Variable 'y' not found...

class sage.ext.fast_callable.ExpressionVariable

Bases: Expression

An Expression that represents a variable.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: type(etb.var(x))
<class 'sage.ext.fast_callable.ExpressionVariable'>

variable_index()

Return the variable index of an ExpressionVariable.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=(x,))
sage: etb(x).variable_index()
0

class sage.ext.fast_callable.FastCallableFloatWrapper(ff, imag_tol)

Bases: object

A class to alter the return types of the fast-callable functions.

When applying numerical routines (including plotting) to symbolic expressions and functions, we generally first convert them to a faster form with fast_callable(). That function takes a domain parameter that forces the end (and all intermediate) results of evaluation to a specific type. Though usually always want the end result to be of type float, correctly choosing the domain presents some problems:

• float is a bad choice because it’s common for real functions to have complex terms in them. Moreover precision issues can produce terms like 1.0 + 1e-12*I that are hard to avoid if calling real() on everything is infeasible.

• complex has essentially the same problem as float. There are several symbolic functions like min_symbolic(), max_symbolic(), and floor() that are unable to operate on complex numbers.

• None leaves the types of the inputs/outputs alone, but due to the lack of a specialized interpreter, slows down evaluation by an unacceptable amount.

• CDF has none of the other issues, because CDF has its own specialized interpreter, a lexicographic ordering (for min/max), and supports floor(). However, most numerical functions cannot handle complex numbers, so using CDF would require us to wrap every evaluation in a CDF-to-float conversion routine. That would slow things down less than a domain of None would, but is unattractive mainly because of how invasive it would be to “fix” the output everywhere.

Creating a new fast-callable interpreter that has different input and output types solves most of the problems with a CDF domain, but fast_callable() and the interpreter classes in sage.ext.interpreters are not really written with that in mind. The domain parameter to fast_callable(), for example, is expecting a single Sage ring that corresponds to one interpreter. You can make it accept, for example, a string like “CDF-to-float”, but the hacks required to make that work feel wrong.

Thus we arrive at this solution: a class to wrap the result of fast_callable(). Whenever we need to support intermediate complex terms in a numerical routine, we can set domain=CDF while creating its fast-callable incarnation, and then wrap the result in this class. The __call__ method of this class then ensures that the CDF output is converted to a float if its imaginary part is within an acceptable tolerance.

EXAMPLES:

An error is thrown if the answer is complex:

sage: from sage.ext.fast_callable import FastCallableFloatWrapper
sage: f = sqrt(x)
sage: ff = fast_callable(f, vars=[x], domain=CDF)
sage: fff = FastCallableFloatWrapper(ff, imag_tol=1e-8)
sage: fff(1)
1.0
sage: fff(-1)
Traceback (most recent call last):
...
ValueError: complex fast-callable function result
1.0*I for arguments (-1,)

class sage.ext.fast_callable.InstructionStream

Bases: object

An InstructionStream takes a sequence of instructions (passed in by a series of method calls) and computes the data structures needed by the interpreter. This is the stage where we switch from operating on Expression trees to a linear representation. If we had a peephole optimizer (we don’t) it would go here.

Currently, this class is not very general; it only works for interpreters with a fixed set of memory chunks (with fixed names). Basically, it only works for stack-based expression interpreters. It should be generalized, so that the interpreter metadata includes a description of the memory chunks involved and the instruction stream can handle any interpreter.

Once you’re done adding instructions, you call get_current() to retrieve the information needed by the interpreter (as a Python dictionary).

current_op_list()

Return the list of instructions that have been added to this InstructionStream so far.

It’s OK to call this, then add more instructions.

EXAMPLES:

sage: from sage.ext.interpreters.wrapper_rdf import metadata
sage: from sage.ext.fast_callable import InstructionStream
sage: instr_stream = InstructionStream(metadata, 1)
sage: instr_stream.instr('load_arg', 0)
sage: instr_stream.instr('py_call', math.sin, 1)
sage: instr_stream.instr('abs')
sage: instr_stream.instr('return')
sage: instr_stream.current_op_list()
[('load_arg', 0), ('py_call', <built-in function sin>, 1), 'abs', 'return']

get_current()

Return the current state of the InstructionStream, as a dictionary suitable for passing to a wrapper class.

NOTE: The dictionary includes internal data structures of the InstructionStream; you must not modify it.

EXAMPLES:

sage: from sage.ext.interpreters.wrapper_rdf import metadata
sage: from sage.ext.fast_callable import InstructionStream
sage: instr_stream = InstructionStream(metadata, 1)
sage: instr_stream.get_current()
{'args': 1,
 'code': [],
 'constants': [],
 'domain': None,
 'py_constants': [],
 'stack': 0}
sage: instr_stream.instr('load_arg', 0)
sage: instr_stream.instr('py_call', math.sin, 1)
sage: instr_stream.instr('abs')
sage: instr_stream.instr('return')
sage: instr_stream.current_op_list()
[('load_arg', 0), ('py_call', <built-in function sin>, 1), 'abs', 'return']
sage: instr_stream.get_current()
{'args': 1,
 'code': [0, 0, 3, 0, 1, 12, 2],
 'constants': [],
 'domain': None,
 'py_constants': [<built-in function sin>],
 'stack': 1}

get_metadata()

Return the interpreter metadata being used by the current InstructionStream.

The code generator sometimes uses this to decide which code to generate.

EXAMPLES:

sage: from sage.ext.interpreters.wrapper_rdf import metadata
sage: from sage.ext.fast_callable import InstructionStream
sage: instr_stream = InstructionStream(metadata, 1)
sage: md = instr_stream.get_metadata()
sage: type(md)
<class 'sage.ext.fast_callable.InterpreterMetadata'>

has_instr(opname)

Check whether this InstructionStream knows how to generate code for a given instruction.

EXAMPLES:

sage: from sage.ext.interpreters.wrapper_rdf import metadata
sage: from sage.ext.fast_callable import InstructionStream
sage: instr_stream = InstructionStream(metadata, 1)
sage: instr_stream.has_instr('return')
True
sage: instr_stream.has_instr('factorial')
False
sage: instr_stream.has_instr('abs')
True

instr(opname, *args)

Generate code in this InstructionStream for the given instruction and arguments.

The opname is used to look up a CompilerInstrSpec; the CompilerInstrSpec describes how to interpret the arguments. (This is documented in the class docstring for CompilerInstrSpec.)

EXAMPLES:

sage: from sage.ext.interpreters.wrapper_rdf import metadata
sage: from sage.ext.fast_callable import InstructionStream
sage: instr_stream = InstructionStream(metadata, 1)
sage: instr_stream.instr('load_arg', 0)
sage: instr_stream.instr('sin')
sage: instr_stream.instr('py_call', math.sin, 1)
sage: instr_stream.instr('abs')
sage: instr_stream.instr('factorial')
Traceback (most recent call last):
...
KeyError: 'factorial'
sage: instr_stream.instr('return')
sage: instr_stream.current_op_list()
[('load_arg', 0), 'sin', ('py_call', <built-in function sin>, 1), 'abs', 'return']

load_arg(n)

Add a ‘load_arg’ instruction to this InstructionStream.

EXAMPLES:

sage: from sage.ext.interpreters.wrapper_rdf import metadata
sage: from sage.ext.fast_callable import InstructionStream
sage: instr_stream = InstructionStream(metadata, 12)
sage: instr_stream.load_arg(5)
sage: instr_stream.current_op_list()
[('load_arg', 5)]
sage: instr_stream.load_arg(3)
sage: instr_stream.current_op_list()
[('load_arg', 5), ('load_arg', 3)]

load_const(c)

Add a ‘load_const’ instruction to this InstructionStream.

EXAMPLES:

sage: from sage.ext.interpreters.wrapper_rdf import metadata
sage: from sage.ext.fast_callable import InstructionStream, op_list
sage: instr_stream = InstructionStream(metadata, 1)
sage: instr_stream.load_const(5)
sage: instr_stream.current_op_list()
[('load_const', 5)]
sage: instr_stream.load_const(7)
sage: instr_stream.load_const(5)
sage: instr_stream.current_op_list()
[('load_const', 5), ('load_const', 7), ('load_const', 5)]

Note that constants are shared: even though we load 5 twice, it only appears once in the constant table.

sage: instr_stream.get_current()['constants']
[5, 7]

class sage.ext.fast_callable.IntegerPowerFunction(n)

Bases: object

This class represents the function x^n for an arbitrary integral power n. That is, IntegerPowerFunction(2) is the squaring function; IntegerPowerFunction(-1) is the reciprocal function.

EXAMPLES:

sage: from sage.ext.fast_callable import IntegerPowerFunction
sage: square = IntegerPowerFunction(2)
sage: square
(^2)
sage: square(pi)
pi^2
sage: square(I)
-1
sage: square(RIF(-1, 1)).str(style='brackets')
'[0.0000000000000000 .. 1.0000000000000000]'
sage: IntegerPowerFunction(-1)
(^(-1))
sage: IntegerPowerFunction(-1)(22/7)
7/22
sage: v = Integers(123456789)(54321)
sage: v^9876543210
79745229
sage: IntegerPowerFunction(9876543210)(v)
79745229

class sage.ext.fast_callable.InterpreterMetadata

Bases: object

The interpreter metadata for a fast_callable interpreter. Currently consists of a dictionary mapping instruction names to (CompilerInstrSpec, opcode) pairs, a list mapping opcodes to (instruction name, CompilerInstrSpec) pairs, and a range of exponents for which the ipow instruction can be used. This range can be False (if the ipow instruction should never be used), a pair of two integers (a,b), if ipow should be used for a<=n<=b, or True, if ipow should always be used. When ipow cannot be used, then we fall back on calling IntegerPowerFunction.

See the class docstring for CompilerInstrSpec for more information.

NOTE: You must not modify the metadata.

by_opcode

by_opname

ipow_range

class sage.ext.fast_callable.Wrapper

Bases: object

The parent class for all fast_callable wrappers. Implements shared behavior (currently only debugging).

get_orig_args()

Get the original arguments used when initializing this wrapper.

(Probably only useful when writing doctests.)

EXAMPLES:

sage: fast_callable(sin(x)/x, vars=[x], domain=RDF).get_orig_args()
{'args': 1,
 'code': [0, 0, 16, 0, 0, 8, 2],
 'constants': [],
 'domain': Real Double Field,
 'py_constants': [],
 'stack': 2}

op_list()

Return the list of instructions in this wrapper.

EXAMPLES:

sage: fast_callable(cos(x)*x, vars=[x], domain=RDF).op_list()
[('load_arg', 0), ('load_arg', 0), 'cos', 'mul', 'return']

python_calls()

List the Python functions that are called in this wrapper.

(Python function calls are slow, so ideally this list would be empty. If it is not empty, then perhaps there is an optimization opportunity where a Sage developer could speed this up by adding a new instruction to the interpreter.)

EXAMPLES:

sage: fast_callable(abs(sin(x)), vars=[x], domain=RDF).python_calls()
[]
sage: fast_callable(abs(sin(factorial(x))), vars=[x]).python_calls()
[factorial, sin]

sage.ext.fast_callable.fast_callable(x, domain=None, vars=None, expect_one_var=False)

Given an expression x, compile it into a form that can be quickly evaluated, given values for the variables in x.

Currently, x can be an expression object, an element of SR, or a (univariate or multivariate) polynomial; this list will probably be extended soon.

By default, x is evaluated the same way that a Python function would evaluate it – addition maps to PyNumber_Add, etc. However, you can specify domain=D where D is some Sage parent or Python type; in this case, all arithmetic is done in that domain. If we have a special-purpose interpreter for that parent (like RDF or float), domain=… will trigger the use of that interpreter.

If vars is None and x is a polynomial, then we will use the generators of parent(x) as the variables; otherwise, vars must be specified (unless x is a symbolic expression with only one variable, and expect_one_var is True, in which case we will use that variable).

EXAMPLES:

sage: var('x')
x
sage: expr = sin(x) + 3*x^2
sage: f = fast_callable(expr, vars=[x])
sage: f(2)
sin(2) + 12
sage: f(2.0)
12.9092974268257

We have special fast interpreters for domain=float and domain=RDF. (Actually it’s the same interpreter; only the return type varies.) Note that the float interpreter is not actually more accurate than the RDF interpreter; elements of RDF just don’t display all their digits. We have special fast interpreter for domain=CDF:

sage: f_float = fast_callable(expr, vars=[x], domain=float)
sage: f_float(2)
12.909297426825681
sage: f_rdf = fast_callable(expr, vars=[x], domain=RDF)
sage: f_rdf(2)
12.909297426825681
sage: f_cdf = fast_callable(expr, vars=[x], domain=CDF)
sage: f_cdf(2)
12.909297426825681
sage: f_cdf(2+I)
10.40311925062204 + 11.510943740958707*I
sage: f = fast_callable(expr, vars=('z','x','y'))
sage: f(1, 2, 3)
sin(2) + 12
sage: K.<x> = QQ[]
sage: p = -1/4*x^6 + 1/2*x^5 - x^4 - 12*x^3 + 1/2*x^2 - 1/95*x - 1/2
sage: fp = fast_callable(p, domain=RDF)
sage: fp.op_list()
[('load_arg', 0), ('load_const', -0.25), 'mul', ('load_const', 0.5), 'add', ('load_arg', 0), 'mul', ('load_const', -1.0), 'add', ('load_arg', 0), 'mul', ('load_const', -12.0), 'add', ('load_arg', 0), 'mul', ('load_const', 0.5), 'add', ('load_arg', 0), 'mul', ('load_const', -0.010526315789473684), 'add', ('load_arg', 0), 'mul', ('load_const', -0.5), 'add', 'return']
sage: fp(3.14159)
-552.4182988917153
sage: K.<x,y,z> = QQ[]
sage: p = x*y^2 + 1/3*y^2 - x*z - y*z
sage: fp = fast_callable(p, domain=RDF)
sage: fp.op_list()
[('load_const', 0.0), ('load_const', 1.0), ('load_arg', 0), ('ipow', 1), ('load_arg', 1), ('ipow', 2), 'mul', 'mul', 'add', ('load_const', 0.3333333333333333), ('load_arg', 1), ('ipow', 2), 'mul', 'add', ('load_const', -1.0), ('load_arg', 0), ('ipow', 1), ('load_arg', 2), ('ipow', 1), 'mul', 'mul', 'add', ('load_const', -1.0), ('load_arg', 1), ('ipow', 1), ('load_arg', 2), ('ipow', 1), 'mul', 'mul', 'add', 'return']
sage: fp(e, pi, sqrt(2))   # abs tol 3e-14
21.831120464939584
sage: symbolic_result = p(e, pi, sqrt(2)); symbolic_result
pi^2*e + 1/3*pi^2 - sqrt(2)*pi - sqrt(2)*e
sage: n(symbolic_result)
21.8311204649396

sage: from sage.ext.fast_callable import ExpressionTreeBuilder
sage: etb = ExpressionTreeBuilder(vars=('x','y'), domain=float)
sage: x = etb.var('x')
sage: y = etb.var('y')
sage: expr = etb.call(sin, x^2 + y); expr
sin(add(ipow(v_0, 2), v_1))
sage: fc = fast_callable(expr, domain=float)
sage: fc(5, 7)
0.5514266812416906

Check that fast_callable also works for symbolic functions with evaluation functions:

sage: def evalf_func(self, x, y, parent): return parent(x*y) if parent is not None else x*y
sage: x,y = var('x,y')
sage: f = function('f', evalf_func=evalf_func)
sage: fc = fast_callable(f(x, y), vars=[x, y])
sage: fc(3, 4)
f(3, 4)

And also when there are complex values involved:

sage: def evalf_func(self, x, y, parent): return parent(I*x*y) if parent is not None else I*x*y
sage: g = function('g', evalf_func=evalf_func)
sage: fc = fast_callable(g(x, y), vars=[x, y])
sage: fc(3, 4)
g(3, 4)
sage: fc2 = fast_callable(g(x, y), domain=complex, vars=[x, y])
sage: fc2(3, 4)
12j
sage: fc3 = fast_callable(g(x, y), domain=float, vars=[x, y])
sage: fc3(3, 4)
Traceback (most recent call last):
    ...
TypeError: unable to simplify to float approximation

sage.ext.fast_callable.function_name(fn)

Given a function, return a string giving a name for the function.

For functions we recognize, we use our standard opcode name for the function (so operator.add becomes ‘add’, and sage.all.sin becomes ‘sin’).

For functions we don’t recognize, we try to come up with a name, but the name will be wrapped in braces; this is a signal that we’ll definitely use a slow Python call to call this function. (We may use a slow Python call even for functions we do recognize, if we’re targeting an interpreter without an opcode for the function.)

Only used when printing Expressions.

EXAMPLES:

sage: from sage.ext.fast_callable import function_name
sage: function_name(operator.pow)
'pow'
sage: function_name(cos)
'cos'
sage: function_name(factorial)
'{factorial}'

sage.ext.fast_callable.generate_code(expr, stream)

Generate code from an Expression tree; write the result into an InstructionStream.

In fast_callable, first we create an Expression, either directly with an ExpressionTreeBuilder or with _fast_callable_ methods. Then we optimize the Expression in tree form. (Unfortunately, this step is currently missing – we do no optimizations.)

Then we linearize the Expression into a sequence of instructions, by walking the Expression and sending the corresponding stack instructions to an InstructionStream.

EXAMPLES:

sage: from sage.ext.fast_callable import ExpressionTreeBuilder, generate_code, InstructionStream
sage: etb = ExpressionTreeBuilder('x')
sage: x = etb.var('x')
sage: expr = ((x+pi)*(x+1))
sage: from sage.ext.interpreters.wrapper_py import metadata, Wrapper_py
sage: instr_stream = InstructionStream(metadata, 1)
sage: generate_code(expr, instr_stream)
sage: instr_stream.instr('return')
sage: v = Wrapper_py(instr_stream.get_current())
sage: type(v)
<class 'sage.ext.interpreters.wrapper_py.Wrapper_py'>
sage: v(7)
8*pi + 56

sage.ext.fast_callable.get_builtin_functions()

To handle ExpressionCall, we need to map from Sage and Python functions to opcode names.

This returns a dictionary which is that map.

We delay building builtin_functions to break a circular import between sage.calculus and this file.

EXAMPLES:

sage: from sage.ext.fast_callable import get_builtin_functions
sage: builtins = get_builtin_functions()
sage: sorted(set(builtins.values()))
['abs', 'acos', 'acosh', 'add', 'asin', 'asinh', 'atan', 'atanh', 'ceil', 'cos', 'cosh', 'cot', 'csc', 'div', 'exp', 'floor', 'floordiv', 'inv', 'log', 'mul', 'neg', 'pow', 'sec', 'sin', 'sinh', 'sqrt', 'sub', 'tan', 'tanh']
sage: builtins[sin]
'sin'
sage: builtins[ln]
'log'

sage.ext.fast_callable.op_list(args, metadata)

Given a dictionary with the result of calling get_current on an InstructionStream, and the corresponding interpreter metadata, return a list of the instructions, in a simple somewhat human-readable format.

For debugging only. (That is, it’s probably not a good idea to try to programmatically manipulate the result of this function; the expected use is just to print the returned list to the screen.)

There’s probably no reason to call this directly; if you have a wrapper object, call op_list on it; if you have an InstructionStream object, call current_op_list on it.

EXAMPLES:

sage: from sage.ext.interpreters.wrapper_rdf import metadata
sage: from sage.ext.fast_callable import InstructionStream, op_list
sage: instr_stream = InstructionStream(metadata, 1)
sage: instr_stream.instr('load_arg', 0)
sage: instr_stream.instr('abs')
sage: instr_stream.instr('return')
sage: instr_stream.current_op_list()
[('load_arg', 0), 'abs', 'return']
sage: op_list(instr_stream.get_current(), metadata)
[('load_arg', 0), 'abs', 'return']