Ticket #7377: trac_7377fastcalculusrebased.patch
File trac_7377fastcalculusrebased.patch, 14.6 KB (added by , 11 years ago) 


sage/calculus/calculus.py
# HG changeset patch # User Nils Bruin <nbruin@sfu.ca> # Date 1265693495 28800 # Node ID b9db51bba0134fddc27b3af64ca95f987abb4423 # Parent d1def36a6fbc238838fca35740574e46fc6f8064 Fast calculus interface diff r d1def36a6fbc r b9db51bba013 sage/calculus/calculus.py
a b 386 386 f1 387 387 """ 388 388 maxima = sage.interfaces.maxima_lib.maxima 389 #maxima = Maxima(init_code = ['display2d:false', 'domain: complex', 390 # 'keepfloat: true', 'load(to_poly_solver)', 'load(simplify_sum)'], 391 # script_subdirectory=None) 389 392 390 393 ######################################################## 391 394 def symbolic_sum(expression, v, a, b, algorithm='maxima'): … … 529 532 raise ValueError, "summation limits must not depend on the summation variable" 530 533 531 534 if algorithm == 'maxima': 532 sum = "'sum(%s, %s, %s, %s)" % tuple([repr(expr._maxima_()) for expr in (expression, v, a, b)])533 535 try: 534 result = maxima.simplify_sum(sum) 535 result = result.ratsimp() 536 return expression.parent()(result) 536 return maxima.sr_sum(expression,v,a,b) 537 537 except TypeError, error: 538 538 s = str(error) 539 539 if "divergent" in s or 'Pole encountered' in s: … … 1101 1101 1102 1102 if algorithm == 'maxima': 1103 1103 if dir is None: 1104 l = ex._maxima_().limit(v, a)1104 l = maxima.sr_limit(ex, v, a) 1105 1105 elif dir in ['plus', '+', 'right', 'above']: 1106 1106 if dir == 'above': 1107 1107 from sage.misc.misc import deprecation 1108 1108 deprecation("the keyword 'above' is deprecated. Please use 'right' or '+' instead.", 'Sage version 4.6') 1109 l = ex._maxima_().limit(v, a, 'plus')1109 l = maxima.sr_limit(ex, v, a, 'plus') 1110 1110 elif dir in ['minus', '', 'left', 'below']: 1111 1111 if dir == 'below': 1112 1112 from sage.misc.misc import deprecation 1113 1113 deprecation("the keyword 'below' is deprecated. Please use 'left' or '' instead.", 'Sage version 4.6') 1114 l = ex._maxima_().limit(v, a, 'minus')1114 l = maxima.sr_limit(ex, v, a, 'minus') 1115 1115 elif algorithm == 'maxima_taylor': 1116 1116 if dir is None: 1117 l = ex._maxima_().tlimit(v, a)1117 l = maxima.sr_tlimit(ex, v, a) 1118 1118 elif dir == 'plus' or dir == 'above' or dir == 'from_right': 1119 l = ex._maxima_().tlimit(v, a, 'plus')1119 l = maxima.sr_tlimit(ex, v, a, 'plus') 1120 1120 elif dir == 'minus' or dir == 'below' or dir == 'from_left': 1121 l = ex._maxima_().tlimit(v, a, 'minus')1121 l = maxima.sr_tlimit(ex, v, a, 'minus') 1122 1122 elif algorithm == 'sympy': 1123 1123 if dir is None: 1124 1124 import sympy 
sage/interfaces/maxima.py
diff r d1def36a6fbc r b9db51bba013 sage/interfaces/maxima.py
a b 1008 1008 """ 1009 1009 return maxima_version() 1010 1010 1011 ##some helper functions to wrap tha calculus use of the maxima interface. 1012 ##these routines expect arguments living in the symbolic ring and return something 1013 ##that is hopefully coercible into the symbolic ring again. 1011 1014 1015 def sr_integral(self,*args): 1016 return args[0]._maxima_().integrate(*args[1:]) 1017 1018 def sr_sum(self,expression,v,a,b): 1019 sum = "'sum(%s, %s, %s, %s)" % tuple([repr(expr._maxima_()) for expr in (expression, v, a, b)]) 1020 result = self.simplify_sum(sum) 1021 result = result.ratsimp() 1022 return expression.parent()(result) 1023 1024 def sr_limit(self,ex,*args): 1025 return ex._maxima_().limit(*args) 1026 1027 def sr_tlimit(self,ex,*args): 1028 return ex._maxima_().tlimit(*args) 1029 1012 1030 ## def display2d(self, flag=True): 1013 1031 ## """ 1014 1032 ## Set the flag that determines whether Maxima objects are … … 1074 1092 def __doctest_cleanup(): 1075 1093 import sage.interfaces.quit 1076 1094 sage.interfaces.quit.expect_quitall() 1077 1078 
sage/interfaces/maxima_lib.py
diff r d1def36a6fbc r b9db51bba013 sage/interfaces/maxima_lib.py
a b 470 470 import sage.server.support 471 471 472 472 import maxima_abstract 473 from maxima_abstract import MaximaFunctionElement, MaximaExpectFunction, Maxima Element, MaximaFunction, maxima_console473 from maxima_abstract import MaximaFunctionElement, MaximaExpectFunction, MaximaFunction, maxima_console 474 474 475 475 # The Maxima "apropos" command, e.g., apropos(det) gives a list 476 476 # of all identifiers that begin in a certain way. This could … … 485 485 ecl_eval("(inpackage :maxima)") 486 486 ecl_eval("(setq $nolabels t))") 487 487 ecl_eval("(defvar *MAXIMALANGSUBDIR* NIL)") 488 ecl_eval("(setlocalesubdir)") 488 489 ecl_eval("(setpathnames)") 489 490 ecl_eval("(defun addlineinfo (x) x)") 490 491 ecl_eval("""(defun retrieve (msg flag &aux (print? nil))(error (concatenate 'string "Maxima asks:" (meval (list '($string) msg)))))""") … … 527 528 528 529 maxprint=EclObject("$STRING") 529 530 meval=EclObject("MEVAL") 531 msetq=EclObject("MSETQ") 532 cadadr=EclObject("CADADR") 533 534 max_integrate=EclObject("$INTEGRATE") 535 max_sum=EclObject("$SUM") 536 max_simplify_sum=EclObject("$SIMPLIFY_SUM") 537 max_ratsimp=EclObject("$RATSIMP") 538 max_limit=EclObject("$LIMIT") 539 max_tlimit=EclObject("$TLIMIT") 540 max_plus=EclObject("$PLUS") 541 max_minus=EclObject("$MINUS") 530 542 531 543 def max_to_string(s): 532 544 return meval(EclObject([[maxprint],s])).python()[1:1] … … 575 587 for l in init_code: 576 588 ecl_eval("#$%s$"%l) 577 589 ecl_eval("(setf *standardoutput* originalstandardoutput)") 590 591 def _coerce_from_special_method(self, x): 592 if isinstance(x, EclObject): 593 return MaximaElement(self,self._create(x)) 594 else: 595 return maxima_abstract.Maxima._coerce_from_special_method(self,x) 596 578 597 579 598 def _function_class(self): 580 599 """ … … 929 948 """ 930 949 s = self._eval_line('%s;'%var) 931 950 return s 932 951 952 def _create(self, value, name=None): 953 name = self._next_var_name() if name is None else name 954 if isinstance(value,EclObject): 955 maxima_eval([[msetq],cadadr("#$%s$#$"%name),value]) 956 else: 957 self.set(name, value) 958 return name 959 933 960 def version(self): 934 961 """ 935 962 Return the version of Maxima that Sage includes. … … 941 968 """ 942 969 return maxima_version() 943 970 971 ##some helper functions to wrap tha calculus use of the maxima interface. 972 ##these routines expect arguments living in the symbolic ring and return something 973 ##that is hopefully coercible into the symbolic ring again. 974 975 def sr_integral(self,*args): 976 try: 977 return max_to_sr(maxima_eval(([max_integrate],[sr_to_max(SR(a)) for a in args]))) 978 except RuntimeError, error: 979 s = str(error) 980 if "Divergent" in s or "divergent" in s: 981 raise ValueError, "Integral is divergent." 982 else: 983 raise error 984 985 def sr_sum(self,*args): 986 try: 987 return max_to_sr(maxima_eval([[max_ratsimp],[[max_simplify_sum],([max_sum],[sr_to_max(SR(a)) for a in args])]])); 988 except RuntimeError, error: 989 s = str(error) 990 if "divergent" in s: 991 raise ValueError, "Sum is divergent." 992 else: 993 raise error 994 995 def sr_limit(self,expr,v,a,dir=None): 996 L=[sr_to_max(SR(a)) for a in [expr,v,a]] 997 if dir == "plus": 998 L.append(max_plus) 999 elif dir == "minus": 1000 L.append(max_minus) 1001 return max_to_sr(maxima_eval(([max_limit],L))) 1002 1003 def sr_tlimit(self,expr,v,a,dir=None): 1004 L=[sr_to_max(SR(a)) for a in [expr,v,a]] 1005 if dir == "plus": 1006 L.append(max_plus) 1007 elif dir == "minus": 1008 L.append(max_minus) 1009 return max_to_sr(maxima_eval(([max_tlimit],L))) 944 1010 945 1011 ## def display2d(self, flag=True): 946 1012 ## """ … … 964 1030 ## """ 965 1031 ## self._display2d = bool(flag) 966 1032 1033 class MaximaElement(maxima_abstract.MaximaElement): 1034 """ 1035 """ 1036 def ecl(self): 1037 try: 1038 return self._ecl 1039 except AttributeError: 1040 self._ecl=maxima_eval("#$%s$"%self._name) 1041 return self._ecl 1042 967 1043 def is_MaximaElement(x): 968 1044 """ 969 1045 Returns True if x is of type MaximaElement. … … 1007 1083 sage.interfaces.quit.expect_quitall() 1008 1084 1009 1085 1086 import sage.symbolic.expression 1087 from sage.symbolic.ring import SR 1088 1089 import sage.symbolic.expression 1090 import sage.functions.trig 1091 import sage.functions.log 1092 import sage.functions.other 1093 1094 car=EclObject("car") 1095 cdr=EclObject("cdr") 1096 caar=EclObject("caar") 1097 cadadr=EclObject("cadadr") 1098 meval=EclObject("meval") 1099 NIL=EclObject("NIL") 1100 ratdisrep=EclObject("ratdisrep") 1101 1102 sage_op_dict={} 1103 1104 sage_op_dict = { 1105 sage.symbolic.expression.operator.abs : "MABS", 1106 sage.symbolic.expression.operator.add : "MPLUS", 1107 sage.symbolic.expression.operator.div : "MQUOTIENT", 1108 sage.symbolic.expression.operator.eq : "MEQUAL", 1109 sage.symbolic.expression.operator.ge : "MGEQP", 1110 sage.symbolic.expression.operator.gt : "MGREATERP", 1111 sage.symbolic.expression.operator.le : "MLEQP", 1112 sage.symbolic.expression.operator.lt : "MLESSP", 1113 sage.symbolic.expression.operator.mul : "MTIMES", 1114 sage.symbolic.expression.operator.ne : "MNOTEQUAL", 1115 sage.symbolic.expression.operator.neg : "MMINUS", 1116 sage.symbolic.expression.operator.pow : "MEXPT", 1117 sage.symbolic.expression.operator.or_ : "MOR", 1118 sage.symbolic.expression.operator.and_ : "MAND", 1119 sage.functions.trig.acos : "%ACOS", 1120 sage.functions.trig.acot : "%ACOT", 1121 sage.functions.trig.acsc : "%ACSC", 1122 sage.functions.trig.asec : "%ASEC", 1123 sage.functions.trig.asin : "%ASIN", 1124 sage.functions.trig.atan : "%ATAN", 1125 sage.functions.trig.cos : "%COS", 1126 sage.functions.trig.cot : "%COT", 1127 sage.functions.trig.csc : "%CSC", 1128 sage.functions.trig.sec : "%SEC", 1129 sage.functions.trig.sin : "%SIN", 1130 sage.functions.trig.tan : "%TAN", 1131 sage.functions.log.exp : "%EXP", 1132 sage.functions.log.ln : "%LOG", 1133 sage.functions.log.log : "%LOG", 1134 sage.functions.other.factorial : "MFACTORIAL", 1135 sage.functions.other.erf : "%ERF", 1136 } 1137 #we compile the dictionary 1138 sage_op_dict = dict([(k,EclObject(sage_op_dict[k])) for k in sage_op_dict]) 1139 max_op_dict = dict([(sage_op_dict[k],k) for k in sage_op_dict]) 1140 def add_vararg(*args): 1141 S=0 1142 for a in args: 1143 S=S+a 1144 return S 1145 1146 def mul_vararg(*args): 1147 P=1 1148 for a in args: 1149 P=P*a 1150 return P 1151 1152 def sage_rat(x,y): 1153 return x/y 1154 1155 mplus=EclObject("MPLUS") 1156 mtimes=EclObject("MTIMES") 1157 rat=EclObject("RAT") 1158 max_i=EclObject("$%I") 1159 max_op_dict[mplus]=add_vararg 1160 max_op_dict[mtimes]=mul_vararg 1161 max_op_dict[rat]=sage_rat 1162 1163 def mrat_to_sage(expr): 1164 r""" 1165 Convert a maxima MRAT expression to Sage SR 1166 1167 Maxima has an optimised representation for multivariate rational expressions. 1168 The easiest way to translate those to SR is by first asking maxima to give 1169 the generic representation of the object. That is what RATDISREP does in 1170 maxima. 1171 """ 1172 return max_to_sage(meval(EclObject([[ratdisrep],expr]))) 1173 1174 special_max_to_sage={ 1175 EclObject("MRAT") : mrat_to_sage 1176 } 1177 1178 sage_sym_dict={} 1179 max_sym_dict={} 1180 1181 def pyobject_to_max(obj): 1182 if isinstance(obj,sage.rings.rational.Rational): 1183 return EclObject(obj) if (obj.denom().is_one()) else EclObject([[rat], obj.numer(),obj.denom()]) 1184 elif isinstance(obj,sage.rings.number_field.number_field_element_quadratic.NumberFieldElement_quadratic) and obj.parent().defining_polynomial().list() == [1,0,1]: 1185 re, im = obj.list() 1186 return EclObject([[mplus], pyobject_to_max(re), [[mtimes], pyobject_to_max(im), max_i]]) 1187 1188 return EclObject(obj) 1189 1190 def sr_to_max(expr): 1191 r""" 1192 """ 1193 global sage_op_dict, max_op_dict 1194 global sage_sym_dict, max_sym_dict 1195 op = expr.operator() 1196 if op: 1197 if not (op in sage_op_dict): 1198 op_max=caar(maxima(expr).ecl()) 1199 sage_op_dict[op]=op_max 1200 max_op_dict[op_max]=op 1201 return EclObject(([sage_op_dict[op]], 1202 [sr_to_max(o) for o in expr.operands()])) 1203 elif expr._is_symbol() or expr._is_constant(): 1204 if not expr in sage_sym_dict: 1205 sym_max=maxima(expr).ecl() 1206 sage_sym_dict[expr]=sym_max 1207 max_sym_dict[sym_max]=expr 1208 return sage_sym_dict[expr] 1209 else: 1210 try: 1211 return pyobject_to_max(expr.pyobject()) 1212 except TypeError: 1213 return maxima(expr).ecl() 1214 1215 def max_to_sr(expr): 1216 if expr.consp(): 1217 op_max=caar(expr) 1218 if op_max in special_max_to_sage: 1219 return special_max_to_sage[op_max](expr) 1220 if not(op_max in max_op_dict): 1221 sage_expr=SR(maxima(expr)) 1222 max_op_dict[op_max]=sage_expr.operator() 1223 sage_op_dict[sage_expr.operator()]=op_max 1224 op=max_op_dict[op_max] 1225 max_args=cdr(expr) 1226 args=[] 1227 while not(max_args.nullp()): 1228 args.append(max_to_sr(car(max_args))) 1229 max_args=cdr(max_args) 1230 return op(*args) 1231 elif expr.symbolp(): 1232 if not(expr in max_sym_dict): 1233 sage_symbol=SR(maxima(expr)) 1234 sage_sym_dict[sage_symbol]=expr 1235 max_sym_dict[expr]=sage_symbol 1236 return max_sym_dict[expr] 1237 else: 1238 return expr.python() 1239 
sage/symbolic/integration/external.py
diff r d1def36a6fbc r b9db51bba013 sage/symbolic/integration/external.py
a b 12 12 sage: maxima_integrator(f(x), x) 13 13 integrate(f(x), x) 14 14 """ 15 from sage.calculus.calculus import maxima 15 16 if not isinstance(expression, Expression): 16 17 expression = SR(expression) 17 18 if a is None: 18 result = expression._maxima_().integrate(v)19 result = maxima.sr_integral(expression,v) 19 20 else: 20 21 try: 21 result = expression._maxima_().integrate(v, a, b)22 result = maxima.sr_integral(expression, v, a, b) 22 23 except TypeError, error: 23 24 s = str(error) 24 25 if "divergent" in s or 'Principal Value' in s: 25 26 raise ValueError, "Integral is divergent." 26 27 else: 27 28 raise 28 return result .sage()29 return result 29 30 30 31 def sympy_integrator(expression, v, a=None, b=None): 31 32 """