fast_eval.pyx

r"""
Fast Numerical 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. 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.

Up until now the solution has been to use lambda expressions, but this
is neither intuitive, Sage-like, nor efficient (compared to operating
on raw C doubles).  This module provides a representation of algebraic
expression in Reverse Polish Notation, and provides an efficient
interpreter on C double values as a callable python object. It does
what it can in C, and will call out to Python if necessary.

Essential to the understanding of this class is the distinction
between symbolic expressions and callable symbolic expressions (where
the latter binds argument names to argument positions). The
\code{*vars} parameter passed around encapsulates this information.

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
constant function), \code{fast_float_arg} (selects the $n$-th value
when called with $\ge n$ arguments), and \code{fast_float_func} which
wraps a callable Python function. These may be combined with the
standard Python arithmetic operators, and support many of the basic
math functions such sqrt, exp, and trig functions.

EXAMPLES: 
    sage: from sage.ext.fast_eval import fast_float
    sage: f = fast_float(sqrt(x^7+1), 'x', old=True)
    sage: f(1)
    1.4142135623730951
    sage: f.op_list()
    ['load 0', 'push 7.0', 'pow', 'push 1.0', 'add', 'call sqrt(1)']
    
    To interpret that last line, we load argument 0 ('x' in this case) onto 
    the stack, push the constant 2.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_eval import fast_float_arg
    sage: g = fast_float_arg(0).sin()
    sage: (f+g).op_list()
    ['load 0', 'push 7.0', 'pow', 'push 1.0', 'add', 'call sqrt(1)', 'load 0', 'call sin(1)', 'add']

TESTS:
This used to segfault because of an assumption that assigning None to a
variable would raise a TypeError:
    sage: from sage.ext.fast_eval import fast_float_arg, fast_float
    sage: fast_float_arg(0)+None
    Traceback (most recent call last):
    ...
    TypeError

AUTHOR: 
    -- Robert Bradshaw (2008-10): Initial version
"""


#*****************************************************************************
#       Copyright (C) 2008 Robert Bradshaw <[email protected]>
#
#  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/
#*****************************************************************************

from sage.ext.fast_callable import fast_callable, Wrapper

include "stdsage.pxi"

cdef extern from "Python.h":
    int PyInt_AS_LONG(PyObject*)
    PyObject* PyTuple_New(Py_ssize_t size)
    PyObject* PyTuple_GET_ITEM(PyObject* t, Py_ssize_t index)
    void PyTuple_SET_ITEM(PyObject* t, Py_ssize_t index, PyObject* item)
    object PyObject_CallObject(PyObject* func, PyObject* args)
    PyObject* PyFloat_FromDouble(double d)
    void Py_DECREF(PyObject *)

cdef extern from "math.h":
    double sqrt(double)
    double pow(double, double)

    double ceil(double)
    double floor(double)

    double sin(double)
    double cos(double)
    double tan(double)
    
    double asin(double)
    double acos(double)
    double atan(double)
    double atan2(double, double)
    
    double sinh(double)
    double cosh(double)
    double tanh(double)
    
    double asinh(double)
    double acosh(double)
    double atanh(double)
    
    double exp(double)
    double log(double)
    double log10(double)
    double log2_ "log2"(double)


# This is only needed on Cygwin since log2 is a macro.
# If we don't do this the cygwin GCC gets very confused.
cdef inline double log2(double x): 
    return log2_(x)

cdef extern from *:
    void* memcpy(void* dst, void* src, size_t len)

cdef inline int max(int a, int b):
    return a if a > b else b

cdef inline int min(int a, int b):
    return a if a < b else b

cdef enum:
# stack
    LOAD_ARG    # push input argument n onto the stack
    PUSH_CONST
    POP
    POP_N
    DUP

# basic arithmetic
    ADD
    SUB
    MUL
    DIV
    NEG
    ABS
    INVERT
    POW
    
# basic comparison
    LT
    LE
    EQ
    NE
    GT
    GE

# functional
    ONE_ARG_FUNC
    TWO_ARG_FUNC
    PY_FUNC
    

# These two dictionaries are for printable and machine independent representation. 

op_names = {
    LOAD_ARG: 'load',
    PUSH_CONST: 'push',
    POP: 'pop',
    POP_N: 'popn',
    DUP: 'dup',
    
    ADD: 'add',
    SUB: 'sub',
    MUL: 'mul',
    DIV: 'div',
    NEG: 'neg',
    ABS: 'abs',
    INVERT: 'invert',
    POW: 'pow',


    LT: 'lt',
    LE: 'le',
    EQ: 'eq',
    NE: 'ne',
    GT: 'gt',
    GE: 'ge',


    ONE_ARG_FUNC: 'call',
    TWO_ARG_FUNC: 'call',
    PY_FUNC: 'py_call',
}

