Ticket #7377: fastcalculus.patch

File fastcalculus.patch, 14.0 KB (added by nbruin, 9 years ago)

Patch to get faster calculus routines

  • sage/calculus/calculus.py

    # HG changeset patch
    # User Nils Bruin <nbruin@sfu.ca>
    # Date 1265693495 28800
    # Node ID f74ac05eae85bb1150ced5c460a6e6bc5acd521a
    # Parent  db5502ca296c8253c940ee9b9450d48234b61e46
    Fast calculus interface
    
    diff -r db5502ca296c -r f74ac05eae85 sage/calculus/calculus.py
    a b  
    329329arc_functions =  ['asin', 'acos', 'atan', 'asinh', 'acosh', 'atanh', 'acoth', 'asech', 'acsch', 'acot', 'acsc', 'asec']
    330330
    331331import sage.interfaces.maxima_lib
     332maxima = sage.interfaces.maxima_lib.maxima
    332333
    333 maxima = sage.interfaces.maxima_lib.maxima
     334#maxima = Maxima(init_code = ['display2d:false', 'domain: complex', 'keepfloat: true', 'load(to_poly_solver)', 'load(simplify_sum)'], script_subdirectory=None)
    334335
    335336########################################################
    336337def symbolic_sum(expression, v, a, b, algorithm='maxima'):
     
    466467        raise ValueError, "summation limits must not depend on the summation variable"
    467468
    468469    if algorithm == 'maxima':
    469         sum  = "'sum(%s, %s, %s, %s)" % tuple([repr(expr._maxima_()) for expr in (expression, v, a, b)])
    470470        try:
    471             result = maxima.simplify_sum(sum)
    472             result = result.ratsimp()
    473             return expression.parent()(result)
     471            return maxima.sr_sum(expression,v,a,b)
    474472        except TypeError, error:
    475473            s = str(error)
    476474            if "divergent" in s or 'Pole encountered' in s:
     
    819817
    820818    if algorithm == 'maxima':
    821819        if a is None:
    822             result = expression._maxima_().integrate(v)
     820            result = maxima.sr_integral(expression,v)
    823821        else:
    824822            try:
    825                 result = expression._maxima_().integrate(v, a, b)
     823                result = maxima.sr_integral(expression,v,a,b)
    826824            except TypeError, error:
    827825                s = str(error)
    828826                if "divergent" in s or 'Principal Value' in s:
     
    13681366
    13691367    if algorithm == 'maxima':
    13701368        if dir is None:
    1371             l = ex._maxima_().limit(v, a)
     1369            l = maxima.sr_limit(ex, v, a)
    13721370        elif dir == 'plus' or dir == 'above':
    1373             l = ex._maxima_().limit(v, a, 'plus')               
     1371            l = maxima.sr_limit(ex, v, a, 'plus')               
    13741372        elif dir == 'minus' or dir == 'below':
    1375             l = ex._maxima_().limit(v, a, 'minus')
     1373            l = maxima.sr_limit(ex, v, a, 'minus')
    13761374    elif algorithm == 'maxima_taylor':
    13771375        if dir is None:
    1378             l = ex._maxima_().tlimit(v, a)
     1376            l = maxima.sr_tlimit(ex, v, a)
    13791377        elif dir == 'plus' or dir == 'above':
    1380             l = ex._maxima_().tlimit(v, a, 'plus')
     1378            l = maxima.sr_tlimit(ex, v, a, 'plus')
    13811379        elif dir == 'minus' or dir == 'below':
    1382             l = ex._maxima_().tlimit(v, a, 'minus')
     1380            l = maxima.sr_tlimit(ex, v, a, 'minus')
    13831381    elif algorithm == 'sympy':
    13841382        if dir is None:
    13851383            import sympy
  • sage/interfaces/maxima.py

    diff -r db5502ca296c -r f74ac05eae85 sage/interfaces/maxima.py
    a b  
    944944        """
    945945        return maxima_version()
    946946
     947##some helper functions to wrap tha calculus use of the maxima interface.
     948##these routines expect arguments living in the symbolic ring and return something
     949##that is hopefully coercible into the symbolic ring again.
    947950
     951    def sr_integral(self,*args):
     952        return args[0]._maxima_().integrate(*args[1:])
     953
     954    def sr_sum(self,expression,v,a,b):
     955        sum  = "'sum(%s, %s, %s, %s)" % tuple([repr(expr._maxima_()) for expr in (expression, v, a, b)])
     956        result = self.simplify_sum(sum)
     957        result = result.ratsimp()
     958        return expression.parent()(result)
     959
     960    def sr_limit(self,ex,*args):
     961        return ex._maxima_().limit(*args)
     962
     963    def sr_tlimit(self,ex,*args):
     964        return ex._maxima_().tlimit(*args)
     965   
    948966##     def display2d(self, flag=True):
    949967##         """
    950968##         Set the flag that determines whether Maxima objects are
     
    10091027def __doctest_cleanup():
    10101028    import sage.interfaces.quit
    10111029    sage.interfaces.quit.expect_quitall()
    1012 
    1013 
  • sage/interfaces/maxima_lib.py

    diff -r db5502ca296c -r f74ac05eae85 sage/interfaces/maxima_lib.py
    a b  
    470470import sage.server.support
    471471
    472472import maxima_abstract
    473 from maxima_abstract import MaximaFunctionElement, MaximaExpectFunction, MaximaElement, MaximaFunction, maxima_console
     473from maxima_abstract import MaximaFunctionElement, MaximaExpectFunction, MaximaFunction, maxima_console
    474474
    475475# The Maxima "apropos" command, e.g., apropos(det) gives a list
    476476# of all identifiers that begin in a certain way.  This could
     
    527527
    528528maxprint=EclObject("$STRING")
    529529meval=EclObject("MEVAL")
     530msetq=EclObject("MSETQ")
     531cadadr=EclObject("CADADR")
     532
     533max_integrate=EclObject("$INTEGRATE")
     534max_sum=EclObject("$SUM")
     535max_simplify_sum=EclObject("$SIMPLIFY_SUM")
     536max_ratsimp=EclObject("$RATSIMP")
     537max_limit=EclObject("$LIMIT")
     538max_tlimit=EclObject("$TLIMIT")
     539max_plus=EclObject("$PLUS")
     540max_minus=EclObject("$MINUS")
    530541
    531542def max_to_string(s):
    532543     return meval(EclObject([[maxprint],s])).python()[1:-1]
     
    575586        for l in init_code:
    576587            ecl_eval("#$%s$"%l)
    577588        ecl_eval("(setf *standard-output* original-standard-output)")
     589 
     590    def _coerce_from_special_method(self, x):
     591        if isinstance(x, EclObject):
     592            return MaximaElement(self,self._create(x))
     593        else:
     594            return maxima_abstract.Maxima._coerce_from_special_method(self,x)
     595       
    578596       
    579597    def _function_class(self):
    580598        """
     
    929947        """
    930948        s = self._eval_line('%s;'%var)
    931949        return s
    932    
     950 
     951    def _create(self, value, name=None):
     952        name = self._next_var_name() if name is None else name
     953        if isinstance(value,EclObject):
     954            maxima_eval([[msetq],cadadr("#$%s$#$"%name),value])
     955        else:
     956            self.set(name, value)
     957        return name
     958               
    933959    def version(self):
    934960        """
    935961        Return the version of Maxima that Sage includes.
     
    941967        """
    942968        return maxima_version()
    943969
     970##some helper functions to wrap tha calculus use of the maxima interface.
     971##these routines expect arguments living in the symbolic ring and return something
     972##that is hopefully coercible into the symbolic ring again.
     973
     974    def sr_integral(self,*args):
     975        try:
     976            return max_to_sr(maxima_eval(([max_integrate],[sr_to_max(SR(a)) for a in args])))
     977        except RuntimeError, error:
     978            s = str(error)
     979            if "Divergent" in s or "divergent" in s:
     980                raise ValueError, "Integral is divergent."
     981            else:
     982                raise error
     983
     984    def sr_sum(self,*args):
     985        try:
     986            return max_to_sr(maxima_eval([[max_ratsimp],[[max_simplify_sum],([max_sum],[sr_to_max(SR(a)) for a in args])]]));
     987        except RuntimeError, error:
     988            s = str(error)
     989            if "divergent" in s:
     990                raise ValueError, "Sum is divergent."
     991            else:
     992                raise error
     993
     994    def sr_limit(self,expr,v,a,dir=None):
     995        L=[sr_to_max(SR(a)) for a in [expr,v,a]]
     996        if dir == "plus":
     997            L.append(max_plus)
     998        elif dir == "minus":
     999            L.append(max_minus)
     1000        return max_to_sr(maxima_eval(([max_limit],L)))
     1001
     1002    def sr_tlimit(self,expr,v,a,dir=None):
     1003        L=[sr_to_max(SR(a)) for a in [expr,v,a]]
     1004        if dir == "plus":
     1005            L.append(max_plus)
     1006        elif dir == "minus":
     1007            L.append(max_minus)
     1008        return max_to_sr(maxima_eval(([max_tlimit],L)))
    9441009
    9451010##     def display2d(self, flag=True):
    9461011##         """
     
    9641029##         """
    9651030##         self._display2d = bool(flag)
    9661031
     1032class MaximaElement(maxima_abstract.MaximaElement):
     1033    """
     1034    """   
     1035    def ecl(self):
     1036        try:
     1037            return self._ecl
     1038        except AttributeError:
     1039            self._ecl=maxima_eval("#$%s$"%self._name)
     1040            return self._ecl
     1041
    9671042def is_MaximaElement(x):
    9681043    """
    9691044    Returns True if x is of type MaximaElement.
     
    10071082    sage.interfaces.quit.expect_quitall()
    10081083
    10091084
     1085import sage.symbolic.expression
     1086from sage.symbolic.ring import SR
     1087
     1088import sage.symbolic.expression
     1089import sage.functions.trig
     1090import sage.functions.log
     1091import sage.functions.other
     1092
     1093car=EclObject("car")
     1094cdr=EclObject("cdr")
     1095caar=EclObject("caar")
     1096cadadr=EclObject("cadadr")
     1097meval=EclObject("meval")
     1098NIL=EclObject("NIL")
     1099ratdisrep=EclObject("ratdisrep")
     1100
     1101sage_op_dict={}
     1102
     1103sage_op_dict = {
     1104    sage.symbolic.expression.operator.abs : "MABS",
     1105    sage.symbolic.expression.operator.add : "MPLUS",
     1106    sage.symbolic.expression.operator.div : "MQUOTIENT",
     1107    sage.symbolic.expression.operator.eq : "MEQUAL",
     1108    sage.symbolic.expression.operator.ge : "MGEQP",
     1109    sage.symbolic.expression.operator.gt : "MGREATERP",
     1110    sage.symbolic.expression.operator.le : "MLEQP",
     1111    sage.symbolic.expression.operator.lt : "MLESSP",
     1112    sage.symbolic.expression.operator.mul : "MTIMES",
     1113    sage.symbolic.expression.operator.ne : "MNOTEQUAL",
     1114    sage.symbolic.expression.operator.neg : "MMINUS",
     1115    sage.symbolic.expression.operator.pow : "MEXPT",
     1116    sage.symbolic.expression.operator.or_ : "MOR",
     1117    sage.symbolic.expression.operator.and_ : "MAND",
     1118    sage.functions.trig.acos : "%ACOS",
     1119    sage.functions.trig.acot : "%ACOT",
     1120    sage.functions.trig.acsc : "%ACSC",
     1121    sage.functions.trig.asec : "%ASEC",
     1122    sage.functions.trig.asin : "%ASIN",
     1123    sage.functions.trig.atan : "%ATAN",
     1124    sage.functions.trig.cos : "%COS",
     1125    sage.functions.trig.cot : "%COT",
     1126    sage.functions.trig.csc : "%CSC",
     1127    sage.functions.trig.sec : "%SEC",
     1128    sage.functions.trig.sin : "%SIN",
     1129    sage.functions.trig.tan : "%TAN",
     1130    sage.functions.log.exp : "%EXP",
     1131    sage.functions.log.ln : "%LOG",
     1132    sage.functions.log.log : "%LOG",
     1133    sage.functions.other.factorial : "MFACTORIAL",
     1134    sage.functions.other.erf : "%ERF",
     1135}
     1136#we compile the dictionary
     1137sage_op_dict = dict([(k,EclObject(sage_op_dict[k])) for k in sage_op_dict])
     1138max_op_dict = dict([(sage_op_dict[k],k) for k in sage_op_dict])
     1139def add_vararg(*args):
     1140    S=0
     1141    for a in args:
     1142        S=S+a
     1143    return S
     1144
     1145def mul_vararg(*args):
     1146    P=1
     1147    for a in args:
     1148        P=P*a
     1149    return P
     1150
     1151def sage_rat(x,y):
     1152    return x/y
     1153
     1154mplus=EclObject("MPLUS")
     1155mtimes=EclObject("MTIMES")
     1156rat=EclObject("RAT")
     1157max_i=EclObject("$%I")
     1158max_op_dict[mplus]=add_vararg
     1159max_op_dict[mtimes]=mul_vararg
     1160max_op_dict[rat]=sage_rat
     1161
     1162def mrat_to_sage(expr):
     1163    r"""
     1164    Convert a maxima MRAT expression to Sage SR
     1165   
     1166    Maxima has an optimised representation for multivariate rational expressions.
     1167    The easiest way to translate those to SR is by first asking maxima to give
     1168    the generic representation of the object. That is what RATDISREP does in
     1169    maxima.
     1170    """
     1171    return max_to_sage(meval(EclObject([[ratdisrep],expr])))
     1172
     1173special_max_to_sage={
     1174    EclObject("MRAT") : mrat_to_sage
     1175}
     1176
     1177sage_sym_dict={}
     1178max_sym_dict={}
     1179
     1180def pyobject_to_max(obj):
     1181    if isinstance(obj,sage.rings.rational.Rational):
     1182        return EclObject(obj) if (obj.denom().is_one()) else EclObject([[rat], obj.numer(),obj.denom()])
     1183    elif isinstance(obj,sage.rings.number_field.number_field_element_quadratic.NumberFieldElement_quadratic) and obj.parent().defining_polynomial().list() == [1,0,1]:
     1184        re, im = obj.list()
     1185        return EclObject([[mplus], pyobject_to_max(re), [[mtimes], pyobject_to_max(im), max_i]])
     1186   
     1187    return EclObject(obj)
     1188
     1189def sr_to_max(expr):
     1190    r"""
     1191    """
     1192    global sage_op_dict, max_op_dict
     1193    global sage_sym_dict, max_sym_dict
     1194    op = expr.operator()
     1195    if op:
     1196        if not (op in sage_op_dict):
     1197            op_max=caar(maxima(expr).ecl())
     1198            sage_op_dict[op]=op_max
     1199            max_op_dict[op_max]=op
     1200        return EclObject(([sage_op_dict[op]],
     1201                     [sr_to_max(o) for o in expr.operands()]))
     1202    elif expr._is_symbol() or expr._is_constant():
     1203        if not expr in sage_sym_dict:
     1204            sym_max=maxima(expr).ecl()
     1205            sage_sym_dict[expr]=sym_max
     1206            max_sym_dict[sym_max]=expr
     1207        return sage_sym_dict[expr]
     1208    else:
     1209        try:
     1210            return pyobject_to_max(expr.pyobject())
     1211        except TypeError:
     1212            return maxima(expr).ecl()
     1213   
     1214def max_to_sr(expr):
     1215    if expr.consp():
     1216        op_max=caar(expr)
     1217        if op_max in special_max_to_sage:
     1218            return special_max_to_sage[op_max](expr)
     1219        if not(op_max in max_op_dict):
     1220            sage_expr=SR(maxima(expr))
     1221            max_op_dict[op_max]=sage_expr.operator()
     1222            sage_op_dict[sage_expr.operator()]=op_max
     1223        op=max_op_dict[op_max]
     1224        max_args=cdr(expr)
     1225        args=[]
     1226        while not(max_args.nullp()):
     1227            args.append(max_to_sr(car(max_args)))
     1228            max_args=cdr(max_args)
     1229        return op(*args)
     1230    elif expr.symbolp():
     1231        if not(expr in max_sym_dict):
     1232            sage_symbol=SR(maxima(expr))
     1233            sage_sym_dict[sage_symbol]=expr
     1234            max_sym_dict[expr]=sage_symbol
     1235        return max_sym_dict[expr]
     1236    else:
     1237        return expr.python()
     1238
  • sage/libs/ecl/ecl.pyx

    diff -r db5502ca296c -r f74ac05eae85 sage/libs/ecl/ecl.pyx
    a b  
    277277            return Ct
    278278        else:
    279279            return Cnil
    280     elif pyobj==None:
     280    elif pyobj is None:
    281281        return Cnil
    282282    elif isinstance(pyobj,int):
    283283        return ecl_make_integer(pyobj)
     
    993993
    994994#input: a cl-object. Output: EclObject wrapping that.
    995995cdef EclObject ecl_wrap(cl_object o):
    996     cdef EclObject obj
    997     obj = EclObject()
     996    cdef EclObject obj = EclObject.__new__(EclObject)
    998997    obj.set_obj(o)
    999998    return obj
    1000999