# HG changeset patch # User Carl Witty # Date 1237396960 25200 # Node ID a97f9d176f07425f7dc6c9965a931152469b7bdf # Parent c8ebd6134a6ff8bde681161097f56c7892502ab8 Fix bugs (and make minor enhancements) brought up by the reviewers/testers diff -r c8ebd6134a6f -r a97f9d176f07 sage/calculus/calculus.py --- a/sage/calculus/calculus.py Sat Feb 28 20:13:00 2009 -0800 +++ b/sage/calculus/calculus.py Wed Mar 18 10:22:40 2009 -0700 @@ -5624,9 +5624,18 @@ 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) + + TESTS:: + + sage: etb = ExpressionTreeBuilder(vars=['x'], domain=RDF) + sage: (x^7)._fast_callable_(etb) + ipow(v_0, 7) + """ + # This used to convert the operands first. Doing it this way + # instead gives a chance to notice powers with an integer + # exponent before the exponent gets (potentially) converted + # to another type. + return etb.call(self._operator, *self._operands) def _convert(self, typ): """ @@ -6232,7 +6241,7 @@ sage: z._fast_callable_(etb) Traceback (most recent call last): ... - ValueError: list.index(x): x not in list + ValueError: Variable 'z' not found """ return etb.var(self) diff -r c8ebd6134a6f -r a97f9d176f07 sage/ext/fast_callable.pyx --- a/sage/ext/fast_callable.pyx Sat Feb 28 20:13:00 2009 -0800 +++ b/sage/ext/fast_callable.pyx Wed Mar 18 10:22:40 2009 -0700 @@ -32,10 +32,10 @@ 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 +625 loops, best of 3: 1.43 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 +625 loops, best of 3: 9.72 us per loop You can specify a particular domain for the evaluation using \code{domain=}: @@ -59,7 +59,7 @@ 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 +625 loops, best of 3: 15.4 us per loop However, for other types, using domain=D will get a large speedup, because we have special-purpose interpreters for those types. One @@ -70,13 +70,23 @@ 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 +625 loops, best of 3: 7 us per loop + +The domain does not need to be a Sage type; for instance, domain=float +also works. (We actually use the same fast interpreter for domain=float +and domain=RDF; the only difference is that when domain=RDF is used, +the return value is an RDF element, and when domain=float is used, +the return value is a Python float.) + +sage: fc_wilk_float = fast_callable(wilk, domain=float) +sage: timeit('fc_wilk_float(30.0)') # random, long time +625 loops, best of 3: 5.04 us per loop We also have support for RR: sage: fc_wilk_rr = fast_callable(wilk, domain=RR) sage: timeit('fc_wilk_rr(30.0)') # random, long time -625 loops, best of 3: 12.9 us per loop +625 loops, best of 3: 13 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; @@ -165,12 +175,12 @@ EXAMPLES: sage: var('x') x - sage: f = fast_callable(sqrt(x^7+1), domain=RDF) + sage: f = fast_callable(sqrt(x^7+1), domain=float) sage: f(1) 1.4142135623730951 sage: f.op_list() - [('load_arg', 0), ('load_const', 7.0), 'pow', ('load_const', 1.0), 'add', 'sqrt', 'return'] + [('load_arg', 0), ('ipow', 7), ('load_const', 1.0), 'add', 'sqrt', 'return'] To interpret that last line, we load argument 0 ('x' in this case) onto the stack, push the constant 7.0 onto the stack, call the pow function @@ -185,7 +195,7 @@ 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'] + [('load_arg', 0), ('ipow', 7), ('load_const', 1), 'add', ('py_call', sqrt, 1), ('load_arg', 0), ('py_call', sin, 1), 'add', 'return'] AUTHOR: @@ -265,6 +275,8 @@ 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 +from sage.rings.integer import Integer +from sage.rings.integer_ring import ZZ include "stdsage.pxi" @@ -277,12 +289,12 @@ 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. + By default, x is evaluated the same way that a Python function + would evaluate it -- addition maps to PyNumber_Add, etc. However, + you can specify domain=D where D is some Sage parent or Python + type; in this case, all arithmetic is done in that domain. If we + have a special-purpose interpreter for that parent (like RDF or float), + domain=... will trigger the use of that interpreter. If vars is None, then we will attempt to determine the set of variables from x; otherwise, we will use the given set. @@ -296,9 +308,19 @@ sin(2) + 12 sage: f(2.0) 12.9092974268257 + + We have special fast interpreters for domain=float and domain=RDF. + (Actually it's the same interpreter; only the return type varies.) + Note that the float interpreter is not actually more accurate than + the RDF interpreter; elements of RDF just don't display all + their digits. + + sage: f_float = fast_callable(expr, domain=float) + sage: f_float(2) + 12.909297426825681 sage: f_rdf = fast_callable(expr, domain=RDF) sage: f_rdf(2) - 12.909297426825681 + 12.9092974268 sage: f = fast_callable(expr, vars=('z','x','y')) sage: f(1, 2, 3) sin(2) + 12 @@ -309,27 +331,27 @@ 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 + -4594.16182364 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'] + [('load_const', 0.0), ('load_const', -3.0), ('load_arg', 0), ('ipow', 1), ('load_arg', 2), ('ipow', 1), 'mul', 'mul', 'add', ('load_const', -1.0), ('load_arg', 0), ('ipow', 1), ('load_arg', 1), ('ipow', 2), 'mul', 'mul', 'add', ('load_const', -6.0), ('load_arg', 0), ('ipow', 2), 'mul', 'add', ('load_const', -1.0), ('load_arg', 1), ('ipow', 2), 'mul', 'add', ('load_const', -1.0), ('load_arg', 0), ('ipow', 1), ('load_arg', 2), ('ipow', 2), 'mul', 'mul', 'add', 'return'] sage: fp(e, pi, sqrt(2)) - -98.001564033629322 + -98.0015640336 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: etb = ExpressionTreeBuilder(vars=('x','y'), domain=float) sage: x = etb.var('x') sage: y = etb.var('y') sage: expr = etb.call(sin, x^2 + y); expr - sin(add(pow(v_0, 2.0), v_1)) - sage: fc = fast_callable(expr, domain=RDF) + sin(add(ipow(v_0, 2), v_1)) + sage: fc = fast_callable(expr, domain=float) sage: fc(5, 7) 0.55142668124169059 """ @@ -357,11 +379,12 @@ str = InstructionStream(sage.ext.interpreters.wrapper_rr.metadata, len(vars), domain) - elif domain == RDF: + elif domain == RDF or domain is float: import sage.ext.interpreters.wrapper_rdf builder = sage.ext.interpreters.wrapper_rdf.Wrapper_rdf str = InstructionStream(sage.ext.interpreters.wrapper_rdf.metadata, - len(vars)) + len(vars), + domain) elif domain is None: import sage.ext.interpreters.wrapper_py builder = sage.ext.interpreters.wrapper_py.Wrapper_py @@ -550,9 +573,13 @@ sage: etb.var('y') Traceback (most recent call last): ... - ValueError: list.index(x): x not in list + ValueError: Variable 'y' not found """ - ind = self._vars.index(self._clean_var(v)) + var_name = self._clean_var(v) + try: + ind = self._vars.index(var_name) + except ValueError: + raise ValueError, "Variable '%s' not found" % var_name return ExpressionVariable(self, ind) def _var_number(self, n): @@ -581,6 +608,10 @@ The arguments will be converted to Expressions using ExpressionTreeBuilder.__call__. + As a special case, notices if the function is operator.pow and + the second argument is integral, and constructs an ExpressionIPow + instead. + EXAMPLES: sage: from sage.ext.fast_callable import ExpressionTreeBuilder sage: etb = ExpressionTreeBuilder(vars=(x,)) @@ -592,8 +623,14 @@ sin(1) sage: etb.call(factorial, x+57) {factorial}(add(v_0, 57)) + sage: etb.call(operator.pow, x, 543) + ipow(v_0, 543) """ - return ExpressionCall(self, fn, map(self, args)) + if fn is operator.pow: + base, exponent = args + return self(base)**exponent + else: + return ExpressionCall(self, fn, map(self, args)) def choice(self, cond, iftrue, iffalse): r""" @@ -791,6 +828,9 @@ r""" Compute a power expression from two Expressions. + If the second Expression is a constant integer, then return + an ExpressionIPow instead of an ExpressionCall. + EXAMPLES: sage: from sage.ext.fast_callable import ExpressionTreeBuilder sage: etb = ExpressionTreeBuilder(vars=(x,)) @@ -798,9 +838,11 @@ sage: x^x pow(v_0, v_0) sage: x^1 - pow(v_0, 1) + ipow(v_0, 1) sage: x.__pow__(1) - pow(v_0, 1) + ipow(v_0, 1) + sage: x.__pow__(1.0) + pow(v_0, 1.00000000000000) sage: x.__rpow__(1) pow(1, v_0) """ @@ -811,7 +853,19 @@ # (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) + + cdef Expression es + if isinstance(o, (int, long, Integer)): + es = s + return ExpressionIPow(es._etb, s, o) + else: + # I really don't like this, but I can't think of a better way + from sage.calculus.calculus import is_SymbolicExpression + if is_SymbolicExpression(o) and o in ZZ: + es = s + return ExpressionIPow(es._etb, s, ZZ(o)) + else: + return _expression_binop_helper(s, o, op_pow) def __neg__(self): r""" @@ -1088,6 +1142,88 @@ fn = function_name(self._function) return '%s(%s)' % (fn, ', '.join(map(repr, self._arguments))) +cdef class ExpressionIPow(Expression): + r""" + A power Expression with an integer exponent. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: type(etb.var('x')^17) + + """ + cdef object _base + cdef object _exponent + + def __init__(self, etb, base, exponent): + r""" + Initialize an ExpressionIPow. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder, ExpressionIPow + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb(x) + sage: x^(-12) + ipow(v_0, -12) + sage: v = ExpressionIPow(etb, x, 55); v + ipow(v_0, 55) + sage: v._get_etb() is etb + True + sage: v.base() + v_0 + sage: v.exponent() + 55 + """ + Expression.__init__(self, etb) + self._base = base + self._exponent = exponent + + def base(self): + r""" + Return the base from this ExpressionIPow. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: (etb(33)^42).base() + 33 + """ + return self._base + + def exponent(self): + r""" + Return the exponent from this ExpressionIPow. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: (etb(x)^(-1)).exponent() + -1 + """ + return self._exponent + + def __repr__(self): + r""" + Give a string representing this ExpressionIPow. + + EXAMPLES: + sage: from sage.ext.fast_callable import ExpressionTreeBuilder + sage: etb = ExpressionTreeBuilder(vars=(x,)) + sage: x = etb.var(x) + sage: x^3 + ipow(v_0, 3) + sage: x^(-2) + ipow(v_0, -2) + sage: v = (x+1)^3 + sage: v + ipow(add(v_0, 1), 3) + sage: repr(v) + 'ipow(add(v_0, 1), 3)' + sage: v.__repr__() + 'ipow(add(v_0, 1), 3)' + """ + return 'ipow(%s, %d)' % (repr(self._base), self._exponent) + cdef class ExpressionChoice(Expression): r""" A conditional expression. @@ -1239,6 +1375,85 @@ return ExpressionCall(self._etb, op, [self, other]) +class IntegerPowerFunction(object): + r""" + This class represents the function x^n for an arbitrary integral + power n. That is, IntegerPowerFunction(2) is the squaring function; + IntegerPowerFunction(-1) is the reciprocal function. + + EXAMPLES: + sage: from sage.ext.fast_callable import IntegerPowerFunction + sage: square = IntegerPowerFunction(2) + sage: square + (^2) + sage: square(pi) + pi^2 + sage: square(I) + -1 + sage: square(RIF(-1, 1)).str(style='brackets') + '[0.00000000000000000 .. 1.0000000000000000]' + sage: IntegerPowerFunction(-1) + (^(-1)) + sage: IntegerPowerFunction(-1)(22/7) + 7/22 + sage: v = Integers(123456789)(54321) + sage: v^9876543210 + 79745229 + sage: IntegerPowerFunction(9876543210)(v) + 79745229 + """ + + def __init__(self, n): + r""" + Initializes an IntegerPowerFunction. + + EXAMPLES: + sage: from sage.ext.fast_callable import IntegerPowerFunction + sage: cube = IntegerPowerFunction(3) + sage: cube + (^3) + sage: cube(AA(7)^(1/3)) + 7.000000000000000? + sage: cube.exponent + 3 + """ + self.exponent = n + + def __repr__(self): + r""" + Return a string representing this IntegerPowerFunction. + + EXAMPLES: + sage: from sage.ext.fast_callable import IntegerPowerFunction + sage: square = IntegerPowerFunction(2) + sage: square + (^2) + sage: repr(square) + '(^2)' + sage: square.__repr__() + '(^2)' + sage: repr(IntegerPowerFunction(-57)) + '(^(-57))' + """ + if self.exponent >= 0: + return "(^%s)" % self.exponent + else: + return "(^(%s))" % self.exponent + + def __call__(self, x): + r""" + Call this IntegerPowerFunction, to compute a power of its argument. + + EXAMPLES: + sage: from sage.ext.fast_callable import IntegerPowerFunction + sage: square = IntegerPowerFunction(2) + sage: square.__call__(5) + 25 + sage: square(5) + 25 + """ + return x**self.exponent + builtin_functions = None cpdef get_builtin_functions(): r""" @@ -1324,14 +1539,14 @@ 8.2831853071795864769252867665590057684 sage: fc = fast_callable(expr, domain=RDF) sage: fc(0) - 3.1415926535897931 + 3.14159265359 sage: fc(1) - 8.2831853071795862 + 8.28318530718 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 + 1.84147098481 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)) @@ -1341,7 +1556,7 @@ [('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 + 0.14112000806 sage: fc = fast_callable(etb.call(my_sin, x), domain=RealField(100)) sage: fc(3) 0.14112000805986722210074480281 @@ -1349,7 +1564,9 @@ [('load_arg', 0), ('py_call', , 1), 'return'] sage: fc = fast_callable(etb.call(my_sqrt, x), domain=RDF) sage: fc(3) - 1.7320508075688772 + 1.73205080757 + sage: parent(fc(3)) + Real Double Field sage: fc(-3) Traceback (most recent call last): ... @@ -1399,6 +1616,50 @@ Traceback (most recent call last): ... TypeError: unable to convert x (=sin(3)) to an integer + + sage: fc = fast_callable(etb(x)^100) + sage: fc(pi) + pi^100 + sage: fc = fast_callable(etb(x)^100, domain=ZZ) + sage: fc(2) + 1267650600228229401496703205376 + sage: fc = fast_callable(etb(x)^100, domain=RIF) + sage: fc(RIF(-2)) + 1.2676506002282295?e30 + sage: fc = fast_callable(etb(x)^100, domain=RDF) + sage: fc.op_list() + [('load_arg', 0), ('ipow', 100), 'return'] + sage: fc(1.1) + 13780.6123398 + sage: fc = fast_callable(etb(x)^100, domain=RR) + sage: fc.op_list() + [('load_arg', 0), ('ipow', 100), 'return'] + sage: fc(1.1) + 13780.6123398224 + sage: fc = fast_callable(etb(x)^(-100), domain=RDF) + sage: fc.op_list() + [('load_arg', 0), ('ipow', -100), 'return'] + sage: fc(1.1) + 7.25657159015e-05 + sage: fc = fast_callable(etb(x)^(-100), domain=RR) + sage: fc(1.1) + 0.0000725657159014814 + sage: expo = 2^32 + sage: base = (1.0).nextabove() + sage: fc = fast_callable(etb(x)^expo, domain=RDF) + sage: fc.op_list() + [('load_arg', 0), ('py_call', (^4294967296), 1), 'return'] + sage: fc(base) + 1.00000095367 + sage: RDF(base)^expo + 1.00000095367 + sage: fc = fast_callable(etb(x)^expo, domain=RR) + sage: fc.op_list() + [('load_arg', 0), ('py_call', (^4294967296), 1), 'return'] + sage: fc(base) + 1.00000095367477 + sage: base^expo + 1.00000095367477 """ cdef ExpressionConstant econst cdef ExpressionVariable evar @@ -1426,6 +1687,23 @@ stream.instr('py_call', fn, len(ecall._arguments)) else: raise ValueError, "Unhandled function %s in generate_code" % fn + elif isinstance(expr, ExpressionIPow): + base = expr.base() + exponent = expr.exponent() + metadata = stream.get_metadata() + ipow_range = metadata.ipow_range + if ipow_range is True: + use_ipow = True + elif isinstance(ipow_range, tuple): + a,b = ipow_range + use_ipow = (a <= exponent <= b) + else: + use_ipow = False + generate_code(base, stream) + if use_ipow: + stream.instr('ipow', exponent) + else: + stream.instr('py_call', IntegerPowerFunction(exponent), 1) else: raise ValueError, "Unhandled expression kind %s in generate_code" % type(expr) @@ -1604,6 +1882,8 @@ self._bytecode.append(loc) elif spec.parameters[i] == 'args': self._bytecode.append(args[i]) + elif spec.parameters[i] == 'code': + self._bytecode.append(args[i]) elif spec.parameters[i] == 'n_inputs': self._bytecode.append(args[i]) n_inputs = args[i] @@ -1631,7 +1911,8 @@ Returns the interpreter metadata being used by the current InstructionStream. - (Probably only useful for writing doctests.) + The code generator sometimes uses this to decide which code + to generate. EXAMPLES: sage: from sage.ext.interpreters.wrapper_rdf import metadata @@ -1684,7 +1965,7 @@ 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': []} + {'domain': None, 'code': [0, 0, 3, 0, 1, 12, 2], 'py_constants': [], 'args': 1, 'stack': 1, 'constants': []} """ d = {'args': self._n_args, 'constants': self._constants, @@ -1697,16 +1978,21 @@ 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. + consists of a dictionary mapping instruction names to + (CompilerInstrSpec, opcode) pairs, a list mapping opcodes to + (instruction name, CompilerInstrSpec) pairs, and a range of exponents + for which the ipow instruction can be used. This range can be + False (if the ipow instruction should never be used), a pair of + two integers (a,b), if ipow should be used for a<=n<=b, or True, + if ipow should always be used. When ipow cannot be used, then + we fall back on calling IntegerPowerFunction. See the class docstring for CompilerInstrSpec for more information. NOTE: You must not modify the metadata. """ - def __init__(self, by_opname, by_opcode): + def __init__(self, by_opname, by_opcode, ipow_range): r""" Initialize an InterpreterMetadata object. @@ -1715,14 +2001,17 @@ 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 = InterpreterMetadata(by_opname='opname dict goes here', by_opcode='opcode list goes here', ipow_range=(2, 57)) sage: metadata.by_opname 'opname dict goes here' sage: metadata.by_opcode 'opcode list goes here' + sage: metadata.ipow_range + (2, 57) """ self.by_opname = by_opname self.by_opcode = by_opcode + self.ipow_range = ipow_range class CompilerInstrSpec(object): r""" @@ -1817,7 +2106,7 @@ for p in instr.parameters: p_loc = code[0] code = code[1:] - if p in ('args', 'n_inputs', 'n_outputs'): + if p in ('args', 'code', 'n_inputs', 'n_outputs'): op.append(p_loc) else: op.append(args[p][p_loc]) @@ -1869,7 +2158,7 @@ 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': []} + {'domain': Real Double Field, 'code': [0, 0, 16, 0, 0, 7, 2], 'py_constants': [], 'args': 1, 'stack': 2, 'constants': []} """ return self._orig_args diff -r c8ebd6134a6f -r a97f9d176f07 sage/ext/fast_eval.pyx --- a/sage/ext/fast_eval.pyx Sat Feb 28 20:13:00 2009 -0800 +++ b/sage/ext/fast_eval.pyx Wed Mar 18 10:22:40 2009 -0700 @@ -88,7 +88,6 @@ #***************************************************************************** from sage.ext.fast_callable import fast_callable, Wrapper -from sage.rings.all import RDF include "stdsage.pxi" @@ -1348,7 +1347,7 @@ if old: return f._fast_float_(*vars) else: - return fast_callable(f, vars=vars, domain=RDF) + return fast_callable(f, vars=vars, domain=float) except AttributeError: pass diff -r c8ebd6134a6f -r a97f9d176f07 sage/ext/gen_interpreters.py --- a/sage/ext/gen_interpreters.py Sat Feb 28 20:13:00 2009 -0800 +++ b/sage/ext/gen_interpreters.py Wed Mar 18 10:22:40 2009 -0700 @@ -104,12 +104,12 @@ from distutils.extension import Extension ############################## -# This module is used during the Sage buld process, so it should not +# This module is used during the Sage build 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, +# interpreter configuration is split out into a separate file, # that will have to be changed. ############################## @@ -1644,7 +1644,8 @@ chunk = chunks[chunk_code] addr = None ch_len = None - if chunk.is_stack(): + # shouldn't hardcode 'code' here + if chunk.is_stack() or chunk.name == 'code': pass else: m = re.match(r'\[(?:([0-9]+)|([a-zA-Z]))\]', s) @@ -1959,23 +1960,33 @@ 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 + extra_class_members -- Class members for the wrapper that + don't correspond to memory chunks + extra_members_initialize -- Code to initialize extra_class_members EXAMPLES: sage: from sage.ext.gen_interpreters import * sage: interp = RDFInterpreter() sage: interp.header - '' + '\n#include \n' sage: interp.pyx_header '' sage: interp.err_return '-1094648009105371' sage: interp.mc_code {MC:code} + sage: interp = RRInterpreter() + sage: interp.extra_class_members + '' + sage: interp.extra_members_initialize + '' """ self.header = '' self.pyx_header = '' self.err_return = 'NULL' self.mc_code = MemoryChunkConstants('code', ty_int) + self.extra_class_members = '' + self.extra_members_initialize = '' def _set_opcodes(self): r""" @@ -2017,6 +2028,11 @@ 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 + ipow_range -- the range of exponents supported by the ipow + instruction (default is False, meaning never use ipow) + adjust_retval -- None, or a string naming a function to call + in the wrapper's __call__ to modify the return + value of the interpreter EXAMPLES: sage: from sage.ext.gen_interpreters import * @@ -2044,11 +2060,16 @@ else: self.return_type = None self.mc_retval = mc_retval + self.ipow_range = False + self.adjust_retval = None class RDFInterpreter(StackInterpreter): r""" A subclass of StackInterpreter, specifying an interpreter over - machine-floating-point values (C doubles). + machine-floating-point values (C doubles). This is used for + both domain=RDF and domain=float; currently the only difference + between the two is the type of the value returned from the + wrapper (they use the same wrapper and interpreter). """ def __init__(self): @@ -2060,6 +2081,12 @@ sage: interp = RDFInterpreter() sage: interp.name 'rdf' + sage: interp.extra_class_members + 'cdef object _domain\n' + sage: interp.extra_members_initialize + "self._domain = args['domain']\n" + sage: interp.adjust_retval + 'self._domain' sage: interp.mc_py_constants {MC:py_constants} sage: interp.chunks @@ -2087,6 +2114,9 @@ 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 +""" instrs = [ InstrSpec('load_arg', pg('A[D]', 'S'), code='o0 = i0;'), @@ -2123,6 +2153,7 @@ ('mul', '*'), ('div', '/')]: instrs.append(instr_infix(name, pg('SS', 'S'), op)) instrs.append(instr_funcall_2args('pow', pg('SS', 'S'), 'pow')) + instrs.append(instr_funcall_2args('ipow', pg('SD', 'S'), 'gsl_pow_int')) for (name, op) in [('neg', '-i0'), ('invert', '1/i0'), ('abs', 'fabs(i0)')]: instrs.append(instr_unary(name, pg('S', 'S'), op)) @@ -2132,6 +2163,11 @@ instrs.append(instr_unary(name, pg('S', 'S'), "%s(i0)" % name)) self.instr_descs = instrs self._set_opcodes() + # supported for exponents that fit in an int + self.ipow_range = (int(-2**31), int(2**31-1)) + self.extra_class_members = "cdef object _domain\n" + self.extra_members_initialize = "self._domain = args['domain']\n" + self.adjust_retval = 'self._domain' class RRInterpreter(StackInterpreter): r""" @@ -2255,6 +2291,7 @@ ('mul', 'mpfr_mul'), ('div', 'mpfr_div'), ('pow', 'mpfr_pow')]: instrs.append(instr_funcall_2args_mpfr(name, pg('SS', 'S'), op)) + instrs.append(instr_funcall_2args_mpfr('ipow', pg('SD', 'S'), 'mpfr_pow_si')) for name in ['neg', 'abs', 'log', 'log2', 'log10', 'exp', 'exp2', 'exp10', @@ -2277,6 +2314,11 @@ code='mpfr_ui_div(o0, 1, i0, GMP_RNDN);')) self.instr_descs = instrs self._set_opcodes() + # Supported for exponents that fit in a long, so we could use + # a much wider range on a 64-bit machine. On the other hand, + # it's easier to write the code this way, and constant integer + # exponents outside this range probably aren't very common anyway. + self.ipow_range = (int(-2**31), int(2**31-1)) class PythonInterpreter(StackInterpreter): r""" @@ -2300,7 +2342,7 @@ 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 + (and ensures) that empty parts of the stack contain NULL, so it doesn't bother to Py_XDECREF before it pushes onto the stack. """ @@ -2368,12 +2410,16 @@ 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);')) + instrs.append(InstrSpec('ipow', pg('SC[D]', '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() + # Always use ipow + self.ipow_range = True class ElementInterpreter(PythonInterpreter): r""" @@ -2542,7 +2588,10 @@ 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: + # Shouldn't hardcode 'code' here + if ch.name == 'code': + w(" %s i%d = %s;\n" % (chst.c_local_type(), i, string_of_addr(ch))) + elif input_len is not None: w(" %s i%d = %s + ai%d;\n" % (chst.c_ptr_type(), i, ch.name, i)) else: @@ -2742,11 +2791,13 @@ the_call = je(""" {% if s.return_type %}return {% endif -%} +{% if s.adjust_retval %}{{ s.adjust_retval }}({% endif %} interp_{{ s.name }}({{ arg_ch.pass_argument() }} {% for ch in s.chunks[1:] %} , {{ ch.pass_argument() }} {% endfor %} - ) + ){% if s.adjust_retval %}){% endif %} + """, s=s, arg_ch=arg_ch) w(je(""" @@ -2783,6 +2834,7 @@ {% for ch in s.chunks %} {% print ch.declare_class_members() %} {% endfor %} +{% print indent_lines(4, s.extra_class_members) %} def __init__(self, args): Wrapper.__init__(self, args, metadata) @@ -2795,6 +2847,7 @@ {% for ch in s.chunks %} {% print ch.init_class_members() %} {% endfor %} +{% print indent_lines(8, s.extra_members_initialize) %} def __dealloc__(self): cdef int i @@ -2840,7 +2893,8 @@ ('{{ instr.name }}', CompilerInstrSpec({{ instr.n_inputs }}, {{ instr.n_outputs }}, {{ instr.parameters }})), {% endfor %} - ]) + ], + ipow_range={{ s.ipow_range }}) """, 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): @@ -2895,7 +2949,7 @@ simplest instructions: sage: print rdf_interp /* ... */ ... - case 9: /* neg */ + case 10: /* neg */ { double i0 = *--stack; double o0; @@ -2913,7 +2967,7 @@ type. sage: print rr_interp /* ... */ ... - case 9: /* neg */ + case 10: /* neg */ { mpfr_ptr i0 = *--stack; mpfr_ptr o0 = *stack++; @@ -2932,7 +2986,7 @@ Python-object element interpreter. sage: print el_interp /* ... */ ... - case 9: /* neg */ + case 10: /* neg */ { PyObject* i0 = *--stack; *stack = NULL; @@ -3153,6 +3207,12 @@ 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. + + We optionally adjust the return value of the interpreter + (currently only the RDF/float interpreter performs this step; + this is the only place where domain=RDF differs than + domain=float): + sage: print rdf_wrapper # ... def __call__(self, *args): @@ -3161,12 +3221,12 @@ cdef int i for i from 0 <= i < len(args): self._args[i] = args[i] - return interp_rdf(c_args + return self._domain(interp_rdf(c_args , self._constants , self._py_constants , self._stack , self._code - ) + )) ... In Python-object based interpreters, the call to the C @@ -3200,12 +3260,13 @@ documented at InterpreterMetadata; for now, we'll just show what it looks like. - Currently, there are two parts to the metadata; the first maps + Currently, there are three 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. + we'll need at runtime. The third part says what range the + ipow instruction is defined over. First the part that maps instruction names to (CompilerInstrSpec, opcode) pairs. @@ -3237,7 +3298,14 @@ ('add', CompilerInstrSpec(2, 1, [])), ... - ]) + ], ...) + + And then the ipow range: + sage: print rdf_wrapper + # ... + metadata = InterpreterMetadata(..., + ipow_range=(-2147483648, 2147483647)) + And that's it for the wrapper. """ @@ -3367,7 +3435,8 @@ modules = [ Extension('sage.ext.interpreters.wrapper_rdf', sources = ['sage/ext/interpreters/wrapper_rdf.pyx', - 'sage/ext/interpreters/interp_rdf.c']), + 'sage/ext/interpreters/interp_rdf.c'], + libraries = ['gsl']), Extension('sage.ext.interpreters.wrapper_rr', sources = ['sage/ext/interpreters/wrapper_rr.pyx', diff -r c8ebd6134a6f -r a97f9d176f07 sage/numerical/optimize.py --- a/sage/numerical/optimize.py Sat Feb 28 20:13:00 2009 -0800 +++ b/sage/numerical/optimize.py Wed Mar 18 10:22:40 2009 -0700 @@ -235,10 +235,10 @@ if isinstance(func,SymbolicExpression): var_list=func.variables() var_names=map(str,var_list) - fast_f=fast_callable(func, vars=var_names, domain=RDF) + fast_f=fast_callable(func, vars=var_names, domain=float) f=lambda p: fast_f(*p) gradient_list=func.gradient() - fast_gradient_functions=[fast_callable(gradient_list[i], vars=var_names, domain=RDF) for i in xrange(len(gradient_list))] + fast_gradient_functions=[fast_callable(gradient_list[i], vars=var_names, domain=float) for i in xrange(len(gradient_list))] gradient=lambda p: scipy.array([ a(*p) for a in fast_gradient_functions]) else: f=func @@ -260,7 +260,7 @@ elif algorithm=="ncg": if isinstance(func,SymbolicExpression): hess=func.hessian() - hess_fast= [ [fast_callable(a, vars=var_names, domain=RDF) for a in row] for row in hess] + hess_fast= [ [fast_callable(a, vars=var_names, domain=float) 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 c8ebd6134a6f -r a97f9d176f07 sage/rings/polynomial/multi_polynomial.pyx --- a/sage/rings/polynomial/multi_polynomial.pyx Sat Feb 28 20:13:00 2009 -0800 +++ b/sage/rings/polynomial/multi_polynomial.pyx Wed Mar 18 10:22:40 2009 -0700 @@ -260,9 +260,15 @@ 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))) + add(add(add(0, mul(-6/5, mul(mul(ipow(v_0, 1), ipow(v_1, 1)), ipow(v_2, 1)))), mul(2, mul(ipow(v_1, 1), ipow(v_2, 2)))), mul(-1, ipow(v_0, 1))) TESTS: + sage: v = K(0) + sage: vf = fast_callable(v) + sage: type(v(0r, 0r, 0r)) + + sage: type(vf(0r, 0r, 0r)) + sage: K. = QQ[] sage: from sage.ext.fast_eval import fast_float sage: fast_float(K(0)).op_list() @@ -270,13 +276,13 @@ 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'] + [('load_const', 0.0), ('load_const', 1.0), ('load_arg', 1), ('ipow', 1), '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) + expr = etb.constant(self.base_ring()(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))