cfunc_names = {
    <size_t>&sqrt: 'sqrt',
    <size_t>&pow: 'pow',
    
    <size_t>&ceil: 'ceil',
    <size_t>&floor: 'floor',
    
    <size_t>&sin: 'sin',
    <size_t>&cos: 'cos',
    <size_t>&tan: 'tan',
    
    <size_t>&asin: 'asin',
    <size_t>&atan: 'atan',
    <size_t>&atan2: 'atan2',
    
    <size_t>&sinh: 'sinh',
    <size_t>&cosh: 'cosh',
    <size_t>&tanh: 'tanh',

    <size_t>&asinh: 'asinh',
    <size_t>&acosh: 'acosh',
    <size_t>&atanh: 'atanh',

    <size_t>&exp: 'exp',
    <size_t>&log: 'log',
    <size_t>&log2: 'log2',
    <size_t>&log10: 'log10',

}

cdef reverse_map(m):
    r = {}
    for key, value in m.iteritems():
        r[value] = key
    return r

# With all the functionality around the op struct, perhaps there should be 
# a wrapper class, though we still wish to operate on pure structs for speed. 

cdef op_to_string(fast_double_op op):
    s = op_names[op.type]
    if op.type in [LOAD_ARG, POP_N]:
        s += " %s" % op.params.n
    elif op.type == PUSH_CONST:
        s += " %s" % op.params.c
    elif op.type in [ONE_ARG_FUNC, TWO_ARG_FUNC]:
        try:
            cname = cfunc_names[<size_t>op.params.func]
        except KeyError:
            cname = "0x%x" % <size_t>op.params.func
        s += " %s(%s)" % (cname, 1 if op.type == ONE_ARG_FUNC else 2)
    elif op.type == PY_FUNC:
        n, func = <object>(op.params.func)
        s += " %s(%s)" % (func, n)
    return s
    
cdef op_to_tuple(fast_double_op op):
    s = op_names[op.type]
    if op.type in [LOAD_ARG, POP_N]:
        param = op.params.n
    elif op.type == PUSH_CONST:
        param = op.params.c
    elif op.type in [ONE_ARG_FUNC, TWO_ARG_FUNC]:
        param_count = 1 if op.type == ONE_ARG_FUNC else 2
        try:
            param = param_count, cfunc_names[<size_t>op.params.func]
        except KeyError:
            raise ValueError, "Unknown C function: 0x%x" % <size_t>op.params.func
    elif op.type == PY_FUNC:
        param = <object>(op.params.func)
    else:
        param = None
    if param is None:
        return (s,)
    else:
        return s, param
        
def _unpickle_FastDoubleFunc(nargs, max_height, op_list):
    cdef FastDoubleFunc self = PY_NEW(FastDoubleFunc)
    self.nops = len(op_list)
    self.nargs = nargs
    self.max_height = max_height
    self.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op) * self.nops)
    self.allocate_stack()
    cfunc_addresses = reverse_map(cfunc_names)
    op_enums = reverse_map(op_names)
    cdef size_t address
    cdef int i = 0, type
    for op in op_list:
        self.ops[i].type = type = op_enums[op[0]]
        if type in [LOAD_ARG, POP_N]:
            self.ops[i].params.n = op[1]
        elif type == PUSH_CONST:
            self.ops[i].params.c = op[1]
        elif type in [ONE_ARG_FUNC, TWO_ARG_FUNC]:
            param_count, cfunc = op[1]
            address = cfunc_addresses[cfunc]
            self.ops[i].params.func = <void *>address
            self.ops[i].type = ['', ONE_ARG_FUNC, TWO_ARG_FUNC][param_count]
        elif type == PY_FUNC:
            if self.py_funcs is None:
                self.py_funcs = op[1]
            else:
                self.py_funcs = self.py_funcs + (op[1],)
            self.ops[i].params.func = <void *>op[1]
        i += 1
    return self
    
    
# This is where we wish we had case statements...
# It looks like gcc might be smart enough to figure it out. 
cdef inline int process_op(fast_double_op op, double* stack, double* argv, int top) except -2:

