# HG changeset patch # User Carl Witty # Date 1235880780 28800 # Node ID 0949e5d7cc4cc9bade2b54973a1fedf2005b49c9 # Parent 5ba8b66a0ae4da7d663e961da9e50f2fcbd2b829 Add fast_callable (rewritten clone of fast_eval) diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c .hgignore --- a/.hgignore Wed Feb 25 22:34:12 2009 -0800 +++ b/.hgignore Sat Feb 28 20:13:00 2009 -0800 @@ -110,6 +110,7 @@ ^c_lib/config.status$ ^c_lib/libtool$ ^c_lib/libcsage.la$ +^sage/ext/interpreters/ ^sage/libs/ntl/misc.h$ ^sage/libs/ntl/ntl.h$ ^sage/rings/integer.h$ diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c module_list.py --- a/module_list.py Wed Feb 25 22:34:12 2009 -0800 +++ b/module_list.py Sat Feb 28 20:13:00 2009 -0800 @@ -151,6 +151,9 @@ ## ################################ + Extension('sage.ext.fast_callable', + sources = ['sage/ext/fast_callable.pyx']), + Extension('sage.ext.fast_eval', sources = ['sage/ext/fast_eval.pyx']), diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/all.py --- a/sage/all.py Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/all.py Sat Feb 28 20:13:00 2009 -0800 @@ -126,6 +126,9 @@ from sage.parallel.all import * +from sage.ext.fast_callable import fast_callable +from sage.ext.fast_eval import fast_float + from copy import copy, deepcopy # The code executed here uses a large amount of Sage components diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/calculus/calculus.py --- a/sage/calculus/calculus.py Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/calculus/calculus.py Sat Feb 28 20:13:00 2009 -0800 @@ -4741,7 +4741,7 @@ return fast_float.fast_float_constant(float(self)) except: raise NotImplementedError # return lambda x: float(self(x)) - + class Symbolic_object(SymbolicExpression): r""" @@ -4768,6 +4768,23 @@ def _fast_float_(self, *vars): return fast_float.fast_float_constant(float(self)) + def _fast_callable_(self, etb): + r""" + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES:: + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['x']) + sage: SR(pi)._fast_callable_(etb) + pi + sage: etb = ExpressionTreeBuilder(vars=['x'], domain=RDF) + sage: SR(pi)._fast_callable_(etb) + 3.14159265359 + """ + return etb.constant(self._obj) + def __float__(self): """ EXAMPLES: @@ -5592,6 +5609,25 @@ fops = [op._fast_float_(*vars) for op in self._operands] return self._operator(*fops) + def _fast_callable_(self, etb): + r""" + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES:: + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['x','y']) + sage: var('x,y') + (x, y) + sage: (x+y)._fast_callable_(etb) + add(v_0, v_1) + sage: (-x)._fast_callable_(etb) + neg(v_0) + """ + fops = [op._fast_callable_(etb) for op in self._operands] + return self._operator(*fops) + def _convert(self, typ): """ Convert self to the given type by converting each of the operands @@ -6168,6 +6204,28 @@ except TypeError: raise ValueError, "free variable: %s" % self._name + def _fast_callable_(self, etb): + r""" + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES:: + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['x','y']) + sage: var('x,y,z') + (x, y, z) + sage: x._fast_callable_(etb) + v_0 + sage: y._fast_callable_(etb) + v_1 + sage: z._fast_callable_(etb) + Traceback (most recent call last): + ... + ValueError: list.index(x): x not in list + """ + return etb.var(self) + def _recursive_sub(self, kwds): # do the replacement if needed if kwds.has_key(self): @@ -6675,6 +6733,21 @@ vars = self.arguments() return self._expr._fast_float_(*vars) + def _fast_callable_(self, etb): + r""" + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES:: + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['x','y']) + sage: g(x) = sin(x) + 2 + sage: g._fast_callable_(etb) + add(sin(v_0), 2) + """ + return self._expr._fast_callable_(etb) + def __float__(self): return float(self._expr) @@ -7225,6 +7298,22 @@ else: return fast_float.fast_float_func(f, g) + def _fast_callable_(self, etb): + r""" + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES:: + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['x','y']) + sage: var('x,y') + (x, y) + sage: sin(sqrt(x+y))._fast_callable_(etb) + sin(sqrt(add(v_0, v_1))) + """ + return etb.call(self._operands[0], etb(self._operands[1])) + def __complex__(self): """ Convert this symbolic composition to a Python complex number. @@ -7506,7 +7595,7 @@ sage: from sage.ext.fast_eval import fast_float sage: fast_float(sin) - + sage: sin._fast_float_() sage: sin._fast_float_()(0) @@ -7537,6 +7626,23 @@ except TypeError: return fast_float.fast_float_func(self, *args) + def _fast_callable_(self, etb): + r""" + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES:: + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['x','y']) + sage: sin._fast_callable_(etb) + sin(v_0) + sage: erf._fast_callable_(etb) + {erf}(v_0) + """ + args = [etb._var_number(n) for n in range(self.number_of_arguments())] + return etb.call(self, *args) + _syms = {} class Function_erf(PrimitiveFunction): diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/ext/fast_callable.pxd --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/ext/fast_callable.pxd Sat Feb 28 20:13:00 2009 -0800 @@ -0,0 +1,3 @@ +cdef class Wrapper: + cdef readonly object _orig_args + cdef readonly object _metadata diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/ext/fast_callable.pyx --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/ext/fast_callable.pyx Sat Feb 28 20:13:00 2009 -0800 @@ -0,0 +1,1907 @@ +r""" +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, \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) +sage: ff(3.5) +36.3992167723104 +sage: ff(RIF(3.5)) +36.39921677231038? + +By default, \function{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 - 20)*(x - 19)*(x - 18)*(x - 17)*(x - 16)*(x - 15)*(x - 14)*(x - 13)*(x - 12)*(x - 11)*(x - 10)*(x - 9)*(x - 8)*(x - 7)*(x - 6)*(x - 5)*(x - 4)*(x - 3)*(x - 2)*(x - 1) +sage: timeit('wilk.subs(x=30)') # random, long time +625 loops, best of 3: 1.46 ms per loop +sage: fc_wilk = fast_callable(wilk) +sage: timeit('fc_wilk(30)') # random, long time +625 loops, best of 3: 10.4 us per loop + +You can specify a particular domain for the evaluation using +\code{domain=}: + +sage: fc_wilk_zz = fast_callable(wilk, 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.6 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, domain=RDF) +sage: timeit('fc_wilk_rdf(30.0)') # random, long time +625 loops, best of 3: 5.13 us per loop + +We also have support for RR: + +sage: fc_wilk_rr = fast_callable(wilk, domain=RR) +sage: timeit('fc_wilk_rr(30.0)') # random, long time +625 loops, best of 3: 12.9 us per loop + +By default, \function{fast_callable} uses the same variable names in the +same order that the \method{__call__} method on its argument would use; +for instance, on symbolic expressions, the variables are used in alphabetical +order. + +sage: var('y,z,x') +(y, z, x) +sage: f = 10*y + 100*z + x +sage: f(1,2,3) +321 +sage: ff = fast_callable(f) +sage: ff(1,2,3) +321 + +However, this can be overridden with \code{vars=}: + +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 \function{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 \function{fast_callable} an expression +with common subexpressions. + +Internally, \function{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 +\function{fast_callable}, using an \class{ExpressionTreeBuilder}. + +sage: from sage.ext.fast_callable import ExpressionTreeBuilder +sage: etb = ExpressionTreeBuilder(vars=('x','y','z')) + +An \class{ExpressionTreeBuilder} has three interesting methods: +\method{constant}, \method{var}, and \method{call}. +All of these methods return \class{ExpressionTree} objects. + +The \method{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 \method{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 \method{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 \function{sin} over \code{RDF}, +the C math library function \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 \function{fast_callable} for your own class (so that +\code{fast_callable(x)} works when \variable{x} is an instance of your +class), implement a method \code{_fast_callable_(self, etb)} for your class. +This method takes an \class{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), domain=RDF) + + sage: f(1) + 1.4142135623730951 + sage: f.op_list() + [('load_arg', 0), ('load_const', 7.0), 'pow', ('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), ('load_const', 7), 'pow', ('load_const', 1), 'add', ('py_call', sqrt, 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) +""" + + +#***************************************************************************** +# Copyright (C) 2008 Robert Bradshaw +# Copyright (C) 2009 Carl Witty +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# This code is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# The full text of the GPL is available at: +# +# http://www.gnu.org/licenses/ +#***************************************************************************** + +# 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 docstring. +#------------------------------ WRONG (for now) docs follow +# The final interesting method of \class{ExpressionTreeBuilder} is +# \method{choice}. This produces conditional expressions, like the C +# \code{COND ? T : F} expression or the Python {T if COND else F}. +# This lets you define piecewise functions using \function{fast_callable}. + +# sage: v4 = etb.choice(v3 >= etb.constant(0), v1, v2) + +# The arguments are \code{(COND, T, F)} (the same order as in C), so the +# above means that if \variable{v3} evaluates to a nonnegative number, +# then \variable{v4} will evaluate to the result of \variable{v1}; +# otherwise, \variable{v4} will evaluate to the result of \variable{v2}. + +# Let's see an example where we see that \function{fast_callable} does not +# evaluate common subexpressions more than once. We'll make a +# \function{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) # not tested + +# Now \variable{ff} does 32 complex arithmetic operations on each call +# (16 additions and 16 multiplications). However, if \code{z*z} produced +# code that evaluated \variable{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 \code{x+1} twice: + +# sage: etb = ExpressionTreeBuilder('x') +# sage: x = etb.var('x') +# sage: (x+1)*(x+1) +# *(+(v_0, 1), +(v_0, 1)) + +# but this code will only evaluate \code{x+1} once: + +# sage: v = x+1; v*v +# *(+(v_0, 1), +(v_0, 1)) +#------------------------------ done with WRONG (for now) docs + + +import operator +from copy import copy +from sage.rings.real_mpfr cimport RealField, RealNumber +from sage.structure.element cimport Element +from sage.rings.all import RDF +from sage.libs.mpfr cimport mpfr_t, mpfr_ptr, mpfr_init2, mpfr_set, GMP_RNDN + +include "stdsage.pxi" + +def fast_callable(x, domain=None, vars=None): + r""" + Given an expression x, compiles it into a form that can be quickly + evaluated, given values for the variables in x. + + x can be an expression object, or anything that has a .variables() + method and a ._fast_callable_() method (this includes SR, univariate + polynomials, and multivariate polynomials). + + 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; in this case, + all arithmetic is done in that parent. If we have a special-purpose + interpreter for that parent (like RDF), domain=RDF will trigger the + use of that interpreter. + + If vars is None, then we will attempt to determine the set of + variables from x; otherwise, we will use the given set. + + EXAMPLES: + sage: var('x') + x + sage: expr = sin(x) + 3*x^2 + sage: f = fast_callable(expr) + sage: f(2) + sin(2) + 12 + sage: f(2.0) + 12.9092974268257 + sage: f_rdf = fast_callable(expr, domain=RDF) + sage: f_rdf(2) + 12.909297426825681 + sage: f = fast_callable(expr, vars=('z','x','y')) + sage: f(1, 2, 3) + sin(2) + 12 + sage: K. = QQ[] + sage: p = K.random_element(6); p + -x^6 - 12*x^5 + 1/2*x^4 - 1/95*x^3 - 1/2*x^2 - 4 + sage: fp = fast_callable(p, domain=RDF) + sage: fp.op_list() + [('load_arg', 0), ('load_const', -1.0), 'mul', ('load_const', -12.0), 'add', ('load_arg', 0), 'mul', ('load_const', 0.5), 'add', ('load_arg', 0), 'mul', ('load_const', -0.0105263157895), 'add', ('load_arg', 0), 'mul', ('load_const', -0.5), 'add', ('load_arg', 0), 'mul', ('load_arg', 0), 'mul', ('load_const', -4.0), 'add', 'return'] + sage: fp(3.14159) + -4594.1618236401764 + sage: K. = QQ[] + sage: p = K.random_element(degree=3, terms=5); p + -x*y^2 - x*z^2 - 6*x^2 - y^2 - 3*x*z + sage: fp = fast_callable(p, domain=RDF) + sage: fp.op_list() + [('load_const', 0.0), ('load_const', -3.0), ('load_arg', 0), ('load_const', 1.0), 'pow', ('load_arg', 2), ('load_const', 1.0), 'pow', 'mul', 'mul', 'add', ('load_const', -1.0), ('load_arg', 0), ('load_const', 1.0), 'pow', ('load_arg', 1), ('load_const', 2.0), 'pow', 'mul', 'mul', 'add', ('load_const', -6.0), ('load_arg', 0), ('load_const', 2.0), 'pow', 'mul', 'add', ('load_const', -1.0), ('load_arg', 1), ('load_const', 2.0), 'pow', 'mul', 'add', ('load_const', -1.0), ('load_arg', 0), ('load_const', 1.0), 'pow', ('load_arg', 2), ('load_const', 2.0), 'pow', 'mul', 'mul', 'add', 'return'] + sage: fp(e, pi, sqrt(2)) + -98.001564033629322 + sage: symbolic_result = p(e, pi, sqrt(2)); symbolic_result + -1*e*pi^2 - pi^2 - 6*e^2 - 3*sqrt(2)*e - 2*e + sage: n(symbolic_result) + -98.0015640336293 + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=('x','y'), domain=RDF) + sage: x = etb.var('x') + sage: y = etb.var('y') + sage: expr = etb.call(sin, x^2 + y); expr + sin(add(pow(v_0, 2.0), v_1)) + sage: fc = fast_callable(expr, domain=RDF) + sage: fc(5, 7) + 0.55142668124169059 + """ + cdef Expression et + if isinstance(x, Expression): + et = x + vars = et._etb._vars + else: + if vars is None or len(vars) == 0: + vars = x.variables() + # XXX This is pretty gross... there should be a "callable_variables" + # method that does all this. + from sage.calculus.calculus import SR + from sage.rings.all import is_PolynomialRing, is_MPolynomialRing + if x.parent() is SR and x.number_of_arguments() > len(vars): + vars = list(vars) + ['EXTRA_VAR%d' % n for n in range(len(vars), x.number_of_arguments())] + if is_PolynomialRing(x.parent()) or is_MPolynomialRing(x.parent()): + vars = x.parent().variable_names() + etb = ExpressionTreeBuilder(vars=vars, domain=domain) + et = x._fast_callable_(etb) + if isinstance(domain, RealField): + import sage.ext.interpreters.wrapper_rr + builder = sage.ext.interpreters.wrapper_rr.Wrapper_rr + + str = InstructionStream(sage.ext.interpreters.wrapper_rr.metadata, + len(vars), + domain) + elif domain == RDF: + import sage.ext.interpreters.wrapper_rdf + builder = sage.ext.interpreters.wrapper_rdf.Wrapper_rdf + str = InstructionStream(sage.ext.interpreters.wrapper_rdf.metadata, + len(vars)) + elif domain is None: + import sage.ext.interpreters.wrapper_py + builder = sage.ext.interpreters.wrapper_py.Wrapper_py + str = InstructionStream(sage.ext.interpreters.wrapper_py.metadata, + len(vars)) + else: + import sage.ext.interpreters.wrapper_el + builder = sage.ext.interpreters.wrapper_el.Wrapper_el + str = InstructionStream(sage.ext.interpreters.wrapper_el.metadata, + len(vars), + domain) + generate_code(et, str) + str.instr('return') + return builder(str.get_current()) + +def function_name(fn): + r""" + Given a function, returns 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}' + """ + builtins = get_builtin_functions() + if fn in builtins: + return builtins[fn] + try: + return "{%s}" % fn.__name__ + except AttributeError: + return "{%r}" % fn + +cdef class ExpressionTreeBuilder: + r""" + 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) + """ + + cdef readonly object _domain + cdef readonly object _vars + + def __init__(self, vars, domain=None): + r""" + Initialize an instance of ExpressionTreeBuilder. Takes + a list or tuple of variable names to use, and also an optional + domain. If a domain is given, then creating an ExpressionConstant + node with the __call__, make, or constant methods will convert + the value into the given domain. + + Note that this is the only effect of the domain parameter. It + is quite possible to use different domains for + ExpressionTreeBuilder and for fast_callable; in that case, + constants will be converted twice (once when building the + Expression, and once when generating code). + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder('x') + sage: etb(3^50) + 717897987691852588770249 + sage: etb = ExpressionTreeBuilder('x', domain=RR) + sage: etb(3^50) + 7.17897987691853e23 + """ + + if isinstance(vars, tuple): + vars = list(vars) + elif not isinstance(vars, list): + vars = [vars] + + vars = map(self._clean_var, vars) + + self._domain = domain + self._vars = vars + + def __call__(self, x): + r""" + Try to convert the given value to an Expression. If it is already + an Expression, just return it. If it has a _fast_callable_ + method, then call the method with self as an argument. Otherwise, + use self.constant() to turn it into a constant. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder('x') + sage: v = etb(3); v, type(v) + (3, ) + sage: v = etb(polygen(QQ)); v, type(v) + (v_0, ) + sage: v is etb(v) + True + """ + if isinstance(x, Expression): + return x + + try: + fc = x._fast_callable_ + except AttributeError: + return self.constant(x) + + return fc(self) + + def _clean_var(self, v): + r""" + Give a variable name, given a variable. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder('x') + sage: var('x') + x + sage: etb._clean_var(x) + 'x' + sage: x = polygen(RR); x + 1.00000000000000*x + sage: etb._clean_var(x) + 'x' + """ + # There should be a better way to do this. (Maybe there is.) + if not PY_TYPE_CHECK(v, str): + v = str(v) + if '*' in v: + v = v[v.index('*')+1:] + return v + + def constant(self, c): + r""" + 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 + """ + if self._domain is not None: + c = self._domain(c) + return ExpressionConstant(self, c) + + def var(self, v): + r""" + Turn the argument into an ExpressionVariable. Looks 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: list.index(x): x not in list + """ + ind = self._vars.index(self._clean_var(v)) + return ExpressionVariable(self, ind) + + def _var_number(self, n): + r""" + Given an integer, return the variable with that index. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=('a','b','c','d')) + sage: etb._var_number(0) + v_0 + sage: etb._var_number(3) + v_3 + sage: etb._var_number(4) + Traceback (most recent call last): + ... + ValueError: Variable number 4 out of range + """ + if 0 <= n < len(self._vars): + return ExpressionVariable(self, n) + raise ValueError, "Variable number %d out of range" % n + + def call(self, fn, *args): + r""" + Construct a call node, given a function and a list of arguments. + The arguments will be converted to Expressions using + ExpressionTreeBuilder.__call__. + + 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)) + """ + return ExpressionCall(self, fn, map(self, args)) + + def choice(self, cond, iftrue, iffalse): + r""" + 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)) + """ + return ExpressionChoice(self, + cond, + self(iftrue), + self(iffalse)) + +# Cache these values, to make expression building a tiny bit faster +# (by skipping the hash-table lookup in the operator module). +cdef op_add = operator.add +cdef op_sub = operator.sub +cdef op_mul = operator.mul +cdef op_div = operator.div +cdef op_floordiv = operator.floordiv +cdef op_pow = operator.pow +cdef op_neg = operator.neg +cdef op_abs = operator.abs +cdef op_inv = operator.inv + +cdef class Expression: + r""" + Represents 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))) + """ + + cdef ExpressionTreeBuilder _etb + + def __init__(self, etb): + r""" + Initialize an Expression. Sets the _etb member. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: v = etb(3); v # indirect doctest + 3 + sage: v._get_etb() is etb + True + """ + self._etb = etb + + def _get_etb(self): + r""" + Returns the ExpressionTreeBuilder used to build a given expression. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: v = etb(3); v + 3 + sage: v._get_etb() is etb + True + """ + return self._etb + + def __add__(s, o): + r""" + Compute a sum of two Expressions. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: x+x + add(v_0, v_0) + sage: x+1 + add(v_0, 1) + sage: 1+x + add(1, v_0) + sage: x.__add__(1) + add(v_0, 1) + sage: x.__radd__(1) + add(1, v_0) + """ + return _expression_binop_helper(s, o, op_add) + + def __sub__(s, o): + r""" + Compute a difference of two Expressions. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: x-x + sub(v_0, v_0) + sage: x-1 + sub(v_0, 1) + sage: 1-x + sub(1, v_0) + sage: x.__sub__(1) + sub(v_0, 1) + sage: x.__rsub__(1) + sub(1, v_0) + """ + return _expression_binop_helper(s, o, op_sub) + + def __mul__(s, o): + r""" + Compute a product of two Expressions. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: x*x + mul(v_0, v_0) + sage: x*1 + mul(v_0, 1) + sage: 1*x + mul(1, v_0) + sage: x.__mul__(1) + mul(v_0, 1) + sage: x.__rmul__(1) + mul(1, v_0) + """ + return _expression_binop_helper(s, o, op_mul) + + def __div__(s, o): + r""" + Compute a quotient of two Expressions. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: x/x + div(v_0, v_0) + sage: x/1 + div(v_0, 1) + sage: 1/x + div(1, v_0) + sage: x.__div__(1) + div(v_0, 1) + sage: x.__rdiv__(1) + div(1, v_0) + """ + return _expression_binop_helper(s, o, op_div) + + def __floordiv__(s, o): + r""" + Compute the floordiv (the floor of the quotient) of two Expressions. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: x//x + floordiv(v_0, v_0) + sage: x//1 + floordiv(v_0, 1) + sage: 1//x + floordiv(1, v_0) + sage: x.__floordiv__(1) + floordiv(v_0, 1) + sage: x.__rfloordiv__(1) + floordiv(1, v_0) + """ + return _expression_binop_helper(s, o, op_floordiv) + + def __pow__(s, o, dummy): + r""" + Compute a power expression from two Expressions. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: x^x + pow(v_0, v_0) + sage: x^1 + pow(v_0, 1) + sage: x.__pow__(1) + pow(v_0, 1) + sage: x.__rpow__(1) + pow(1, v_0) + """ + # XXX There is a performance regression from the original + # fast_float here; it would replace small integer powers with + # multiplication. We can't do this safely until we support + # common subexpression elimination (or at least the dup instruction). + # (Plus, we should consider how strict a semantics we want; + # probably this sort of optimization should be controlled by a + # flag.) + return _expression_binop_helper(s, o, op_pow) + + def __neg__(self): + r""" + Compute the negation of an Expression. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: -x + neg(v_0) + sage: x.__neg__() + neg(v_0) + """ + return ExpressionCall(self._etb, op_neg, [self]) + + def __abs__(self): + r""" + 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) + """ + return ExpressionCall(self._etb, op_abs, [self]) + + def abs(self): + r""" + 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) + """ + return ExpressionCall(self._etb, op_abs, [self]) + + def __invert__(self): + r""" + Compute the inverse of an Expression. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: ~x + inv(v_0) + sage: x.__invert__() + inv(v_0) + """ + return ExpressionCall(self._etb, op_inv, [self]) + + +cdef class ExpressionConstant(Expression): + r""" + 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)) + + """ + + cdef object _value + + def __init__(self, etb, c): + r""" + Initialize an ExpressionConstant. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder, ExpressionConstant + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: etb(3) + 3 + sage: v = ExpressionConstant(etb, 3); v + 3 + sage: v._get_etb() is etb + True + sage: v.value() + 3 + sage: v.value() == 3 + True + """ + Expression.__init__(self, etb) + self._value = c + + def value(self): + r""" + 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 + """ + return self._value + + def __repr__(self): + r""" + Give a string representing this ExpressionConstant. + (We use the repr of its value.) + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: v = etb.constant(pi) + sage: v + pi + sage: repr(v) + 'pi' + sage: v.__repr__() + 'pi' + """ + return repr(self._value) + +cdef class ExpressionVariable(Expression): + r""" + 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)) + + """ + cdef int _variable_index + + def __init__(self, etb, int n): + r""" + Initialize an ExpressionVariable. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder, ExpressionVariable + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: etb(x) + v_0 + sage: v = ExpressionVariable(etb, 0); v + v_0 + sage: v._get_etb() is etb + True + sage: v.variable_index() + 0 + """ + Expression.__init__(self, etb) + self._variable_index = n + + def variable_index(self): + r""" + 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 + """ + return self._variable_index + + def __repr__(self): + r""" + Give a string representing this ExpressionVariable. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: v = etb._var_number(0) + sage: v + v_0 + sage: repr(v) + 'v_0' + sage: v.__repr__() + 'v_0' + """ + # Should we look up the variable name in self._etb, instead? + # I think not.. I like the emphasis that we're totally removed + # from the original expression when we have an Expression. + return "v_%d" % self._variable_index + +cdef class ExpressionCall(Expression): + r""" + 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)) + + """ + cdef object _function + cdef object _arguments + + def __init__(self, etb, fn, args): + r""" + Initialize an ExpressionCall. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder, ExpressionCall + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: etb.call(factorial, x) + {factorial}(v_0) + sage: v = ExpressionCall(etb, factorial, [x]); v + {factorial}(v_0) + sage: v._get_etb() is etb + True + sage: v.function() + factorial + sage: v.arguments() + [v_0] + """ + Expression.__init__(self, etb) + self._function = fn + self._arguments = args + + def function(self): + r""" + 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 + """ + return self._function + + def arguments(self): + r""" + 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] + """ + return copy(self._arguments) + + def __repr__(self): + r""" + Give a string representing this ExpressionCall. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb.var(x) + sage: etb.call(operator.add, x, 1) + add(v_0, 1) + sage: etb.call(factorial, x) + {factorial}(v_0) + sage: v = etb.call(sin, x) + sage: v + sin(v_0) + sage: repr(v) + 'sin(v_0)' + sage: v.__repr__() + 'sin(v_0)' + """ + fn = function_name(self._function) + return '%s(%s)' % (fn, ', '.join(map(repr, self._arguments))) + +cdef class ExpressionChoice(Expression): + r""" + 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)) + """ + + cdef object _cond + cdef object _iftrue + cdef object _iffalse + + def __init__(self, etb, cond, iftrue, iffalse): + r""" + Initialize an ExpressionChoice. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder, ExpressionChoice + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: etb.choice(x, ~x, 0) + (inv(v_0) if v_0 else 0) + sage: v = ExpressionChoice(etb, x, ~x, etb(0)); v + (inv(v_0) if v_0 else 0) + sage: v._get_etb() is etb + True + sage: v.condition() + v_0 + sage: v.if_true() + inv(v_0) + sage: v.if_false() + 0 + """ + Expression.__init__(self, etb) + self._cond = cond + self._iftrue = iftrue + self._iffalse = iffalse + + def condition(self): + r""" + 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 + """ + return self._cond + + def if_true(self): + r""" + 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) + """ + return self._iftrue + + def if_false(self): + r""" + 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 + """ + return self._iffalse + + def __repr__(self): + r""" + Give a string representation for this ExpressionChoice. + (Based on the syntax for Python conditional expressions.) + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: v = etb.choice(x, ~x, 0) + sage: v + (inv(v_0) if v_0 else 0) + sage: repr(v) + '(inv(v_0) if v_0 else 0)' + sage: v.__repr__() + '(inv(v_0) if v_0 else 0)' + """ + return '(%s if %s else %s)' % (repr(self._iftrue), + repr(self._cond), + repr(self._iffalse)) + +cpdef _expression_binop_helper(s, o, op): + r""" + Makes an Expression for (s op o). Either s or o (or both) must already + be an expression. + + EXAMPLES: + sage: from sage.ext.fast_callable import _expression_binop_helper, ExpressionTreeBuilder + sage: var('x,y') + (x, y) + sage: etb = ExpressionTreeBuilder(vars=(x,y)) + sage: x = etb(x) + + Now x is an Expression, but y is not. Still, all the following + cases work. + sage: _expression_binop_helper(x, x, operator.add) + add(v_0, v_0) + sage: _expression_binop_helper(x, y, operator.add) + add(v_0, v_1) + sage: _expression_binop_helper(y, x, operator.add) + add(v_1, v_0) + """ + # The Cython way of handling operator overloading on cdef classes + # (which is inherited from Python) is quite annoying. Inside the + # code for a binary operator, you know that either the first or + # second argument (or both) is a member of your class, but you + # don't know which. + + # If there is an arithmetic operator between an Expression and + # a non-Expression, I want to convert the non-Expression into + # an Expression. But to do that, I need the ExpressionTreeBuilder + # from the Expression. + + cdef Expression self + cdef Expression other + + if not isinstance(o, Expression): + self = s + other = self._etb(o) + elif not isinstance(s, Expression): + other = o + self = other._etb(s) + else: + self = s + other = o + assert self._etb is other._etb + + return ExpressionCall(self._etb, op, [self, other]) + +builtin_functions = None +cpdef get_builtin_functions(): + r""" + 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(list(builtins.values())) + ['abs', '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' + """ + # We delay building builtin_functions to break a circular import + # between sage.calculus and this file. + global builtin_functions + if builtin_functions is not None: + return builtin_functions + builtin_functions = { + operator.add: 'add', + operator.sub: 'sub', + operator.mul: 'mul', + operator.div: 'div', + operator.floordiv: 'floordiv', + operator.abs: 'abs', + operator.neg: 'neg', + operator.inv: 'inv', + operator.pow: 'pow', + } + # not handled: atan2, log2, log10 + import sage.calculus.all as calc_all + for fn in ('sqrt', 'ceil', 'floor', + 'sin', 'cos', 'tan', 'sec', 'csc', 'cot', + 'asin', 'acos', 'atan', 'sinh', 'cosh', 'tanh', + 'asinh', 'acosh', 'atanh', 'exp', 'log'): + builtin_functions[getattr(calc_all, fn)] = fn + builtin_functions[calc_all.abs_symbolic] = 'abs' + return builtin_functions + +cpdef generate_code(Expression expr, stream): + r""" + 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) + + sage: v(7) + 8*(pi + 7) + + TESTS: + sage: def my_sin(x): return sin(x) + sage: def my_norm(x, y): return x*x + y*y + sage: def my_sqrt(x): return sqrt(x, extend=False) + sage: fc = fast_callable(expr, domain=RealField(130)) + sage: fc(0) + 3.1415926535897932384626433832795028842 + sage: fc(1) + 8.2831853071795864769252867665590057684 + sage: fc = fast_callable(expr, domain=RDF) + sage: fc(0) + 3.1415926535897931 + sage: fc(1) + 8.2831853071795862 + sage: fc.op_list() + [('load_arg', 0), ('load_const', pi), 'add', ('load_arg', 0), ('load_const', 1), 'add', 'mul', 'return'] + sage: fc = fast_callable(etb.call(sin, x) + etb.call(sqrt, x), domain=RDF) + sage: fc(1) + 1.8414709848078965 + sage: fc.op_list() + [('load_arg', 0), 'sin', ('load_arg', 0), 'sqrt', 'add', 'return'] + sage: fc = fast_callable(etb.call(sin, x) + etb.call(sqrt, x)) + sage: fc(1) + sin(1) + 1 + sage: fc.op_list() + [('load_arg', 0), ('py_call', sin, 1), ('load_arg', 0), ('py_call', sqrt, 1), 'add', 'return'] + sage: fc = fast_callable(etb.call(my_sin, x), domain=RDF) + sage: fc(3) + 0.14112000805986721 + sage: fc = fast_callable(etb.call(my_sin, x), domain=RealField(100)) + sage: fc(3) + 0.14112000805986722210074480281 + sage: fc.op_list() + [('load_arg', 0), ('py_call', , 1), 'return'] + sage: fc = fast_callable(etb.call(my_sqrt, x), domain=RDF) + sage: fc(3) + 1.7320508075688772 + sage: fc(-3) + Traceback (most recent call last): + ... + ValueError: math domain error + sage: fc = fast_callable(etb.call(my_sqrt, x), domain=RR) + sage: fc(3) + 1.73205080756888 + sage: fc(-3) + Traceback (most recent call last): + ... + ValueError: negative number -3.00000000000000 does not have a square root in the real field + sage: etb2 = ExpressionTreeBuilder(('y','z')) + sage: y = etb2.var('y') + sage: z = etb2.var('z') + sage: fc = fast_callable(etb2.call(sqrt, etb2.call(my_norm, y, z)), domain=RDF) + sage: fc(3, 4) + 5.0 + sage: fc.op_list() + [('load_arg', 0), ('load_arg', 1), ('py_call', , 2), 'sqrt', 'return'] + sage: fc.python_calls() + [] + sage: fc = fast_callable(etb2.call(sqrt, etb2.call(my_norm, y, z)), domain=RR) + sage: fc(3, 4) + 5.00000000000000 + sage: fc = fast_callable(etb2.call(my_norm, y, z), domain=ZZ) + sage: fc(3, 4) + 25 + sage: fc.op_list() + [('load_arg', 0), ('load_arg', 1), ('py_call', , 2), 'return'] + sage: fc = fast_callable(expr) + sage: fc(3.0r) + 4.0*(pi + 3.0) + sage: fc = fast_callable(x+3, domain=ZZ) + sage: fc(4) + 7 + sage: fc = fast_callable(x/3, domain=ZZ) + sage: fc(4) + Traceback (most recent call last): + ... + TypeError: no conversion of this rational to integer + sage: fc(6) + 2 + sage: fc = fast_callable(etb.call(sin, x), domain=ZZ) + sage: fc(0) + 0 + sage: fc(3) + Traceback (most recent call last): + ... + TypeError: unable to convert x (=sin(3)) to an integer + """ + cdef ExpressionConstant econst + cdef ExpressionVariable evar + cdef ExpressionCall ecall + cdef ExpressionChoice echoice + + if isinstance(expr, ExpressionConstant): + econst = expr + stream.load_const(econst._value) + elif isinstance(expr, ExpressionVariable): + evar = expr + stream.load_arg(evar._variable_index) + elif isinstance(expr, ExpressionCall): + ecall = expr + fn = ecall._function + for arg in ecall._arguments: + generate_code(arg, stream) + builtins = get_builtin_functions() + if fn in builtins: + opname = builtins[fn] + if stream.has_instr(opname): + stream.instr(builtin_functions[fn]) + return + if stream.has_instr('py_call'): + stream.instr('py_call', fn, len(ecall._arguments)) + else: + raise ValueError, "Unhandled function %s in generate_code" % fn + else: + raise ValueError, "Unhandled expression kind %s in generate_code" % type(expr) + + +cdef class InstructionStream: + r""" + 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). + """ + + cdef object _metadata + cdef object _instrs + cdef object _bytecode + cdef object _constants + cdef object _constant_locs + cdef object _py_constants + cdef object _py_constant_locs + cdef int _stack_cur_size + cdef int _stack_max_size + cdef int _n_args + cdef object _domain + + def __init__(self, metadata, n_args, domain=None): + r""" + Initialize an InstructionStream. + + INPUTS: + metadata - The metadata_by_opname from a wrapper module + n_args - The number of arguments accessible by the generated code + (this is just passed to the wrapper class) + domain - The domain of interpretation (this is just passed to the + wrapper class) + + 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() + {'domain': None, 'code': [], 'py_constants': [], 'args': 1, 'stack': 0, 'constants': []} + sage: md = instr_stream.get_metadata() + sage: type(md) + + sage: md.by_opname['py_call'] + (CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs']), 3) + sage: md.by_opcode[3] + ('py_call', CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs'])) + """ + self._metadata = metadata + self._instrs = [] + self._bytecode = [] + self._constants = [] + self._constant_locs = {} + self._py_constants = [] + self._py_constant_locs = {} + self._stack_cur_size = 0 + self._stack_max_size = 0 + self._domain = domain + self._n_args = n_args + + def load_const(self, c): + r""" + 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] + """ + self.instr('load_const', c) + + def load_arg(self, n): + r""" + 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)] + """ + self.instr('load_arg', n) + + def has_instr(self, opname): + r""" + 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 + """ + return (opname in self._metadata.by_opname) + + def instr(self, opname, *args): + r""" + 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', , 1), 'abs', 'return'] + """ + cdef int i + + spec, opcode = self._metadata.by_opname[opname] + assert len(spec.parameters) == len(args) + + n_inputs = spec.n_inputs + n_outputs = spec.n_outputs + + self._bytecode.append(opcode) + for i in range(len(args)): + if spec.parameters[i] == 'constants': + # XXX bad for strict-mode floating-point constants + # (doesn't handle signed 0, NaN) + arg = args[i] + if arg in self._constant_locs: + self._bytecode.append(self._constant_locs[arg]) + else: + loc = len(self._constants) + self._constants.append(arg) + self._constant_locs[arg] = loc + self._bytecode.append(loc) + elif spec.parameters[i] == 'args': + self._bytecode.append(args[i]) + elif spec.parameters[i] == 'n_inputs': + self._bytecode.append(args[i]) + n_inputs = args[i] + elif spec.parameters[i] == 'n_outputs': + self._bytecode.append(args[i]) + n_outputs = args[i] + elif spec.parameters[i] == 'py_constants': + arg = args[i] + if arg in self._py_constant_locs: + self._bytecode.append(self._py_constant_locs[arg]) + else: + loc = len(self._py_constants) + self._py_constants.append(arg) + self._py_constant_locs[arg] = loc + self._bytecode.append(loc) + else: + raise ValueError + + self._stack_cur_size -= n_inputs + self._stack_cur_size += n_outputs + self._stack_max_size = max(self._stack_max_size, self._stack_cur_size) + + def get_metadata(self): + r""" + Returns the interpreter metadata being used by the current + InstructionStream. + + (Probably only useful for writing doctests.) + + 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) + + """ + return self._metadata + + def current_op_list(self): + r""" + Returns the list of intructions 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', , 1), 'abs', 'return'] + """ + return op_list(self.get_current(), self._metadata) + + def get_current(self): + r""" + 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() + {'domain': None, 'code': [], 'py_constants': [], 'args': 1, 'stack': 0, 'constants': []} + 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', , 1), 'abs', 'return'] + sage: instr_stream.get_current() + {'domain': None, 'code': [0, 0, 3, 0, 1, 11, 2], 'py_constants': [], 'args': 1, 'stack': 1, 'constants': []} + """ + d = {'args': self._n_args, + 'constants': self._constants, + 'py_constants': self._py_constants, + 'stack': self._stack_max_size, + 'code': self._bytecode, + 'domain': self._domain} + return d + +class InterpreterMetadata(object): + r""" + The interpreter metadata for a fast_callable interpreter. Currently + only consists of a dictionary mapping instruction names to + (CompilerInstrSpec, opcode) pairs, and a list mapping opcodes to + (instruction name, CompilerInstrSpec) pairs. + + See the class docstring for CompilerInstrSpec for more information. + + NOTE: You must not modify the metadata. + """ + + def __init__(self, by_opname, by_opcode): + r""" + Initialize an InterpreterMetadata object. + + EXAMPLES: + sage: from sage.ext.fast_callable import InterpreterMetadata + + Currently we do no error checking or processing, so we can + use this simple test: + sage: metadata = InterpreterMetadata(by_opname='opname dict goes here', by_opcode='opcode list goes here') + sage: metadata.by_opname + 'opname dict goes here' + sage: metadata.by_opcode + 'opcode list goes here' + """ + self.by_opname = by_opname + self.by_opcode = by_opcode + +class CompilerInstrSpec(object): + r""" + Describes 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. + """ + + def __init__(self, n_inputs, n_outputs, parameters): + r""" + Initialize a CompilerInstrSpec. + + EXAMPLES: + sage: from sage.ext.fast_callable import CompilerInstrSpec + sage: CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs']) + CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs']) + """ + self.n_inputs = n_inputs + self.n_outputs = n_outputs + self.parameters = parameters + + def __repr__(self): + r""" + Give a string representation for this CompilerInstrSpec. + + EXAMPLES: + sage: from sage.ext.fast_callable import CompilerInstrSpec + sage: v = CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs']) + sage: v + CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs']) + sage: repr(v) + "CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs'])" + sage: v.__repr__() + "CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs'])" + """ + return "CompilerInstrSpec(%d, %d, %s)" % (self.n_inputs, self.n_outputs, self.parameters) + +def op_list(args, metadata): + r""" + 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'] + """ + ops = [] + code = args['code'] + while len(code): + opcode = code[0] + code = code[1:] + (opname, instr) = metadata.by_opcode[opcode] + if len(instr.parameters): + op = [opname] + for p in instr.parameters: + p_loc = code[0] + code = code[1:] + if p in ('args', 'n_inputs', 'n_outputs'): + op.append(p_loc) + else: + op.append(args[p][p_loc]) + ops.append(tuple(op)) + else: + ops.append(opname) + return ops + + +cdef class Wrapper: + r""" + The parent class for all fast_callable wrappers. Implements shared + behavior (currently only debugging). + """ + + def __init__(self, args, metadata): + r""" + Initialize a Wrapper object. + + 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: v.get_orig_args() + {'domain': None, 'code': [0, 0, 1, 0, 4, 0, 0, 1, 1, 4, 6, 2], 'py_constants': [], 'args': 1, 'stack': 3, 'constants': [pi, 1]} + sage: v.op_list() + [('load_arg', 0), ('load_const', pi), 'add', ('load_arg', 0), ('load_const', 1), 'add', 'mul', 'return'] + """ + + # We only keep the original arguments for debugging (op_list(), etc.); + # is it worth the memory cost? (Note that we may be holding on to + # large objects that could otherwise be garbage collected, for + # instance.) + self._orig_args = args + self._metadata = metadata + + def get_orig_args(self): + r""" + Get the original arguments used when initializing this + wrapper. + + (Probably only useful when writing doctests.) + + EXAMPLES: + sage: fast_callable(sin(x)/x, domain=RDF).get_orig_args() + {'domain': None, 'code': [0, 0, 15, 0, 0, 7, 2], 'py_constants': [], 'args': 1, 'stack': 2, 'constants': []} + """ + return self._orig_args + + def op_list(self): + r""" + Return the list of instructions in this wrapper. + + EXAMPLES: + sage: fast_callable(cos(x)*x, domain=RDF).op_list() + [('load_arg', 0), 'cos', ('load_arg', 0), 'mul', 'return'] + """ + return op_list(self._orig_args, self._metadata) + + def python_calls(self): + r""" + 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)), domain=RDF).python_calls() + [] + sage: fast_callable(abs(sin(factorial(x)))).python_calls() + [factorial, sin] + """ + ops = self.op_list() + py_calls = [] + for op in ops: + if isinstance(op, tuple) and op[0] == 'py_call': + py_calls.append(op[1]) + return py_calls + diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/ext/fast_eval.pyx --- a/sage/ext/fast_eval.pyx Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/ext/fast_eval.pyx Sat Feb 28 20:13:00 2009 -0800 @@ -24,6 +24,12 @@ See the function \code{fast_float(f, *vars)} to create a fast-callable version of f. +NOTE: Sage temporarily has two implementations of this functionality; +one in this file, which will probably be deprecated soon, and one in +fast_callable.pyx. The following instructions are for the old +implementation; you probably want to be looking at fast_callable.pyx +instead. + To provide this interface for a class, implement \code{_fast_float_(self, *vars)}. The basic building blocks are provided by the functions \code{fast_float_constant} (returns a @@ -35,7 +41,7 @@ EXAMPLES: sage: from sage.ext.fast_eval import fast_float - sage: f = fast_float(sqrt(x^7+1), 'x') + sage: f = fast_float(sqrt(x^7+1), 'x', old=True) sage: f(1) 1.4142135623730951 sage: f.op_list() @@ -81,6 +87,9 @@ # http://www.gnu.org/licenses/ #***************************************************************************** +from sage.ext.fast_callable import fast_callable, Wrapper +from sage.rings.all import RDF + include "stdsage.pxi" cdef extern from "Python.h": @@ -1288,16 +1297,21 @@ return FastDoubleFunc('callable', f, *args) -def fast_float(f, *vars): +new_fast_float=True + +def fast_float(f, *vars, old=None): """ - Tries to create a fast float function out of the - input, if possible. + Tries to create a function that evaluates f quickly using + floating-point numbers, if possible. There are two implementations + of fast_float in Sage; by default we use the newer, which is + slightly faster on most tests. On failure, returns the input unchanged. INPUT: f -- an expression vars -- the names of the arguments + old -- use the original algorithm for fast_float EXAMPLES: sage: from sage.ext.fast_eval import fast_float @@ -1315,6 +1329,9 @@ sage: f(1,2) 1.0 """ + if old is None: + old = not new_fast_float + if isinstance(f, (tuple, list)): return tuple([fast_float(x, *vars) for x in f]) @@ -1328,7 +1345,10 @@ vars = vars[:i] + (v,) + vars[i+1:] try: - return f._fast_float_(*vars) + if old: + return f._fast_float_(*vars) + else: + return fast_callable(f, vars=vars, domain=RDF) except AttributeError: pass @@ -1339,7 +1359,7 @@ try: from sage.calculus.calculus import SR - return SR(f)._fast_float_(*vars) + return fast_float(SR(f), *vars) except TypeError: pass @@ -1349,4 +1369,5 @@ return f def is_fast_float(x): - return PY_TYPE_CHECK(x, FastDoubleFunc) + return PY_TYPE_CHECK(x, FastDoubleFunc) or PY_TYPE_CHECK(x, Wrapper) + diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/ext/gen_interpreters.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/ext/gen_interpreters.py Sat Feb 28 20:13:00 2009 -0800 @@ -0,0 +1,3385 @@ +r""" +Generate interpreters for fast_callable + +AUTHORS: + +- Carl Witty + +This file is part of the Sage support for "planned" computations; +that is, computations that are separated into a planning stage and +a plan-execution stage. Here, we generate fast interpreters for plan +executions. + +There are at least two kinds of computations that are often planned in +this fashion. First is arithmetic expression evaluation, where we +take an arbitrary user-specified arithmetic expression and compile it +into a bytecode form for fast interpretation. Second is things like +FFTs and large multiplications, where large problems are split into +multiple smaller problems... we can do the logical "splitting" for a +given size only once, producing a plan which can be reused as often as +we want for different problems of the same size. Currently only +arithmetic expression evaluation is implemented, but other kinds of +planned computations should be easy to add. + +Typically, for arithmetic expressions, we want the storage of +intermediate results to be handled automatically (on a stack); for +FFTs/multiplications/etc., the planner will keep track of intermediate +results itself. + +For arithmetic expression evaluation, we want to have lots of +interpreters (at least one, and possibly several, per +specially-handled type). Also, for any given type, we have many +possible variants of instruction encoding, etc.; some of these could +be handled with conditional compilation, but some are more +complicated. So we end up writing an interpreter generator. + +We want to share as much code as possible across all of these +interpreters, while still maintaining the freedom to make drastic +changes in the interpretation strategy (which may change the +generated code, the calling convention for the interpreter, etc.) + +To make this work, the interpreter back-end is divided into three +parts: + + 1. The interpreter itself, in C or C++. + + 2. The wrapper, which is a Cython object holding the + constants, code, etc., and which actually calls the interpreter. + + 3. The code generator. + +We generate parts 1 and 2. The code generator is table-driven, +and we generate the tables for the code generator. + +There are a lot of techniques for fast interpreters that we do not yet +use; hopefully at least some of these will eventually be implemented: + +- using gcc's "labels as values" extension where available + +- top-of-stack caching + +- superinstructions and/or superoperators + +- static stack caching + +- context threading/subrouting threading + +- selective inlining/dynamic superinstructions + +- automatic replication + +Interpreters may be stack-based or register-based. Recent research +suggests that register-based interpreters are better, but the +researchers are investigating interpreters for entire programming +languages, rather than interpreters for expressions. I suspect +that stack-based expression interpreters may be better. However, +we'll implement both varieties and see what's best. + +The relative costs of stack- and register-based interpreters will +depend on the costs of moving values. For complicated types (like +mpz_t), a register-based interpreter will quite likely be better, +since it will avoid moving values. + +We will NOT support any sort of storage of bytecode; instead, the +code must be re-generated from expression trees in every Sage run. +This means that we can trivially experiment with different styles of +interpreter, or even use quite different interpreters depending on +the architecture, without having to worry about forward and backward +compatibility. +""" +#***************************************************************************** +# Copyright (C) 2009 Carl Witty +# +# Distributed under the terms of the GNU General Public License (GPL) +# http://www.gnu.org/licenses/ +#***************************************************************************** + +from __future__ import with_statement + +import os +import re +from jinja import Environment +from jinja.datastructure import ComplainingUndefined +from collections import defaultdict +from distutils.extension import Extension + +############################## +# This module is used during the Sage buld process, so it should not +# use any other Sage modules. (In particular, it MUST NOT use any +# Cython modules -- they won't be built yet!) +# Also, we have some trivial dependency tracking, where we don't +# rebuild the interpreters if this file hasn't changed; if +# interpreter configuation is split out into a separate file, +# that will have to be changed. +############################## + + +# We share a single jinja environment among all templating in this file. +# We use trim_blocks=True (which means that we ignore white space after +# "%}" jinja command endings), and set undefined_singleton to complain +# if we use an undefined variable. +jinja_env = Environment(trim_blocks=True, + undefined_singleton=ComplainingUndefined) +# Allow 'i' as a shorter alias for the built-in 'indent' filter. +jinja_env.filters['i'] = jinja_env.filters['indent'] + +def indent_lines(n, text): + r""" + INPUTS: + n -- indentation amount + text -- text to indent + + Indents each line in text by n spaces. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import indent_lines + sage: indent_lines(3, "foo") + ' foo' + sage: indent_lines(3, "foo\nbar") + ' foo\n bar' + sage: indent_lines(3, "foo\nbar\n") + ' foo\n bar\n' + """ + lines = text.splitlines(True) + spaces = ' ' * n + return ''.join(spaces + line for line in lines) + +def je(template, **kwargs): + r""" + A convenience method for creating strings with Jinja templates. + The name je stands for "Jinja evaluate". + + The first argument is the template string; remaining keyword + arguments define Jinja variables. + + If the first character in the template string is a newline, it is + removed (this feature is useful when using multi-line templates defined + with triple-quoted strings -- the first line doesn't have to be on + the same line as the quotes, which would screw up the indentation). + + (This is very inefficient, because it recompiles the Jinja + template on each call; don't use it in situations where + performance is important.) + + EXAMPLES: + sage: from sage.ext.gen_interpreters import je + sage: je("{{ a }} > {{ b }} * {{ c }}", a='"a suffusion of yellow"', b=3, c=7) + u'"a suffusion of yellow" > 3 * 7' + """ + if len(template) > 0 and template[0] == '\n': + template = template[1:] + + # It looks like Jinja automatically removes one trailing newline? + if len(template) > 0 and template[-1] == '\n': + template = template + '\n' + + tmpl = jinja_env.from_string(template) + return tmpl.render(kwargs) + +class StorageType(object): + r""" + A StorageType specifies the C types used to deal with values of a + given type. + + We currently support three categories of types. + + First are the "simple" types. These are types where: the + representation is small, functions expect arguments to be passed + by value, and the C/C++ assignment operator works. This would + include built-in C types (long, float, etc.) and small structs + (like gsl_complex). + + Second is 'PyObject*'. This is just like a simple type, except + that we have to incref/decref at appropriate places. + + Third is "auto-reference" types. This is how + GMP/MPIR/MPFR/MPFI/FLINT types work. For these types, functions + expect arguments to be passed by reference, and the C assignment + operator does not do what we want. In addition, they take + advantage of a quirk in C (where arrays are automatically + converted to pointers) to automatically pass arguments by + reference. + + Support for further categories would not be difficult to add (such + as reference-counted types other than PyObject*, or + pass-by-reference types that don't use the GMP auto-reference + trick), if we ever run across a use for them. + """ + + def __init__(self): + r""" + Initialize an instance of StorageType. + + This sets several properties: + + class_member_declarations: + A string giving variable declarations that must be members of any + wrapper class using this type. + + class_member_initializations: + A string initializing the class_member_declarations; will be + inserted into the __init__ method of any wrapper class using this + type. + + local_declarations: + A string giving variable declarations that must be local variables + in Cython methods using this storage type. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.class_member_declarations + '' + sage: ty_double.class_member_initializations + '' + sage: ty_double.local_declarations + '' + sage: ty_mpfr.class_member_declarations + 'cdef RealField domain\n' + sage: ty_mpfr.class_member_initializations + "self.domain = args['domain']\n" + sage: ty_mpfr.local_declarations + 'cdef RealNumber rn\n' + """ + self.class_member_declarations = '' + self.class_member_initializations = '' + self.local_declarations = '' + + def cheap_copies(self): + r""" + Returns True or False, depending on whether this StorageType + supports cheap copies -- whether it is cheap to copy values of + this type from one location to another. This is true for + primitive types, and for types like PyObject* (where you're only + copying a pointer, and possibly changing some reference counts). + It is false for types like mpz_t and mpfr_t, where copying values + can involve arbitrarily much work (including memory allocation). + + The practical effect is that if cheap_copies is True, + instructions with outputs of this type write the results into + local variables, and the results are then copied to their + final locations. If cheap_copies is False, then the addresses + of output locations are passed into the instruction and the + instruction writes outputs directly in the final location. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.cheap_copies() + True + sage: ty_python.cheap_copies() + True + sage: ty_mpfr.cheap_copies() + False + """ + return False + + def python_refcounted(self): + r""" + Says whether this storage type is a Python type, so we need to + use INCREF/DECREF. + + (If we needed to support any non-Python refcounted types, it + might be better to make this object-oriented and have methods + like "generate an incref" and "generate a decref". But as + long as we only support Python, this way is probably simpler.) + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.python_refcounted() + False + sage: ty_python.python_refcounted() + True + """ + return False + + def cython_decl_type(self): + r""" + Gives the Cython type for a single value of this type (as a string). + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.cython_decl_type() + 'double' + sage: ty_python.cython_decl_type() + 'object' + sage: ty_mpfr.cython_decl_type() + 'mpfr_t' + """ + return self.c_decl_type() + + def cython_array_type(self): + r""" + Gives the Cython type for referring to an array of values of + this type (as a string). + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.cython_array_type() + 'double*' + sage: ty_python.cython_array_type() + 'PyObject**' + sage: ty_mpfr.cython_array_type() + 'mpfr_t*' + """ + return self.c_ptr_type() + + def needs_cython_init_clear(self): + r""" + Says whether values/arrays of this type need to be initialized + before use and cleared before the underlying memory is freed. + + (We could remove this method, always call .cython_init() to + generate initialization code, and just let .cython_init() + generate empty code if no initialization is required; that would + generate empty loops, which are ugly and potentially might not + be optimized away.) + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.needs_cython_init_clear() + False + sage: ty_mpfr.needs_cython_init_clear() + True + sage: ty_python.needs_cython_init_clear() + True + """ + return False + + def c_decl_type(self): + r""" + Gives the C type for a single value of this type (as a string). + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.c_decl_type() + 'double' + sage: ty_python.c_decl_type() + 'PyObject*' + sage: ty_mpfr.c_decl_type() + 'mpfr_t' + """ + raise NotImplementedError + + def c_ptr_type(self): + r""" + Gives the C type for a pointer to this type (as a reference to + either a single value or an array) (as a string). + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.c_ptr_type() + 'double*' + sage: ty_python.c_ptr_type() + 'PyObject**' + sage: ty_mpfr.c_ptr_type() + 'mpfr_t*' + """ + return self.c_decl_type() + '*' + + def c_local_type(self): + r""" + Gives the C type used for a value of this type inside an + instruction. For assignable/cheap_copy types, this is the + same as c_decl_type; for auto-reference types, this is the + pointer type. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.c_local_type() + 'double' + sage: ty_python.c_local_type() + 'PyObject*' + sage: ty_mpfr.c_local_type() + 'mpfr_ptr' + """ + raise NotImplementedError + + def assign_c_from_py(self, c, py): + r""" + Given a Cython variable/array reference/etc. of this storage type, + and a Python expression, generate code to assign to the Cython + variable from the Python expression. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.assign_c_from_py('foo', 'bar') + u'foo = bar' + sage: ty_python.assign_c_from_py('foo[i]', 'bar[j]') + u'foo[i] = bar[j]; Py_INCREF(foo[i])' + sage: ty_mpfr.assign_c_from_py('foo', 'bar') + u'rn = self.domain(bar)\nmpfr_set(foo, rn.value, GMP_RNDN)' + """ + return je("{{ c }} = {{ py }}", c=c, py=py) + + def declare_chunk_class_members(self, name): + r""" + Returns a string giving the declarations of the class members + in a wrapper class for a memory chunk with this storage type + and the given name. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.declare_chunk_class_members('args') + u' cdef int _n_args\n cdef mpfr_t* _args\n' + """ + return je(""" +{# XXX Variables here (and everywhere, really) should actually be Py_ssize_t #} + cdef int _n_{{ name }} + cdef {{ self.cython_array_type() }} _{{ name }} +""", self=self, name=name) + + def alloc_chunk_data(self, name, len): + r""" + Returns a string allocating the memory for the class members for + a memory chunk with this storage type and the given name. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: print ty_mpfr.alloc_chunk_data('args', 'MY_LENGTH') + self._n_args = MY_LENGTH + self._args = sage_malloc(sizeof(mpfr_t) * MY_LENGTH) + if self._args == NULL: raise MemoryError + for i in range(MY_LENGTH): + mpfr_init2(self._args[i], self.domain.prec()) + + """ + return je(""" + self._n_{{ name }} = {{ len }} + self._{{ name }} = <{{ self.c_ptr_type() }}>sage_malloc(sizeof({{ self.c_decl_type() }}) * {{ len }}) + if self._{{ name }} == NULL: raise MemoryError +{% if self.needs_cython_init_clear() %} + for i in range({{ len }}): + {{ self.cython_init('self._%s[i]' % name) }} +{% endif %} +""", self=self, name=name, len=len) + + def dealloc_chunk_data(self, name): + r""" + Returns a string to be put in the __dealloc__ method of a + wrapper class using a memory chunk with this storage type, to + deallocate the corresponding class members. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: print ty_double.dealloc_chunk_data('args') + if self._args: + sage_free(self._args) + + sage: print ty_mpfr.dealloc_chunk_data('constants') + if self._constants: + for i in range(self._n_constants): + mpfr_clear(self._constants[i]) + sage_free(self._constants) + + """ + return je(""" + if self._{{ name }}: +{% if self.needs_cython_init_clear() %} + for i in range(self._n_{{ name }}): + {{ self.cython_clear('self._%s[i]' % name) }} +{% endif %} + sage_free(self._{{ name }}) +""", self=self, name=name) + +class StorageTypeAssignable(StorageType): + r""" + StorageTypeAssignable is a subtype of StorageType that deals with + types with cheap copies, like primitive types and PyObject*. + """ + + def __init__(self, ty): + r""" + Initializes the property type (the C/Cython name for this type), + as well as the properties described in the documentation for + StorageType.__init__. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.class_member_declarations + '' + sage: ty_double.class_member_initializations + '' + sage: ty_double.local_declarations + '' + sage: ty_double.type + 'double' + sage: ty_python.type + 'PyObject*' + """ + StorageType.__init__(self) + self.type = ty + + def cheap_copies(self): + r""" + Returns True or False, depending on whether this StorageType + supports cheap copies -- whether it is cheap to copy values of + this type from one location to another. (See StorageType.cheap_copies + for more on this property.) + + Since having cheap copies is essentially the definition of + StorageTypeAssignable, this always returns True. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.cheap_copies() + True + sage: ty_python.cheap_copies() + True + """ + return True + + def c_decl_type(self): + r""" + Gives the C type for a single value of this type (as a string). + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.c_decl_type() + 'double' + sage: ty_python.c_decl_type() + 'PyObject*' + """ + return self.type + + def c_local_type(self): + r""" + Gives the C type used for a value of this type inside an + instruction. For assignable/cheap_copy types, this is the + same as c_decl_type; for auto-reference types, this is the + pointer type. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_double.c_local_type() + 'double' + sage: ty_python.c_local_type() + 'PyObject*' + """ + return self.type + +class StorageTypeSimple(StorageTypeAssignable): + r""" + StorageTypeSimple is a subtype of StorageTypeAssignable that deals + with non-reference-counted types with cheap copies, like primitive + types. As of yet, it has no functionality differences from + StorageTypeAssignable. + """ + pass + +ty_int = StorageTypeSimple('int') +ty_double = StorageTypeSimple('double') + +class StorageTypePython(StorageTypeAssignable): + r""" + StorageTypePython is a subtype of StorageTypeAssignable that deals + with Python objects. + + Just allocating an array full of PyObject* leads to problems, + because the Python garbage collector must be able to get to every + Python object, and it wouldn't know how to get to these arrays. + So we allocate the array as a Python list, but then we immediately + pull the ob_item out of it and deal only with that from then on. + + We often leave these lists with NULL entries. This is safe for + the garbage collector and the deallocator, which is all we care + about; but it would be unsafe to provide Python-level access to + these lists. + + There is one special thing about StorageTypePython: memory that is + used by the interpreter as scratch space (for example, the stack) + must be cleared after each call (so we don't hold on to + potentially-large objects and waste memory). Since we have to do + this anyway, the interpreter gains a tiny bit of speed by assuming + that the scratch space is cleared on entry; for example, when + pushing a value onto the stack, it doesn't bother to XDECREF the + previous value because it's always NULL. + """ + + def __init__(self): + r""" + Initializes the properties described in the documentation + for StorageTypeAssignable.__init__. The type is always + 'PyObject*'. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.class_member_declarations + '' + sage: ty_python.class_member_initializations + '' + sage: ty_python.local_declarations + '' + sage: ty_python.type + 'PyObject*' + """ + StorageTypeAssignable.__init__(self, 'PyObject*') + + def python_refcounted(self): + r""" + Says whether this storage type is a Python type, so we need to + use INCREF/DECREF. + + Returns True. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.python_refcounted() + True + """ + return True + + def cython_decl_type(self): + r""" + Gives the Cython type for a single value of this type (as a string). + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.cython_decl_type() + 'object' + """ + return 'object' + + def declare_chunk_class_members(self, name): + r""" + Returns a string giving the declarations of the class members + in a wrapper class for a memory chunk with this storage type + and the given name. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.declare_chunk_class_members('args') + u' cdef object _list_args\n cdef int _n_args\n cdef PyObject** _args\n' + """ + return je(""" + cdef object _list_{{ name }} + cdef int _n_{{ name }} + cdef {{ self.cython_array_type() }} _{{ name }} +""", self=self, name=name) + + def alloc_chunk_data(self, name, len): + r""" + Returns a string allocating the memory for the class members for + a memory chunk with this storage type and the given name. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: print ty_python.alloc_chunk_data('args', 'MY_LENGTH') + self._n_args = MY_LENGTH + self._list_args = PyList_New(self._n_args) + self._args = (self._list_args).ob_item + + """ + return je(""" + self._n_{{ name }} = {{ len }} + self._list_{{ name }} = PyList_New(self._n_{{ name }}) + self._{{ name }} = (self._list_{{ name }}).ob_item +""", self=self, name=name, len=len) + + def dealloc_chunk_data(self, name): + r""" + Returns a string to be put in the __dealloc__ method of a + wrapper class using a memory chunk with this storage type, to + deallocate the corresponding class members. + + Our array was allocated as a Python list; this means we actually + don't need to do anything to deallocate it. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.dealloc_chunk_data('args') + '' + """ + return '' + + def needs_cython_init_clear(self): + r""" + Says whether values/arrays of this type need to be initialized + before use and cleared before the underlying memory is freed. + + Returns True. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.needs_cython_init_clear() + True + """ + return True + + def assign_c_from_py(self, c, py): + r""" + Given a Cython variable/array reference/etc. of this storage type, + and a Python expression, generate code to assign to the Cython + variable from the Python expression. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.assign_c_from_py('foo[i]', 'bar[j]') + u'foo[i] = bar[j]; Py_INCREF(foo[i])' + """ + return je("""{{ c }} = {{ py }}; Py_INCREF({{ c }})""", + c=c, py=py) + + def cython_init(self, loc): + r""" + Generates code to initialize a variable (or array reference) + holding a PyObject*. Sets it to NULL. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.cython_init('foo[i]') + u'foo[i] = NULL' + """ + return je("{{ loc }} = NULL", loc=loc) + + def cython_clear(self, loc): + r""" + Generates code to clear a variable (or array reference) holding + a PyObject*. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_python.cython_clear('foo[i]') + u'Py_CLEAR(foo[i])' + """ + return je("Py_CLEAR({{ loc }})", loc=loc) + +ty_python = StorageTypePython() + +class StorageTypeAutoReference(StorageType): + r""" + StorageTypeAutoReference is a subtype of StorageType that deals with + types in the style of GMP/MPIR/MPFR/MPFI/FLINT, where copies are + not cheap, functions expect arguments to be passed by reference, + and the API takes advantage of the C quirk where arrays are + automatically converted to pointers to automatically pass + arguments by reference. + """ + + def __init__(self, decl_ty, ref_ty): + r""" + Initializes the properties decl_type and ref_type (the C type + names used when declaring variables and function parameters, + respectively), as well as the properties described in + the documentation for StorageType.__init__. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.class_member_declarations + 'cdef RealField domain\n' + sage: ty_mpfr.class_member_initializations + "self.domain = args['domain']\n" + sage: ty_mpfr.local_declarations + 'cdef RealNumber rn\n' + sage: ty_mpfr.decl_type + 'mpfr_t' + sage: ty_mpfr.ref_type + 'mpfr_ptr' + """ + StorageType.__init__(self) + self.decl_type = decl_ty + self.ref_type = ref_ty + + def c_decl_type(self): + r""" + Gives the C type for a single value of this type (as a string). + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.c_decl_type() + 'mpfr_t' + """ + return self.decl_type + + def c_local_type(self): + r""" + Gives the C type used for a value of this type inside an + instruction. For assignable/cheap_copy types, this is the + same as c_decl_type; for auto-reference types, this is the + pointer type. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.c_local_type() + 'mpfr_ptr' + """ + return self.ref_type + + def needs_cython_init_clear(self): + r""" + Says whether values/arrays of this type need to be initialized + before use and cleared before the underlying memory is freed. + + All known examples of auto-reference types do need a special + initialization call, so this always returns True. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.needs_cython_init_clear() + True + """ + return True + +class StorageTypeMPFR(StorageTypeAutoReference): + r""" + StorageTypeMPFR is a subtype of StorageTypeAutoReference that deals + the MPFR's mpfr_t type. + + For any given program that we're interpreting, ty_mpfr can only + refer to a single precision. An interpreter that needs to use + two precisions of mpfr_t in the same program should instantiate two + separate instances of StorageTypeMPFR. (Interpreters that need + to handle arbitrarily many precisions in the same program are not + handled at all.) + """ + + def __init__(self, id=''): + r""" + Initializes the id property, as well as the properties described + in the documentation for StorageTypeAutoReference.__init__. + + The id property is used if you want to have an interpreter + that handles two instances of StorageTypeMPFR (that is, + handles mpfr_t variables at two different precisions + simultaneously). It's a string that's used to generate + variable names that don't conflict. (The id system has + never actually been used, so bugs probably remain.) + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.class_member_declarations + 'cdef RealField domain\n' + sage: ty_mpfr.class_member_initializations + "self.domain = args['domain']\n" + sage: ty_mpfr.local_declarations + 'cdef RealNumber rn\n' + sage: ty_mpfr.decl_type + 'mpfr_t' + sage: ty_mpfr.ref_type + 'mpfr_ptr' + + TESTS: + sage: ty_mpfr2 = StorageTypeMPFR(id='_the_second') + sage: ty_mpfr2.class_member_declarations + 'cdef RealField domain_the_second\n' + sage: ty_mpfr2.class_member_initializations + "self.domain_the_second = args['domain_the_second']\n" + sage: ty_mpfr2.local_declarations + 'cdef RealNumber rn_the_second\n' + """ + StorageTypeAutoReference.__init__(self, 'mpfr_t', 'mpfr_ptr') + self.id = id + self.class_member_declarations = "cdef RealField domain%s\n" % self.id + self.class_member_initializations = \ + "self.domain%s = args['domain%s']\n" % (self.id, self.id) + self.local_declarations = "cdef RealNumber rn%s\n" % self.id + + def cython_init(self, loc): + r""" + Generates code to initialize an mpfr_t reference (a variable, an + array reference, etc.) + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.cython_init('foo[i]') + u'mpfr_init2(foo[i], self.domain.prec())' + """ + return je("mpfr_init2({{ loc }}, self.domain{{ self.id }}.prec())", + self=self, loc=loc) + + def cython_clear(self, loc): + r""" + Generates code to clear an mpfr_t reference (a variable, an + array reference, etc.) + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.cython_clear('foo[i]') + 'mpfr_clear(foo[i])' + """ + return 'mpfr_clear(%s)' % loc + + def assign_c_from_py(self, c, py): + r""" + Given a Cython variable/array reference/etc. of this storage type, + and a Python expression, generate code to assign to the Cython + variable from the Python expression. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: ty_mpfr.assign_c_from_py('foo[i]', 'bar[j]') + u'rn = self.domain(bar[j])\nmpfr_set(foo[i], rn.value, GMP_RNDN)' + """ + return je(""" +rn{{ self.id }} = self.domain({{ py }}) +mpfr_set({{ c }}, rn.value, GMP_RNDN)""", self=self, c=c, py=py) + +ty_mpfr = StorageTypeMPFR() + +class MemoryChunk(object): + r""" + Memory chunks control allocation, deallocation, iniialization, + etc. of the vectors and objects in the interpreter. Basically, + there is one memory chunk per argument to the C interpreter. + + There are three "generic" varieties of memory chunk: "constants", + "arguments", and "scratch". These are named after their most + common use, but they could be used for other things in some + interpreters. + + All three kinds of chunks are allocated in the wrapper class. + Constants are initialized when the wrapper is constructed; + arguments are initialized in the __call__ method, from the + caller's arguments. "scratch" chunks are not initialized at all; + they are used for scratch storage (often, but not necessarily, for + a stack) in the interpreter. + + Interpreters which need memory chunks that don't fit into these + categories can create new subclasses of MemoryChunk. + """ + + def __init__(self, name, storage_type): + r""" + Initialize an instance of MemoryChunk. + + This sets the properties "name" (the name of this memory chunk; + used in generated variable names, etc.) and "storage_type", + which is a StorageType object. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: mc.name + 'args' + sage: mc.storage_type is ty_mpfr + True + """ + self.name = name + self.storage_type = storage_type + + def __repr__(self): + r""" + Give a string representation of this memory chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: mc + {MC:args} + sage: mc.__repr__() + '{MC:args}' + """ + return '{MC:%s}' % self.name + + def declare_class_members(self): + r""" + Returns a string giving the declarations of the class members + in a wrapper class for this memory chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: mc.declare_class_members() + u' cdef int _n_args\n cdef mpfr_t* _args\n' + """ + return self.storage_type.declare_chunk_class_members(self.name) + + def init_class_members(self): + r""" + Returns a string to be put in the __init__ method of a wrapper + class using this memory chunk, to initialize the corresponding + class members. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: print mc.init_class_members() + count = args['args'] + self._n_args = count + self._args = sage_malloc(sizeof(mpfr_t) * count) + if self._args == NULL: raise MemoryError + for i in range(count): + mpfr_init2(self._args[i], self.domain.prec()) + + """ + return "" + + def dealloc_class_members(self): + r""" + Returns a string to be put in the __dealloc__ method of a wrapper + class using this memory chunk, to deallocate the corresponding + class members. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: print mc.dealloc_class_members() + if self._args: + for i in range(self._n_args): + mpfr_clear(self._args[i]) + sage_free(self._args) + + """ + return "" + + def declare_parameter(self): + r""" + Returns the string to use to declare the interpreter parameter + corresponding to this memory chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: mc.declare_parameter() + 'mpfr_t* args' + """ + return '%s %s' % (self.storage_type.c_ptr_type(), self.name) + + def declare_call_locals(self): + r""" + Returns a string to put in the __call__ method of a wrapper + class using this memory chunk, to allocate local variables. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkRRRetval('retval', ty_mpfr) + sage: mc.declare_call_locals() + u' cdef RealNumber retval = (self.domain)()\n' + """ + return "" + + def pass_argument(self): + r""" + Returns the string to pass the argument corresponding to this + memory chunk to the interpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkConstants('constants', ty_mpfr) + sage: mc.pass_argument() + 'self._constants' + """ + raise NotImplementedError + + def needs_cleanup_on_error(self): + r""" + In an interpreter that can terminate prematurely (due to an + exception from calling Python code, or divide by zero, or + whatever) it will just return at the end of the current instruction, + skipping the rest of the program. Thus, it may still have + values pushed on the stack, etc. + + This method returns True if this memory chunk is modified by the + interpreter and needs some sort of cleanup when an error happens. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkConstants('constants', ty_mpfr) + sage: mc.needs_cleanup_on_error() + False + """ + return False + + def is_stack(self): + r""" + Says whether this memory chunk is a stack. This affects code + generation for instructions using this memory chunk. + + It would be nicer to make this object-oriented somehow, so + that the code generator called MemoryChunk methods instead of + using + if ch.is_stack(): + ... hardcoded stack code + else: + ... hardcoded non-stack code + but that hasn't been done yet. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkScratch('scratch', ty_mpfr) + sage: mc.is_stack() + False + sage: mc = MemoryChunkScratch('stack', ty_mpfr, is_stack=True) + sage: mc.is_stack() + True + """ + return False + + def is_python_refcounted_stack(self): + r""" + Says whether this memory chunk refers to a stack where the entries + need to be INCREF/DECREF'ed. + + It would be nice to make this object-oriented, so that the + code generator called MemoryChunk methods to do the potential + INCREF/DECREF and didn't have to explicitly test + is_python_refcounted_stack. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkScratch('args', ty_python) + sage: mc.is_python_refcounted_stack() + False + sage: mc = MemoryChunkScratch('args', ty_python, is_stack=True) + sage: mc.is_python_refcounted_stack() + True + sage: mc = MemoryChunkScratch('args', ty_mpfr, is_stack=True) + sage: mc.is_python_refcounted_stack() + False + """ + return self.is_stack() and self.storage_type.python_refcounted() + +class MemoryChunkLonglivedArray(MemoryChunk): + r""" + MemoryChunkLonglivedArray is a subtype of MemoryChunk that deals + with memory chunks that are both 1) allocated as class members (rather + than being allocated in __call__) and 2) are arrays. + """ + + def init_class_members(self): + r""" + Returns a string to be put in the __init__ method of a wrapper + class using this memory chunk, to initialize the corresponding + class members. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_double) + sage: print mc.init_class_members() + count = args['args'] + self._n_args = count + self._args = sage_malloc(sizeof(double) * count) + if self._args == NULL: raise MemoryError + + """ + return je(""" + count = args['{{ self.name }}'] +{% print self.storage_type.alloc_chunk_data(self.name, 'count') %} +""", self=self) + + def dealloc_class_members(self): + r""" + Returns a string to be put in the __dealloc__ method of a wrapper + class using this memory chunk, to deallocate the corresponding + class members. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: print mc.dealloc_class_members() + if self._args: + for i in range(self._n_args): + mpfr_clear(self._args[i]) + sage_free(self._args) + + """ + return self.storage_type.dealloc_chunk_data(self.name) + + def pass_argument(self): + r""" + Returns the string to pass the argument corresponding to this + memory chunk to the interpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkConstants('constants', ty_mpfr) + sage: mc.pass_argument() + 'self._constants' + """ + return 'self._%s' % self.name + +class MemoryChunkConstants(MemoryChunkLonglivedArray): + r""" + MemoryChunkConstants is a subtype of MemoryChunkLonglivedArray. + + MemoryChunkConstants chunks have their contents set in the + wrapper's __init__ method (and not changed afterward). + """ + + def init_class_members(self): + r""" + Returns a string to be put in the __init__ method of a wrapper + class using this memory chunk, to initialize the corresponding + class members. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkConstants('constants', ty_mpfr) + sage: print mc.init_class_members() + val = args['constants'] + self._n_constants = len(val) + self._constants = sage_malloc(sizeof(mpfr_t) * len(val)) + if self._constants == NULL: raise MemoryError + for i in range(len(val)): + mpfr_init2(self._constants[i], self.domain.prec()) + for i in range(len(val)): + rn = self.domain(val[i]) + mpfr_set(self._constants[i], rn.value, GMP_RNDN) + + """ + return je(""" + val = args['{{ self.name }}'] +{% print self.storage_type.alloc_chunk_data(self.name, 'len(val)') %} + for i in range(len(val)): + {{ self.storage_type.assign_c_from_py('self._%s[i]' % self.name, 'val[i]') | i(12) }} +""", self=self) + +class MemoryChunkArguments(MemoryChunkLonglivedArray): + r""" + MemoryChunkArguments is a subtype of MemoryChunkLonglivedArray, + for dealing with arguments to the wrapper's __call__ method. + + Currently the __call__ method is declared to take a varargs + *args argument tuple. We assume that the MemoryChunk named 'args' + will deal with that tuple. + """ + + def setup_args(self): + r""" + Handle the arguments of __call__ -- copy them into a pre-allocated + array, ready to pass to the interpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: print mc.setup_args() + cdef mpfr_t* c_args = self._args + cdef int i + for i from 0 <= i < len(args): + rn = self.domain(args[i]) + mpfr_set(self._args[i], rn.value, GMP_RNDN) + + """ + return je(""" +cdef {{ self.storage_type.c_ptr_type() }} c_args = self._args +cdef int i +for i from 0 <= i < len(args): + {{ self.storage_type.assign_c_from_py('self._args[i]', 'args[i]') | i(4) }} +""", self=self) + + def pass_argument(self): + r""" + Returns the string to pass the argument corresponding to this + memory chunk to the interpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkArguments('args', ty_mpfr) + sage: mc.pass_argument() + 'c_args' + """ + return 'c_args' + +class MemoryChunkScratch(MemoryChunkLonglivedArray): + r""" + MemoryChunkScratch is a subtype of MemoryChunkLonglivedArray + for dealing with memory chunks that are allocated in the wrapper, + but only used in the interpreter -- stacks, scratch registers, etc. + + (Currently these are only used as stacks.) + """ + + def __init__(self, name, storage_type, is_stack=False): + r""" + Initialize an instance of MemoryChunkScratch. + + Initializes the _is_stack property, as well as + the properties described in the documentation for + MemoryChunk.__init__. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkScratch('stack', ty_double, is_stack=True) + sage: mc.name + 'stack' + sage: mc.storage_type is ty_double + True + sage: mc._is_stack + True + """ + MemoryChunkLonglivedArray.__init__(self, name, storage_type) + self._is_stack = is_stack + + def is_stack(self): + r""" + Says whether this memory chunk is a stack. This affects code + generation for instructions using this memory chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkScratch('stack', ty_mpfr, is_stack=True) + sage: mc.is_stack() + True + """ + return self._is_stack + + def needs_cleanup_on_error(self): + r""" + In an interpreter that can terminate prematurely (due to an + exception from calling Python code, or divide by zero, or + whatever) it will just return at the end of the current instruction, + skipping the rest of the program. Thus, it may still have + values pushed on the stack, etc. + + This method returns True if this memory chunk is modified by the + interpreter and needs some sort of cleanup when an error happens. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkScratch('registers', ty_python) + sage: mc.needs_cleanup_on_error() + True + """ + return self.storage_type.python_refcounted() + + def handle_cleanup(self): + r""" + Handle the cleanup if the interpreter exits with an error. + + For scratch/stack chunks that hold Python-refcounted values, + we assume that they are filled with NULL on every entry to the + interpreter. If the interpreter exited with an error, it may + have left values in the chunk, so we need to go through + the chunk and Py_CLEAR it. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkScratch('registers', ty_python) + sage: print mc.handle_cleanup() + for i in range(self._n_registers): + Py_CLEAR(self._registers[i]) + """ + # XXX This is a lot slower than it needs to be, because + # we don't have a "cdef int i" in scope here. + return je(""" +for i in range(self._n_{{ self.name }}): + Py_CLEAR(self._{{ self.name }}[i]) +""", self=self) + +class MemoryChunkRRRetval(MemoryChunk): + r""" + A special-purpose memory chunk, for dealing with the return value + of the RR-based interpreter. + """ + + def declare_class_members(self): + r""" + Returns a string giving the declarations of the class members + in a wrapper class for this memory chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkRRRetval('retval', ty_mpfr) + sage: mc.declare_class_members() + '' + """ + return "" + + def declare_call_locals(self): + r""" + Returns a string to put in the __call__ method of a wrapper + class using this memory chunk, to allocate local variables. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkRRRetval('retval', ty_mpfr) + sage: mc.declare_call_locals() + u' cdef RealNumber retval = (self.domain)()\n' + """ + return je(""" + cdef RealNumber {{ self.name }} = (self.domain)() +""", self=self) + + def pass_argument(self): + r""" + Returns the string to pass the argument corresponding to this + memory chunk to the interpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkRRRetval('retval', ty_mpfr) + sage: mc.pass_argument() + u'&retval.value' + """ + return je("""&{{ self.name }}.value""", self=self) + +class MemoryChunkPythonArguments(MemoryChunk): + r""" + A special-purpose memory chunk, for the generic Python-object based + interpreter. Rather than copy the arguments into an array allocated + in the wrapper, we use the PyTupleObject internals and pass the array + that's inside the argument tuple. + """ + + def declare_class_members(self): + r""" + Returns a string giving the declarations of the class members + in a wrapper class for this memory chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPythonArguments('args', ty_python) + """ + return " cdef int _n_%s\n" % self.name + + def init_class_members(self): + r""" + Returns a string to be put in the __init__ method of a wrapper + class using this memory chunk, to initialize the corresponding + class members. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPythonArguments('args', ty_python) + sage: mc.init_class_members() + u" count = args['args']\n self._n_args = count\n" + """ + return je(""" + count = args['{{ self.name }}'] + self._n_args = count +""", self=self) + + def setup_args(self): + r""" + Handle the arguments of __call__. Nothing to do. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPythonArguments('args', ty_python) + sage: mc.setup_args() + '' + """ + return '' + + def pass_argument(self): + r""" + Pass the innards of the argument tuple to the interpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPythonArguments('args', ty_python) + sage: mc.pass_argument() + '(args).ob_item' + """ + return "(args).ob_item" + +class MemoryChunkElementArguments(MemoryChunkPythonArguments): + r""" + A special-purpose memory chunk, for the Python-object based + interpreters that want to process (and perhaps modify) the data. + + We allocate a new list (via the map function) on every call to + hold the modified arguments. That's not strictly necessary -- + we could pre-allocate a list and map into it -- but this lets us + use simpler code for a very-likely-negligible efficiency cost. + (The Element interpreter is going to allocate lots of objects + as it runs, anyway.) + """ + + def setup_args(self): + r""" + Handle the arguments of __call__. Note: This hardcodes + "self._domain". + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkElementArguments('args', ty_python) + sage: mc.setup_args() + 'mapped_args = map(self._domain, args)\n' + """ + return "mapped_args = map(self._domain, args)\n" + + def pass_argument(self): + r""" + Pass the innards of the argument tuple to the interpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkElementArguments('args', ty_python) + sage: mc.pass_argument() + '(mapped_args).ob_item' + """ + return "(mapped_args).ob_item" + +class MemoryChunkPyConstant(MemoryChunk): + r""" + A special-purpose memory chunk, for holding a single Python constant + and passing it to the interpreter as a PyObject*. + """ + + def __init__(self, name): + r""" + Initialize an instance of MemoryChunkPyConstant. + + Always uses the type ty_python. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPyConstant('domain') + sage: mc.name + 'domain' + sage: mc.storage_type is ty_python + True + """ + MemoryChunk.__init__(self, name, ty_python) + + def declare_class_members(self): + r""" + Returns a string giving the declarations of the class members + in a wrapper class for this memory chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPyConstant('domain') + sage: mc.declare_class_members() + u' cdef object _domain\n' + """ + return je(""" + cdef object _{{ self.name }} +""", self=self) + + def init_class_members(self): + r""" + Returns a string to be put in the __init__ method of a wrapper + class using this memory chunk, to initialize the corresponding + class members. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPyConstant('domain') + sage: mc.init_class_members() + u" self._domain = args['domain']\n" + """ + return je(""" + self._{{ self.name }} = args['{{ self.name }}'] +""", self=self) + + def declare_parameter(self): + r""" + Returns the string to use to declare the interpreter parameter + corresponding to this memory chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPyConstant('domain') + sage: mc.declare_parameter() + 'PyObject* domain' + """ + return 'PyObject* %s' % self.name + + def pass_argument(self): + r""" + Returns the string to pass the argument corresponding to this + memory chunk to the interpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc = MemoryChunkPyConstant('domain') + sage: mc.pass_argument() + 'self._domain' + """ + return 'self._%s' % self.name + +def params_gen(**chunks): + r""" + Instructions have a parameter specification that says where they get + their inputs and where their outputs go. Each parameter has + the same form: it is a triple (chunk, addr, len). The chunk says + where the parameter is read from/written to. The addr says which + value in the chunk is used. If the chunk is a stack chunk, then + addr must be null; the value will be read from/written to the top + of the stack. Otherwise, addr must be an integer, or another chunk; + if addr is another chunk, then the next value is read from that chunk + to be the address. + + The len says how many values to read/write. It can be either None + (meaning to read/write only a single value), an integer, or + another chunk; if it is a chunk, then the next value is read from that + chunk to be the len. Note that specifying len changes the types + given to the instruction, so len==None is different than len==1 even + though both mean to use a single value. + + These parameter specifications are cumbersome to write by hand, so + there's also a simple string format for them. This (curried) + function parses the simple string format and produces parameter + specifications. The params_gen function takes keyword arguments + mapping single-character names to memory chunks. The string format + uses these names. The params_gen function returns another function, + that takes two strings and returns a pair of lists of parameter + specifications. + + Each string is the concatenation of arbitrarily many specifications. + Each specification consists of an address and a length. The + address is either a single character naming a stack chunk, + or a string of the form 'A[B]' where A names a non-stack chunk + and B names the code chunk. The length is either empty, or '@n' + for a number n (meaning to use that many arguments), or '@C', where + C is the code chunk. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc_stack = MemoryChunkScratch('stack', ty_double, is_stack=True) + sage: mc_args = MemoryChunkArguments('args', ty_double) + sage: mc_code = MemoryChunkConstants('code', ty_int) + + sage: pg = params_gen(D=mc_code, A=mc_args, S=mc_stack) + sage: pg('S', '') + ([({MC:stack}, None, None)], []) + sage: pg('A[D]', '') + ([({MC:args}, {MC:code}, None)], []) + sage: pg('S@5', '') + ([({MC:stack}, None, 5)], []) + sage: pg('S@D', '') + ([({MC:stack}, None, {MC:code})], []) + sage: pg('A[D]@D', '') + ([({MC:args}, {MC:code}, {MC:code})], []) + sage: pg('SSS@D', 'A[D]S@D') + ([({MC:stack}, None, None), ({MC:stack}, None, None), ({MC:stack}, None, {MC:code})], [({MC:args}, {MC:code}, None), ({MC:stack}, None, {MC:code})]) + """ + + def make_params(s): + p = [] + s = s.strip() + while s: + chunk_code = s[0] + s = s[1:] + chunk = chunks[chunk_code] + addr = None + ch_len = None + if chunk.is_stack(): + pass + else: + m = re.match(r'\[(?:([0-9]+)|([a-zA-Z]))\]', s) + if m.group(1): + addr = int(m.group(1)) + else: + ch = chunks[m.group(2)] + assert ch.storage_type is ty_int + addr = ch + s = s[m.end():].strip() + if len(s) and s[0] == '@': + m = re.match(r'@(?:([0-9]+)|([a-zA-Z]))', s) + if m.group(1): + ch_len = int(m.group(1)) + else: + ch = chunks[m.group(2)] + assert ch.storage_type is ty_int + ch_len = ch + s = s[m.end():].strip() + p.append((chunk, addr, ch_len)) + return p + + def params(s_ins, s_outs): + ins = make_params(s_ins) + outs = make_params(s_outs) + return (ins, outs) + + return params + +def string_of_addr(a): + r""" + An address or a length from a parameter specification may be + either None, an integer, or a MemoryChunk. If the address or + length is an integer or a MemoryChunk, this function will convert + it to a string giving an expression that will evaluate to the correct + address or length. (See the docstring for params_gen for more + information on parameter specifications.) + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: mc_code = MemoryChunkConstants('code', ty_int) + sage: string_of_addr(mc_code) + '*code++' + sage: string_of_addr(42r) + '42' + """ + if isinstance(a, (int, long)): + return str(a) + assert(isinstance(a, MemoryChunk)) + return '*%s++' % a.name + +class InstrSpec(object): + r""" + Each instruction in an interpreter is represented as an InstrSpec. + This contains all the information that we need to generate code + to interpret the instruction; it also is used to build the tables + that fast_callable uses, so this is the nexus point between + users of the interpreter (possibly pure Python) and the + generated C interpreter. + + The underlying instructions are matched to the caller by name. + For instance, fast_callable assumes that if the interpreter has an + instruction named 'cos', then it will take a single argument, + return a single result, and implement the cos() function. + + The print representation of an instruction (which will probably + only be used when doctesting this file) consists of the name, + a simplified stack effect, and the code (truncated if it's long). + The stack effect has two parts, the input and the output, separated + by '->'; the input shows what will be popped from the stack, + the output what will be placed on the stack. Each consists of + a sequence of 'S' and '*' characters, where 'S' refers to a single + argument and '*' refers to a variable number of arguments. + + The code for an instruction is a small snippet of C code. It has + available variables 'i0', 'i1', ..., 'o0', 'o1', ...; one variable + for each input and output; its job is to assign values to the output + variables, based on the values of the input variables. + + Normally, in an interpreter that uses doubles, each of the input + and output variables will be a double. If i0 actually represents + a variable number of arguments, then it will be a pointer to + double instead, and there will be another variable n_i0 giving + the actual number of arguments. + + When instructions refer to auto-reference types, they actually + get a pointer to the data in its original location; it is + not copied into a local variable. Mostly, this makes no difference, + but there is one potential problem to be aware of. It is possible + for an output variable to point to the same object as an input + variable; in fact, this usually will happen when you're working + with the stack. If the instruction maps to a single function call, + then this is fine; the standard auto-reference implementations + (GMP, MPFR, etc.) are careful to allow having the input and output + be the same. But if the instruction maps to multiple function + calls, you may need to use a temporary variable. + + Here's an example of this issue. Suppose you want to make an + instruction that does ``out = a+b*c``. You write code like this: + out = b*c + out = a+out + But out will actually share the same storage as a; so the first line + modifies a, and you actually end up computing 2*(b+c). The fix + is to only write to the output once, at the very end of your + instruction. + + Instructions are also allowed to access memory chunks (other than + the stack and code) directly. They are available as C variables + with the same name as the chunk. This is useful if some type of + memory chunk doesn't fit well with the params_gen interface. + + There are additional reference-counting rules that must be + followed if your interpreter operates on Python objects; these + rules are described in the docstring of the PythonInterpreter + class. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: pg = RDFInterpreter().pg + sage: InstrSpec('add', pg('SS','S'), code='o0 = i0+i1;') + add: SS->S = 'o0 = i0+i1;' + """ + + def __init__(self, name, io, code=None, uses_error_handler=False, handles_own_decref=False): + r""" + Initialize an InstrSpec. + + INPUTS: + name -- the name of the instruction + io -- a pair of lists of parameter specifications for I/O of the + instruction + code -- a string containing a snippet of C code to read + from the input variables and write to the output variables + uses_error_handler -- True if the instruction calls Python + and jumps to error: on a Python error + handles_own_decref -- True if the instruction handles Python + objects and includes its own + reference-counting + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + + sage: pg = RDFInterpreter().pg + sage: InstrSpec('add', pg('SS','S'), code='o0 = i0+i1;') + add: SS->S = 'o0 = i0+i1;' + sage: instr = InstrSpec('py_call', pg('P[D]S@D', 'S'), code=('This is very complicated. ' + 'blah ' * 30)); instr + py_call: *->S = 'This is very compli... blah blah blah ' + sage: instr.name + 'py_call' + sage: instr.inputs + [({MC:py_constants}, {MC:code}, None), ({MC:stack}, None, {MC:code})] + sage: instr.outputs + [({MC:stack}, None, None)] + sage: instr.code + 'This is very complicated. blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah blah ' + sage: instr.parameters + ['py_constants', 'n_inputs'] + sage: instr.n_inputs + 0 + sage: instr.n_outputs + 1 + """ + self.name = name + self.inputs = io[0] + self.outputs = io[1] + self.uses_error_handler = uses_error_handler + self.handles_own_decref = handles_own_decref + if code is not None: + self.code = code + # XXX We assume that there is only one stack + n_inputs = 0 + n_outputs = 0 + in_effect = '' + out_effect = '' + p = [] + for (ch, addr, len) in self.inputs: + if ch.is_stack(): + if len is None: + n_inputs += 1 + in_effect += 'S' + elif isinstance(len, (int, long)): + n_inputs += len + in_effect += 'S%d' % len + else: + p.append('n_inputs') + in_effect += '*' + else: + p.append(ch.name) + for (ch, addr, len) in self.outputs: + if ch.is_stack(): + if len is None: + n_outputs += 1 + out_effect += 'S' + elif isinstance(len, (int, long)): + n_outputs += len + out_effect += 'S%d' % len + else: + p.append('n_outputs') + out_effect += '*' + else: + p.append(ch.name) + self.parameters = p + self.n_inputs = n_inputs + self.n_outputs = n_outputs + self.in_effect = in_effect + self.out_effect = out_effect + + def __repr__(self): + r""" + Produce a string representing a given instruction, consisting + of its name, a brief stack specification, and its code + (possibly abbreviated). + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: pg = RDFInterpreter().pg + sage: InstrSpec('add', pg('SS','S'), code='o0 = i0+i1;') + add: SS->S = 'o0 = i0+i1;' + """ + rcode = repr(self.code) + if len(rcode) > 40: + rcode = rcode[:20] + '...' + rcode[-17:] + return '%s: %s->%s = %s' % \ + (self.name, self.in_effect, self.out_effect, rcode) + +# Now we have a series of helper functions that make it slightly easier +# to create instructions. + +def instr_infix(name, io, op): + r""" + A helper function for creating instructions implemented by + a single infix binary operator. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: pg = RDFInterpreter().pg + sage: instr_infix('mul', pg('SS', 'S'), '*') + mul: SS->S = 'o0 = i0 * i1;' + """ + return InstrSpec(name, io, code='o0 = i0 %s i1;' % op) + +def instr_funcall_2args(name, io, op): + r""" + A helper function for creating instructions implemented by + a two-argument function call. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: pg = RDFInterpreter().pg + sage: instr_funcall_2args('atan2', pg('SS', 'S'), 'atan2') + atan2: SS->S = 'o0 = atan2(i0, i1);' + """ + return InstrSpec(name, io, code='o0 = %s(i0, i1);' % op) + +def instr_unary(name, io, op): + r""" + A helper function for creating instructions with one input + and one output. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: pg = RDFInterpreter().pg + sage: instr_unary('sin', pg('S','S'), 'sin(i0)') + sin: S->S = 'o0 = sin(i0);' + sage: instr_unary('neg', pg('S','S'), '-i0') + neg: S->S = 'o0 = -i0;' + """ + return InstrSpec(name, io, code='o0 = ' + op + ';') + +def instr_funcall_2args_mpfr(name, io, op): + r""" + A helper function for creating MPFR instructions with two inputs + and one output. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: pg = RRInterpreter().pg + sage: instr_funcall_2args_mpfr('add', pg('SS','S'), 'mpfr_add') + add: SS->S = 'mpfr_add(o0, i0, i1, GMP_RNDN);' + """ + return InstrSpec(name, io, code='%s(o0, i0, i1, GMP_RNDN);' % op) + +def instr_funcall_1arg_mpfr(name, io, op): + r""" + A helper function for creating MPFR instructions with one input + and one output. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: pg = RRInterpreter().pg + sage: instr_funcall_1arg_mpfr('exp', pg('S','S'), 'mpfr_exp') + exp: S->S = 'mpfr_exp(o0, i0, GMP_RNDN);' + """ + return InstrSpec(name, io, code='%s(o0, i0, GMP_RNDN);' % op) + +class InterpreterSpec(object): + r""" + Each interpreter to be generated by this module is represented + by an InterpreterSpec. + """ + + def __init__(self): + r""" + Initialize an InterpreterSpec. + + Initializes the following fields: + + header -- a code snippet to go at the top of the C interpreter + source file + pyx_header -- a code snippet to go at the top of the wrapper + class source file + err_return -- a string indicating the value to be returned + in case of a Python exception + mc_code -- a memory chunk to use for the interpreted code + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = RDFInterpreter() + sage: interp.header + '' + sage: interp.pyx_header + '' + sage: interp.err_return + '-1094648009105371' + sage: interp.mc_code + {MC:code} + """ + self.header = '' + self.pyx_header = '' + self.err_return = 'NULL' + self.mc_code = MemoryChunkConstants('code', ty_int) + + def _set_opcodes(self): + r""" + Assign opcodes to the instructions in this interpreter. + + Must be called at the end of __init__ by any subclass of + InterpreterSpec. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = RDFInterpreter() + sage: interp.instr_descs[5].opcode + 5 + """ + for i in range(len(self.instr_descs)): + self.instr_descs[i].opcode = i + + +class StackInterpreter(InterpreterSpec): + r""" + A subclass of InterpreterSpec, specialized for stack-based + interpreters. (Currently all interpreters are stack-based.) + """ + + def __init__(self, type, mc_retval=None): + r""" + Initialize a StackInterpreter. + + INPUTS: + type -- A StorageType; the basic type that this interpreter + operates on + mc_retval -- default None; if not None, a special-purpose + MemoryChunk to use as a return value + + Initializes the fields described in the documentation for + InterpreterSpec.__init__, as well as the following: + + mc_args, mc_constants, mc_stack -- MemoryChunk values + return_type -- the type returned by the C interpreter (None for int, + where 1 means success and 0 means error) + mc_retval -- None, or the MemoryChunk to use as a return value + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: rdf = RDFInterpreter() + sage: rr = RRInterpreter() + sage: rdf.mc_args + {MC:args} + sage: rdf.mc_constants + {MC:constants} + sage: rdf.mc_stack + {MC:stack} + sage: rr.mc_retval + {MC:retval} + sage: rr.return_type is None + True + sage: rdf.return_type.type + 'double' + """ + InterpreterSpec.__init__(self) + self.mc_args = MemoryChunkArguments('args', type) + self.mc_constants = MemoryChunkConstants('constants', type) + self.mc_stack = MemoryChunkScratch('stack', type, is_stack=True) + if isinstance(type, StorageTypeAssignable): + self.return_type = type + else: + self.return_type = None + self.mc_retval = mc_retval + +class RDFInterpreter(StackInterpreter): + r""" + A subclass of StackInterpreter, specifying an interpreter over + machine-floating-point values (C doubles). + """ + + def __init__(self): + r""" + Initialize an RDFInterpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = RDFInterpreter() + sage: interp.name + 'rdf' + sage: interp.mc_py_constants + {MC:py_constants} + sage: interp.chunks + [{MC:args}, {MC:constants}, {MC:py_constants}, {MC:stack}, {MC:code}] + sage: interp.pg('A[D]', 'S') + ([({MC:args}, {MC:code}, None)], [({MC:stack}, None, None)]) + sage: instrs = dict([(ins.name, ins) for ins in interp.instr_descs]) + sage: instrs['add'] + add: SS->S = 'o0 = i0 + i1;' + sage: instrs['py_call'] + py_call: *->S = ' \nPyObject *py_arg...goto error;\n}\n' + """ + + StackInterpreter.__init__(self, ty_double) + self.name = 'rdf' + self.mc_py_constants = MemoryChunkConstants('py_constants', ty_python) + # This is a randomly chosen number. Whenever this number is + # returned, the wrapper has to check whether an exception actually + # happened, so if an expression evaluates to this number execution + # is slightly slower. Hopefully that won't happen too often :) + self.err_return = '-1094648009105371' + self.chunks = [self.mc_args, self.mc_constants, self.mc_py_constants, + self.mc_stack, + self.mc_code] + pg = params_gen(A=self.mc_args, C=self.mc_constants, D=self.mc_code, + S=self.mc_stack, P=self.mc_py_constants) + self.pg = pg + instrs = [ + InstrSpec('load_arg', pg('A[D]', 'S'), + code='o0 = i0;'), + InstrSpec('load_const', pg('C[D]', 'S'), + code='o0 = i0;'), + InstrSpec('return', pg('S', ''), + code='return i0;'), + InstrSpec('py_call', pg('P[D]S@D', 'S'), + uses_error_handler=True, + code=""" +PyObject *py_args = PyTuple_New(n_i1); +if (py_args == NULL) goto error; +int i; +for (i = 0; i < n_i1; i++) { + PyObject *arg = PyFloat_FromDouble(i1[i]); + if (arg == NULL) { + Py_DECREF(py_args); + goto error; + } + PyTuple_SET_ITEM(py_args, i, arg); +} +PyObject *result = PyObject_CallObject(i0, py_args); +Py_DECREF(py_args); +if (result == NULL) goto error; +/* If result is not a float, then this will turn it into a float first. */ +o0 = PyFloat_AsDouble(result); +Py_DECREF(result); +if (o0 == -1 && PyErr_Occurred()) { + goto error; +} +""") + ] + for (name, op) in [('add', '+'), ('sub', '-'), + ('mul', '*'), ('div', '/')]: + instrs.append(instr_infix(name, pg('SS', 'S'), op)) + instrs.append(instr_funcall_2args('pow', pg('SS', 'S'), 'pow')) + for (name, op) in [('neg', '-i0'), ('invert', '1/i0'), + ('abs', 'fabs(i0)')]: + instrs.append(instr_unary(name, pg('S', 'S'), op)) + for name in ['sqrt', 'ceil', 'floor', 'sin', 'cos', 'tan', + 'asin', 'acos', 'atan', 'sinh', 'cosh', 'tanh', + 'asinh', 'acosh', 'atanh', 'exp', 'log']: + instrs.append(instr_unary(name, pg('S', 'S'), "%s(i0)" % name)) + self.instr_descs = instrs + self._set_opcodes() + +class RRInterpreter(StackInterpreter): + r""" + A subclass of StackInterpreter, specifying an interpreter over + MPFR arbitrary-precision floating-point numbers. + """ + + def __init__(self): + r""" + Initialize an RDFInterpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = RRInterpreter() + sage: interp.name + 'rr' + sage: interp.mc_py_constants + {MC:py_constants} + sage: interp.chunks + [{MC:args}, {MC:retval}, {MC:constants}, {MC:py_constants}, {MC:stack}, {MC:code}, {MC:domain}] + sage: interp.pg('A[D]', 'S') + ([({MC:args}, {MC:code}, None)], [({MC:stack}, None, None)]) + sage: instrs = dict([(ins.name, ins) for ins in interp.instr_descs]) + sage: instrs['add'] + add: SS->S = 'mpfr_add(o0, i0, i1, GMP_RNDN);' + sage: instrs['py_call'] + py_call: *->S = '\nif (!rr_py_call_h...goto error;\n}\n' + + That py_call instruction is particularly interesting, and + demonstrates a useful technique to let you use Cython code + in an interpreter. Let's look more closely: + + sage: print instrs['py_call'].code + if (!rr_py_call_helper(domain, i0, n_i1, i1, o0)) { + goto error; + } + + This instruction makes use of the function rr_py_call_helper, + which is declared... + + sage: print interp.header + #include + #include "wrapper_rr.h" + + in wrapper_rr.h. This file is automatically generated by Cython + from wrapper_rr.pyx, which is automatically generated from + this spec; in particular, rr_py_call_helper comes from: + + sage: print interp.pyx_header + from ... + cdef public bint rr_py_call_helper(object domain, object fn, + int n_args, + mpfr_t* args, mpfr_t* retval) except 0: + py_args = [] + cdef int i + cdef RealNumber rn + for i from 0 <= i < n_args: + rn = domain() + mpfr_set(rn.value, args[i], GMP_RNDN) + py_args.append(rn) + cdef RealNumber result = domain(fn(*py_args)) + mpfr_set(retval[0], result.value, GMP_RNDN) + return 1 + + + So instructions where you need to interact with Python can + call back into Cython code fairly easily. + """ + + StackInterpreter.__init__(self, ty_mpfr, mc_retval= MemoryChunkRRRetval('retval', ty_mpfr)) + self.name = 'rr' + self.err_return = '0' + self.mc_py_constants = MemoryChunkConstants('py_constants', ty_python) + self.mc_domain = MemoryChunkPyConstant('domain') + self.chunks = [self.mc_args, self.mc_retval, self.mc_constants, + self.mc_py_constants, + self.mc_stack, self.mc_code, self.mc_domain] + pg = params_gen(A=self.mc_args, C=self.mc_constants, D=self.mc_code, + S=self.mc_stack, + P=self.mc_py_constants) + self.pg = pg + self.header = """ +#include +#include "wrapper_rr.h" +""" + self.pyx_header = """ +from sage.rings.real_mpfr cimport RealField, RealNumber +from sage.libs.mpfr cimport * + +cdef public bint rr_py_call_helper(object domain, object fn, + int n_args, + mpfr_t* args, mpfr_t* retval) except 0: + py_args = [] + cdef int i + cdef RealNumber rn + for i from 0 <= i < n_args: + rn = domain() + mpfr_set(rn.value, args[i], GMP_RNDN) + py_args.append(rn) + cdef RealNumber result = domain(fn(*py_args)) + mpfr_set(retval[0], result.value, GMP_RNDN) + return 1 + +"""[1:] + instrs = [ + InstrSpec('load_arg', pg('A[D]', 'S'), + code='mpfr_set(o0, i0, GMP_RNDN);'), + InstrSpec('load_const', pg('C[D]', 'S'), + code='mpfr_set(o0, i0, GMP_RNDN);'), + InstrSpec('return', pg('S', ''), + code='mpfr_set(retval[0], i0, GMP_RNDN);\nreturn 1;\n'), + InstrSpec('py_call', pg('P[D]S@D', 'S'), + uses_error_handler=True, + code=""" +if (!rr_py_call_helper(domain, i0, n_i1, i1, o0)) { + goto error; +} +""") + ] + for (name, op) in [('add', 'mpfr_add'), ('sub', 'mpfr_sub'), + ('mul', 'mpfr_mul'), ('div', 'mpfr_div'), + ('pow', 'mpfr_pow')]: + instrs.append(instr_funcall_2args_mpfr(name, pg('SS', 'S'), op)) + for name in ['neg', 'abs', + 'log', 'log2', 'log10', + 'exp', 'exp2', 'exp10', + 'cos', 'sin', 'tan', + 'sec', 'csc', 'cot', + 'acos', 'asin', 'atan', + 'cosh', 'sinh', 'tanh', + 'sech', 'csch', 'coth', + 'acosh', 'asinh', 'atanh', + 'log1p', 'expm1', 'eint', + 'gamma', 'lngamma', + 'zeta', 'erf', 'erfc', + 'j0', 'j1', 'y0', 'y1']: + instrs.append(instr_funcall_1arg_mpfr(name, pg('S', 'S'), 'mpfr_' + name)) + # mpfr_ui_div constructs a temporary mpfr_t and then calls mpfr_div; + # it would probably be (slightly) faster to use a permanent copy + # of "one" (on the other hand, the constructed temporary copy is + # on the stack, so it's very likely to be in the cache). + instrs.append(InstrSpec('invert', pg('S', 'S'), + code='mpfr_ui_div(o0, 1, i0, GMP_RNDN);')) + self.instr_descs = instrs + self._set_opcodes() + +class PythonInterpreter(StackInterpreter): + r""" + A subclass of StackInterpreter, specifying an interpreter over + Python objects. + + Let's discuss how the reference-counting works in Python-object + based interpreters. + + There is a simple rule to remember: when executing the code + snippets, the input variables contain borrowed references; + you must fill in the output variables with references you own. + + As an optimization, an instruction may set .handles_own_decref; in + that case, it must decref any input variables that came from the + stack. (Input variables that came from arguments/constants chunks + must NOT be decref'ed!) In addition, with .handles_own_decref, if + any of your input variables are arbitrary-count, then you must + NULL out these variables as you decref them. (Use Py_CLEAR to do + this, unless you understand the documentation of Py_CLEAR and why + it's different than Py_XDECREF followed by assigning NULL.) + + Note that as a tiny optimization, the interpreter always assumes + (and assures) that empty parts of the stack contain NULL, so + it doesn't bother to Py_XDECREF before it pushes onto the stack. + """ + + def __init__(self): + r""" + Initialize a PythonInterpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = PythonInterpreter() + sage: interp.name + 'py' + sage: interp.mc_args + {MC:args} + sage: interp.chunks + [{MC:args}, {MC:constants}, {MC:stack}, {MC:code}] + sage: instrs = dict([(ins.name, ins) for ins in interp.instr_descs]) + sage: instrs['add'] + add: SS->S = 'o0 = PyNumber_Add(i0, i1);' + sage: instrs['py_call'] + py_call: *->S = '\nPyObject *py_args...CREF(py_args);\n' + """ + + StackInterpreter.__init__(self, ty_python) + self.name = 'py' + # StackInterpreter.__init__ gave us a MemoryChunkArguments. + # Override with MemoryChunkPythonArguments. + self.mc_args = MemoryChunkPythonArguments('args', ty_python) + self.chunks = [self.mc_args, self.mc_constants, self.mc_stack, + self.mc_code] + pg = params_gen(A=self.mc_args, C=self.mc_constants, D=self.mc_code, + S=self.mc_stack) + self.pg = pg + self.header = """ +#include +#define CHECK(x) (x != NULL) +""" + instrs = [ + InstrSpec('load_arg', pg('A[D]', 'S'), + code='o0 = i0; Py_INCREF(o0);'), + InstrSpec('load_const', pg('C[D]', 'S'), + code='o0 = i0; Py_INCREF(o0);'), + InstrSpec('return', pg('S', ''), + code='return i0;', + handles_own_decref=True), + InstrSpec('py_call', pg('C[D]S@D', 'S'), + handles_own_decref=True, + code=""" +PyObject *py_args = PyTuple_New(n_i1); +if (py_args == NULL) goto error; +int i; +for (i = 0; i < n_i1; i++) { + PyObject *arg = i1[i]; + PyTuple_SET_ITEM(py_args, i, arg); + i1[i] = NULL; +} +o0 = PyObject_CallObject(i0, py_args); +Py_DECREF(py_args); +""") + ] + for (name, op) in [('add', 'PyNumber_Add'), + ('sub', 'PyNumber_Subtract'), + ('mul', 'PyNumber_Multiply'), + ('div', 'PyNumber_Divide')]: + instrs.append(instr_funcall_2args(name, pg('SS', 'S'), op)) + instrs.append(InstrSpec('pow', pg('SS', 'S'), + code='o0 = PyNumber_Power(i0, i1, Py_None);')) + for (name, op) in [('neg', 'PyNumber_Negative'), + ('invert', 'PyNumber_Invert'), + ('abs', 'PyNumber_Absolute')]: + instrs.append(instr_unary(name, pg('S', 'S'), '%s(i0)'%op)) + self.instr_descs = instrs + self._set_opcodes() + +class ElementInterpreter(PythonInterpreter): + r""" + A subclass of PythonInterpreter, specifying an interpreter over + Sage elements with a particular parent. + + This is very similar to the PythonInterpreter, but after every + instruction, the result is checked to make sure it actually an + element with the correct parent; if not, we attempt to convert it. + + Uses the same instructions (with the same implementation) as + PythonInterpreter. + """ + + def __init__(self): + r""" + Initialize an ElementInterpreter. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = ElementInterpreter() + sage: interp.name + 'el' + sage: interp.mc_args + {MC:args} + sage: interp.chunks + [{MC:args}, {MC:constants}, {MC:stack}, {MC:domain}, {MC:code}] + sage: instrs = dict([(ins.name, ins) for ins in interp.instr_descs]) + sage: instrs['add'] + add: SS->S = 'o0 = PyNumber_Add(i0, i1);' + sage: instrs['py_call'] + py_call: *->S = '\nPyObject *py_args...CREF(py_args);\n' + """ + + PythonInterpreter.__init__(self) + self.name = 'el' + # PythonInterpreter.__init__ gave us a MemoryChunkPythonArguments. + # Override with MemoryChunkElementArguments. + self.mc_args = MemoryChunkElementArguments('args', ty_python) + self.mc_domain_info = MemoryChunkPyConstant('domain') + self.chunks = [self.mc_args, self.mc_constants, self.mc_stack, + self.mc_domain_info, self.mc_code] + self.header = """ +#include +#include "wrapper_el.h" + +#define CHECK(x) do_check(&(x), domain) + +static inline int do_check(PyObject **x, PyObject *domain) { + if (*x == NULL) return 0; + PyObject *new_x = el_check_element(*x, domain); + Py_DECREF(*x); + *x = new_x; + if (*x == NULL) return 0; + return 1; +} +""" + self.pyx_header = """ +from sage.structure.element cimport Element + +cdef public object el_check_element(object v, parent): + cdef Element v_el + + if PY_TYPE_CHECK(v, Element): + v_el = v + if v_el._parent is parent: + return v_el + + return parent(v) + +"""[1:] + +class InterpreterGenerator(object): + r""" + This class takes an InterpreterSpec and generates the corresponding + C interpreter and Cython wrapper. + + See the documentation for methods get_wrapper and get_interpreter + for more information. + """ + + def __init__(self, spec): + r""" + Initialize an InterpreterGenerator. + + INPUTS: + spec -- an InterpreterSpec + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = RDFInterpreter() + sage: gen = InterpreterGenerator(interp) + sage: gen._spec is interp + True + sage: gen.uses_error_handler + False + """ + + self._spec = spec + self.uses_error_handler = False + + def gen_code(self, instr_desc, write): + r""" + Generates code for a single instruction. + + INPUTS: + instr_desc -- an InstrSpec + write -- a Python callable + + This function calls its write parameter successively with + strings; when these strings are concatenated, the result is + the code for the given instruction. + + See the documentation for the get_interpreter method for more + information. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = RDFInterpreter() + sage: gen = InterpreterGenerator(interp) + sage: import cStringIO + sage: buff = cStringIO.StringIO() + sage: instrs = dict([(ins.name, ins) for ins in interp.instr_descs]) + sage: gen.gen_code(instrs['div'], buff.write) + sage: print buff.getvalue() + case 7: /* div */ + { + double i1 = *--stack; + double i0 = *--stack; + double o0; + o0 = i0 / i1; + *stack++ = o0; + } + break; + + """ + + d = instr_desc + w = write + s = self._spec + + if d.uses_error_handler: + self.uses_error_handler = True + + w(je(""" + case {{ d.opcode }}: /* {{ d.name }} */ + { +""", d=d)) + + # If the inputs to an instruction come from the stack, + # then we want to generate code for the inputs in reverse order: + # for instance, the divide instruction, which takes inputs A and B + # and generates A/B, needs to pop B off the stack first. + # On the other hand, if the inputs come from the constant pool, + # then we want to generate code for the inputs in normal order, + # because the addresses in the code stream will be in that order. + # We handle this by running through the inputs in two passes: + # first a forward pass, where we handle non-stack inputs + # (and lengths for stack inputs), and then a reverse pass, + # where we handle stack inputs. + for i in range(len(d.inputs)): + (ch, addr, input_len) = d.inputs[i] + chst = ch.storage_type + if addr is not None: + w(" int ai%d = %s;\n" % (i, string_of_addr(addr))) + if input_len is not None: + w(" int n_i%d = %s;\n" % (i, string_of_addr(input_len))) + if not ch.is_stack(): + if input_len is not None: + w(" %s i%d = %s + ai%d;\n" % + (chst.c_ptr_type(), i, ch.name, i)) + else: + w(" %s i%d = %s[ai%d];\n" % + (chst.c_local_type(), i, ch.name, i)) + + for i in reversed(range(len(d.inputs))): + (ch, addr, input_len) = d.inputs[i] + chst = ch.storage_type + if ch.is_stack(): + if input_len is not None: + w(" %s -= n_i%d;\n" % (ch.name, i)) + w(" %s i%d = %s;\n" % (chst.c_ptr_type(), i, ch.name)) + else: + w(" %s i%d = *--%s;\n" % (chst.c_local_type(), i, ch.name)) + if ch.is_python_refcounted_stack(): + w(" *%s = NULL;\n" % ch.name) + + for i in range(len(d.outputs)): + (ch, addr, output_len) = d.outputs[i] + chst = ch.storage_type + if addr is not None: + w(" int ao%d = %s;\n" % (i, string_of_addr(addr))) + if output_len is not None: + w(" int n_o%d = %s;\n" % (i, string_of_addr(output_len))) + if ch.is_stack(): + w(" %s o%d = %s;\n" % + (chst.c_ptr_type(), i, ch.name)) + w(" %s += n_o%d;\n" % (ch.name, i)) + else: + w(" %s o%d = %s + ao%d;\n" % + (chst.c_ptr_type(), i, ch.name, i)) + + else: + if not chst.cheap_copies(): + if ch.is_stack(): + w(" %s o%d = *%s++;\n" % + (chst.c_local_type(), i, ch.name)) + else: + w(" %s o%d = %s[ao%d];\n" % + (chst.c_local_type(), i, ch.name, i)) + else: + w(" %s o%d;\n" % (chst.c_local_type(), i)) + + w(indent_lines(8, d.code.rstrip('\n') + '\n')) + + stack_offsets = defaultdict(int) + for i in range(len(d.inputs)): + (ch, addr, input_len) = d.inputs[i] + chst = ch.storage_type + if ch.is_python_refcounted_stack() and not d.handles_own_decref: + if input_len is None: + w(" Py_DECREF(i%d);\n" % i) + stack_offsets[ch] += 1 + else: + w(je(""" + int {{ iter }}; + for ({{ iter }} = 0; {{ iter }} < n_i{{ i }}; {{ iter }}++) { + Py_CLEAR(i{{ i }}[{{ iter }}]); + } +""", iter='_interp_iter_%d' % i, i=i)) + + for i in range(len(d.outputs)): + ch = d.outputs[i][0] + chst = ch.storage_type + if chst.python_refcounted(): + # We don't yet support code chunks + # that produce multiple Python values, because of + # the way it complicates error handling. + assert i == 0 + w(" if (!CHECK(o%d)) {\n" % i) + w(" Py_XDECREF(o%d);\n" % i) + w(" goto error;\n") + w(" }\n") + self.uses_error_handler = True + if chst.cheap_copies(): + if ch.is_stack(): + w(" *%s++ = o%d;\n" % (ch.name, i)) + else: + w(" %s[ao%d] = o%d;\n" % (ch.name, i, i)) + + w(je(""" + } + break; +""")) + + def func_header(self, cython=False): + r""" + Generates the function header for the declaration (in the Cython + wrapper) or the definition (in the C interpreter) of the interpreter + function. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = ElementInterpreter() + sage: gen = InterpreterGenerator(interp) + sage: print gen.func_header() + PyObject* interp_el(PyObject** args, + PyObject** constants, + PyObject** stack, + PyObject* domain, + int* code) + sage: print gen.func_header(cython=True) + object interp_el(PyObject** args, + PyObject** constants, + PyObject** stack, + PyObject* domain, + int* code) + """ + s = self._spec + ret_ty = 'bint' if cython else 'int' + if s.return_type: + ret_ty = s.return_type.c_decl_type() + if cython: + ret_ty = s.return_type.cython_decl_type() + return je("""{{ ret_ty }} interp_{{ s.name }}( +{%- for ch in s.chunks %} +{% if not loop.first %}, + {% endif %}{{ ch.declare_parameter() }} +{%- endfor %})""", ret_ty=ret_ty, s=s) + + def write_interpreter(self, write): + r""" + Generate the code for the C interpreter. + + This function calls its write parameter successively with + strings; when these strings are concatenated, the result is + the code for the interpreter. + + See the documentation for the get_interpreter method for more + information. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = RDFInterpreter() + sage: gen = InterpreterGenerator(interp) + sage: import cStringIO + sage: buff = cStringIO.StringIO() + sage: gen.write_interpreter(buff.write) + sage: print buff.getvalue() + /* Automatically generated. Do not edit! */ ... + """ + s = self._spec + w = write + w(je(""" +/* Automatically generated. Do not edit! */ +#include +{% print s.header %} +{{ self.func_header() }} { + while (1) { + switch (*code++) { +""", s=s, self=self, i=indent_lines)) + for instr_desc in s.instr_descs: + self.gen_code(instr_desc, w) + w(je(""" + } + } +{% if self.uses_error_handler %} +error: + return {{ s.err_return }}; +{% endif %} +} + +""", s=s, i=indent_lines, self=self)) + + def write_wrapper(self, write): + r""" + Generate the code for the Cython wrapper. + This function calls its write parameter successively with + strings; when these strings are concatenated, the result is + the code for the interpreter. + + See the documentation for the get_wrapper method for more + information. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: interp = RDFInterpreter() + sage: gen = InterpreterGenerator(interp) + sage: import cStringIO + sage: buff = cStringIO.StringIO() + sage: gen.write_wrapper(buff.write) + sage: print buff.getvalue() + # Automatically generated. Do not edit! ... + """ + s = self._spec + w = write + types = set() + do_cleanup = False + for ch in s.chunks: + if ch.storage_type is not None: + types.add(ch.storage_type) + do_cleanup = do_cleanup or ch.needs_cleanup_on_error() + for ch in s.chunks: + if ch.name == 'args': + arg_ch = ch + + the_call = je(""" + {% if s.return_type %}return {% endif -%} +interp_{{ s.name }}({{ arg_ch.pass_argument() }} +{% for ch in s.chunks[1:] %} + , {{ ch.pass_argument() }} +{% endfor %} + ) +""", s=s, arg_ch=arg_ch) + + w(je(""" +# Automatically generated. Do not edit! + +include "../stdsage.pxi" +from python_object cimport PyObject +cdef extern from "Python.h": + void Py_DECREF(PyObject *o) + void Py_INCREF(PyObject *o) + void Py_CLEAR(PyObject *o) + +cdef extern from "listobject.h": + object PyList_New(Py_ssize_t len) + ctypedef struct PyListObject: + PyObject **ob_item + +cdef extern from "tupleobject.h": + ctypedef struct PyTupleObject: + PyObject **ob_item + +from sage.ext.fast_callable cimport Wrapper +{% print s.pyx_header %} + +cdef extern {{ self.func_header(cython=true) -}} +{% if s.err_return != 'NULL' %} + except? {{ s.err_return -}} +{% endif %} + +cdef class Wrapper_{{ s.name }}(Wrapper): +{% for ty in types %} +{% print indent_lines(4, ty.class_member_declarations) %} +{% endfor %} +{% for ch in s.chunks %} +{% print ch.declare_class_members() %} +{% endfor %} + + def __init__(self, args): + Wrapper.__init__(self, args, metadata) + cdef int i + cdef int count +{% for ty in types %} +{% print indent_lines(8, ty.local_declarations) %} +{% print indent_lines(8, ty.class_member_initializations) %} +{% endfor %} +{% for ch in s.chunks %} +{% print ch.init_class_members() %} +{% endfor %} + + def __dealloc__(self): + cdef int i +{% for ch in s.chunks %} +{% print ch.dealloc_class_members() %} +{% endfor %} + + def __call__(self, *args): + if self._n_args != len(args): raise ValueError +{% for ty in types %} +{% print indent_lines(8, ty.local_declarations) %} +{% endfor %} +{% print indent_lines(8, arg_ch.setup_args()) %} +{% for ch in s.chunks %} +{% print ch.declare_call_locals() %} +{% endfor %} +{% if do_cleanup %} + try: +{% print indent_lines(4, the_call) %} + except: +{% for ch in s.chunks %} +{% if ch.needs_cleanup_on_error() %} +{% print indent_lines(12, ch.handle_cleanup()) %} +{% endif %} +{% endfor %} + raise +{% else %} +{% print the_call %} +{% endif %} +{% if not s.return_type %} + return retval +{% endif %} + +from sage.ext.fast_callable import CompilerInstrSpec, InterpreterMetadata +metadata = InterpreterMetadata(by_opname={ +{% for instr in s.instr_descs %} + '{{ instr.name }}': + (CompilerInstrSpec({{ instr.n_inputs }}, {{ instr.n_outputs }}, {{ instr.parameters }}), {{ instr.opcode }}), +{% endfor %} + }, + by_opcode=[ +{% for instr in s.instr_descs %} + ('{{ instr.name }}', + CompilerInstrSpec({{ instr.n_inputs }}, {{ instr.n_outputs }}, {{ instr.parameters }})), +{% endfor %} + ]) +""", s=s, self=self, types=types, arg_ch=arg_ch, indent_lines=indent_lines, the_call=the_call, do_cleanup=do_cleanup)) + + def get_interpreter(self): + r""" + Returns the code for the C interpreter. + + EXAMPLES: + + First we get the InterpreterSpec for several interpreters: + sage: from sage.ext.gen_interpreters import * + sage: rdf_spec = RDFInterpreter() + sage: rr_spec = RRInterpreter() + sage: el_spec = ElementInterpreter() + + Then we get the actual interpreter code: + sage: rdf_interp = InterpreterGenerator(rdf_spec).get_interpreter() + sage: rr_interp = InterpreterGenerator(rr_spec).get_interpreter() + sage: el_interp = InterpreterGenerator(el_spec).get_interpreter() + + Now we can look through these interpreters. + + Each interpreter starts with a file header; this can be + customized on a per-interpreter basis: + sage: print rr_interp + /* Automatically generated. Do not edit! */ + #include + #include + #include "wrapper_rr.h" + ... + + Next is the function header, with one argument per memory chunk + in the interpreter spec. + sage: print el_interp + /* ... */ ... + PyObject* interp_el(PyObject** args, + PyObject** constants, + PyObject** stack, + PyObject* domain, + int* code) { + ... + + Currently, the interpreters have a very simple structure; just + grab the next instruction and execute it, in a switch + statement. + sage: print rdf_interp + /* ... */ ... + while (1) { + switch (*code++) { + ... + + Then comes the code for each instruction. Here is one of the + simplest instructions: + sage: print rdf_interp + /* ... */ ... + case 9: /* neg */ + { + double i0 = *--stack; + double o0; + o0 = -i0; + *stack++ = o0; + } + break; + ... + + We simply pull the top of the stack into a variable, negate it, + and write the result back onto the stack. + + Let's look at the MPFR-based version of this instruction. + This is an example of an interpreter with an auto-reference + type. + sage: print rr_interp + /* ... */ ... + case 9: /* neg */ + { + mpfr_ptr i0 = *--stack; + mpfr_ptr o0 = *stack++; + mpfr_neg(o0, i0, GMP_RNDN); + } + break; + ... + + Here we see that the input and output variables are actually + just pointers into the stack. But due to the auto-reference + trick, the actual code snippet, ``mpfr_net(o0, i0, GMP_RNDN);``, + is exactly the same as if i0 and o0 were declared as local + mpfr_t variables. + + For completeness, let's look at this instruction in the + Python-object element interpreter. + sage: print el_interp + /* ... */ ... + case 9: /* neg */ + { + PyObject* i0 = *--stack; + *stack = NULL; + PyObject* o0; + o0 = PyNumber_Negative(i0); + Py_DECREF(i0); + if (!CHECK(o0)) { + Py_XDECREF(o0); + goto error; + } + *stack++ = o0; + } + break; + ... + + The original code snippet was only ``o0 = PyNumber_Negative(i0);``; + all the rest is automatically generated. For ElementInterpreter, + the CHECK macro actually checks for an exception (makes sure that + o0 is not NULL), tests if the o0 is an element with the correct + parent, and if not converts it into the correct parent. (That is, + it can potentially modify the variable o0.) + """ + import cStringIO + buff = cStringIO.StringIO() + self.write_interpreter(buff.write) + return buff.getvalue() + + def get_wrapper(self): + r""" + Returns the code for the Cython wrapper. + + EXAMPLES: + + First we get the InterpreterSpec for several interpreters: + sage: from sage.ext.gen_interpreters import * + sage: rdf_spec = RDFInterpreter() + sage: rr_spec = RRInterpreter() + sage: el_spec = ElementInterpreter() + + Then we get the actual wrapper code: + sage: rdf_wrapper = InterpreterGenerator(rdf_spec).get_wrapper() + sage: rr_wrapper = InterpreterGenerator(rr_spec).get_wrapper() + sage: el_wrapper = InterpreterGenerator(el_spec).get_wrapper() + + Now we can look through these wrappers. + + Each wrapper starts with a file header; this can be + customized on a per-interpreter basis (some blank lines have been + elided below): + sage: print rdf_wrapper + # Automatically generated. Do not edit! + include "../stdsage.pxi" + from python_object cimport PyObject + cdef extern from "Python.h": + void Py_DECREF(PyObject *o) + void Py_INCREF(PyObject *o) + void Py_CLEAR(PyObject *o) + cdef extern from "listobject.h": + object PyList_New(Py_ssize_t len) + ctypedef struct PyListObject: + PyObject **ob_item + cdef extern from "tupleobject.h": + ctypedef struct PyTupleObject: + PyObject **ob_item + from sage.ext.fast_callable cimport Wrapper + ... + + Next is the declaration of the C interpreter function. + sage: print rdf_wrapper + # ... + cdef extern double interp_rdf(double* args, + double* constants, + PyObject** py_constants, + double* stack, + int* code) except? -1094648009105371 + ... + + We need a way to propagate exceptions back to the wrapper, + even though we only return a double from interp_rdf. The + ``except? -1094648009105371`` (that's a randomly chosen + number) means that we will return that number if there's an + exception, but the wrapper still has to check whether that's a + legitimate return or an exception. (Cython does this + automatically.) + + Next comes the actual wrapper class, which starts off with + a list of member declarations. + sage: print rdf_wrapper + # ... + cdef class Wrapper_rdf(Wrapper): + cdef int _n_args + cdef double* _args + cdef int _n_constants + cdef double* _constants + cdef object _list_py_constants + cdef int _n_py_constants + cdef PyObject** _py_constants + cdef int _n_stack + cdef double* _stack + cdef int _n_code + cdef int* _code + ... + + Contrast the declaration of ``_stack`` here with the + ElementInterpreter version. To simplify our handling of + reference counting and garbage collection, in a Python-object + based interpreter, we allocate arrays as Python lists, + and then pull the array out of the innards of the list. + sage: print el_wrapper + # ... + cdef object _list_stack + cdef int _n_stack + cdef PyObject** _stack + ... + + Next is the __init__ method, which starts like this: + sage: print rdf_wrapper + # ... + def __init__(self, args): + Wrapper.__init__(self, args, metadata) + cdef int i + cdef int count + ... + + To make it possible to generate code for all expression + interpreters with a single code generator, all wrappers + have the same API. The __init__ method takes a single + argument (here called *args*), which is a dictionary holding + all the information needed to initialize this wrapper. + + We call Wrapper.__init__, which saves a copy of this arguments + object and of the interpreter metadata in the wrapper. (This is + only used for debugging.) + + Now we allocate memory for each memory chunk. (We allocate + the memory here, and reuse it on each call of the + wrapper/interpreter. This is for speed reasons; in a fast + interpreter like RDFInterpreter, there are no memory allocations + involved in a call of the wrapper, except for the ones that + are required by the Python calling convention. Eventually + we will support alternate Cython-only entry points that do + absolutely no memory allocation.) + + Basically the same code is repeated, with minor variations, for + each memory chunk; for brevity, we'll only show the code + for 'constants'. + + sage: print rdf_wrapper + # ... + val = args['constants'] + self._n_constants = len(val) + self._constants = sage_malloc(sizeof(double) * len(val)) + if self._constants == NULL: raise MemoryError + for i in range(len(val)): + self._constants[i] = val[i] + ... + + Recall that _n_constants is an int, and _constants is a + double*. + + The RRInterpreter version is more complicated, because it has to + call mpfr_init. + sage: print rr_wrapper + # ... + cdef RealNumber rn + ... + val = args['constants'] + self._n_constants = len(val) + self._constants = sage_malloc(sizeof(mpfr_t) * len(val)) + if self._constants == NULL: raise MemoryError + for i in range(len(val)): + mpfr_init2(self._constants[i], self.domain.prec()) + for i in range(len(val)): + rn = self.domain(val[i]) + mpfr_set(self._constants[i], rn.value, GMP_RNDN) + ... + + And as we mentioned above, in Python-object based interpreters + we actually allocate the memory as a Python list. + sage: print el_wrapper + # ... + val = args['constants'] + self._n_constants = len(val) + self._list_constants = PyList_New(self._n_constants) + self._constants = (self._list_constants).ob_item + for i in range(len(val)): + self._constants[i] = val[i]; Py_INCREF(self._constants[i]) + ... + + Of course, once we've allocated the memory, we eventually have + to free it. (Again, we'll only look at 'constants'.) + sage: print rdf_wrapper + # ... + def __dealloc__(self): + ... + if self._constants: + sage_free(self._constants) + ... + + The RRInterpreter code is more complicated again because it has + to call mpfr_clear. + sage: print rr_wrapper + # ... + def __dealloc__(self): + cdef int i + ... + if self._constants: + for i in range(self._n_constants): + mpfr_clear(self._constants[i]) + sage_free(self._constants) + ... + + But the ElementInterpreter code is extremely simple -- + it doesn't have to do anything to deallocate constants! + (Since the memory for constants is actually allocated as a + Python list, and Cython knows how to deallocate Python lists.) + + Finally we get to the __call__ method. We grab the arguments + passed by the caller, stuff them in our pre-allocated + argument array, and then call the C interpreter. + sage: print rdf_wrapper + # ... + def __call__(self, *args): + if self._n_args != len(args): raise ValueError + cdef double* c_args = self._args + cdef int i + for i from 0 <= i < len(args): + self._args[i] = args[i] + return interp_rdf(c_args + , self._constants + , self._py_constants + , self._stack + , self._code + ) + ... + + In Python-object based interpreters, the call to the C + interpreter has to be a little more complicated. We don't + want to hold on to Python objects from an old computation by + leaving them referenced from the stack. In normal operation, + the C interpreter clears out the stack as it runs, leaving the + stack totally clear when the interpreter finishes. However, + this doesn't happen if the C interpreter raises an exception. + In that case, we have to clear out any remnants from the stack + in the wrapper. + sage: print el_wrapper + # ... + try: + return interp_el((mapped_args).ob_item + , self._constants + , self._stack + , self._domain + , self._code + ) + except: + for i in range(self._n_stack): + Py_CLEAR(self._stack[i]) + raise + ... + + That's it for the wrapper class. The only thing remaining is + the interpreter metadata. This is the information necessary + for the code generator to map instruction names to opcodes; it + also gives information about stack usage, etc. This is fully + documented at InterpreterMetadata; for now, we'll just show + what it looks like. + + Currently, there are two parts to the metadata; the first maps + instruction names to instruction descriptions. The second one + maps opcodes to instruction descriptions. Note that we don't + use InstrSpec objects here; instead, we use CompilerInstrSpec + objects, which are much simpler and contain only the information + we'll need at runtime. + + First the part that maps instruction names to + (CompilerInstrSpec, opcode) pairs. + + sage: print rdf_wrapper + # ... + from sage.ext.fast_callable import CompilerInstrSpec, InterpreterMetadata + metadata = InterpreterMetadata(by_opname={ + ... + 'return': + (CompilerInstrSpec(1, 0, []), 2), + 'py_call': + (CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs']), 3), + 'add': + (CompilerInstrSpec(2, 1, []), 4), + ... + }, ...) + + There's also a table that maps opcodes to (instruction name, + CompilerInstrSpec) pairs: + sage: print rdf_wrapper + # ... + metadata = InterpreterMetadata(..., by_opcode=[ + ... + ('return', + CompilerInstrSpec(1, 0, [])), + ('py_call', + CompilerInstrSpec(0, 1, ['py_constants', 'n_inputs'])), + ('add', + CompilerInstrSpec(2, 1, [])), + ... + ]) + + And that's it for the wrapper. + """ + import cStringIO + buff = cStringIO.StringIO() + self.write_wrapper(buff.write) + return buff.getvalue() + +def write_if_changed(fn, value): + r""" + Writes value to the file named fn, if value is different than + the current contents. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: def last_modification(fn): return os.stat(fn).st_mtime + sage: fn = tmp_filename('gen_interp') + sage: write_if_changed(fn, 'Hello, world') + sage: t1 = last_modification(fn) + sage: open(fn).read() + 'Hello, world' + sage: sleep(2) # long time + sage: write_if_changed(fn, 'Goodbye, world') + sage: t2 = last_modification(fn) + sage: open(fn).read() + 'Goodbye, world' + sage: sleep(2) # long time + sage: write_if_changed(fn, 'Goodbye, world') + sage: t3 = last_modification(fn) + sage: open(fn).read() + 'Goodbye, world' + sage: t1 == t2 # long time + False + sage: t2 == t3 + True + """ + old_value = None + try: + with open(fn) as file: + old_value = file.read() + except IOError: + pass + + if value != old_value: + # We try to remove the file, in case it exists. This is to + # automatically break hardlinks... see #5350 for motivation. + try: + os.remove(fn) + except OSError: + pass + + with open(fn, 'w') as file: + file.write(value) + +def build_interp(interp_spec, dir): + r""" + Given an InterpreterSpec, writes the C interpreter and the Cython + wrapper. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: testdir = tmp_filename() + sage: os.mkdir(testdir) + sage: rdf_interp = RDFInterpreter() + sage: build_interp(rdf_interp, testdir) + sage: open(testdir + '/interp_rdf.c').readline() + '/* Automatically generated. Do not edit! */\n' + """ + ig = InterpreterGenerator(interp_spec) + interp_fn = '%s/interp_%s.c' % (dir, interp_spec.name) + wrapper_fn = '%s/wrapper_%s.pyx' % (dir, interp_spec.name) + interp = ig.get_interpreter() + wrapper = ig.get_wrapper() + write_if_changed(interp_fn, interp) + write_if_changed(wrapper_fn, wrapper) + +def rebuild(dir): + r""" + Check whether the interpreter and wrapper sources have been written + since the last time this module was changed. If not, write them. + + EXAMPLES: + sage: from sage.ext.gen_interpreters import * + sage: testdir = tmp_filename() + sage: os.mkdir(testdir) + sage: rebuild(testdir) + Building interpreters for fast_callable + sage: rebuild(testdir) + sage: open(testdir + '/wrapper_el.pyx').readline() + '# Automatically generated. Do not edit!\n' + """ + module_mtime = os.stat(__file__).st_mtime + try: + if os.stat(dir + '/timestamp').st_mtime > module_mtime: + # No need to rebuild. + return + except OSError: + pass + + # This line will show up "sage -b" (once per upgrade, not every time + # you run it). + print "Building interpreters for fast_callable" + + interp = RDFInterpreter() + build_interp(interp, dir) + + interp = RRInterpreter() + build_interp(interp, dir) + + interp = PythonInterpreter() + build_interp(interp, dir) + + interp = ElementInterpreter() + build_interp(interp, dir) + + # Do this last, so we don't do it if there's an error above. + with open(dir + '/timestamp', 'w'): + pass + +# This list of modules gets added to the list in module_list.py. +# For now, that's not important -- we could have just put this +# list in module_list.py directly. But eventually, we'll have +# interpreters that are conditionally built (for example, +# interpreters that rely on SSE so they only work on x86), so +# it makes sense to keep the decisions about which interpreters +# to write and which interpreters to build in the same place. +modules = [ + Extension('sage.ext.interpreters.wrapper_rdf', + sources = ['sage/ext/interpreters/wrapper_rdf.pyx', + 'sage/ext/interpreters/interp_rdf.c']), + + Extension('sage.ext.interpreters.wrapper_rr', + sources = ['sage/ext/interpreters/wrapper_rr.pyx', + 'sage/ext/interpreters/interp_rr.c'], + libraries=['mpfr']), + + Extension('sage.ext.interpreters.wrapper_py', + sources = ['sage/ext/interpreters/wrapper_py.pyx', + 'sage/ext/interpreters/interp_py.c']), + + Extension('sage.ext.interpreters.wrapper_el', + sources = ['sage/ext/interpreters/wrapper_el.pyx', + 'sage/ext/interpreters/interp_el.c']), + +] diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/ext/interpreters/__init__.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/ext/interpreters/__init__.py Sat Feb 28 20:13:00 2009 -0800 @@ -0,0 +1,1 @@ + diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/functions/constants.py --- a/sage/functions/constants.py Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/functions/constants.py Sat Feb 28 20:13:00 2009 -0800 @@ -434,6 +434,35 @@ from sage.ext.fast_eval import fast_float_constant return fast_float_constant(self) + def _fast_callable_(self, etb): + r""" + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES:: + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['x']) + sage: pi._fast_callable_(etb) + pi + sage: e._fast_callable_(etb) + e + sage: golden_ratio._fast_callable_(etb) + golden_ratio + sage: log2._fast_callable_(etb) + log2 + sage: etb = ExpressionTreeBuilder(vars=['x'], domain=RealField(150)) + sage: pi._fast_callable_(etb) + 3.1415926535897932384626433832795028841971694 + sage: e._fast_callable_(etb) + 2.7182818284590452353602874713526624977572471 + sage: golden_ratio._fast_callable_(etb) + 1.6180339887498948482045868343656381177203092 + sage: log2._fast_callable_(etb) + 0.69314718055994530941723212145817656807550013 + """ + return etb.constant(self) + def floor(self): """ Returns the floor of self. diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/misc/log.py --- a/sage/misc/log.py Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/misc/log.py Sat Feb 28 20:13:00 2009 -0800 @@ -65,8 +65,6 @@ import latex import misc -import sage.plot.all - from sage.misc.viewer import browser, dvi_viewer offset = 0 diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/misc/mathml.py --- a/sage/misc/mathml.py Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/misc/mathml.py Sat Feb 28 20:13:00 2009 -0800 @@ -27,8 +27,6 @@ import os.path -import sage.plot.all - from misc import tmp_dir def list_function(x): diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/numerical/optimize.py --- a/sage/numerical/optimize.py Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/numerical/optimize.py Sat Feb 28 20:13:00 2009 -0800 @@ -229,15 +229,16 @@ """ from sage.calculus.calculus import SymbolicExpression + from sage.ext.fast_eval import fast_callable import scipy from scipy import optimize if isinstance(func,SymbolicExpression): var_list=func.variables() var_names=map(str,var_list) - fast_f=func._fast_float_(*var_names) + fast_f=fast_callable(func, vars=var_names, domain=RDF) f=lambda p: fast_f(*p) gradient_list=func.gradient() - fast_gradient_functions=[gradient_list[i]._fast_float_(*var_names) for i in xrange(len(gradient_list))] + fast_gradient_functions=[fast_callable(gradient_list[i], vars=var_names, domain=RDF) for i in xrange(len(gradient_list))] gradient=lambda p: scipy.array([ a(*p) for a in fast_gradient_functions]) else: f=func @@ -259,7 +260,7 @@ elif algorithm=="ncg": if isinstance(func,SymbolicExpression): hess=func.hessian() - hess_fast= [ [a._fast_float_(*var_names) for a in row] for row in hess] + hess_fast= [ [fast_callable(a, vars=var_names, domain=RDF) for a in row] for row in hess] hessian=lambda p: [[a(*p) for a in row] for row in hess_fast] hessian_p=lambda p,v: scipy.dot(scipy.array(hessian(p)),v) min= optimize.fmin_ncg(f,map(float,x0),fprime=gradient,fhess=hessian,fhess_p=hessian_p,**args) diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/plot/plot.py --- a/sage/plot/plot.py Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/plot/plot.py Sat Feb 28 20:13:00 2009 -0800 @@ -2792,7 +2792,7 @@ sage: x,y = var('x,y') sage: sage.plot.plot.setup_for_eval_on_grid([x^2 + y^2], (x,0,5), (y,0,pi), 11) - ([], + ([], 0.5, 0.31415926535897931, (0.0, 5.0), @@ -2803,7 +2803,7 @@ :: sage: sage.plot.plot.setup_for_eval_on_grid([x^2+y^2], (x,0,1), (y,-1,1), 1) - ([], + ([], 1.0, 2.0, (0.0, 1.0), diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/rings/polynomial/multi_polynomial.pyx --- a/sage/rings/polynomial/multi_polynomial.pyx Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/rings/polynomial/multi_polynomial.pyx Sat Feb 28 20:13:00 2009 -0800 @@ -225,11 +225,11 @@ TESTS: sage: from sage.ext.fast_eval import fast_float - sage: list(fast_float(K(0))) + sage: list(fast_float(K(0), old=True)) ['push 0.0'] - sage: list(fast_float(K(17))) + sage: list(fast_float(K(17), old=True)) ['push 0.0', 'push 17.0', 'add'] - sage: list(fast_float(y)) + sage: list(fast_float(y, old=True)) ['push 0.0', 'push 1.0', 'load 1', 'mul', 'add'] """ from sage.ext.fast_eval import fast_float_arg, fast_float_constant @@ -248,6 +248,40 @@ expr = expr + monom return expr + def _fast_callable_(self, etb): + """ + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['x','y','z']) + sage: K. = QQ[] + sage: v = K.random_element(degree=3, terms=4); v + -6/5*x*y*z + 2*y*z^2 - x + sage: v._fast_callable_(etb) + add(add(add(0, mul(-6/5, mul(mul(pow(v_0, 1), pow(v_1, 1)), pow(v_2, 1)))), mul(2, mul(pow(v_1, 1), pow(v_2, 2)))), mul(-1, pow(v_0, 1))) + + TESTS: + sage: K. = QQ[] + sage: from sage.ext.fast_eval import fast_float + sage: fast_float(K(0)).op_list() + [('load_const', 0.0), 'return'] + sage: fast_float(K(17)).op_list() + [('load_const', 0.0), ('load_const', 17.0), 'add', 'return'] + sage: fast_float(y).op_list() + [('load_const', 0.0), ('load_const', 1.0), ('load_arg', 1), ('load_const', 1.0), 'pow', 'mul', 'add', 'return'] + """ + my_vars = self.parent().variable_names() + x = [etb.var(v) for v in my_vars] + n = len(x) + + expr = etb.constant(0) + for (m, c) in self.dict().iteritems(): + monom = misc.mul([ x[i]**m[i] for i in range(n) if m[i] != 0], + etb.constant(c)) + expr = expr + monom + return expr def derivative(self, *args): r""" diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c sage/rings/polynomial/polynomial_element.pyx --- a/sage/rings/polynomial/polynomial_element.pyx Wed Feb 25 22:34:12 2009 -0800 +++ b/sage/rings/polynomial/polynomial_element.pyx Sat Feb 28 20:13:00 2009 -0800 @@ -618,6 +618,43 @@ expr *= x return expr + def _fast_callable_(self, etb): + r""" + Given an ExpressionTreeBuilder, return an Expression representing + this value. + + EXAMPLES:: + + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=['t']) + sage: R. = QQ[] + sage: v = R.random_element(6); v + -t^6 - 12*t^5 + 1/2*t^4 - 1/95*t^3 - 1/2*t^2 - 4 + sage: v._fast_callable_(etb) + add(mul(mul(add(mul(add(mul(add(mul(add(mul(v_0, -1), -12), v_0), 1/2), v_0), -1/95), v_0), -1/2), v_0), v_0), -4) + """ + x = etb.var(self.variable_name()) + expr = x + cdef int i, d = self.degree() + coeff = self[d] + # We handle polynomial rings like QQ['x']['y']; that gives us some + # slowdown. Optimize away some of that: + if len(etb._vars) == 1: + # OK, we're in the (very common) univariate case. + coeff_maker = etb.constant + else: + # There may be variables in our coefficients... + coeff_maker = etb.make + if coeff != 1: + expr *= coeff_maker(coeff) + for i from d > i >= 0: + coeff = self[i] + if coeff: + expr += coeff_maker(coeff) + if i > 0: + expr *= x + return expr + cdef int _cmp_c_impl(self, Element other) except -2: """ Compare the two polynomials self and other. diff -r 5ba8b66a0ae4 -r 0949e5d7cc4c setup.py --- a/setup.py Wed Feb 25 22:34:12 2009 -0800 +++ b/setup.py Sat Feb 28 20:13:00 2009 -0800 @@ -8,10 +8,13 @@ ### List of Extensions ### ### Since Sage 3.2 the list of extensions resides in -### module_list.py in the same directory as this file. +### module_list.py in the same directory as this file +### (augmented by the list of interpreters +### generated by sage/ext/gen_interpreters.py) ######################################################### from module_list import ext_modules +import sage.ext.gen_interpreters ######################################################### ### Configuration @@ -72,6 +75,10 @@ if DEVEL: extra_compile_args.append('-ggdb') +# Generate interpreters + +sage.ext.gen_interpreters.rebuild(SAGE_DEVEL + 'sage/sage/ext/interpreters') +ext_modules = ext_modules + sage.ext.gen_interpreters.modules ###################################################################### @@ -562,6 +569,7 @@ 'sage.databases', 'sage.ext', + 'sage.ext.interpreters', 'sage.finance',