#    print [stack[i] for i from 0 <= i <= top], ':', op.type

    cdef int i, n
    cdef PyObject* py_args
    
    # We have to do some trickery because Pyrex disallows function pointer casts
    # This will be removed in a future version of Cython.
    cdef double (*f)(double)
    cdef void** fp = <void **>&f
    cdef double (*ff)(double, double)
    cdef void** ffp = <void **>&ff

    if op.type == LOAD_ARG:
        stack[top+1] = argv[op.params.n]
        return top+1
        
    elif op.type == PUSH_CONST:
        stack[top+1] = op.params.c
        return top+1
        
    elif op.type == POP:
        return top-1
        
    elif op.type == POP_N:
        return top-op.params.n
        
    elif op.type == DUP:
        stack[top+1] = stack[top]
        return top+1

    elif op.type == ADD:
        stack[top-1] += stack[top]
        return top-1

    elif op.type == SUB:
        stack[top-1] -= stack[top]
        return top-1

    elif op.type == MUL:
        stack[top-1] *= stack[top]
        return top-1

    elif op.type == DIV:
        stack[top-1] /= stack[top]
        return top-1
        
    elif op.type == NEG:
        stack[top] = -stack[top]
        return top
        
    elif op.type == ABS:
        if stack[top] < 0:
            stack[top] = -stack[top]
        return top
        
    elif op.type == INVERT:
        stack[top] = 1/stack[top]
        return top
        
    elif op.type == POW:
        if stack[top-1] < 0 and stack[top] != floor(stack[top]):
            raise ValueError, "negative number to a fractional power not real"
        stack[top-1] = pow(stack[top-1], stack[top])
        return top-1
        
    elif op.type == LT:
        stack[top-1] = 1.0 if stack[top-1] < stack[top] else 0.0
        return top-1

    elif op.type == LE:
        stack[top-1] = 1.0 if stack[top-1] <= stack[top] else 0.0
        return top-1

    elif op.type == EQ:
        stack[top-1] = 1.0 if stack[top-1] == stack[top] else 0.0
        return top-1

    elif op.type == NE:
        stack[top-1] = 1.0 if stack[top-1] != stack[top] else 0.0
        return top-1

    elif op.type == GT:
        stack[top-1] = 1.0 if stack[top-1] > stack[top] else 0.0
        return top-1

    elif op.type == GE:
        stack[top-1] = 1.0 if stack[top-1] >= stack[top] else 0.0
        return top-1

    elif op.type == ONE_ARG_FUNC:
        fp[0] = op.params.func
        stack[top] = f(stack[top])
        return top

    elif op.type == TWO_ARG_FUNC:
        ffp[0] = op.params.func
        stack[top-1] = ff(stack[top-1], stack[top])
        return top-1

    elif op.type == PY_FUNC:
        # Even though it's python, optimize this because it'll be used often...
        # We also want to avoid cluttering the header and footer
        # of this function Python variables bring. 
        n = PyInt_AS_LONG(PyTuple_GET_ITEM(op.params.func, 0))
        top = top - n + 1
        py_args = PyTuple_New(n)
        for i from 0 <= i < n:
            PyTuple_SET_ITEM(py_args, i, PyFloat_FromDouble(stack[top+i]))
        stack[top] = PyObject_CallObject(PyTuple_GET_ITEM(op.params.func, 1), py_args)
        Py_DECREF(py_args)
        return top
        
    raise RuntimeError, "Bad op code %s" % op.type


cdef class FastDoubleFunc:
    """
    This class is for fast evaluation of algebraic expressions over 
    the real numbers (e.g. for plotting). It represents an expression 
    as a stack-based series of operations. 
    
    EXAMPLES: 
        sage: from sage.ext.fast_eval import FastDoubleFunc
        sage: f = FastDoubleFunc('const', 1.5) # the constant function
        sage: f()
        1.5
        sage: g = FastDoubleFunc('arg', 0) # the first argument
        sage: g(5)
        5.0
        sage: h = f+g
        sage: h(17)
        18.5
        sage: h = h.sin()
        sage: h(pi/2-1.5)
        1.0
        sage: h.is_pure_c()
        True
        sage: list(h)
        ['push 1.5', 'load 0', 'add', 'call sin(1)']
        
    We can wrap Python functions too:
        sage: h = FastDoubleFunc('callable', lambda x,y: x*x*x - y, g, f)
        sage: h(10)
        998.5
        sage: h.is_pure_c()
        False
        sage: list(h)        # random address
        ['load 0', 'push 1.5', 'py_call <function <lambda> at 0x9fedf70>(2)']
        
    Here's a more complicated expression:
        sage: from sage.ext.fast_eval import fast_float_constant, fast_float_arg
        sage: a = fast_float_constant(1.5)
        sage: b = fast_float_constant(3.14)
        sage: c = fast_float_constant(7)
        sage: x = fast_float_arg(0)
        sage: y = fast_float_arg(1)
        sage: f = a*x^2 + b*x + c - y/sqrt(sin(y)^2+a)
        sage: f(2,3)
        16.846610528508116
        sage: f.max_height
        4
        sage: f.is_pure_c()
        True
        sage: list(f)
        ['push 1.5', 'load 0', 'dup', 'mul', 'mul', 'push 3.14', 'load 0', 'mul', 'add', 'push 7.0', 'add', 'load 1', 'load 1', 'call sin(1)', 'dup', 'mul', 'push 1.5', 'add', 'call sqrt(1)', 'div', 'sub']
        
    
    AUTHOR: 
        -- Robert Bradshaw
    """
    def __init__(self, type, param, *args):

        cdef FastDoubleFunc arg
        cdef int i
        
        if type == 'arg':
            self.nargs = param+1
            self.nops = 1
            self.max_height = 1
            self.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op))
            self.ops[0].type = LOAD_ARG
            self.ops[0].params.n = param
            
        elif type == 'const':
            self.nargs = 0
            self.nops = 1
            self.max_height = 1
            self.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op))
            self.ops[0].type = PUSH_CONST
            self.ops[0].params.c = param
            
        elif type == 'callable':
            py_func = len(args), param
            self.py_funcs = (py_func,) # just so it doesn't get garbage collected
            self.nops = 1
            self.nargs = 0
            for i from 0 <= i < len(args):
                a = args[i]
                if not isinstance(a, FastDoubleFunc):
                     a = FastDoubleFunc('const', a)
                     args = args[:i] + (a,) + args[i+1:]
                arg = a
                self.nops += arg.nops
                if arg.py_funcs is not None:
                    self.py_funcs += arg.py_funcs
                self.nargs = max(self.nargs, arg.nargs)
                self.max_height = max(self.max_height, arg.max_height+i)
            self.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op) * self.nops)
            if self.ops == NULL:
                raise MemoryError
            i = 0
            for arg in args:
                memcpy(self.ops + i, arg.ops, sizeof(fast_double_op) * arg.nops)
                i += arg.nops
            self.ops[self.nops-1].type = PY_FUNC
            self.ops[self.nops-1].params.func = <void *>py_func
            
        else:
            raise ValueError, "Unknown operation: %s" % type
            
        self.allocate_stack()
            
            
    cdef int allocate_stack(FastDoubleFunc self) except -1:
        self.argv = <double*>sage_malloc(sizeof(double) * self.nargs)
        if self.argv == NULL:
            raise MemoryError
        self.stack = <double*>sage_malloc(sizeof(double) * self.max_height)
        if self.stack == NULL:
            raise MemoryError
            
    def __dealloc__(self):
        if self.ops:
            sage_free(self.ops)
        if self.stack:
            sage_free(self.stack)
        if self.argv:
            sage_free(self.argv)
            
    def __reduce__(self):
        """
        TESTS:
            sage: from sage.ext.fast_eval import fast_float_arg, fast_float_func
            sage: f = fast_float_arg(0).sin() * 10 + fast_float_func(hash, fast_float_arg(1))
            sage: loads(dumps(f)) == f
            True
        """
        L = [op_to_tuple(self.ops[i]) for i from 0 <= i < self.nops] 
        return _unpickle_FastDoubleFunc, (self.nargs, self.max_height, L)
        
    def __cmp__(self, other):
        """
        Two functions are considered equal if they represent the same
        exact sequence of operations. 
        
        TESTS:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: fast_float_arg(0) == fast_float_arg(0)
            True
            sage: fast_float_arg(0) == fast_float_arg(1)
            False
            sage: fast_float_arg(0) == fast_float_arg(0).sin()
            False
        """
        cdef int c, i
        cdef FastDoubleFunc left, right
        try:
            left, right = self, other
            c = cmp((left.nargs, left.nops, left.max_height), 
                    (right.nargs, right.nops, right.max_height))
            if c != 0:
                return c
            for i from 0 <= i < self.nops:
                if left.ops[i].type != right.ops[i].type:
                    return cmp(left.ops[i].type, right.ops[i].type)
            for i from 0 <= i < self.nops:
                c = cmp(op_to_tuple(left.ops[i]), op_to_tuple(right.ops[i]))
                if c != 0:
                    return c
            return c
        except TypeError:
            return cmp(type(self), type(other))
    
    def __call__(FastDoubleFunc self, *args):
        """
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(2)
            sage: f(0,1,2,3)
            2.0
            sage: f(10)
            Traceback (most recent call last):
            ...
            TypeError: Wrong number of arguments (need at least 3, got 1)
            sage: f('blah', 1, 2, 3)
            Traceback (most recent call last):
            ...
            TypeError: a float is required
        """
        if len(args) < self.nargs:
            raise TypeError, "Wrong number of arguments (need at least %s, got %s)" % (self.nargs, len(args))
        cdef int i = 0
        for i from 0 <= i < self.nargs:
            self.argv[i] = args[i]
        res = self._call_c(self.argv)
        return res
        
    cdef double _call_c(FastDoubleFunc self, double* argv) except? -2:
        # The caller must assure that argv has length at least self.nargs
        # The bulk of this function is in the (inlined) function process_op. 
        cdef int i, top = -1
        for i from 0 <= i < self.nops:
            top = process_op(self.ops[i], self.stack, argv, top)
        cdef double res = self.stack[0]
        return res
        
    def _fast_float_(self, *vars):
        r"""
        Returns \code{self} if there are enough arguments, otherwise raises a TypeError. 
        
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(1)
            sage: f._fast_float_('x','y') is f
            True
            sage: f._fast_float_('x') is f
            Traceback (most recent call last):
            ...
            TypeError: Needs at least 2 arguments (1 provided)
        """
        if self.nargs > len(vars):
            raise TypeError, "Needs at least %s arguments (%s provided)" % (self.nargs, len(vars))
        return self
        
    def op_list(self):
        """
        Returns a list of string representations of the 
        operations that make up this expression. 
        
        Python and C function calls may be only available by function pointer addresses. 
        
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_constant, fast_float_arg
            sage: a = fast_float_constant(17)
            sage: x = fast_float_arg(0)
            sage: a.op_list()
            ['push 17.0']
            sage: x.op_list()
            ['load 0']
            sage: (a*x).op_list()
            ['push 17.0', 'load 0', 'mul']
            sage: (a+a*x^2).sqrt().op_list()
            ['push 17.0', 'push 17.0', 'load 0', 'dup', 'mul', 'mul', 'add', 'call sqrt(1)']
        """
        cdef int i
        return [op_to_string(self.ops[i]) for i from 0 <= i < self.nops]
        
    def __iter__(self):
        """
        Returns the list of operations of self. 
        
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0)*2 + 3
            sage: list(f)
            ['load 0', 'push 2.0', 'mul', 'push 3.0', 'add']
        """
        return iter(self.op_list())
        
    cpdef bint is_pure_c(self):
        """
        Returns True if this function can be evaluated without
        any python calls (at any level). 
        
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_constant, fast_float_arg, fast_float_func
            sage: fast_float_constant(2).is_pure_c()
            True
            sage: fast_float_arg(2).sqrt().sin().is_pure_c()
            True
            sage: fast_float_func(lambda _: 2).is_pure_c()
            False
        """
        cdef int i
        for i from 0 <= i < self.nops:
            if self.ops[i].type == PY_FUNC:
                return 0
        return 1
        
    def python_calls(self):
        """
        Returns a list of all python calls used by function. 
        
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_func, fast_float_arg
            sage: x = fast_float_arg(0)
            sage: f = fast_float_func(hash, sqrt(x))
            sage: f.op_list()
            ['load 0', 'call sqrt(1)', 'py_call <built-in function hash>(1)']
            sage: f.python_calls()
            [<built-in function hash>]
        """
        L = []
        cdef int i
        for i from 0 <= i < self.nops:
            if self.ops[i].type == PY_FUNC:
                L.append((<object>self.ops[i].params.func)[1])
        return L
        
    ###################################################################
    #   Basic Arithmetic
    ###################################################################
        
    def __add__(left, right):
        """
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0) + fast_float_arg(1)
            sage: f(3,4)
            7.0
        """
        return binop(left, right, ADD)
        
    def __sub__(left, right):
        """
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0) - fast_float_arg(2)
            sage: f(3,4,5)
            -2.0
        """
        return binop(left, right, SUB)
        
    def __mul__(left, right):
        """
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0) * 2
            sage: f(17)
            34.0
        """
        return binop(left, right, MUL)
        
    def __div__(left, right):
        """
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0) / 7
            sage: f(14)
            2.0
        """
        return binop(left, right, DIV)
        
    def __pow__(FastDoubleFunc left, right, dummy):
        """
        EXAMPLES: 
            sage: from sage.ext.fast_eval import FastDoubleFunc
            sage: f = FastDoubleFunc('arg', 0)^2
            sage: f(2)
            4.0
            sage: f = FastDoubleFunc('arg', 0)^4
            sage: f(2)
            16.0
            sage: f = FastDoubleFunc('arg', 0)^-3
            sage: f(2)
            0.125
            sage: f = FastDoubleFunc('arg', 0)^FastDoubleFunc('arg', 1)
            sage: f(5,3)
            125.0
        
        TESTS:
            sage: var('a,b')
            (a, b)
            sage: ff = (a^b)._fast_float_(a,b)
            sage: ff(2, 9)
            512.0
            sage: ff(-2, 9)
            -512.0
            sage: ff(-2, 9.1)
            Traceback (most recent call last):
            ...
            ValueError: negative number to a fractional power not real
        """
        if isinstance(right, FastDoubleFunc) and right.nargs == 0:
            right = float(right)
        if not isinstance(right, FastDoubleFunc):
            if right == int(float(right)):
                if right == 1:
                    return left
                elif right == 2:
                    return left.unop(DUP).unop(MUL)
                elif right == 3:
                    return left.unop(DUP).unop(DUP).unop(MUL).unop(MUL)
                elif right == 4:
                    return left.unop(DUP).unop(MUL).unop(DUP).unop(MUL)
                elif right < 0:
                    return (~left)**(-right)
            right = FastDoubleFunc('const', right)
        cdef FastDoubleFunc feval = binop(left, right, POW)
        return feval
        
    def __neg__(FastDoubleFunc self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = -fast_float_arg(0)
            sage: f(3.5)
            -3.5
        """
        return self.unop(NEG)
        
    def __abs__(FastDoubleFunc self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = abs(fast_float_arg(0))
            sage: f(-3)
            3.0
        """
        return self.unop(ABS)
        
    def __float__(self):
        """
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_constant, fast_float_arg
            sage: ff = fast_float_constant(17)
            sage: float(ff)
            17.0
            sage: ff = fast_float_constant(17) - fast_float_constant(2)^2
            sage: float(ff)
            13.0
            sage: ff = fast_float_constant(17) - fast_float_constant(2)^2 + fast_float_arg(1)
            sage: float(ff)
            Traceback (most recent call last):
            ...
            TypeError: Not a constant.
        """
        if self.nargs == 0:
            return self._call_c(NULL)
        else:
            raise TypeError, "Not a constant."
        
    def abs(FastDoubleFunc self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).abs()
            sage: f(3)
            3.0
        """
        return self.unop(ABS)
        
    def __invert__(FastDoubleFunc self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = ~fast_float_arg(0)
            sage: f(4)
            0.25
        """
        return self.unop(INVERT)
        
    def sqrt(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).sqrt()
            sage: f(4)
            2.0
        """
        return self.cfunc(&sqrt)
        
    ###################################################################
    #   Basic Comparison
    ###################################################################

    def _richcmp_(left, right, op):
        """
        Compare left and right. 
        EXAMPLES: 
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: import operator
            sage: f = fast_float_arg(0)._richcmp_(2, operator.lt)
            sage: [f(i) for i in (1..3)]
            [1.0, 0.0, 0.0]
            sage: f = fast_float_arg(0)._richcmp_(2, operator.le)
            sage: [f(i) for i in (1..3)]
            [1.0, 1.0, 0.0]
            sage: f = fast_float_arg(0)._richcmp_(2, operator.eq)
            sage: [f(i) for i in (1..3)]
            [0.0, 1.0, 0.0]
            sage: f = fast_float_arg(0)._richcmp_(2, operator.ne)
            sage: [f(i) for i in (1..3)]
            [1.0, 0.0, 1.0]
            sage: f = fast_float_arg(0)._richcmp_(2, operator.ge)
            sage: [f(i) for i in (1..3)]
            [0.0, 1.0, 1.0]
            sage: f = fast_float_arg(0)._richcmp_(2, operator.gt)
            sage: [f(i) for i in (1..3)]
            [0.0, 0.0, 1.0]
        """
        import operator
        if op == operator.lt:  #<
            return binop(left, right, LT)
        elif op == operator.eq: #==
            return binop(left, right, EQ)
        elif op == operator.gt: #>
            return binop(left, right, GT)
        elif op == operator.le: #<=
            return binop(left, right, LE)
        elif op == operator.ne: #!=
            return binop(left, right, NE)
        elif op == operator.ge: #>=
            return binop(left, right, GE)


    ###################################################################
    #   Exponential and log
    ###################################################################
        
    def log(self, base=None):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).log()
            sage: f(2)
            0.693147180559945...
            sage: f = fast_float_arg(0).log(2)
            sage: f(2)
            1.0
            sage: f = fast_float_arg(0).log(3)
            sage: f(9)
            2.0...
        """
        if base is None:
            return self.cfunc(&log)
        elif base == 2:
            return self.cfunc(&log2)
        elif base == 10:
            return self.cfunc(&log10)
        else:
            try:
                base = fast_float_constant(log(float(base)))
            except TypeError, e:
                base = fast_float(base.log())
            return binop(self.cfunc(&log), base, DIV)
        
    def exp(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).exp()
            sage: f(1)
            2.718281828459045...
            sage: f(100)
            2.6881171418161356e+43
        """
        return self.cfunc(&exp)
        
    ###################################################################
    #   Rounding
    ###################################################################

    def ceil(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).ceil()
            sage: f(1.5)
            2.0
            sage: f(-1.5)
            -1.0
        """
        return self.cfunc(&ceil)
        
    def floor(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).floor()
            sage: f(11.5)
            11.0
            sage: f(-11.5)
            -12.0
        """
        return self.cfunc(&floor)
        
    ###################################################################
    #   Trigonometric
    ###################################################################

    def sin(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).sin()
            sage: f(pi/2)
            1.0
        """
        return self.cfunc(&sin)
        
    def cos(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).cos()
            sage: f(0)
            1.0
        """
        return self.cfunc(&cos)
        
    def tan(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).tan()
            sage: f(pi/3)
            1.73205080756887...
        """
        return self.cfunc(&tan)
        
    def csc(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).csc()
            sage: f(pi/2)
            1.0
        """
        return ~self.sin()
        
    def sec(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).sec()
            sage: f(pi)
            -1.0
        """
        return ~self.cos()
        
    def cot(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).cot()
            sage: f(pi/4)
            1.0...
        """
        return ~self.tan()
        
    def arcsin(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).arcsin()
            sage: f(0.5)
            0.523598775598298...
        """
        return self.cfunc(&asin)
        
    def arccos(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).arccos()
            sage: f(sqrt(3)/2)
            0.5235987755982989...
        """
        return self.cfunc(&acos)
        
    def arctan(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).arctan()
            sage: f(1)
            0.785398163397448...
        """
        return self.cfunc(&atan)
        
    ###################################################################
    #   Hyperbolic
    ###################################################################

    def sinh(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).sinh()
            sage: f(log(2))
            0.75
        """
        return self.cfunc(&sinh)
        
    def cosh(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).cosh()
            sage: f(log(2))
            1.25
        """
        return self.cfunc(&cosh)
        
    def tanh(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).tanh()
            sage: f(0)
            0.0
        """
        return self.cfunc(&tanh)
        
    def arcsinh(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).arcsinh()
            sage: f(sinh(5))
            5.0
        """
        return self.cfunc(&asinh)
        
    def arccosh(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).arccosh()
            sage: f(cosh(5))
            5.0
        """
        return self.cfunc(&acosh)
        
    def arctanh(self):
        """
        EXAMPLE:
            sage: from sage.ext.fast_eval import fast_float_arg
            sage: f = fast_float_arg(0).arctanh()
            sage: abs(f(tanh(0.5)) - 0.5) < 0.0000001
            True
        """
        return self.cfunc(&atanh)
        
    cdef FastDoubleFunc cfunc(FastDoubleFunc self, void* func):
        cdef FastDoubleFunc feval = self.unop(ONE_ARG_FUNC)
        feval.ops[feval.nops - 1].params.func = func
        feval.allocate_stack()
        return feval
        
    ###################################################################
    #   Utility functions
    ###################################################################

    cdef FastDoubleFunc unop(FastDoubleFunc self, char type):
        cdef FastDoubleFunc feval = PY_NEW(FastDoubleFunc)
        feval.nargs = self.nargs
        feval.nops = self.nops + 1
        feval.max_height = self.max_height
        if type == DUP:
            feval.max_height += 1
        feval.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op) * feval.nops)
        memcpy(feval.ops, self.ops, sizeof(fast_double_op) * self.nops)
        feval.ops[feval.nops - 1].type = type
        feval.py_funcs = self.py_funcs
        feval.allocate_stack()
        return feval

cdef FastDoubleFunc binop(_left, _right, char type):
    r"""
    Returns a function that calculates left and right on the stack, leaving 
    their results on the top, and then calls operation \code{type}. 
    
    EXAMPLES:
        sage: from sage.ext.fast_eval import fast_float_arg
        sage: f = fast_float_arg(1)
        sage: g = fast_float_arg(2) * 11
        sage: f.op_list()
        ['load 1']
        sage: g.op_list()
        ['load 2', 'push 11.0', 'mul']
        sage: (f+g).op_list()
        ['load 1', 'load 2', 'push 11.0', 'mul', 'add']
        
    Correctly calculates the maximum stack heights and number of arguments:
        sage: f.max_height
        1
        sage: g.max_height
        2
        sage: (f+g).max_height
        3
        sage: (g+f).max_height
        2
        
        sage: f.nargs
        2
        sage: g.nargs
        3
        sage: (f+g).nargs
        3
    """
    cdef FastDoubleFunc left, right
    try:
        left = _left
    except TypeError:
        left = fast_float(_left)
    try:
        right = _right
    except TypeError:
        right = fast_float(_right)

    # In Cython assigning None does NOT raise a TypeError above. 
    if left is None or right is None:
        raise TypeError
    
    cdef FastDoubleFunc feval = PY_NEW(FastDoubleFunc)
    feval.nargs = max(left.nargs, right.nargs)
    feval.nops = left.nops + right.nops + 1
    feval.max_height = max(left.max_height, right.max_height+1)
    feval.ops = <fast_double_op *>sage_malloc(sizeof(fast_double_op) * feval.nops)
    memcpy(feval.ops, left.ops, sizeof(fast_double_op) * left.nops)
    memcpy(feval.ops + left.nops, right.ops, sizeof(fast_double_op) * right.nops)
    feval.ops[feval.nops - 1].type = type
    if left.py_funcs is None:
        feval.py_funcs = right.py_funcs
    elif right.py_funcs is None:
        feval.py_funcs = left.py_funcs
    else:
        feval.py_funcs = left.py_funcs + right.py_funcs
    feval.allocate_stack()
    return feval
    
    
def fast_float_constant(x):
    """
    Return a fast-to-evaluate constant function. 
    
    EXAMPLES: 
        sage: from sage.ext.fast_eval import fast_float_constant
        sage: f = fast_float_constant(-2.75)
        sage: f()
        -2.75

    This is all that goes on under the hood: 
        sage: fast_float_constant(pi).op_list()
        ['push 3.14159265359']
    """
    return FastDoubleFunc('const', x)
    
def fast_float_arg(n):
    """
    Return a fast-to-evaluate argument selector. 
    
    INPUT: 
        n -- the (zero-indexed) argument to select
    
    EXAMPLES: 
        sage: from sage.ext.fast_eval import fast_float_arg
        sage: f = fast_float_arg(0)
        sage: f(1,2)
        1.0
        sage: f = fast_float_arg(1)
        sage: f(1,2)
        2.0
        
    This is all that goes on under the hood: 
        sage: fast_float_arg(10).op_list()
        ['load 10']
    """
    return FastDoubleFunc('arg', n)
    
def fast_float_func(f, *args):
    """
    Returns a wrapper around a python function. 
    
    INPUT: 
        f -- a callable python object
        args -- a list of FastDoubleFunc inputs
    
    EXAMPLES: 
        sage: from sage.ext.fast_eval import fast_float_func, fast_float_arg
        sage: f = fast_float_arg(0)
        sage: g = fast_float_arg(1)
        sage: h = fast_float_func(lambda x,y: x-y, f, g)
        sage: h(5, 10)
        -5.0

    This is all that goes on under the hood: 
        sage: h.op_list() # random memory address
        ['load 0', 'load 1', 'py_call <function <lambda> at 0xb62b230>(2)']
    """
    return FastDoubleFunc('callable', f, *args)


new_fast_float=True

def fast_float(f, *vars, old=None, expect_one_var=False):
    """
    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
        expect_one_var -- don't give deprecation warning if vars is
                          omitted, as long as expression has only one var
        
    EXAMPLES: 
        sage: from sage.ext.fast_eval import fast_float
        sage: x,y = var('x,y')
        sage: f = fast_float(sqrt(x^2+y^2), 'x', 'y')
        sage: f(3,4)
        5.0

    Specifying the argument names is essential, as fast_float objects
    only distinguish between arguments by order. 
        sage: f = fast_float(x-y, 'x','y')
        sage: f(1,2)
        -1.0
        sage: f = fast_float(x-y, 'y','x')
        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, expect_one_var=expect_one_var) for x in f])
    
    cdef int i
    for i from 0 <= i < len(vars):
        if not PY_TYPE_CHECK(vars[i], str):
            v = str(vars[i])
            # inexact generators display as 1.00..0*x
            if '*' in v:
                v = v[v.index('*')+1:]
            vars = vars[:i] + (v,) + vars[i+1:]

    try:
        if old:
            return f._fast_float_(*vars)
        else:
            return fast_callable(f, vars=vars, domain=float, _autocompute_vars_for_backward_compatibility_with_deprecated_fast_float_functionality=True, expect_one_var=expect_one_var)
    except AttributeError:
        pass

    try:
        return FastDoubleFunc('const', float(f))
    except TypeError:
        pass

    try:
        from sage.symbolic.ring import SR
        return fast_float(SR(f), *vars)
    except TypeError:
        pass

    if f is None:
        raise TypeError, "no way to make fast_float from None"
        
    return f

def is_fast_float(x):
    return PY_TYPE_CHECK(x, FastDoubleFunc) or PY_TYPE_CHECK(x, Wrapper)