Context Navigation

source: sage/calculus/calculus.py @ 7429:868c3a8c5462

Visit:

Revision 7429:868c3a8c5462, 170.6 KB checked in by Mike Hansen <mhansen@…>, 5 years ago (diff)
Fixed #847.

Line
1	r"""
2	Symbolic Computation.
3
4	AUTHORS:
5	Bobby Moretti and William Stein: 2006--2007
6
7	The \sage calculus module is loosely based on the \sage Enhahcement Proposal
8	found at: http://www.sagemath.org:9001/CalculusSEP.
9
10	EXAMPLES:
11
12	The basic units of the calculus package are symbolic expressions
13	which are elements of the symbolic expression ring (SR). There are
14	many subclasses of SymbolicExpression. The most basic of these is
15	the formal indeterminate class, SymbolicVariable. To create a
16	SymbolicVariable object in \sage, use the var() method, whose
17	argument is the text of that variable. Note that \sage is
18	intelligent about {\LaTeX}ing variable names.
19
20	sage: x1 = var('x1'); x1
21	x1
22	sage: latex(x1)
23	x_{1}
24	sage: theta = var('theta'); theta
25	theta
26	sage: latex(theta)
27	\theta
28
29	\sage predefines upper and lowercase letters as global
30	indeterminates. Thus the following works:
31	sage: x^2
32	x^2
33	sage: type(x)
34	<class 'sage.calculus.calculus.SymbolicVariable'>
35
36	More complicated expressions in SAGE can be built up using
37	ordinary arithmetic. The following are valid, and follow the rules
38	of Python arithmetic: (The '=' operator represents assignment, and
39	not equality)
40	sage: var('x,y,z')
41	(x, y, z)
42	sage: f = x + y + z/(2sin(yz/55))
43	sage: g = f^f; g
44	(z/(2sin(yz/55)) + y + x)^(z/(2sin(yz/55)) + y + x)
45
46	Differentiation and integration are available, but behind the
47	scenes through maxima:
48
49	sage: f = sin(x)/cos(2*y)
50	sage: f.derivative(y)
51	2sin(x)sin(2y)/cos(2y)^2
52	sage: g = f.integral(x); g
53	-cos(x)/cos(2*y)
54
55	Note that these methods require an explicit variable name. If none
56	is given, \sage will try to find one for you.
57	sage: f = sin(x); f.derivative()
58	cos(x)
59
60	However when this is ambiguous, \sage will raise an exception:
61	sage: f = sin(x+y); f.derivative()
62	Traceback (most recent call last):
63	...
64	ValueError: must supply an explicit variable for an expression containing more than one variable
65
66	Substitution works similarly. We can substitute with a python dict:
67	sage: f = sin(x*y - z)
68	sage: f({x: var('t'), y: z})
69	sin(t*z - z)
70
71	Also we can substitute with keywords:
72	sage: f = sin(x*y - z)
73	sage: f(x = t, y = z)
74	sin(t*z - z)
75
76	If there is no ambiguity of variable names, we don't have to specify them:
77	sage: f = sin(x)
78	sage: f(y)
79	sin(y)
80	sage: f(pi)
81	0
82
83	However if there is ambiguity, we must explicitly state what
84	variables we're substituting for:
85
86	sage: f = sin(2pix/y)
87	sage: f(4)
88	sin(8*pi/y)
89
90	We can also make CallableSymbolicExpressions, which is a SymbolicExpression
91	that are functions of variables in a fixed order. Each
92	SymbolicExpression has a function() method used to create a
93	CallableSymbolicExpression.
94
95	sage: u = log((2-x)/(y+5))
96	sage: f = u.function(x, y); f
97	(x, y) \|--> log((2 - x)/(y + 5))
98
99	There is an easier way of creating a CallableSymbolicExpression, which
100	relies on the \sage preparser.
101
102	sage: f(x,y) = log(x)*cos(y); f
103	(x, y) \|--> log(x)*cos(y)
104
105	Then we have fixed an order of variables and there is no ambiguity
106	substituting or evaluating:
107
108	sage: f(x,y) = log((2-x)/(y+5))
109	sage: f(7,t)
110	log(-5/(t + 5))
111
112	Some further examples:
113
114	sage: f = 5*sin(x)
115	sage: f
116	5*sin(x)
117	sage: f(2)
118	5*sin(2)
119	sage: f(pi)
120	0
121	sage: float(f(pi)) # random low order bits
122	6.1232339957367663e-16
123
124	Another example:
125	sage: f = integrate(1/sqrt(9+x^2), x); f
126	asinh(x/3)
127	sage: f(3)
128	asinh(1)
129	sage: f.diff(x)
130	1/(3*sqrt(x^2/9 + 1))
131
132	We compute the length of the parabola from 0 to 2:
133	sage: x = var('x')
134	sage: y = x^2
135	sage: dy = diff(y,x)
136	sage: z = integral(sqrt(1 + dy^2), x, 0, 2)
137	sage: print z
138	asinh(4) + 4 sqrt(17)
139	---------------------
140	4
141	sage: n(z,200)
142	4.6467837624329358733826155674904591885104869874232887508703
143	sage: float(z)
144	4.6467837624329356
145
146	COERCION EXAMPLES:
147
148	We coerce various symbolic expressions into the complex numbers:
149
150	sage: CC(I)
151	1.00000000000000*I
152	sage: CC(2*I)
153	2.00000000000000*I
154	sage: ComplexField(200)(2*I)
155	2.0000000000000000000000000000000000000000000000000000000000*I
156	sage: ComplexField(200)(sin(I))
157	1.1752011936438014568823818505956008151557179813340958702296*I
158	sage: f = sin(I) + cos(I/2); f
159	I*sinh(1) + cosh(1/2)
160	sage: CC(f)
161	1.12762596520638 + 1.17520119364380*I
162	sage: ComplexField(200)(f)
163	1.1276259652063807852262251614026720125478471180986674836290 + 1.1752011936438014568823818505956008151557179813340958702296*I
164	sage: ComplexField(100)(f)
165	1.1276259652063807852262251614 + 1.1752011936438014568823818506*I
166
167	We illustrate construction of an inverse sum where each denominator
168	has a new variable name:
169	sage: f = sum(1/var('n%s'%i)^i for i in range(10))
170	sage: print f
171	1 1 1 1 1 1 1 1 1
172	--- + --- + --- + --- + --- + --- + --- + --- + -- + 1
173	9 8 7 6 5 4 3 2 n1
174	n9 n8 n7 n6 n5 n4 n3 n2
175
176	Note that after calling var, the variables are immediately available for use:
177	sage: (n1 + n2)^5
178	(n2 + n1)^5
179
180	We can, of course, substitute:
181	sage: print f(n9=9,n7=n6)
182	1 1 1 1 1 1 1 1 387420490
183	--- + --- + --- + --- + --- + --- + --- + -- + ---------
184	8 6 7 5 4 3 2 n1 387420489
185	n8 n6 n6 n5 n4 n3 n2
186
187	TESTS:
188	We test pickling:
189	sage: f = -sqrt(pi)*(x^3 + sin(x/cos(y)))
190	sage: bool(loads(dumps(f)) == f)
191	True
192
193	Substitution:
194	sage: f = x
195	sage: f(x=5)
196	5
197
198	The symbolic Calculus package uses its own copy of Maxima for
199	simplification, etc., which is separate from the default system-wide
200	version:
201	sage: maxima.eval('[x,y]: [1,2]')
202	'[1,2]'
203	sage: maxima.eval('expand((x+y)^3)')
204	'27'
205
206	If the copy of maxima used by the symbolic calculus package were
207	the same as the default one, then the following would return 27,
208	which would be very confusing indeed!
209	sage: x, y = var('x,y')
210	sage: expand((x+y)^3)
211	y^3 + 3xy^2 + 3x^2y + x^3
212
213	Set x to be 5 in maxima:
214	sage: maxima('x: 5')
215	5
216	sage: maxima('x + x + %pi')
217	%pi+10
218
219	This simplification is done using maxima (behind the scenes):
220	sage: x + x + pi
221	2*x + pi
222
223	Note that x is still x, since the maxima used by the calculus package
224	is different than the one in the interactive interpreter.
225
226	"""
227
228	import weakref
229
230	from sage.rings.all import (CommutativeRing, RealField, is_Polynomial,
231	is_MPolynomial, is_MPolynomialRing, is_FractionFieldElement,
232	is_RealNumber, is_ComplexNumber, RR,
233	Integer, Rational, CC,
234	QuadDoubleElement,
235	PolynomialRing, ComplexField,
236	algdep, Integer, RealNumber)
237
238	from sage.structure.element import RingElement, is_Element
239	from sage.structure.parent_base import ParentWithBase
240
241	import operator
242	from sage.misc.latex import latex, latex_variable_name
243	from sage.structure.sage_object import SageObject
244
245	from sage.interfaces.maxima import MaximaElement, Maxima
246
247	# The calculus package uses its own copy of maxima, which is
248	# separate from the default system-wide version.
249	maxima = Maxima()
250
251	from sage.misc.sage_eval import sage_eval
252
253	from sage.calculus.equations import SymbolicEquation
254	from sage.rings.real_mpfr import RealNumber
255	from sage.rings.complex_number import ComplexNumber
256	from sage.rings.real_double import RealDoubleElement
257	from sage.rings.complex_double import ComplexDoubleElement
258	from sage.rings.real_mpfi import RealIntervalFieldElement
259	from sage.rings.infinity import InfinityElement
260
261	from sage.libs.pari.gen import pari, PariError, gen as PariGen
262
263	from sage.rings.complex_double import ComplexDoubleElement
264
265	import sage.functions.constants
266
267	import math
268	import sage.functions.functions
269
270	is_simplified = False
271
272	infixops = {operator.add: '+',
273	operator.sub: '-',
274	operator.mul: '*',
275	operator.div: '/',
276	operator.pow: '^'}
277
278
279	def is_SymbolicExpression(x):
280	"""
281	EXAMPLES:
282	sage: is_SymbolicExpression(sin(x))
283	True
284	sage: is_SymbolicExpression(2/3)
285	False
286	sage: is_SymbolicExpression(sqrt(2))
287	True
288	"""
289	return isinstance(x, SymbolicExpression)
290
291	def is_SymbolicExpressionRing(x):
292	"""
293	EXAMPLES:
294	sage: is_SymbolicExpressionRing(QQ)
295	False
296	sage: is_SymbolicExpressionRing(SR)
297	True
298	"""
299	return isinstance(x, SymbolicExpressionRing_class)
300
301	cache = {}
302	class uniq(object):
303	def __new__(cls):
304	global cache
305	if cache.has_key(cls):
306	return cache[cls]
307	O = object.__new__(cls)
308	cache[cls] = O
309	return O
310
311	class SymbolicExpressionRing_class(CommutativeRing):
312	"""
313	The ring of all formal symbolic expressions.
314
315	EXAMPLES:
316	sage: SR
317	Symbolic Ring
318	sage: type(SR)
319	<class 'sage.calculus.calculus.SymbolicExpressionRing_class'>
320	"""
321	def __init__(self):
322	self._default_precision = 53 # default precision bits
323	ParentWithBase.__init__(self, RR)
324
325	def __cmp__(self, other):
326	return cmp(type(self), type(other))
327
328	def __reduce__(self):
329	return SymbolicExpressionRing, tuple([])
330
331	def __call__(self, x):
332	"""
333	Coerce x into the symbolic expression ring SR.
334
335	EXAMPLES:
336	sage: a = SR(-3/4); a
337	-3/4
338	sage: type(a)
339	<class 'sage.calculus.calculus.SymbolicConstant'>
340	sage: a.parent()
341	Symbolic Ring
342
343	If a is already in the symblic expression ring, coercing returns
344	a itself (not a copy):
345	sage: SR(a) is a
346	True
347
348	A Python complex number:
349	sage: SR(complex(2,-3))
350	2.00000000000000 - 3.00000000000000*I
351	"""
352	if is_Element(x) and x.parent() is self:
353	return x
354	elif hasattr(x, '_symbolic_'):
355	return x._symbolic_(self)
356	return self._coerce_impl(x)
357
358	def _coerce_impl(self, x):
359	"""
360	EXAMPLES:
361	sage: x=var('x'); y0,y1=PolynomialRing(ZZ,2,'y').gens()
362	sage: x+y0/y1
363	y0/y1 + x
364	sage: x.subs(x=y0/y1)
365	y0/y1
366	"""
367	if isinstance(x, CallableSymbolicExpression):
368	return x._expr
369	elif isinstance(x, SymbolicExpression):
370	return x
371	elif isinstance(x, MaximaElement):
372	return symbolic_expression_from_maxima_element(x)
373	elif is_Polynomial(x) or is_MPolynomial(x):
374	if x.base_ring() != self: # would want coercion to go the other way
375	return SymbolicPolynomial(x)
376	else:
377	raise TypeError, "Basering is Symbolic Ring, please coerce in the other direction."
378	elif is_FractionFieldElement(x) and (is_Polynomial(x.numerator()) or is_MPolynomial(x.numerator())):
379	if x.base_ring() != self: # would want coercion to go the other way
380	return SymbolicPolynomial(x.numerator()) / SymbolicPolynomial(x.denominator())
381	else:
382	raise TypeError, "Basering is Symbolic Ring, please coerce in the other direction."
383	elif isinstance(x, (RealNumber,
384	RealDoubleElement,
385	RealIntervalFieldElement,
386	float,
387	sage.functions.constants.Constant,
388	Integer,
389	int,
390	Rational,
391	PariGen,
392	ComplexNumber,
393	ComplexDoubleElement,
394	QuadDoubleElement,
395	InfinityElement
396	)):
397	return SymbolicConstant(x)
398	elif isinstance(x, complex):
399	return evaled_symbolic_expression_from_maxima_string('%s+%%i*%s'%(x.real,x.imag))
400	else:
401	raise TypeError, "cannot coerce type '%s' into a SymbolicExpression."%type(x)
402
403	def _repr_(self):
404	return 'Symbolic Ring'
405
406	def _latex_(self):
407	return 'SymbolicExpressionRing'
408
409	def var(self, x):
410	return var(x)
411
412	def characteristic(self):
413	return Integer(0)
414
415	def _an_element_impl(self):
416	return SymbolicVariable('_generic_variable_name_')
417
418	def is_field(self):
419	return True
420
421	def is_exact(self):
422	return False
423
424	# Define the unique symbolic expression ring.
425	SR = SymbolicExpressionRing_class()
426
427	# The factory function that returns the unique SR.
428	def SymbolicExpressionRing():
429	"""
430	Return the symbolic expression ring.
431
432	EXAMPLES:
433	sage: SymbolicExpressionRing()
434	Symbolic Ring
435	sage: SymbolicExpressionRing() is SR
436	True
437	"""
438	return SR
439
440	class SymbolicExpression(RingElement):
441	"""
442	A Symbolic Expression.
443
444	EXAMPLES:
445	Some types of SymbolicExpressions:
446
447	sage: a = SR(2+2); a
448	4
449	sage: type(a)
450	<class 'sage.calculus.calculus.SymbolicConstant'>
451
452	"""
453	def __init__(self):
454	RingElement.__init__(self, SR)
455	if is_simplified:
456	self._simp = self
457
458	def __hash__(self):
459	return hash(self._repr_(simplify=False))
460
461	def __nonzero__(self):
462	"""
463	EXAMPLES:
464	sage: k = var('k')
465	sage: pol = 1/(k-1) - 1/k -1/k/(k-1);
466	sage: pol.is_zero()
467	True
468
469	sage: f = sin(x)^2 + cos(x)^2 - 1
470	sage: f.is_zero()
471	True
472	"""
473
474	try:
475	return self.__nonzero
476	except AttributeError:
477	ans = not bool(self == SR.zero_element())
478	self.__nonzero = ans
479	return ans
480
481	def __str__(self):
482	"""
483	Printing an object explicitly gives ASCII art:
484
485	EXAMPLES:
486	sage: var('x y')
487	(x, y)
488	sage: f = y^2/(y+1)^3 + x/(x-1)^3
489	sage: f
490	y^2/(y + 1)^3 + x/(x - 1)^3
491	sage: print f
492	2
493	y x
494	-------- + --------
495	3 3
496	(y + 1) (x - 1)
497
498	"""
499	return self.display2d(onscreen=False)
500
501	def show(self):
502	from sage.misc.functional import _do_show
503	return _do_show(self)
504
505	def display2d(self, onscreen=True):
506	"""
507	Display self using ASCII art.
508
509	INPUT:
510	onscreen -- string (optional, default True) If True,
511	displays; if False, returns string.
512
513	EXAMPLES:
514	We display a fraction:
515	sage: var('x,y')
516	(x, y)
517	sage: f = (x^3+y)/(x+3*y^2+1); f
518	(y + x^3)/(3*y^2 + x + 1)
519	sage: print f
520	3
521	y + x
522	------------
523	2
524	3 y + x + 1
525
526	Use onscreen=False to get the 2d string:
527	sage: f.display2d(onscreen=False)
528	' 3\r\n y + x\r\n ------------\r\n 2\r\n 3 y + x + 1'
529
530	ASCII art is really helps for the following integral:
531	sage: f = integral(sin(x^2)); f
532	sqrt(pi)((sqrt(2)I + sqrt(2))erf((sqrt(2)I + sqrt(2))x/2) + (sqrt(2)I - sqrt(2))erf((sqrt(2)I - sqrt(2))*x/2))/8
533	sage: print f
534	(sqrt(2) I + sqrt(2)) x
535	sqrt( pi) ((sqrt(2) I + sqrt(2)) erf(------------------------)
536	2
537	(sqrt(2) I - sqrt(2)) x
538	+ (sqrt(2) I - sqrt(2)) erf(------------------------))/8
539	2
540	"""
541	if not self.is_simplified():
542	self = self.simplify()
543	s = self._maxima_().display2d(onscreen=False)
544	s = s.replace('%pi',' pi').replace('%i',' I').replace('%e', ' e')
545	if onscreen:
546	print s
547	else:
548	return s
549
550	def is_simplified(self):
551	return hasattr(self, '_simp') and self._simp is self
552
553	def _declare_simplified(self):
554	self._simp = self
555
556	def hash(self):
557	return hash(self._repr_(simplify=False))
558
559	def plot(self, args, *kwds):
560	from sage.plot.plot import plot
561
562	# see if the user passed a variable in.
563	if kwds.has_key('param'):
564	param = kwds['param']
565	else:
566	param = None
567	for i in range(len(args)):
568	if isinstance(args[i], SymbolicVariable):
569	param = args[i]
570	args = args[:i] + args[i+1:]
571	break
572
573	F = self.simplify()
574	if isinstance(F, Symbolic_object):
575	if hasattr(F._obj, '__call__'):
576	f = lambda x: F._obj(x)
577	else:
578	y = float(F._obj)
579	f = lambda x: y
580
581	elif param is None:
582	if isinstance(self, CallableSymbolicExpression):
583	A = self.arguments()
584	if len(A) == 0:
585	raise ValueError, "function has no input arguments"
586	else:
587	param = A[0]
588	f = lambda x: self(x)
589	else:
590	A = self.variables()
591	if len(A) == 0:
592	y = float(self)
593	f = lambda x: y
594	else:
595	param = A[0]
596	f = self.function(param)
597	else:
598	f = self.function(param)
599	return plot(f, args, *kwds)
600
601	def __lt__(self, right):
602	try:
603	return SymbolicEquation(self, SR(right), operator.lt)
604	except TypeError:
605	return False
606
607	def __le__(self, right):
608	try:
609	return SymbolicEquation(self, SR(right), operator.le)
610	except TypeError:
611	return False
612
613	def __eq__(self, right):
614	try:
615	return SymbolicEquation(self, SR(right), operator.eq)
616	except TypeError:
617	return False
618
619	def __ne__(self, right):
620	try:
621	return SymbolicEquation(self, SR(right), operator.ne)
622	except TypeError:
623	return False
624
625	def __ge__(self, right):
626	try:
627	return SymbolicEquation(self, SR(right), operator.ge)
628	except TypeError:
629	return False
630
631	def __gt__(self, right):
632	try:
633	return SymbolicEquation(self, SR(right), operator.gt)
634	except TypeError:
635	return False
636
637	def __cmp__(self, right):
638	"""
639	Compares self and right.
640
641	This is by definition the comparison of the underlying Maxima
642	objects.
643
644	EXAMPLES:
645	These two are equal:
646	sage: cmp(e+e, e*2)
647	0
648	"""
649	return cmp(maxima(self), maxima(right))
650
651	def _richcmp_(left, right, op):
652	"""
653	TESTS:
654	sage: 3 < x
655	3 < x
656	sage: 3 <= x
657	3 <= x
658	sage: 3 == x
659	3 == x
660	sage: 3 >= x
661	3 >= x
662	sage: 3 > x
663	3 > x
664	"""
665	if op == 0: #<
666	return left < right
667	elif op == 2: #==
668	return left == right
669	elif op == 4: #>
670	return left > right
671	elif op == 1: #<=
672	return left <= right
673	elif op == 3: #!=
674	return left != right
675	elif op == 5: #>=
676	return left >= right
677
678
679	def _neg_(self):
680	"""
681	Return the formal negative of self.
682
683	EXAMPLES:
684	sage: var('a,x,y')
685	(a, x, y)
686	sage: -a
687	-a
688	sage: -(x+y)
689	-y - x
690	"""
691	return SymbolicArithmetic([self], operator.neg)
692
693	##################################################################
694	# Coercions to interfaces
695	##################################################################
696	# The maxima one is special:
697	def _maxima_(self, session=None):
698	if session is None:
699	return RingElement._maxima_(self, maxima)
700	else:
701	return RingElement._maxima_(self, session)
702
703	def _maxima_init_(self):
704	return self._repr_(simplify=False)
705
706
707	# The others all go through _sys_init_, which is defined below,
708	# and does all interfaces in a unified manner.
709
710	def _axiom_init_(self):
711	return self._sys_init_('axiom')
712
713	def _gp_init_(self):
714	return self._sys_init_('pari') # yes, gp goes through pari
715
716	def _maple_init_(self):
717	return self._sys_init_('maple')
718
719	def _magma_init_(self):
720	return '"%s"'%self.str()
721
722	def _kash_init_(self):
723	return self._sys_init_('kash')
724
725	def _macaulay2_init_(self):
726	return self._sys_init_('macaulay2')
727
728	def _mathematica_init_(self):
729	return self._sys_init_('mathematica')
730
731	def _octave_init_(self):
732	return self._sys_init_('octave')
733
734	def _pari_init_(self):
735	return self._sys_init_('pari')
736
737	def _sys_init_(self, system):
738	return repr(self)
739
740	##################################################################
741	# These systems have no symbolic or numerical capabilities at all,
742	# really, so we always just coerce to a string
743	##################################################################
744	def _gap_init_(self):
745	"""
746	Conversion of symbolic object to GAP always results in a GAP string.
747
748	EXAMPLES:
749	sage: gap(e+pi^2 + x^3)
750	"x^3 + pi^2 + e"
751	"""
752	return '"%s"'%repr(self)
753
754	def _singular_init_(self):
755	"""
756	Conversion of symbolic object to Singular always results in a Singular string.
757
758	EXAMPLES:
759	sage: singular(e+pi^2 + x^3)
760	x^3 + pi^2 + e
761	"""
762	return '"%s"'%repr(self)
763
764	##################################################################
765	# Non-canonical coercions to compiled built-in rings and fields
766	##################################################################
767	def __int__(self):
768	"""
769	EXAMPLES:
770	sage: int(sin(2)*100)
771	90
772	"""
773	try:
774	return int(repr(self))
775	except (ValueError, TypeError):
776	return int(float(self))
777
778	def __long__(self):
779	"""
780	EXAMPLES:
781	sage: long(sin(2)*100)
782	90L
783	"""
784	return long(int(self))
785
786	def numerical_approx(self, prec=None, digits=None):
787	r"""
788	Return a numerical approximation of self as either a real or
789	complex number with at least the requested number of bits or
790	digits of precision.
791
792	NOTE: You can use \code{foo.n()} as a shortcut for
793	\code{foo.numerical_approx()}.
794
795	INPUT:
796	prec -- an integer: the number of bits of precision
797	digits -- an integer: digits of precision
798
799	OUTPUT:
800	A RealNumber or ComplexNumber approximation of self with
801	prec bits of precision.
802
803	EXAMPLES:
804	sage: cos(3).numerical_approx()
805	-0.989992496600445
806
807	Use the n() shortcut:
808	sage: cos(3).n()
809	-0.989992496600445
810
811	Higher precision:
812	sage: cos(3).numerical_approx(200)
813	-0.98999249660044545727157279473126130239367909661558832881409
814	sage: numerical_approx(cos(3), digits=10)
815	-0.9899924966
816	sage: (i + 1).numerical_approx(32)
817	1.00000000 + 1.00000000*I
818	sage: (pi + e + sqrt(2)).numerical_approx(100)
819	7.2740880444219335226246195788
820	"""
821	if prec is None:
822	if digits is None:
823	prec = 53
824	else:
825	prec = int(digits * 3.4) + 2
826
827	# make sure the field is of the right precision
828	field = RealField(prec)
829
830	try:
831	approx = self._mpfr_(field)
832	except TypeError:
833	# try to return a complex result
834	approx = self._complex_mpfr_field_(ComplexField(prec))
835
836	return approx
837
838	n = numerical_approx
839
840	def minpoly(self, bits=None, degree=None, epsilon=0):
841	"""
842	Return the minimal polynomial of self, if possible.
843
844	INPUT:
845	bits -- the number of bits to use in numerical approx
846	degree -- the expected algebraic degree
847	epsilon -- return without error as long as f(self) < epsilon,
848	in the case that the result cannot be proven.
849
850	All of the above parameters are optional, with epsilon=0,
851	bits and degree tested up to 1000 and 24 by default respectively.
852	If these are known, it will be faster to give them explicitly.
853
854	OUTPUT:
855	The minimal polynomial of self. This is proved symbolically if
856	epsilon=0 (default).
857
858	If the minimal polynomial could not be found, two distinct kinds
859	of errors are raised. If no reasonable candidate was found with
860	the given bit/degree parameters, a ValueError will be raised.
861	If a reasonable candidate was found but (perhaps due to limits
862	in the underlying symbolic package) was unable to be proved
863	correct, a NotImplementedError will be raised.
864
865	ALGORITHM:
866	Use the PARI algdep command on a numerical approximation of self
867	to get a candidate minpoly f. Approximate f(self) to higher
868	precision and if the result is still close enough to 0 then
869	evaluate f(self) symbolically, attempting to prove vanishing.
870	If this fails, and epsilon is non-zero, return f as long as
871	f(self) < epsilon. Otherwise raise an error.
872
873
874	NOTE: Failure of this function does not prove self is
875	not algebraic.
876
877
878	EXAMPLES:
879
880	First some simple examples:
881	sage: sqrt(2).minpoly()
882	x^2 - 2
883	sage: a = 2^(1/3)
884	sage: a.minpoly()
885	x^3 - 2
886	sage: (sqrt(2)-3^(1/3)).minpoly()
887	x^6 - 6x^4 + 6x^3 + 12x^2 + 36x + 1
888
889	Sometimes it fails.
890	sage: sin(1).minpoly()
891	Traceback (most recent call last):
892	...
893	ValueError: Could not find minimal polynomial (1000 bits, degree 24).
894
895	Note that simplification may be necessary.
896	sage: a = sqrt(2)+sqrt(3)+sqrt(5)
897	sage: f = a.minpoly(); f
898	x^8 - 40x^6 + 352x^4 - 960*x^2 + 576
899	sage: f(a)
900	(sqrt(5) + sqrt(3) + sqrt(2))^2((sqrt(5) + sqrt(3) + sqrt(2))^2((sqrt(5) + sqrt(3) + sqrt(2))^2*((sqrt(5) + sqrt(3) + sqrt(2))^2 - 40) + 352) - 960) + 576
901	sage: f(a).simplify_radical()
902	0
903
904	Here we verify it gives the same result as the abstract number field.
905	sage: (sqrt(2) + sqrt(3) + sqrt(6)).minpoly()
906	x^4 - 22x^2 - 48x - 23
907	sage: K.<a,b> = NumberField([x^2-2, x^2-3])
908	sage: (a+b+a*b).absolute_minpoly()
909	x^4 - 22x^2 - 48x - 23
910
911	Works with trig functions too.
912	sage: sin(pi/3).minpoly()
913	x^2 - 3/4
914
915	Here we show use of the epsilon parameter. That this result is
916	actually exact can be shown using the addition formula for sin,
917	but maxima is unable to see that.
918
919	sage: a = sin(pi/5)
920	sage: a.minpoly()
921	Traceback (most recent call last):
922	...
923	NotImplementedError: Could not prove minimal polynomial x^4 - 5/4*x^2 + 5/16 (epsilon 0.00000000000000e-1)
924	sage: f = a.minpoly(epsilon=1e-100); f
925	x^4 - 5/4*x^2 + 5/16
926	sage: f(a).numerical_approx(100)
927	0.00000000000000000000000000000
928
929	The degree must be high enough (default tops out at 24).
930	sage: a = sqrt(3) + sqrt(2)
931	sage: a.minpoly(bits=100, degree=3)
932	Traceback (most recent call last):
933	...
934	ValueError: Could not find minimal polynomial (100 bits, degree 3).
935	sage: a.minpoly(bits=100, degree=10)
936	x^4 - 10*x^2 + 1
937
938	Here we solve a cubic and then recover it from its complicated radical expansion.
939	sage: f = x^3 - x + 1
940	sage: a = f.solve(x)[0].rhs(); a
941	(sqrt(3)I/2 - 1/2)/(3(sqrt(23)/(6sqrt(3)) - 1/2)^(1/3)) + (sqrt(23)/(6sqrt(3)) - 1/2)^(1/3)(-sqrt(3)I/2 - 1/2)
942	sage: a.minpoly()
943	x^3 - x + 1
944	"""
945	bits_list = [bits] if bits else [100,200,500,1000]
946	degree_list = [degree] if degree else [2,4,8,12,24]
947
948	for bits in bits_list:
949	a = self.numerical_approx(bits)
950	check_bits = int(1.25 * bits + 80)
951	aa = self.numerical_approx(check_bits)
952
953	for degree in degree_list:
954
955	f = algdep(a, degree) # TODO: use the known_bits parameter?
956	# If indeed we have found a minimal polynomial,
957	# it should be accurate to a much higher precision.
958	error = abs(f(aa))
959	dx = ~RR(Integer(1) << (check_bits - degree - 2))
960	expected_error = abs(f.derivative()(CC(aa))) * dx
961
962	if error < expected_error:
963	# Degree might have been an over-estimate, factor because we want (irreducible) minpoly.
964	ff = f.factor()
965	for g, e in ff:
966	lead = g.leading_coefficient()
967	if lead != 1:
968	g = g / lead
969	expected_error = abs(g.derivative()(CC(aa))) * dx
970	error = abs(g(aa))
971	if error < expected_error:
972	# See if we can prove equality exactly
973	if g(self).simplify_trig().simplify_radical() == 0:
974	return g
975	# Otherwise fall back to numerical guess
976	elif epsilon and error < epsilon:
977	return g
978	else:
979	raise NotImplementedError, "Could not prove minimal polynomial %s (epsilon %s)" % (g, RR(error).str(no_sci=False))
980
981	raise ValueError, "Could not find minimal polynomial (%s bits, degree %s)." % (bits, degree)
982
983	def _mpfr_(self, field):
984	raise TypeError
985
986	def _complex_mpfr_field_(self, field):
987	raise TypeError
988
989	def _complex_double_(self, C):
990	raise TypeError
991
992	def _real_double_(self, R):
993	raise TypeError
994
995	def _real_rqdf_(self, R):
996	raise TypeError
997
998	def _rational_(self):
999	return Rational(repr(self))
1000
1001	def __abs__(self):
1002	return abs_symbolic(self)
1003
1004	def _integer_(self):
1005	"""
1006	EXAMPLES:
1007
1008	"""
1009	return Integer(repr(self))
1010
1011	def _add_(self, right):
1012	"""
1013	EXAMPLES:
1014	sage: var('x,y')
1015	(x, y)
1016	sage: x + y
1017	y + x
1018	sage: x._add_(y)
1019	y + x
1020	"""
1021	return SymbolicArithmetic([self, right], operator.add)
1022
1023	def _sub_(self, right):
1024	"""
1025	EXAMPLES:
1026	sage: var('x,y')
1027	(x, y)
1028	sage: x - y
1029	x - y
1030	"""
1031	return SymbolicArithmetic([self, right], operator.sub)
1032
1033	def _mul_(self, right):
1034	"""
1035	EXAMPLES:
1036	sage: var('x,y')
1037	(x, y)
1038	sage: x * y
1039	x*y
1040	"""
1041	return SymbolicArithmetic([self, right], operator.mul)
1042
1043	def _div_(self, right):
1044	"""
1045	EXAMPLES:
1046	sage: var('x,y')
1047	(x, y)
1048	sage: x / y
1049	x/y
1050	"""
1051	return SymbolicArithmetic([self, right], operator.div)
1052
1053	def __pow__(self, right):
1054	"""
1055	EXAMPLES:
1056	sage: var('x,n')
1057	(x, n)
1058	sage: x^(n+1)
1059	x^(n + 1)
1060	"""
1061	right = self.parent()(right)
1062	return SymbolicArithmetic([self, right], operator.pow)
1063
1064	def variables(self, vars=tuple([])):
1065	"""
1066	Return sorted list of variables that occur in the simplified
1067	form of self.
1068
1069	OUTPUT:
1070	a Python set
1071
1072	EXAMPLES:
1073	sage: var('x,n')
1074	(x, n)
1075	sage: f = x^(n+1) + sin(pi/19); f
1076	x^(n + 1) + sin(pi/19)
1077	sage: f.variables()
1078	(n, x)
1079
1080	sage: a = e^x
1081	sage: a.variables()
1082	(x,)
1083	"""
1084	return vars
1085
1086	def _first_variable(self):
1087	try:
1088	return self.__first_variable
1089	except AttributeError:
1090	pass
1091	v = self.variables()
1092	if len(v) == 0:
1093	ans = var('x')
1094	else:
1095	ans = v[0]
1096	self.__first_variable = ans
1097	return ans
1098
1099	def _has_op(self, operator):
1100	"""
1101	Recursively searches for the given operator in a SymbolicExpression
1102	object.
1103
1104	INPUT:
1105	operator: the operator to search for
1106
1107	OUTPUT:
1108	True or False
1109
1110	EXAMPLES:
1111	sage: f = 4*(x^2 - 3)
1112	sage: f._has_op(operator.sub)
1113	True
1114	sage: f._has_op(operator.div)
1115	False
1116	"""
1117
1118	# if we are the operator, then return true right away
1119	try:
1120	if operator is self._operator:
1121	return True
1122	except AttributeError:
1123	pass
1124
1125	# now try to look at this guy's operands
1126	try:
1127	ops = self._operands
1128	# if we don't have an operands, then we can return false
1129	except AttributeError:
1130	return False
1131	for oprnd in ops:
1132	if oprnd._has_op(operator): return True
1133	else: pass
1134	# if we get to this point, neither of the operands have the required
1135	# operator
1136	return False
1137
1138
1139	def __call__(self, dict=None, **kwds):
1140	return self.substitute(dict, **kwds)
1141
1142	def power_series(self, base_ring):
1143	"""
1144	Return algebraic power series associated to this symbolic
1145	expression, which must be a polynomial in one variable, with
1146	coefficients coercible to the base ring.
1147
1148	The power series is truncated one more than the degree.
1149
1150	EXAMPLES:
1151	sage: theta = var('theta')
1152	sage: f = theta^3 + (1/3)*theta - 17/3
1153	sage: g = f.power_series(QQ); g
1154	-17/3 + 1/3*theta + theta^3 + O(theta^4)
1155	sage: g^3
1156	-4913/27 + 289/9theta - 17/9theta^2 + 2602/27*theta^3 + O(theta^4)
1157	sage: g.parent()
1158	Power Series Ring in theta over Rational Field
1159	"""
1160	v = self.variables()
1161	if len(v) != 1:
1162	raise ValueError, "self must be a polynomial in one variable but it is in the variables %s"%tuple([v])
1163	f = self.polynomial(base_ring)
1164	from sage.rings.all import PowerSeriesRing
1165	R = PowerSeriesRing(base_ring, names=f.parent().variable_names())
1166	return R(f, f.degree()+1)
1167
1168	def polynomial(self, base_ring):
1169	r"""
1170	Return self as an algebraic polynomial over the given base ring, if
1171	possible.
1172
1173	The point of this function is that it converts purely symbolic
1174	polynomials into optimized algebraic polynomials over a given
1175	base ring.
1176
1177	WARNING: This is different from \code{self.poly(x)} which is used
1178	to rewrite self as a polynomial in x.
1179
1180	INPUT:
1181	base_ring -- a ring
1182
1183	EXAMPLES:
1184	sage: f = x^2 -2/3*x + 1
1185	sage: f.polynomial(QQ)
1186	x^2 - 2/3*x + 1
1187	sage: f.polynomial(GF(19))
1188	x^2 + 12*x + 1
1189
1190	Polynomials can be useful for getting the coefficients
1191	of an expression:
1192	sage: g = 6*x^2 - 5
1193	sage: g.coeffs()
1194	[[-5, 0], [6, 2]]
1195	sage: g.polynomial(QQ).list()
1196	[-5, 0, 6]
1197	sage: g.polynomial(QQ).dict()
1198	{0: -5, 2: 6}
1199
1200	sage: f = x^2*e + x + pi/e
1201	sage: f.polynomial(RDF)
1202	2.71828182846x^2 + 1.0x + 1.15572734979
1203	sage: g = f.polynomial(RR); g
1204	2.71828182845905x^2 + 1.00000000000000x + 1.15572734979092
1205	sage: g.parent()
1206	Univariate Polynomial Ring in x over Real Field with 53 bits of precision
1207	sage: f.polynomial(RealField(100))
1208	2.7182818284590452353602874714x^2 + 1.0000000000000000000000000000x + 1.1557273497909217179100931833
1209	sage: f.polynomial(CDF)
1210	2.71828182846x^2 + 1.0x + 1.15572734979
1211	sage: f.polynomial(CC)
1212	2.71828182845905x^2 + 1.00000000000000x + 1.15572734979092
1213
1214	We coerce a multivariate polynomial with complex symbolic coefficients:
1215	sage: x, y, n = var('x, y, n')
1216	sage: f = pi^3x - y^2e - I; f
1217	-1ey^2 + pi^3*x - I
1218	sage: f.polynomial(CDF)
1219	(-2.71828182846)y^2 + 31.0062766803x - 1.0*I
1220	sage: f.polynomial(CC)
1221	(-2.71828182845905)y^2 + 31.0062766802998x - 1.00000000000000*I
1222	sage: f.polynomial(ComplexField(70))
1223	(-2.7182818284590452354)y^2 + 31.006276680299820175x - 1.0000000000000000000*I
1224
1225	Another polynomial:
1226	sage: f = sum((eI)^nx^n for n in range(5)); f
1227	e^4x^4 - e^3Ix^3 - e^2x^2 + eIx + 1
1228	sage: f.polynomial(CDF)
1229	54.5981500331x^4 + (-20.0855369232I)x^3 + (-7.38905609893)x^2 + 2.71828182846Ix + 1.0
1230	sage: f.polynomial(CC)
1231	54.5981500331442x^4 + (-20.0855369231877I)x^3 + (-7.38905609893065)x^2 + 2.71828182845905Ix + 1.00000000000000
1232
1233	A multivariate polynomial over a finite field:
1234	sage: f = (3x^5 - 5y^5)^7; f
1235	(3x^5 - 5y^5)^7
1236	sage: g = f.polynomial(GF(7)); g
1237	3x^35 + 2y^35
1238	sage: parent(g)
1239	Multivariate Polynomial Ring in x, y over Finite Field of size 7
1240	"""
1241	vars = self.variables()
1242	if len(vars) == 0:
1243	vars = ['x']
1244	R = PolynomialRing(base_ring, names=vars)
1245	G = R.gens()
1246	V = R.variable_names()
1247	return self.substitute_over_ring(
1248	dict([(var(V[i]),G[i]) for i in range(len(G))]), ring=R)
1249
1250	def _polynomial_(self, R):
1251	"""
1252	Coerce this symbolic expression to a polynomial in R.
1253
1254	EXAMPLES:
1255	sage: var('x,y,z,w')
1256	(x, y, z, w)
1257
1258	sage: R = QQ[x,y,z]
1259	sage: R(x^2 + y)
1260	x^2 + y
1261	sage: R = QQ[w]
1262	sage: R(w^3 + w + 1)
1263	w^3 + w + 1
1264	sage: R = GF(7)[z]
1265	sage: R(z^3 + 10*z)
1266	z^3 + 3*z
1267
1268	NOTE: If the base ring of the polynomial ring is the symbolic
1269	ring, then a constant polynomial is always returned.
1270	sage: R = SR[x]
1271	sage: a = R(sqrt(2) + x^3 + y)
1272	sage: a
1273	y + x^3 + sqrt(2)
1274	sage: type(a)
1275	<type 'sage.rings.polynomial.polynomial_element.Polynomial_generic_dense'>
1276	sage: a.degree()
1277	0
1278
1279	We coerce to a double precision complex polynomial ring:
1280	sage: f = ex^3 + piy^3 + sqrt(2) + I; f
1281	piy^3 + ex^3 + I + sqrt(2)
1282	sage: R = CDF[x,y]
1283	sage: R(f)
1284	2.71828182846x^3 + 3.14159265359y^3 + 1.41421356237 + 1.0*I
1285
1286	We coerce to a higher-precision polynomial ring
1287	sage: R = ComplexField(100)[x,y]
1288	sage: R(f)
1289	2.7182818284590452353602874714x^3 + 3.1415926535897932384626433833y^3 + 1.4142135623730950488016887242 + 1.0000000000000000000000000000*I
1290	"""
1291	vars = self.variables()
1292	B = R.base_ring()
1293	if B == SR:
1294	if is_MPolynomialRing(R):
1295	return R({tuple([0]*R.ngens()):self})
1296	else:
1297	return R([self])
1298	G = R.gens()
1299	sub = []
1300	for v in vars:
1301	r = repr(v)
1302	for g in G:
1303	if repr(g) == r:
1304	sub.append((v,g))
1305	if len(sub) == 0:
1306	try:
1307	return R(B(self))
1308	except TypeError:
1309	if len(vars) == 1:
1310	sub = [(vars[0], G[0])]
1311	else:
1312	raise
1313	return self.substitute_over_ring(dict(sub), ring=R)
1314
1315	def function(self, *args):
1316	"""
1317	Return a CallableSymbolicExpression, fixing a variable order
1318	to be the order of args.
1319
1320	EXAMPLES:
1321	We will use several symbolic variables in the examples below:
1322	sage: var('x, y, z, t, a, w, n')
1323	(x, y, z, t, a, w, n)
1324
1325	sage: u = sin(x) + x*cos(y)
1326	sage: g = u.function(x,y)
1327	sage: g(x,y)
1328	x*cos(y) + sin(x)
1329	sage: g(t,z)
1330	t*cos(z) + sin(t)
1331	sage: g(x^2, x^y)
1332	x^2*cos(x^y) + sin(x^2)
1333
1334	sage: f = (x^2 + sin(a*w)).function(a,x,w); f
1335	(a, x, w) \|--> x^2 + sin(a*w)
1336	sage: f(1,2,3)
1337	sin(3) + 4
1338
1339	Using the function method we can obtain the above function f,
1340	but viewed as a function of different variables:
1341	sage: h = f.function(w,a); h
1342	(w, a) \|--> x^2 + sin(a*w)
1343
1344	This notation also works:
1345	sage: h(w,a) = f
1346	sage: h
1347	(w, a) \|--> x^2 + sin(a*w)
1348
1349	You can even make a symbolic expression $f$ into a function by
1350	writing \code{f(x,y) = f}:
1351	sage: f = x^n + y^n; f
1352	y^n + x^n
1353	sage: f(x,y) = f
1354	sage: f
1355	(x, y) \|--> y^n + x^n
1356	sage: f(2,3)
1357	3^n + 2^n
1358	"""
1359	R = CallableSymbolicExpressionRing(args)
1360	return R(self)
1361
1362
1363	###################################################################
1364	# derivative
1365	###################################################################
1366	def derivative(self, *args):
1367	"""
1368	Returns the derivative of itself. If self has exactly one variable, then
1369	it differentiates with respect to that variable. If there is more than one
1370	variable in the expression, then you must explicitly supply a variable.
1371	If you supply a variable $x$ followed by a number $n$, then it will
1372	differentiate with respect to $n$ times with respect to $n$.
1373
1374	You may supply more than one variable. Each variable may optionally be
1375	followed by a positive integer. Then SAGE will differentiate with
1376	respect to the first variable $n$ times, where $n$ is the number
1377	immediately following the variable in the parameter list. If the
1378	variable is not followed by an integer, then SAGE will differentiate
1379	once. Then SAGE will differentiate by the second variables, and if that
1380	is followed by a number $m$, it will differentiate $m$ times, and so on.
1381
1382	EXAMPLES:
1383	sage: h = sin(x)/cos(x)
1384	sage: diff(h,x,x,x)
1385	6sin(x)^4/cos(x)^4 + 8sin(x)^2/cos(x)^2 + 2
1386	sage: diff(h,x,3)
1387	6sin(x)^4/cos(x)^4 + 8sin(x)^2/cos(x)^2 + 2
1388
1389	sage: var('x, y')
1390	(x, y)
1391	sage: u = (sin(x) + cos(y))*(cos(x) - sin(y))
1392	sage: diff(u,x,y)
1393	sin(x)sin(y) - cos(x)cos(y)
1394	sage: f = ((x^2+1)/(x^2-1))^(1/4)
1395	sage: g = diff(f, x); g # this is a complex expression
1396	x/(2(x^2 - 1)^(1/4)(x^2 + 1)^(3/4)) - x(x^2 + 1)^(1/4)/(2(x^2 - 1)^(5/4))
1397	sage: g.simplify_rational()
1398	-x/((x^2 - 1)^(5/4)*(x^2 + 1)^(3/4))
1399
1400	sage: f = y^(sin(x))
1401	sage: diff(f, x)
1402	cos(x)y^sin(x)log(y)
1403
1404	sage: g(x) = sqrt(5-2*x)
1405	sage: g_3 = diff(g, x, 3); g_3(2)
1406	-3
1407
1408	sage: f = x*e^(-x)
1409	sage: diff(f, 100)
1410	xe^(-x) - 100e^(-x)
1411
1412	sage: g = 1/(sqrt((x^2-1)*(x+5)^6))
1413	sage: diff(g, x)
1414	-3/((x + 5)^3sqrt(x^2 - 1)abs(x + 5)) - x/((x^2 - 1)^(3/2)*abs(x + 5)^3)
1415	"""
1416	# check each time
1417	s = ""
1418	# see if we can implicitly supply a variable name
1419	try:
1420	a = args[0]
1421	except IndexError:
1422	# if there were NO arguments, try assuming
1423	a = 1
1424	if a is None or isinstance(a, (int, long, Integer)):
1425	vars = self.variables()
1426	if len(vars) == 1:
1427	s = "%s, %s" % (vars[0], a)
1428	else:
1429	raise ValueError, "must supply an explicit variable for an " +\
1430	"expression containing more than one variable"
1431	for i in range(len(args)):
1432	if isinstance(args[i], SymbolicVariable):
1433	s = s + '%s, ' %repr(args[i])
1434	# check to see if this is followed by an integer
1435	try:
1436	if isinstance(args[i+1], (int, long, Integer)):
1437	s = s + '%s, ' %repr(args[i+1])
1438	else:
1439	s = s + '1, '
1440	except IndexError:
1441	s = s + '1'
1442	elif isinstance(args[i], (int, long, Integer)):
1443	if args[i] == 0:
1444	return self
1445	if args[i] < 0:
1446	raise ValueError, "cannot take negative derivative"
1447	else:
1448	raise TypeError, "arguments must be integers or " +\
1449	"SymbolicVariable objects"
1450
1451	try:
1452	if s[-2] == ',':
1453	s = s[:-2]
1454	except IndexError:
1455	pass
1456	t = maxima('diff(%s, %s)'%(self._maxima_().name(), s))
1457	f = self.parent()(t)
1458	return f
1459
1460	differentiate = derivative
1461	diff = derivative
1462
1463
1464	###################################################################
1465	# Taylor series
1466	###################################################################
1467	def taylor(self, v, a, n):
1468	"""
1469	Expands self in a truncated Taylor or Laurent series in the
1470	variable v around the point a, containing terms through $(x - a)^n$.
1471
1472	INPUT:
1473	v -- variable
1474	a -- number
1475	n -- integer
1476
1477	EXAMPLES:
1478	sage: var('a, x, z')
1479	(a, x, z)
1480	sage: taylor(a*log(z), z, 2, 3)
1481	log(2)a + a(z - 2)/2 - a(z - 2)^2/8 + a(z - 2)^3/24
1482	sage: taylor(sqrt (sin(x) + a*x + 1), x, 0, 3)
1483	1 + (a + 1)x/2 - (a^2 + 2a + 1)x^2/8 + (3a^3 + 9a^2 + 9a - 1)*x^3/48
1484	sage: taylor (sqrt (x + 1), x, 0, 5)
1485	1 + x/2 - x^2/8 + x^3/16 - 5x^4/128 + 7x^5/256
1486	sage: taylor (1/log (x + 1), x, 0, 3)
1487	1/x + 1/2 - x/12 + x^2/24 - 19*x^3/720
1488	sage: taylor (cos(x) - sec(x), x, 0, 5)
1489	-x^2 - x^4/6
1490	sage: taylor ((cos(x) - sec(x))^3, x, 0, 9)
1491	-x^6 - x^8/2
1492	sage: taylor (1/(cos(x) - sec(x))^3, x, 0, 5)
1493	-1/x^6 + 1/(2x^4) + 11/(120x^2) - 347/15120 - 6767x^2/604800 - 15377x^4/7983360
1494	"""
1495	v = var(v)
1496	l = self._maxima_().taylor(v, SR(a), Integer(n))
1497	return self.parent()(l)
1498
1499	###################################################################
1500	# limits
1501	###################################################################
1502	def limit(self, dir=None, taylor=False, **argv):
1503	r"""
1504	Return the limit as the variable v approaches a from the
1505	given direction.
1506
1507	\begin{verbatim}
1508	expr.limit(x = a)
1509	expr.limit(x = a, dir='above')
1510	\end{verbatim}
1511
1512	INPUT:
1513	dir -- (default: None); dir may have the value `plus' (or 'above')
1514	for a limit from above, `minus' (or 'below') for a limit from
1515	below, or may be omitted (implying a two-sided
1516	limit is to be computed).
1517	taylor -- (default: False); if True, use Taylor series, which
1518	allows more integrals to be computed (but may also crash
1519	in some obscure cases due to bugs in Maxima).
1520	**argv -- 1 named parameter
1521
1522	NOTE: Output it may also use `und' (undefined), `ind'
1523	(indefinite but bounded), and `infinity' (complex infinity).
1524
1525	EXAMPLES:
1526	sage: f = (1+1/x)^x
1527	sage: f.limit(x = oo)
1528	e
1529	sage: f.limit(x = 5)
1530	7776/3125
1531	sage: f.limit(x = 1.2)
1532	2.0696157546720...
1533	sage: f.limit(x = I, taylor=True)
1534	(1 - I)^I
1535	sage: f(1.2)
1536	2.0696157546720...
1537	sage: f(I)
1538	(1 - I)^I
1539	sage: CDF(f(I))
1540	2.06287223508 + 0.74500706218*I
1541	sage: CDF(f.limit(x = I))
1542	2.06287223508 + 0.74500706218*I
1543
1544	More examples:
1545	sage: limit(x*log(x), x = 0, dir='above')
1546	0
1547	sage: lim((x+1)^(1/x),x = 0)
1548	e
1549	sage: lim(e^x/x, x = oo)
1550	+Infinity
1551	sage: lim(e^x/x, x = -oo)
1552	0
1553	sage: lim(-e^x/x, x = oo)
1554	-Infinity
1555	sage: lim((cos(x))/(x^2), x = 0)
1556	+Infinity
1557	sage: lim(sqrt(x^2+1) - x, x = oo)
1558	0
1559	sage: lim(x^2/(sec(x)-1), x=0)
1560	2
1561	sage: lim(cos(x)/(cos(x)-1), x=0)
1562	-Infinity
1563	sage: lim(x*sin(1/x), x=0)
1564	0
1565
1566	Traceback (most recent call last):
1567	...
1568	TypeError: Computation failed since Maxima requested additional constraints (use assume):
1569	Is x positive or negative?
1570
1571	sage: f = log(log(x))/log(x)
1572	sage: forget(); assume(x<-2); lim(f, x=0, taylor=True)
1573	und
1574
1575	Here ind means "indefinite but bounded":
1576	sage: lim(sin(1/x), x = 0)
1577	ind
1578	"""
1579	if len(argv) != 1:
1580	raise ValueError, "call the limit function like this, e.g. limit(expr, x=2)."
1581	else:
1582	k = argv.keys()[0]
1583	v = var(k, create=False)
1584	a = argv[k]
1585	if dir is None:
1586	if taylor:
1587	l = self._maxima_().tlimit(v, a)
1588	else:
1589	l = self._maxima_().limit(v, a)
1590	elif dir == 'plus' or dir == 'above':
1591	if taylor:
1592	l = self._maxima_().tlimit(v, a, 'plus')
1593	else:
1594	l = self._maxima_().limit(v, a, 'plus')
1595	elif dir == 'minus' or dir == 'below':
1596	if taylor:
1597	l = self._maxima_().tlimit(v, a, 'minus')
1598	else:
1599	l = self._maxima_().limit(v, a, 'minus')
1600	else:
1601	raise ValueError, "dir must be one of 'plus' or 'minus'"
1602	return self.parent()(l)
1603
1604	###################################################################
1605	# Laplace transform
1606	###################################################################
1607	def laplace(self, t, s):
1608	r"""
1609	Attempts to compute and return the Laplace transform of self
1610	with respect to the variable t and transform parameter s. If
1611	Laplace cannot find a solution, a formal function is returned.
1612
1613	The function that is returned maybe be viewed as a function of s.
1614
1615	DEFINITION:
1616	The Laplace transform of a function $f(t)$, defined for all
1617	real numbers $t \geq 0$, is the function $F(s)$ defined by
1618	$$
1619	F(s) = \int_{0}^{\infty} e^{-st} f(t) dt.
1620	$$
1621
1622	EXAMPLES:
1623	We compute a few Laplace transforms:
1624	sage: var('x, s, z, t, t0')
1625	(x, s, z, t, t0)
1626	sage: sin(x).laplace(x, s)
1627	1/(s^2 + 1)
1628	sage: (z + exp(x)).laplace(x, s)
1629	z/s + 1/(s - 1)
1630	sage: log(t/t0).laplace(t, s)
1631	(-log(t0) - log(s) - euler_gamma)/s
1632
1633	We do a formal calculation:
1634	sage: f = function('f', x)
1635	sage: g = f.diff(x); g
1636	diff(f(x), x, 1)
1637	sage: g.laplace(x, s)
1638	s*laplace(f(x), x, s) - f(0)
1639
1640	EXAMPLE: A BATTLE BETWEEN the X-women and the Y-men (by David Joyner):
1641	Solve
1642	$$
1643	x' = -16y, x(0)=270, y' = -x + 1, y(0) = 90.
1644	$$
1645	This models a fight between two sides, the "X-women"
1646	and the "Y-men", where the X-women have 270 initially and
1647	the Y-men have 90, but the Y-men are better at fighting,
1648	because of the higher factor of "-16" vs "-1", and also get
1649	an occasional reinforcement, because of the "+1" term.
1650
1651	sage: var('t')
1652	t
1653	sage: t = var('t')
1654	sage: x = function('x', t)
1655	sage: y = function('y', t)
1656	sage: de1 = x.diff(t) + 16*y
1657	sage: de2 = y.diff(t) + x - 1
1658	sage: de1.laplace(t, s)
1659	16laplace(y(t), t, s) + slaplace(x(t), t, s) - x(0)
1660	sage: de2.laplace(t, s)
1661	s*laplace(y(t), t, s) + laplace(x(t), t, s) - 1/s - y(0)
1662
1663	Next we form the augmented matrix of the above system:
1664	sage: A = matrix([[s, 16, 270],[1, s, 90+1/s]])
1665	sage: E = A.echelon_form()
1666	sage: xt = E[0,2].inverse_laplace(s,t)
1667	sage: yt = E[1,2].inverse_laplace(s,t)
1668	sage: print xt
1669	4 t - 4 t
1670	91 e 629 e
1671	- -------- + ----------- + 1
1672	2 2
1673	sage: print yt
1674	4 t - 4 t
1675	91 e 629 e
1676	-------- + -----------
1677	8 8
1678	sage: p1 = plot(xt,0,1/2,rgbcolor=(1,0,0))
1679	sage: p2 = plot(yt,0,1/2,rgbcolor=(0,1,0))
1680	sage: (p1+p2).save()
1681	"""
1682	return self.parent()(self._maxima_().laplace(var(t), var(s)))
1683
1684	def inverse_laplace(self, t, s):
1685	r"""
1686	Attempts to compute the inverse Laplace transform of self with
1687	respect to the variable t and transform parameter s. If
1688	Laplace cannot find a solution, a formal function is returned.
1689
1690	The function that is returned maybe be viewed as a function of s.
1691
1692
1693	DEFINITION:
1694	The inverse Laplace transform of a function $F(s)$,
1695	is the function $f(t)$ defined by
1696	$$
1697	F(s) = \frac{1}{2\pi i} \int_{\gamma-i\infty}^{\gamma + i\infty} e^{st} F(s) dt,
1698	$$
1699	where $\gamma$ is chosen so that the contour path of
1700	integration is in the region of convergence of $F(s)$.
1701
1702	EXAMPLES:
1703	sage: var('w, m')
1704	(w, m)
1705	sage: f = (1/(w^2+10)).inverse_laplace(w, m); f
1706	sin(sqrt(10)*m)/sqrt(10)
1707	sage: laplace(f, m, w)
1708	1/(w^2 + 10)
1709	"""
1710	return self.parent()(self._maxima_().ilt(var(t), var(s)))
1711
1712
1713	###################################################################
1714	# a default variable, for convenience.
1715	###################################################################
1716	def default_variable(self):
1717	vars = self.variables()
1718	if len(vars) < 1:
1719	return var('x')
1720	else:
1721	return vars[0]
1722
1723	###################################################################
1724	# integration
1725	###################################################################
1726	def integral(self, v=None, a=None, b=None):
1727	"""
1728	Returns the indefinite integral with respect to the variable
1729	$v$, ignoring the constant of integration. Or, if endpoints
1730	$a$ and $b$ are specified, returns the definite integral over
1731	the interval $[a, b]$.
1732
1733	If \code{self} has only one variable, then it returns the
1734	integral with respect to that variable.
1735
1736	INPUT:
1737	v -- (optional) a variable or variable name
1738	a -- (optional) lower endpoint of definite integral
1739	b -- (optional) upper endpoint of definite integral
1740
1741
1742	EXAMPLES:
1743	sage: h = sin(x)/(cos(x))^2
1744	sage: h.integral(x)
1745	1/cos(x)
1746
1747	sage: f = x^2/(x+1)^3
1748	sage: f.integral()
1749	log(x + 1) + (4x + 3)/(2x^2 + 4*x + 2)
1750
1751	sage: f = x*cos(x^2)
1752	sage: f.integral(x, 0, sqrt(pi))
1753	0
1754	sage: f.integral(a=-pi, b=pi)
1755	0
1756
1757	sage: f(x) = sin(x)
1758	sage: f.integral(x, 0, pi/2)
1759	1
1760
1761	Constraints are sometimes needed:
1762	sage: var('x, n')
1763	(x, n)
1764	sage: integral(x^n,x)
1765	Traceback (most recent call last):
1766	...
1767	TypeError: Computation failed since Maxima requested additional constraints (use assume):
1768	Is n+1 zero or nonzero?
1769	sage: assume(n > 0)
1770	sage: integral(x^n,x)
1771	x^(n + 1)/(n + 1)
1772	sage: forget()
1773
1774
1775	Note that an exception is raised when a definite integral is divergent.
1776	sage: integrate(1/x^3,x,0,1)
1777	Traceback (most recent call last):
1778	...
1779	ValueError: Integral is divergent.
1780	sage: integrate(1/x^3,x,-1,3)
1781	Traceback (most recent call last):
1782	...
1783	ValueError: Integral is divergent.
1784
1785	NOTE: Above, putting assume(n == -1) does not yield the right behavior.
1786	Directly in maxima, doing
1787
1788	The examples in the Maxima documentation:
1789	sage: var('x, y, z, b')
1790	(x, y, z, b)
1791	sage: integral(sin(x)^3)
1792	cos(x)^3/3 - cos(x)
1793	sage: integral(x/sqrt(b^2-x^2))
1794	xlog(2sqrt(b^2 - x^2) + 2*b)
1795	sage: integral(x/sqrt(b^2-x^2), x)
1796	-sqrt(b^2 - x^2)
1797	sage: integral(cos(x)^2 * exp(x), x, 0, pi)
1798	3*e^pi/5 - 3/5
1799	sage: integral(x^2 * exp(-x^2), x, -oo, oo)
1800	sqrt(pi)/2
1801
1802	We integrate the same function in both Mathematica and SAGE (via Maxima):
1803	sage: f = sin(x^2) + y^z
1804	sage: g = mathematica(f) # optional -- requires mathematica
1805	sage: print g # optional
1806	z 2
1807	y + Sin[x ]
1808	sage: print g.Integrate(x) # optional
1809	z Pi 2
1810	x y + Sqrt[--] FresnelS[Sqrt[--] x]
1811	2 Pi
1812	sage: print f.integral(x)
1813	z (sqrt(2) I + sqrt(2)) x
1814	x y + sqrt( pi) ((sqrt(2) I + sqrt(2)) erf(------------------------)
1815	2
1816	(sqrt(2) I - sqrt(2)) x
1817	+ (sqrt(2) I - sqrt(2)) erf(------------------------))/8
1818	2
1819
1820	We integrate the above function in maple now:
1821	sage: g = maple(f); g # optional -- requires maple
1822	sin(x^2)+y^z
1823	sage: g.integrate(x) # optional -- requires maple
1824	1/22^(1/2)Pi^(1/2)FresnelS(2^(1/2)/Pi^(1/2)x)+y^z*x
1825
1826	We next integrate a function with no closed form integral. Notice that
1827	the answer comes back as an expression that contains an integral itself.
1828	sage: A = integral(1/ ((x-4) * (x^3+2*x+1)), x); A
1829	log(x - 4)/73 - integrate((x^2 + 4x + 18)/(x^3 + 2x + 1), x)/73
1830	sage: print A
1831	/ 2
1832	[ x + 4 x + 18
1833	I ------------- dx
1834	] 3
1835	log(x - 4) / x + 2 x + 1
1836	---------- - ------------------
1837	73 73
1838
1839	ALIASES:
1840	integral() and integrate() are the same.
1841
1842	"""
1843
1844	if v is None:
1845	v = self.default_variable()
1846
1847	if not isinstance(v, SymbolicVariable):
1848	v = var(repr(v))
1849	#raise TypeError, 'must integrate with respect to a variable'
1850	if (a is None and (not b is None)) or (b is None and (not a is None)):
1851	raise TypeError, 'only one endpoint given'
1852	if a is None:
1853	return self.parent()(self._maxima_().integrate(v))
1854	else:
1855	try:
1856	return self.parent()(self._maxima_().integrate(v, a, b))
1857	except TypeError, error:
1858	s = str(error)
1859	if "divergent" in s or 'Principal Value' in s:
1860	raise ValueError, "Integral is divergent."
1861	else:
1862	raise TypeError, error
1863
1864
1865	integrate = integral
1866
1867	def nintegral(self, x, a, b,
1868	desired_relative_error='1e-8',
1869	maximum_num_subintervals=200):
1870	r"""
1871	Return a numerical approximation to the integral of self from
1872	a to b.
1873
1874	INPUT:
1875	x -- variable to integrate with respect to
1876	a -- lower endpoint of integration
1877	b -- upper endpoint of integration
1878	desired_relative_error -- (default: '1e-8') the desired
1879	relative error
1880	maximum_num_subintervals -- (default: 200) maxima number
1881	of subintervals
1882
1883	OUTPUT:
1884	-- float: approximation to the integral
1885	-- float: estimated absolute error of the approximation
1886	-- the number of integrand evaluations
1887	-- an error code:
1888	0 -- no problems were encountered
1889	1 -- too many subintervals were done
1890	2 -- excessive roundoff error
1891	3 -- extremely bad integrand behavior
1892	4 -- failed to converge
1893	5 -- integral is probably divergent or slowly convergent
1894	6 -- the input is invalid
1895
1896	ALIAS:
1897	nintegrate is the same as nintegral
1898
1899	REMARK:
1900	There is also a function \code{numerical_integral} that implements
1901	numerical integration using the GSL C library. It is potentially
1902	much faster and applies to arbitrary user defined functions.
1903
1904	EXAMPLES:
1905	sage: f(x) = exp(-sqrt(x))
1906	sage: f.nintegral(x, 0, 1)
1907	(0.52848223531423055, 4.1633141378838452e-11, 231, 0)
1908
1909	We can also use the \code{numerical_integral} function, which calls
1910	the GSL C library.
1911	sage: numerical_integral(f, 0, 1) # random low-order bits
1912	(0.52848223225314706, 6.8392846084921134e-07)
1913	"""
1914	v = self._maxima_().quad_qags(var(x),
1915	a, b, desired_relative_error,
1916	maximum_num_subintervals)
1917	return float(v[0]), float(v[1]), Integer(v[2]), Integer(v[3])
1918
1919	nintegrate = nintegral
1920
1921
1922	###################################################################
1923	# Manipulating epxressions
1924	###################################################################
1925	def coeff(self, x, n=1):
1926	"""
1927	Returns the coefficient of $x^n$ in self.
1928
1929	INPUT:
1930	x -- variable, function, expression, etc.
1931	n -- integer, default 1.
1932
1933	Sometimes it may be necessary to expand or factor first, since
1934	this is not done automatically.
1935
1936	EXAMPLES:
1937	sage: var('a, x, y, z')
1938	(a, x, y, z)
1939	sage: f = (asqrt(2))x^2 + sin(y)*x^(1/2) + z^z
1940	sage: f.coeff(sin(y))
1941	sqrt(x)
1942	sage: f.coeff(x^2)
1943	sqrt(2)*a
1944	sage: f.coeff(x^(1/2))
1945	sin(y)
1946	sage: f.coeff(1)
1947	0
1948	sage: f.coeff(x, 0)
1949	z^z
1950	"""
1951	return self.parent()(self._maxima_().coeff(x, n))
1952
1953	def coeffs(self, x=None):
1954	"""
1955	Coefficients of self as a polynomial in x.
1956
1957	INPUT:
1958	x -- optional variable
1959	OUTPUT:
1960	list of pairs [expr, n], where expr is a symbolic expression and n is a power.
1961
1962	EXAMPLES:
1963	sage: var('x, y, a')
1964	(x, y, a)
1965	sage: p = x^3 - (x-3)*(x^2+x) + 1
1966	sage: p.coeffs()
1967	[[1, 0], [3, 1], [2, 2]]
1968	sage: p = expand((x-a*sqrt(2))^2 + x + 1); p
1969	x^2 - 2sqrt(2)ax + x + 2a^2 + 1
1970	sage: p.coeffs(a)
1971	[[x^2 + x + 1, 0], [-2sqrt(2)x, 1], [2, 2]]
1972	sage: p.coeffs(x)
1973	[[2a^2 + 1, 0], [1 - 2sqrt(2)*a, 1], [1, 2]]
1974
1975	A polynomial with wacky exponents:
1976	sage: p = (17/3a)x^(3/2) + x*y + 1/x + x^x
1977	sage: p.coeffs(x)
1978	[[1, -1], [x^x, 0], [y, 1], [17*a/3, 3/2]]
1979	"""
1980	f = self._maxima_()
1981	P = f.parent()
1982	P._eval_line('load(coeflist)')
1983	if x is None:
1984	x = self.default_variable()
1985	x = var(x)
1986	G = f.coeffs(x)
1987	S = symbolic_expression_from_maxima_string(repr(G))
1988	return S[1:]
1989
1990	def poly(self, x=None):
1991	r"""
1992	Express self as a polynomial in x. Is self is not a polynomial
1993	in x, then some coefficients may be functions of x.
1994
1995	WARNING: This is different from \code{self.polynomial()} which
1996	returns a SAGE polynomial over a given base ring.
1997
1998	EXAMPLES:
1999	sage: var('a, x')
2000	(a, x)
2001	sage: p = expand((x-a*sqrt(2))^2 + x + 1); p
2002	x^2 - 2sqrt(2)ax + x + 2a^2 + 1
2003	sage: p.poly(a)
2004	(x^2 + x + 1)1 + -2sqrt(2)xa + 2*a^2
2005	sage: bool(expand(p.poly(a)) == p)
2006	True
2007	sage: p.poly(x)
2008	(2a^2 + 1)1 + (1 - 2sqrt(2)a)x + 1x^2
2009	"""
2010	f = self._maxima_()
2011	P = f.parent()
2012	P._eval_line('load(coeflist)')
2013	if x is None:
2014	x = self.default_variable()
2015	x = var(x)
2016	G = f.coeffs(x)
2017	ans = None
2018	for i in range(1, len(G)):
2019	Z = G[i]
2020	coeff = SR(Z[0])
2021	n = SR(Z[1])
2022	if repr(coeff) != '0':
2023	if repr(n) == '0':
2024	xpow = SR(1)
2025	elif repr(n) == '1':
2026	xpow = x
2027	else:
2028	xpow = x**n
2029	if ans is None:
2030	ans = coeff*xpow
2031	else:
2032	ans += coeff*xpow
2033	ans._declare_simplified()
2034	return ans
2035
2036	def combine(self):
2037	"""
2038	Simplifies self by combining all terms with the same
2039	denominator into a single term.
2040
2041	EXAMPLES:
2042	sage: var('x, y, a, b, c')
2043	(x, y, a, b, c)
2044	sage: f = x*(x-1)/(x^2 - 7) + y^2/(x^2-7) + 1/(x+1) + b/a + c/a
2045	sage: print f
2046	2
2047	y (x - 1) x 1 c b
2048	------ + --------- + ----- + - + -
2049	2 2 x + 1 a a
2050	x - 7 x - 7
2051	sage: print f.combine()
2052	2
2053	y + (x - 1) x 1 c + b
2054	-------------- + ----- + -----
2055	2 x + 1 a
2056	x - 7
2057	"""
2058	return self.parent()(self._maxima_().combine())
2059
2060	def numerator(self):
2061	"""
2062	EXAMPLES:
2063	sage: var('a,x,y')
2064	(a, x, y)
2065	sage: f = x(x-a)/((x^2 - y)(x-a))
2066	sage: print f
2067	x
2068	------
2069	2
2070	x - y
2071	sage: f.numerator()
2072	x
2073	sage: f.denominator()
2074	x^2 - y
2075	"""
2076	return self.parent()(self._maxima_().num())
2077
2078	def denominator(self):
2079	"""
2080	EXAMPLES:
2081	sage: var('x, y, z, theta')
2082	(x, y, z, theta)
2083	sage: f = (sqrt(x) + sqrt(y) + sqrt(z))/(x^10 - y^10 - sqrt(theta))
2084	sage: print f
2085	sqrt(z) + sqrt(y) + sqrt(x)
2086	---------------------------
2087	10 10
2088	- y + x - sqrt(theta)
2089	sage: f.denominator()
2090	-y^10 + x^10 - sqrt(theta)
2091	"""
2092	return self.parent()(self._maxima_().denom())
2093
2094	def factor_list(self, dontfactor=[]):
2095	"""
2096	Returns a list of the factors of self, as computed by the
2097	factor command.
2098
2099	INPUT:
2100	self -- a symbolic expression
2101	dontfactor -- see docs for self.factor.
2102
2103	REMARK: If you already have a factored expression and just
2104	want to get at the individual factors, use self._factor_list()
2105	instead.
2106
2107	EXAMPLES:
2108	sage: var('x, y, z')
2109	(x, y, z)
2110	sage: f = x^3-y^3
2111	sage: f.factor()
2112	-(y - x)(y^2 + xy + x^2)
2113
2114	Notice that the -1 factor is separated out:
2115	sage: f.factor_list()
2116	[(-1, 1), (y - x, 1), (y^2 + x*y + x^2, 1)]
2117
2118	We factor a fairly straightforward expression:
2119	sage: factor(-8y - 4x + z^2(2y + x)).factor_list()
2120	[(2*y + x, 1), (z - 2, 1), (z + 2, 1)]
2121
2122	This function also works for quotients:
2123	sage: f = -1 - 2x - x^2 + y^2 + 2xy^2 + x^2y^2
2124	sage: g = f/(36(1 + 2y + y^2)); g
2125	(x^2y^2 + 2xy^2 + y^2 - x^2 - 2x - 1)/(36(y^2 + 2y + 1))
2126	sage: g.factor(dontfactor=[x])
2127	(x^2 + 2x + 1)(y - 1)/(36*(y + 1))
2128	sage: g.factor_list(dontfactor=[x])
2129	[(x^2 + 2*x + 1, 1), (y - 1, 1), (36, -1), (y + 1, -1)]
2130
2131	An example, where one of the exponents is not an integer.
2132	sage: var('x, u, v')
2133	(x, u, v)
2134	sage: f = expand((2uv^2-v^2-4u^3)^2 (-u)^3 * (x-sin(x))^3)
2135	sage: f.factor()
2136	u^3(2uv^2 - v^2 - 4u^3)^2*(sin(x) - x)^3
2137	sage: g = f.factor_list(); g
2138	[(u, 3), (2uv^2 - v^2 - 4*u^3, 2), (sin(x) - x, 3)]
2139
2140	This example also illustrates that the exponents do not have
2141	to be integers.
2142	sage: f = x^(2sin(x)) (x-1)^(sqrt(2)*x); f
2143	(x - 1)^(sqrt(2)x)x^(2*sin(x))
2144	sage: f.factor_list()
2145	[(x - 1, sqrt(2)x), (x, 2sin(x))]
2146	"""
2147	return self.factor(dontfactor=dontfactor)._factor_list()
2148
2149	def _factor_list(self):
2150	if isinstance(self, SymbolicArithmetic):
2151	if self._operator == operator.mul:
2152	left, right = self._operands
2153	return left._factor_list() + right._factor_list()
2154	elif self._operator == operator.pow:
2155	left, right = self._operands
2156	return [(left, right)]
2157	elif self._operator == operator.div:
2158	left, right = self._operands
2159	return left._factor_list() + \
2160	[(x,-y) for x, y in right._factor_list()]
2161	elif self._operator == operator.neg:
2162	expr = self._operands[0]
2163	v = expr._factor_list()
2164	return [(SR(-1),SR(1))] + v
2165	return [(self, 1)]
2166
2167
2168	###################################################################
2169	# solve
2170	###################################################################
2171	def roots(self, x=None, explicit_solutions=True):
2172	r"""
2173	Returns the roots of self, with multiplicities.
2174
2175	INPUT:
2176	x -- variable to view the function in terms of
2177	(use default variable if not given)
2178	explicit_solutions -- bool (default True); require that
2179	roots be explicit rather than implicit
2180	OUTPUT:
2181	list of pairs (root, multiplicity)
2182
2183	If there are infinitely many roots, e.g., a function
2184	like sin(x), only one is returned.
2185
2186	EXAMPLES:
2187	sage: var('x, a')
2188	(x, a)
2189
2190	A simple example:
2191	sage: ((x^2-1)^2).roots()
2192	[(-1, 2), (1, 2)]
2193
2194	A complicated example.
2195	sage: f = expand((x^2 - 1)^3(x^2 + 1)(x-a)); f
2196	x^9 - ax^8 - 2x^7 + 2ax^6 + 2x^3 - 2a*x^2 - x + a
2197
2198	The default variable is a, since it is the first in alphabetical order:
2199	sage: f.roots()
2200	[(x, 1)]
2201
2202	As a polynomial in a, x is indeed a root:
2203	sage: f.poly(a)
2204	(x^9 - 2x^7 + 2x^3 - x)1 + (-x^8 + 2x^6 - 2x^2 + 1)a
2205	sage: f(a=x)
2206	0
2207
2208	The roots in terms of x are what we expect:
2209	sage: f.roots(x)
2210	[(a, 1), (-1*I, 1), (I, 1), (1, 3), (-1, 3)]
2211
2212	Only one root of $\sin(x) = 0$ is given:
2213	sage: f = sin(x)
2214	sage: f.roots(x)
2215	[(0, 1)]
2216
2217	We derive the roots of a general quadratic polynomial:
2218	sage: var('a,b,c,x')
2219	(a, b, c, x)
2220	sage: (ax^2 + bx + c).roots(x)
2221	[((-sqrt(b^2 - 4ac) - b)/(2a), 1), ((sqrt(b^2 - 4ac) - b)/(2a), 1)]
2222
2223
2224	By default, all the roots are required to be explicit rather than
2225	implicit. To get implicit roots, pass explicit_solutions=False
2226	to .roots()
2227	sage: var('x')
2228	x
2229	sage: f = x^(1/9) + (2^(8/9) - 2^(1/9))*(x - 1) - x^(8/9)
2230	sage: f.roots()
2231	Traceback (most recent call last):
2232	...
2233	RuntimeError: no explicit roots found
2234	sage: f.roots(explicit_solutions=False)
2235	[((x^(8/9) - x^(1/9) + 2^(8/9) - 2^(1/9))/(2^(8/9) - 2^(1/9)), 1)]
2236
2237	Another example, but involving a degree 5 poly whose roots
2238	don't get computed explicitly:
2239	sage: f = x^5 + x^3 + 17*x + 1
2240	sage: f.roots()
2241	Traceback (most recent call last):
2242	...
2243	RuntimeError: no explicit roots found
2244	sage: f.roots(explicit_solutions=False)
2245	[(x^5 + x^3 + 17*x + 1, 1)]
2246	"""
2247	if x is None:
2248	x = self.default_variable()
2249	S, mul = self.solve(x, multiplicities=True, explicit_solutions=explicit_solutions)
2250	if len(mul) == 0 and explicit_solutions:
2251	raise RuntimeError, "no explicit roots found"
2252	else:
2253	return [(S[i].rhs(), mul[i]) for i in range(len(mul))]
2254
2255	def solve(self, x, multiplicities=False, explicit_solutions=False):
2256	r"""
2257	Solve the equation \code{self == 0} for the variable x.
2258
2259	INPUT:
2260	x -- variable to solve for
2261	multiplicities -- bool (default: False); if True, return corresponding multiplicities.
2262	explicit_solutions -- bool (default:False); if True, require that all solutions
2263	returned be explicit (rather than implicit)
2264
2265	EXAMPLES:
2266	sage: z = var('z')
2267	sage: (z^5 - 1).solve(z)
2268	[z == e^(2Ipi/5), z == e^(4Ipi/5), z == e^(-(4Ipi/5)), z == e^(-(2Ipi/5)), z == 1]
2269	"""
2270	x = var(x)
2271	return (self == 0).solve(x, multiplicities=multiplicities, explicit_solutions=explicit_solutions)
2272
2273	###################################################################
2274	# simplify
2275	###################################################################
2276	def simplify(self):
2277	try:
2278	return self._simp
2279	except AttributeError:
2280	S = evaled_symbolic_expression_from_maxima_string(self._maxima_init_())
2281	S._simp = S
2282	self._simp = S
2283	return S
2284
2285	def simplify_trig(self):
2286	r"""
2287	First expands using trig_expand, then employs the identities
2288	$\sin(x)^2 + \cos(x)^2 = 1$ and $\cosh(x)^2 - \sin(x)^2 = 1$
2289	to simplify expressions containing tan, sec, etc., to sin,
2290	cos, sinh, cosh.
2291
2292	ALIAS: trig_simplify and simplify_trig are the same
2293
2294	EXAMPLES:
2295	sage: f = sin(x)^2 + cos(x)^2; f
2296	sin(x)^2 + cos(x)^2
2297	sage: f.simplify()
2298	sin(x)^2 + cos(x)^2
2299	sage: f.simplify_trig()
2300	1
2301	"""
2302	# much better to expand first, since it often doesn't work
2303	# right otherwise!
2304	return self.parent()(self._maxima_().trigexpand().trigsimp())
2305
2306	trig_simplify = simplify_trig
2307
2308	def simplify_rational(self):
2309	"""
2310	Simplify by expanding repeatedly rational expressions.
2311
2312	ALIAS: rational_simplify and simplify_rational are the same
2313
2314	EXAMPLES:
2315	sage: f = sin(x/(x^2 + x))
2316	sage: f
2317	sin(x/(x^2 + x))
2318	sage: f.simplify_rational()
2319	sin(1/(x + 1))
2320
2321	sage: f = ((x - 1)^(3/2) - (x + 1)sqrt(x - 1))/sqrt((x - 1)(x + 1))
2322	sage: print f
2323	3/2
2324	(x - 1) - sqrt(x - 1) (x + 1)
2325	--------------------------------
2326	sqrt((x - 1) (x + 1))
2327	sage: print f.simplify_rational()
2328	2 sqrt(x - 1)
2329	- -------------
2330	2
2331	sqrt(x - 1)
2332	"""
2333	return self.parent()(self._maxima_().fullratsimp())
2334
2335	rational_simplify = simplify_rational
2336
2337	# TODO: come up with a way to intelligently wrap Maxima's way of
2338	# fine-tuning all simplificationsrational
2339
2340	def simplify_radical(self):
2341	r"""
2342	Wraps the Maxima radcan() command. From the Maxima documentation:
2343
2344	Simplifies this expression, which can contain logs, exponentials,
2345	and radicals, by converting it into a form which is canonical over a
2346	large class of expressions and a given ordering of variables; that
2347	is, all functionally equivalent forms are mapped into a unique form.
2348	For a somewhat larger class of expressions, produces a regular form.
2349	Two equivalent expressions in this class do not necessarily have the
2350	same appearance, but their difference can be simplified by radcan to
2351	zero.
2352
2353	For some expressions radcan is quite time consuming. This is the
2354	cost of exploring certain relationships among the components of the
2355	expression for simplifications based on factoring and partial
2356	fraction expansions of exponents.
2357
2358	ALIAS: radical_simplify, simplify_log, log_simplify, exp_simplify,
2359	simplify_exp are all the same
2360
2361	EXAMPLES:
2362
2363	sage: var('x,y,a')
2364	(x, y, a)
2365
2366	sage: f = log(x*y)
2367	sage: f.simplify_radical()
2368	log(y) + log(x)
2369
2370	sage: f = (log(x+x^2)-log(x))^a/log(1+x)^(a/2)
2371	sage: f.simplify_radical()
2372	log(x + 1)^(a/2)
2373
2374	sage: f = (e^x-1)/(1+e^(x/2))
2375	sage: f.simplify_exp()
2376	e^(x/2) - 1
2377	"""
2378	return self.parent()(self._maxima_().radcan())
2379
2380	radical_simplify = simplify_log = log_simplify = simplify_radical
2381	simplify_exp = exp_simplify = simplify_radical
2382
2383	###################################################################
2384	# factor
2385	###################################################################
2386	def factor(self, dontfactor=[]):
2387	"""
2388	Factors self, containing any number of variables or functions,
2389	into factors irreducible over the integers.
2390
2391	INPUT:
2392	self -- a symbolic expression
2393	dontfactor -- list (default: []), a list of variables with
2394	respect to which factoring is not to occur.
2395	Factoring also will not take place with
2396	respect to any variables which are less
2397	important (using the variable ordering
2398	assumed for CRE form) than those on the
2399	`dontfactor' list.
2400
2401	EXAMPLES:
2402	sage: var('x, y, z')
2403	(x, y, z)
2404
2405	sage: (x^3-y^3).factor()
2406	-(y - x)(y^2 + xy + x^2)
2407	sage: factor(-8y - 4x + z^2(2y + x))
2408	(2y + x)(z - 2)*(z + 2)
2409	sage: f = -1 - 2x - x^2 + y^2 + 2xy^2 + x^2y^2
2410	sage: F = factor(f/(36(1 + 2y + y^2)), dontfactor=[x])
2411	sage: print F
2412	2
2413	(x + 2 x + 1) (y - 1)
2414	----------------------
2415	36 (y + 1)
2416	"""
2417	if len(dontfactor) > 0:
2418	m = self._maxima_()
2419	name = m.name()
2420	cmd = 'block([dontfactor:%s],factor(%s))'%(dontfactor, name)
2421	return evaled_symbolic_expression_from_maxima_string(cmd)
2422	else:
2423	return self.parent()(self._maxima_().factor())
2424
2425	###################################################################
2426	# expand
2427	###################################################################
2428	def expand(self):
2429	"""
2430	"""
2431	try:
2432	return self.__expand
2433	except AttributeError:
2434	e = self.parent()(self._maxima_().expand())
2435	self.__expand = e
2436	return e
2437
2438	def expand_trig(self):
2439	"""
2440	EXAMPLES:
2441	sage: sin(5*x).expand_trig()
2442	sin(x)^5 - 10cos(x)^2sin(x)^3 + 5cos(x)^4sin(x)
2443
2444	sage: cos(2*x + var('y')).trig_expand()
2445	cos(2x)cos(y) - sin(2x)sin(y)
2446
2447	ALIAS: trig_expand and expand_trig are the same
2448	"""
2449	return self.parent()(self._maxima_().trigexpand())
2450
2451	trig_expand = expand_trig
2452
2453
2454	###################################################################
2455	# substitute
2456	###################################################################
2457	def substitute(self, in_dict=None, **kwds):
2458	"""
2459	Takes the symbolic variables given as dict keys or as keywords and
2460	replaces them with the symbolic expressions given as dict values or as
2461	keyword values. Also run when you call a SymbolicExpression.
2462
2463	INPUT:
2464	in_dict -- (optional) dictionary of inputs
2465	**kwds -- named parameters
2466
2467	EXAMPLES:
2468	sage: x,y,t = var('x,y,t')
2469	sage: u = 3*y
2470	sage: u.substitute(y=t)
2471	3*t
2472
2473	sage: u = x^3 - 3y + 4t
2474	sage: u.substitute(x=y, y=t)
2475	y^3 + t
2476
2477	sage: f = sin(x)^2 + 32*x^(y/2)
2478	sage: f(x=2, y = 10)
2479	sin(2)^2 + 1024
2480
2481	sage: f(x=pi, y=t)
2482	32*pi^(t/2)
2483
2484	sage: f = x
2485	sage: f.variables()
2486	(x,)
2487
2488	sage: f = 2*x
2489	sage: f.variables()
2490	(x,)
2491
2492	sage: f = 2*x^2 - sin(x)
2493	sage: f.variables()
2494	(x,)
2495	sage: f(pi)
2496	2*pi^2
2497	sage: f(x=pi)
2498	2*pi^2
2499
2500	AUTHORS:
2501	-- Bobby Moretti: Initial version
2502	"""
2503	X = self.simplify()
2504	kwds = self.__parse_in_dict(in_dict, kwds)
2505	kwds = self.__varify_kwds(kwds)
2506	return X._recursive_sub(kwds)
2507
2508	def subs(self, args, *kwds):
2509	return self.substitute(args, *kwds)
2510
2511	def _recursive_sub(self, kwds):
2512	raise NotImplementedError, "implement _recursive_sub for type '%s'!"%(type(self))
2513
2514	def _recursive_sub_over_ring(self, kwds, ring):
2515	raise NotImplementedError, "implement _recursive_sub_over_ring for type '%s'!"%(type(self))
2516
2517	def __parse_in_dict(self, in_dict, kwds):
2518	if in_dict is None:
2519	return kwds
2520
2521	if not isinstance(in_dict, dict):
2522	if len(kwds) > 0:
2523	raise ValueError, "you must not both give the variable and specify it explicitly when doing a substitution."
2524	in_dict = SR(in_dict)
2525	vars = self.variables()
2526	#if len(vars) > 1:
2527	# raise ValueError, "you must specify the variable when doing a substitution, e.g., f(x=5)"
2528	if len(vars) == 0:
2529	return {}
2530	else:
2531	return {vars[0]: in_dict}
2532
2533	# merge dictionaries
2534	kwds.update(in_dict)
2535	return kwds
2536
2537	def __varify_kwds(self, kwds):
2538	return dict([(var(k),w) for k,w in kwds.iteritems()])
2539
2540	def substitute_over_ring(self, in_dict=None, ring=None, **kwds):
2541	X = self.simplify()
2542	kwds = self.__parse_in_dict(in_dict, kwds)
2543	kwds = self.__varify_kwds(kwds)
2544	if ring is None:
2545	return X._recursive_sub(kwds)
2546	else:
2547	return X._recursive_sub_over_ring(kwds, ring)
2548
2549
2550	###################################################################
2551	# Real and imaginary parts
2552	###################################################################
2553	def real(self):
2554	"""
2555	Return the real part of self.
2556
2557	EXAMPLES:
2558	sage: a = log(3+4*I)
2559	sage: print a
2560	log(4 I + 3)
2561	sage: print a.real()
2562	log(5)
2563	sage: print a.imag()
2564	4
2565	atan(-)
2566	3
2567
2568	Now make a and b symbolic and compute the general real part:
2569	sage: var('a,b')
2570	(a, b)
2571	sage: f = log(a + b*I)
2572	sage: f.real()
2573	log(b^2 + a^2)/2
2574	"""
2575	return self.parent()(self._maxima_().real())
2576
2577	def imag(self):
2578	"""
2579	Return the imaginary part of self.
2580
2581	EXAMPLES:
2582	sage: sqrt(-2).imag()
2583	sqrt(2)
2584
2585	We simplify Ln(Exp(z)) to z for -Pi<Im(z)<=Pi:
2586
2587	sage: z = var('z')
2588	sage: f = log(exp(z))
2589	sage: assume(-pi < imag(z))
2590	sage: assume(imag(z) <= pi)
2591	sage: print f
2592	z
2593	sage: forget()
2594
2595	A more symbolic example:
2596	sage: var('a, b')
2597	(a, b)
2598	sage: f = log(a + b*I)
2599	sage: f.imag()
2600	atan(b/a)
2601	"""
2602	return self.parent()(self._maxima_().imag())
2603
2604	def conjugate(self):
2605	"""
2606	The complex conjugate of self.
2607
2608	EXAMPLES:
2609	sage: a = 1 + 2*I
2610	sage: a.conjugate()
2611	1 - 2*I
2612	sage: a = sqrt(2) + 3^(1/3)*I; a
2613	3^(1/3)*I + sqrt(2)
2614	sage: a.conjugate()
2615	sqrt(2) - 3^(1/3)*I
2616	"""
2617	return self.parent()(self._maxima_().conjugate())
2618
2619	def norm(self):
2620	"""
2621	The complex norm of self, i.e., self times its complex conjugate.
2622
2623	EXAMPLES:
2624	sage: a = 1 + 2*I
2625	sage: a.norm()
2626	5
2627	sage: a = sqrt(2) + 3^(1/3)*I; a
2628	3^(1/3)*I + sqrt(2)
2629	sage: a.norm()
2630	3^(2/3) + 2
2631	sage: CDF(a).norm()
2632	4.08008382305
2633	sage: CDF(a.norm())
2634	4.08008382305
2635	"""
2636	m = self._maxima_()
2637	return self.parent()((m*m.conjugate()).expand())
2638
2639	###################################################################
2640	# Partial fractions
2641	###################################################################
2642	def partial_fraction(self, var=None):
2643	"""
2644	Return the partial fraction expansion of self with respect to
2645	the given variable.
2646
2647	INPUT:
2648	var -- variable name or string (default: first variable)
2649
2650	OUTPUT:
2651	Symbolic expression
2652
2653	EXAMPLES:
2654	sage: var('x')
2655	x
2656	sage: f = x^2/(x+1)^3
2657	sage: f.partial_fraction()
2658	1/(x + 1) - 2/(x + 1)^2 + 1/(x + 1)^3
2659	sage: print f.partial_fraction()
2660	1 2 1
2661	----- - -------- + --------
2662	x + 1 2 3
2663	(x + 1) (x + 1)
2664
2665	Notice that the first variable in the expression is used by default:
2666	sage: var('y')
2667	y
2668	sage: f = y^2/(y+1)^3
2669	sage: f.partial_fraction()
2670	1/(y + 1) - 2/(y + 1)^2 + 1/(y + 1)^3
2671
2672	sage: f = y^2/(y+1)^3 + x/(x-1)^3
2673	sage: f.partial_fraction()
2674	y^2/(y^3 + 3y^2 + 3y + 1) + 1/(x - 1)^2 + 1/(x - 1)^3
2675
2676	You can explicitly specify which variable is used.
2677	sage: f.partial_fraction(y)
2678	1/(y + 1) - 2/(y + 1)^2 + 1/(y + 1)^3 + x/(x^3 - 3x^2 + 3x - 1)
2679	"""
2680	if var is None:
2681	var = self._first_variable()
2682	return self.parent()(self._maxima_().partfrac(var))
2683
2684
2685
2686
2687	class Symbolic_object(SymbolicExpression):
2688	r"""
2689	A class representing a symbolic expression in terms of a SageObject (not
2690	necessarily a 'constant').
2691	"""
2692	def __init__(self, obj):
2693	SymbolicExpression.__init__(self)
2694	self._obj = obj
2695
2696	def __hash__(self):
2697	return hash(self._obj)
2698
2699	#def derivative(self, *args):
2700	# TODO: remove
2701	# return self.parent().zero_element()
2702
2703	def obj(self):
2704	"""
2705	EXAMPLES:
2706	"""
2707	return self._obj
2708
2709	def __float__(self):
2710	"""
2711	EXAMPLES:
2712	"""
2713	return float(self._obj)
2714
2715	def __complex__(self):
2716	"""
2717	EXAMPLES:
2718	"""
2719	return complex(self._obj)
2720
2721	def _mpfr_(self, field):
2722	"""
2723	EXAMPLES:
2724	"""
2725	return field(self._obj)
2726
2727	def _complex_mpfr_field_(self, C):
2728	"""
2729	EXAMPLES:
2730	"""
2731	return C(self._obj)
2732
2733	def _complex_double_(self, C):
2734	"""
2735	EXAMPLES:
2736	"""
2737	return C(self._obj)
2738
2739	def _real_double_(self, R):
2740	"""
2741	EXAMPLES:
2742	"""
2743	return R(self._obj)
2744
2745	def _real_rqdf_(self, R):
2746	"""
2747	EXAMPLES:
2748	"""
2749	return R(self._obj)
2750
2751	def _repr_(self, simplify=True):
2752	"""
2753	EXAMPLES:
2754	"""
2755	return repr(self._obj)
2756
2757	def _latex_(self):
2758	"""
2759	EXAMPLES:
2760	"""
2761	return latex(self._obj)
2762
2763	def str(self, bits=None):
2764	if bits is None:
2765	return str(self._obj)
2766	else:
2767	R = sage.rings.all.RealField(53)
2768	return str(R(self._obj))
2769
2770	def _maxima_init_(self):
2771	return maxima_init(self._obj)
2772
2773	def _sys_init_(self, system):
2774	return sys_init(self._obj, system)
2775
2776	def maxima_init(x):
2777	try:
2778	return x._maxima_init_()
2779	except AttributeError:
2780	return repr(x)
2781
2782	def sys_init(x, system):
2783	try:
2784	return x.__getattribute__('_%s_init_'%system)()
2785	except AttributeError:
2786	try:
2787	return x._system_init_(system)
2788	except AttributeError:
2789	return repr(x)
2790
2791	class SymbolicConstant(Symbolic_object):
2792	def __init__(self, x):
2793	from sage.rings.rational import Rational
2794	if isinstance(x, Rational):
2795	if x.is_integral() and x >= 0:
2796	self._precedence = 10**6
2797	else:
2798	self._precedence = 2000
2799	Symbolic_object.__init__(self, x)
2800	## if hasattr(x, 'parent'):
2801	## Symbolic_object.__init__(self, x)
2802	## elif isinstance(x, int):
2803	## Symbolic_object.__init__(self, Integer(x))
2804	## elif isinstance(x, float):
2805	## Symbolic_object.__init__(self, RealNumber(x))
2806	## else:
2807	## raise ValueError, "%s needs parent"%type(x)
2808
2809	#def _is_atomic(self):
2810	# try:
2811	# return self._atomic
2812	# except AttributeError:
2813	# if isinstance(self, Rational):
2814	# self._atomic = False
2815	# else:
2816	# self._atomic = True
2817	# return self._atomic
2818	def _is_atomic(self):
2819	try:
2820	return self._atomic
2821	except AttributeError:
2822	try:
2823	return self._obj._is_atomic()
2824	except AttributeError:
2825	if isinstance(self._obj, int):
2826	return True
2827
2828	def _recursive_sub(self, kwds):
2829	"""
2830	EXAMPLES:
2831	sage: a = SR(5/6)
2832	sage: type(a)
2833	<class 'sage.calculus.calculus.SymbolicConstant'>
2834	sage: a(x=3)
2835	5/6
2836	"""
2837	return self
2838
2839	def _recursive_sub_over_ring(self, kwds, ring):
2840	return ring(self)
2841
2842	def _algebraic_(self, field):
2843	"""
2844	EXAMPLES:
2845	sage: a = SR(5/6)
2846	sage: AA(a)
2847	5/6
2848	sage: type(AA(a))
2849	<class 'sage.rings.qqbar.AlgebraicReal'>
2850	sage: QQbar(a)
2851	5/6
2852	sage: type(QQbar(a))
2853	<class 'sage.rings.qqbar.AlgebraicNumber'>
2854	sage: from sage.calculus.calculus import SymbolicConstant
2855	sage: i = SymbolicConstant(I)
2856	sage: AA(i)
2857	Traceback (most recent call last):
2858	...
2859	TypeError: Cannot coerce algebraic number with non-zero imaginary part to algebraic real
2860	sage: QQbar(i)
2861	1*I
2862	sage: phi = SymbolicConstant(golden_ratio)
2863	sage: AA(phi)
2864	[1.6180339887498946 .. 1.6180339887498950]
2865	sage: QQbar(phi)
2866	[1.6180339887498946 .. 1.6180339887498950]
2867	"""
2868
2869	# Out of the many kinds of things that can be in a SymbolicConstant,
2870	# we accept only rational numbers (or things that can be coerced
2871	# to rational) and a few instances of Constant.
2872
2873	if isinstance(self._obj, sage.functions.constants.Constant) \
2874	and hasattr(self._obj, '_algebraic_'):
2875	return self._obj._algebraic_(field)
2876	return field(Rational(self._obj))
2877
2878	def _add_(self, right):
2879	"""
2880	EXAMPLES:
2881	sage: SR = SymbolicExpressionRing()
2882	sage: a = SR(2)
2883	sage: b = a+2; b
2884	4
2885	sage: type(b)
2886	<class 'sage.calculus.calculus.SymbolicConstant'>
2887	sage: b = sum([a for i in range(1000)]); b
2888	2000
2889	sage: type(_)
2890	<class 'sage.calculus.calculus.SymbolicConstant'>
2891	"""
2892	if isinstance(right, SymbolicConstant):
2893	try:
2894	self_parent = self._obj.parent()
2895	right_parent = right._obj.parent()
2896	if self_parent != SR and right_parent != SR and ( self_parent.has_coerce_map_from(right_parent) or right_parent.has_coerce_map_from(self_parent) ):
2897	return SymbolicConstant( operator.add(self._obj, right._obj) )
2898
2899	except AttributeError:
2900	#Either self._obj or right._obj doesn't have a
2901	#parent method (like 'int' or 'float')
2902	return SymbolicConstant( operator.add(self._obj, right._obj) )
2903
2904
2905
2906	return SymbolicArithmetic([self, right], operator.add)
2907
2908	def _sub_(self, right):
2909	"""
2910	EXAMPLES:
2911	sage: SR = SymbolicExpressionRing()
2912	sage: a = SR(2)
2913	sage: b = a-2; b
2914	0
2915	sage: type(b)
2916	<class 'sage.calculus.calculus.SymbolicConstant'>
2917	sage: b = SR(2000)
2918	sage: for i in range(1000): b -= a;
2919	sage: b
2920	0
2921	sage: type(b)
2922	<class 'sage.calculus.calculus.SymbolicConstant'>
2923
2924	"""
2925	if isinstance(right, SymbolicConstant):
2926	try:
2927	self_parent = self._obj.parent()
2928	right_parent = right._obj.parent()
2929	if self_parent != SR and right_parent != SR and ( self_parent.has_coerce_map_from(right_parent) or right_parent.has_coerce_map_from(self_parent) ):
2930	return SymbolicConstant( operator.sub(self._obj, right._obj) )
2931
2932	except AttributeError:
2933	#Either self._obj or right._obj doesn't have a
2934	#parent method (like 'int' or 'float')
2935	return SymbolicConstant( operator.sub(self._obj, right._obj) )
2936
2937
2938	return SymbolicArithmetic([self, right], operator.sub)
2939
2940
2941
2942	def _mul_(self, right):
2943	"""
2944	EXAMPLES:
2945	sage: SR = SymbolicExpressionRing()
2946	sage: a = SR(2)
2947	sage: b = a*2; b
2948	4
2949	sage: type(b)
2950	<class 'sage.calculus.calculus.SymbolicConstant'>
2951	sage: prod([a for i in range(1000)])
2952	10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376
2953	sage: type(_)
2954	<class 'sage.calculus.calculus.SymbolicConstant'>
2955	"""
2956	if isinstance(right, SymbolicConstant):
2957	try:
2958	self_parent = self._obj.parent()
2959	right_parent = right._obj.parent()
2960	if self_parent != SR and right_parent != SR and ( self_parent.has_coerce_map_from(right_parent) or right_parent.has_coerce_map_from(self_parent) ):
2961	return SymbolicConstant( operator.mul(self._obj, right._obj) )
2962
2963	except AttributeError:
2964	#Either self._obj or right._obj doesn't have a
2965	#parent method (like 'int' or 'float')
2966	return SymbolicConstant( operator.mul(self._obj, right._obj) )
2967
2968	return SymbolicArithmetic([self, right], operator.mul)
2969
2970	def _div_(self, right):
2971	"""
2972	EXAMPLES:
2973	sage: SR = SymbolicExpressionRing()
2974	sage: a = SR(2)
2975	sage: b = a/2; b
2976	1
2977	sage: type(b)
2978	<class 'sage.calculus.calculus.SymbolicConstant'>
2979	sage: b = SR(2^1000)
2980	sage: for i in range(1000): b /= a;
2981	sage: b
2982	1
2983	sage: type(b)
2984	<class 'sage.calculus.calculus.SymbolicConstant'>
2985
2986	"""
2987	if isinstance(right, SymbolicConstant):
2988	try:
2989	self_parent = self._obj.parent()
2990	right_parent = right._obj.parent()
2991	if self_parent != SR and right_parent != SR and ( self_parent.has_coerce_map_from(right_parent) or right_parent.has_coerce_map_from(self_parent) ):
2992	return SymbolicConstant( operator.div(self._obj, right._obj) )
2993
2994	except AttributeError:
2995	#Either self._obj or right._obj doesn't have a
2996	#parent method (like 'int' or 'float')
2997	return SymbolicConstant( operator.div(self._obj, right._obj) )
2998
2999
3000	return SymbolicArithmetic([self, right], operator.div)
3001
3002
3003
3004	def __pow__(self, right):
3005	"""
3006	EXAMPLES:
3007	sage: SR = SymbolicExpressionRing()
3008	sage: a = SR(2)
3009	sage: b = a^2; b
3010	4
3011	sage: type(b)
3012	<class 'sage.calculus.calculus.SymbolicArithmetic'>
3013	"""
3014	right = self.parent()(right)
3015	return SymbolicArithmetic([self, right], operator.pow)
3016
3017	class SymbolicPolynomial(Symbolic_object):
3018	"""
3019	An element of a polynomial ring as a formal symbolic expression.
3020
3021	EXAMPLES:
3022	A single variate polynomial:
3023	sage: R.<x> = QQ[]
3024	sage: f = SR(x^3 + x)
3025	sage: f(y=7)
3026	x^3 + x
3027	sage: f(x=5)
3028	130
3029	sage: f.integral(x)
3030	x^4/4 + x^2/2
3031	sage: f(x=var('y'))
3032	y^3 + y
3033
3034	A multivariate polynomial:
3035
3036	sage: R.<x,y,theta> = ZZ[]
3037	sage: f = SR(x^3 + x + y + theta^2); f
3038	x^3 + theta^2 + x + y
3039	sage: f(x=y, theta=y)
3040	y^3 + y^2 + 2*y
3041	sage: f(x=5)
3042	y + theta^2 + 130
3043
3044	The polynomial must be over a field of characteristic 0.
3045	sage: R.<w> = GF(7)[]
3046	sage: f = SR(w^3 + 1)
3047	Traceback (most recent call last):
3048	...
3049	TypeError: polynomial must be over a field of characteristic 0.
3050	"""
3051	# for now we do nothing except pass the info on to the superconstructor. It's
3052	# not clear to me why we need anything else in this class -Bobby
3053	def __init__(self, p):
3054	if p.parent().base_ring().characteristic() != 0:
3055	raise TypeError, "polynomial must be over a field of characteristic 0."
3056	Symbolic_object.__init__(self, p)
3057
3058	def _recursive_sub(self, kwds):
3059	return self._recursive_sub_over_ring(kwds, ring)
3060
3061	def _recursive_sub_over_ring(self, kwds, ring):
3062	f = self._obj
3063	if is_Polynomial(f):
3064	# Single variable case
3065	v = f.parent().variable_name()
3066	if kwds.has_key(v):
3067	return ring(f(kwds[v]))
3068	else:
3069	if not ring is SR:
3070	return ring(self)
3071	else:
3072	# Multivariable case
3073	t = []
3074	for g in f.parent().gens():
3075	s = repr(g)
3076	if kwds.has_key(s):
3077	t.append(kwds[s])
3078	else:
3079	t.append(g)
3080	return ring(f(*t))
3081
3082	def polynomial(self, base_ring):
3083	"""
3084	Return self as a polynomial over the given base ring, if possible.
3085
3086	INPUT:
3087	base_ring -- a ring
3088
3089	EXAMPLES:
3090	sage: R.<x> = QQ[]
3091	sage: f = SR(x^2 -2/3*x + 1)
3092	sage: f.polynomial(QQ)
3093	x^2 - 2/3*x + 1
3094	sage: f.polynomial(GF(19))
3095	x^2 + 12*x + 1
3096	"""
3097	f = self._obj
3098	if base_ring is f.base_ring():
3099	return f
3100	return f.change_ring(base_ring)
3101
3102	##################################################################
3103
3104
3105	zero_constant = SymbolicConstant(Integer(0))
3106
3107	class SymbolicOperation(SymbolicExpression):
3108	r"""
3109	A parent class representing any operation on SymbolicExpression objects.
3110	"""
3111	def __init__(self, operands):
3112	SymbolicExpression.__init__(self)
3113	self._operands = operands # don't even make a copy -- ok, since immutable.
3114
3115	def variables(self, vars=tuple([])):
3116	"""
3117	Return sorted list of variables that occur in the simplified
3118	form of self. The ordering is alphabetic.
3119
3120	EXAMPLES:
3121	sage: var('x,y,z,w,a,b,c')
3122	(x, y, z, w, a, b, c)
3123	sage: f = (x - x) + y^2 - z/z + (w^2-1)/(w+1); f
3124	y^2 + (w^2 - 1)/(w + 1) - 1
3125	sage: f.variables()
3126	(w, y)
3127
3128	sage: (x + y + z + a + b + c).variables()
3129	(a, b, c, x, y, z)
3130
3131	sage: (x^2 + x).variables()
3132	(x,)
3133	"""
3134	if not self.is_simplified():
3135	return self.simplify().variables(vars)
3136
3137	try:
3138	return self.__variables
3139	except AttributeError:
3140	pass
3141	vars = list(set(sum([list(op.variables()) for op in self._operands], list(vars))))
3142
3143	vars.sort(var_cmp)
3144	vars = tuple(vars)
3145	self.__variables = vars
3146	return vars
3147
3148	def var_cmp(x,y):
3149	return cmp(repr(x), repr(y))
3150
3151	symbols = {operator.add:' + ', operator.sub:' - ', operator.mul:'*',
3152	operator.div:'/', operator.pow:'^', operator.neg:'-'}
3153
3154
3155	class SymbolicArithmetic(SymbolicOperation):
3156	r"""
3157	Represents the result of an arithemtic operation on
3158	$f$ and $g$.
3159	"""
3160	def __init__(self, operands, op):
3161	SymbolicOperation.__init__(self, operands)
3162	self._operator = op
3163	# assume a really low precedence by default
3164	self._precedence = -1
3165	# set up associativity and precedence rules
3166	if op is operator.neg:
3167	self._binary = False
3168	self._unary = True
3169	self._precedence = 2000
3170	else:
3171	self._binary = True
3172	self._unary = False
3173	if op is operator.pow:
3174	self._precedence = 3000
3175	self._l_assoc = False
3176	self._r_assoc = True
3177	elif op is operator.mul:
3178	self._precedence = 2000
3179	self._l_assoc = True
3180	self._r_assoc = True
3181	elif op is operator.div:
3182	self._precedence = 2000
3183	self._l_assoc = True
3184	self._r_assoc = False
3185	elif op is operator.sub:
3186	self._precedence = 1000
3187	self._l_assoc = True
3188	self._r_assoc = False
3189	elif op is operator.add:
3190	self._precedence = 1000
3191	self._l_assoc = True
3192	self._r_assoc = True
3193
3194	def _recursive_sub(self, kwds):
3195	"""
3196	EXAMPLES:
3197	sage: var('x, y, z, w')
3198	(x, y, z, w)
3199	sage: f = (x - x) + y^2 - z/z + (w^2-1)/(w+1); f
3200	y^2 + (w^2 - 1)/(w + 1) - 1
3201	sage: f(y=10)
3202	(w^2 - 1)/(w + 1) + 99
3203	sage: f(w=1,y=10)
3204	99
3205	sage: f(y=w,w=y)
3206	(y^2 - 1)/(y + 1) + w^2 - 1
3207
3208	sage: f = y^5 - sqrt(2)
3209	sage: f(10)
3210	100000 - sqrt(2)
3211
3212	sage: a = x^2; b = a(2); b
3213	4
3214	sage: type(b)
3215	<class 'sage.calculus.calculus.SymbolicConstant'>
3216	"""
3217	ops = self._operands
3218	new_ops = [SR(op._recursive_sub(kwds)) for op in ops]
3219
3220	#Check to see if all of the new_ops are symbolic constants
3221	#If so, then we should return a symbolic constant.
3222	is_constant = all(map(lambda x: isinstance(x, SymbolicConstant), new_ops))
3223	if is_constant:
3224	return SymbolicConstant( self._operator(*map(lambda x: x._obj, new_ops)) )
3225	else:
3226	return self._operator(*new_ops)
3227
3228	def _recursive_sub_over_ring(self, kwds, ring):
3229	ops = self._operands
3230	if self._operator == operator.pow:
3231	new_ops = [ops[0]._recursive_sub_over_ring(kwds, ring=ring), Integer(ops[1])]
3232	else:
3233	new_ops = [op._recursive_sub_over_ring(kwds, ring=ring) for op in ops]
3234	return ring(self._operator(*new_ops))
3235
3236	def __float__(self):
3237	fops = [float(op) for op in self._operands]
3238	return self._operator(*fops)
3239
3240	def __complex__(self):
3241	fops = [complex(op) for op in self._operands]
3242	return self._operator(*fops)
3243
3244	def _mpfr_(self, field):
3245	rops = [op._mpfr_(field) for op in self._operands]
3246	return self._operator(*rops)
3247
3248	def _complex_mpfr_field_(self, field):
3249	rops = [op._complex_mpfr_field_(field) for op in self._operands]
3250	return self._operator(*rops)
3251
3252	def _complex_double_(self, field):
3253	if not self.is_simplified():
3254	return self.simplify()._complex_double_(field)
3255	rops = [op._complex_double_(field) for op in self._operands]
3256	return self._operator(*rops)
3257
3258	def _real_double_(self, field):
3259	if not self.is_simplified():
3260	return self.simplify()._real_double_(field)
3261	rops = [op._real_double_(field) for op in self._operands]
3262	return self._operator(*rops)
3263
3264	def _real_rqdf_(self, field):
3265	if not self.is_simplified():
3266	return self.simplify()._real_rqdf_(field)
3267	rops = [op._real_rqdf_(field) for op in self._operands]
3268	return self._operator(*rops)
3269
3270	def _algebraic_(self, field):
3271	"""
3272	Convert a symbolic expression to an algebraic number.
3273
3274	EXAMPLES:
3275	sage: QQbar(sqrt(2) + sqrt(8))
3276	[4.2426406871192847 .. 4.2426406871192857]
3277	sage: AA(sqrt(2) ^ 4) == 4
3278	True
3279	sage: AA(-golden_ratio)
3280	[-1.6180339887498950 .. -1.6180339887498946]
3281	sage: QQbar((2*I)^(1/2))
3282	[1.0000000000000000 .. 1.0000000000000000] + [1.0000000000000000 .. 1.0000000000000000]*I
3283	"""
3284
3285	# We try to avoid simplifying, because maxima's simplify command
3286	# can change the value of a radical expression (by changing which
3287	# root is selected).
3288	try:
3289	if self._operator is operator.pow:
3290	base = field(self._operands[0])
3291	expt = Rational(self._operands[1])
3292	return field(field(base) ** expt)
3293	else:
3294	rops = [field(op) for op in self._operands]
3295	return self._operator(*rops)
3296	except TypeError:
3297	if self.is_simplified():
3298	raise
3299	else:
3300	return self.simplify()._algebraic_(field)
3301
3302	def _is_atomic(self):
3303	try:
3304	return self._atomic
3305	except AttributeError:
3306	op = self._operator
3307	# todo: optimize with a dictionary
3308	if op is operator.add:
3309	self._atomic = False
3310	elif op is operator.sub:
3311	self._atomic = False
3312	elif op is operator.mul:
3313	self._atomic = True
3314	elif op is operator.div:
3315	self._atomic = False
3316	elif op is operator.pow:
3317	self._atomic = True
3318	elif op is operator.neg:
3319	self._atomic = False
3320	return self._atomic
3321
3322	def _repr_(self, simplify=True):
3323	"""
3324	TESTS:
3325	sage: var('r')
3326	r
3327	sage: a = (1-1/r)^(-1); a
3328	1/(1 - 1/r)
3329	sage: a.derivative(r)
3330	-1/((1 - 1/r)^2*r^2)
3331
3332	sage: var('a,b')
3333	(a, b)
3334	sage: s = 0(1/a) + -b(1/a)(1 + -10(1/a))(1/(ab + -1b*(1/a)))
3335	sage: s
3336	-b/(a(ab - b/a))
3337	sage: s(a=2,b=3)
3338	-1/3
3339	sage: -3/(2(23-(3/2)))
3340	-1/3
3341	sage: (-1)^(1/4)
3342	(-1)^(1/4)
3343	"""
3344	if simplify:
3345	if hasattr(self, '_simp'):
3346	return self._simp._repr_(simplify=False)
3347	else:
3348	return self.simplify()._repr_(simplify=False)
3349
3350	ops = self._operands
3351	op = self._operator
3352	s = [x._repr_(simplify=simplify) for x in ops]
3353
3354	# if an operand is a rational number, trick SAGE into thinking it's an
3355	# operation
3356	li = []
3357	for o in ops:
3358	try:
3359	obj = o._obj
3360	if isinstance(obj, Rational):
3361	temp = SymbolicConstant(obj)
3362	if not temp._obj.is_integral():
3363	temp._operator = operator.div
3364	temp._l_assoc = True
3365	temp._r_assoc = False
3366	temp._precedence = 2000
3367	temp._binary = True
3368	temp._unary = False
3369	li.append(temp)
3370	else:
3371	li.append(o)
3372	except AttributeError:
3373	li.append(o)
3374
3375	ops = li
3376
3377	rop = ops[0]
3378	if self._binary:
3379	lop = rop
3380	rop = ops[1]
3381
3382	lparens = True
3383	rparens = True
3384
3385	if self._binary:
3386	try:
3387	l_operator = lop._operator
3388	except AttributeError:
3389	# if it's not arithmetic on the left, see if it's atomic
3390	try:
3391	prec = lop._precedence
3392	except AttributeError:
3393	if lop._is_atomic():
3394	# if it has no concept of precedence, leave the parens
3395	lparens = False
3396	else:
3397	# if it a higher precedence, don't draw parens
3398	if self._precedence < lop._precedence:
3399	lparens = False
3400	else:
3401	# if the left op is the same is this operator
3402	if op is l_operator:
3403	# if it's left associative, get rid of the left parens
3404	if self._l_assoc:
3405	lparens = False
3406	# different operators, same precedence, get rid of the left parens
3407	elif self._precedence == lop._precedence:
3408	if self._l_assoc:
3409	lparens = False
3410	# if we have a lower precedence than the left, get rid of the parens
3411	elif self._precedence < lop._precedence:
3412	lparens = False
3413
3414	try:
3415	r_operator = rop._operator
3416	except AttributeError:
3417	try:
3418	prec = rop._precedence
3419	except AttributeError:
3420	if rop._is_atomic():
3421	rparens = False
3422	else:
3423	if self._precedence < rop._precedence:
3424	rparens = False
3425	else:
3426	if rop._binary:
3427	if op is r_operator:
3428	try:
3429	if self._r_assoc:
3430	rparens = False
3431	except AttributeError:
3432	pass
3433	elif self._precedence == rop._precedence:
3434	try:
3435	if self._r_assoc:
3436	rparens = False
3437	except AttributeError:
3438	pass
3439	# if the RHS has higher precedence, it comes first and parens are
3440	# redundant
3441	elif self._precedence < rop._precedence:
3442	rparens = False
3443	if self._binary:
3444	if lparens:
3445	s[0] = '(%s)'% s[0]
3446	if rparens:
3447	s[1] = '(%s)'% s[1]
3448
3449	return '%s%s%s' % (s[0], symbols[op], s[1])
3450
3451	elif self._unary:
3452	if rparens:
3453	s[0] = '(%s)'%s[0]
3454	return '%s%s' % (symbols[op], s[0])
3455
3456	def _latex_(self, simplify=True):
3457	# if we are not simplified, return the latex of a simplified version
3458	if simplify and not self.is_simplified():
3459	return self.simplify()._latex_()
3460	op = self._operator
3461	ops = self._operands
3462	s = []
3463	for x in self._operands:
3464	try:
3465	s.append(x._latex_(simplify=simplify))
3466	except TypeError:
3467	s.append(x._latex_())
3468	## what it used to be --> s = [x._latex_() for x in self._operands]
3469
3470	if op is operator.add:
3471	return '%s + %s' % (s[0], s[1])
3472	elif op is operator.sub:
3473	return '%s - %s' % (s[0], s[1])
3474	elif op is operator.mul:
3475	if ops[0]._has_op(operator.add) or ops[0]._has_op(operator.sub):
3476	s[0] = r'\left( %s \right)' %s[0]
3477	if ops[1]._has_op(operator.add) or ops[1]._has_op(operator.sub):
3478	s[1] = r'\left( %s \right)' %s[1]
3479	return '{%s \\cdot %s}' % (s[0], s[1])
3480	elif op is operator.div:
3481	return '\\frac{%s}{%s}' % (s[0], s[1])
3482	elif op is operator.pow:
3483	if ops[0]._has_op(operator.add) or ops[0]._has_op(operator.sub) \
3484	or ops[0]._has_op(operator.mul) or ops[0]._has_op(operator.div) \
3485	or ops[0]._has_op(operator.pow):
3486	s[0] = r'\left( %s \right)' % s[0]
3487	return '{%s}^{%s} ' % (s[0], s[1])
3488	elif op is operator.neg:
3489	return '-%s' % s[0]
3490
3491	def _maxima_init_(self):
3492	ops = self._operands
3493	if self._operator is operator.neg:
3494	return '-(%s)' % ops[0]._maxima_init_()
3495	else:
3496	return '(%s) %s (%s)' % (ops[0]._maxima_init_(),
3497	infixops[self._operator],
3498	ops[1]._maxima_init_())
3499
3500	def _sys_init_(self, system):
3501	ops = self._operands
3502	if self._operator is operator.neg:
3503	return '-(%s)' % sys_init(ops[0], system)
3504	else:
3505	return '(%s) %s (%s)' % (sys_init(ops[0], system),
3506	infixops[self._operator],
3507	sys_init(ops[1], system))
3508
3509	import re
3510
3511	def is_SymbolicVariable(x):
3512	return isinstance(x, SymbolicVariable)
3513
3514	class SymbolicVariable(SymbolicExpression):
3515	def __init__(self, name):
3516	SymbolicExpression.__init__(self)
3517	self._name = name
3518	if len(name) == 0:
3519	raise ValueError, "variable name must be nonempty"
3520
3521	def __hash__(self):
3522	return hash(self._name)
3523
3524	def _recursive_sub(self, kwds):
3525	# do the replacement if needed
3526	if kwds.has_key(self):
3527	return kwds[self]
3528	else:
3529	return self
3530
3531	def _recursive_sub_over_ring(self, kwds, ring):
3532	if kwds.has_key(self):
3533	return ring(kwds[self])
3534	else:
3535	return ring(self)
3536
3537	def variables(self, vars=tuple([])):
3538	"""
3539	Return sorted list of variables that occur in the simplified
3540	form of self.
3541	"""
3542	return (self, )
3543
3544	def __cmp__(self, right):
3545	if isinstance(right, SymbolicVariable):
3546	return cmp(self._repr_(), right._repr_())
3547	else:
3548	return SymbolicExpression.__cmp__(self, right)
3549
3550	def _repr_(self, simplify=True):
3551	return self._name
3552
3553	def __str__(self):
3554	return self._name
3555
3556	def _latex_(self):
3557	try:
3558	return self.__latex
3559	except AttributeError:
3560	pass
3561	self.__latex = latex_variable_name(self._name)
3562	return self.__latex
3563
3564	def _maxima_init_(self):
3565	return self._name
3566
3567	def _sys_init_(self, system):
3568	return self._name
3569
3570
3571	_vars = {}
3572	def var(s, create=True):
3573	r"""
3574	Create a symbolic variable with the name \emph{s}.
3575
3576	EXAMPLES:
3577	sage: var('xx')
3578	xx
3579	"""
3580	if isinstance(s, SymbolicVariable):
3581	return s
3582	s = str(s)
3583	if ',' in s:
3584	return tuple([var(x.strip()) for x in s.split(',')])
3585	elif ' ' in s:
3586	return tuple([var(x.strip()) for x in s.split()])
3587	try:
3588	v = _vars[s]
3589	_syms[s] = v
3590	return v
3591	except KeyError:
3592	if not create:
3593	raise ValueError, "the variable '%s' has not been defined"%var
3594	pass
3595	v = SymbolicVariable(s)
3596	_vars[s] = v
3597	_syms[s] = v
3598	return v
3599
3600
3601	#########################################################################################
3602	# Callable functions
3603	#########################################################################################
3604	def is_CallableSymbolicExpressionRing(x):
3605	"""
3606	EXAMPLES:
3607	sage: is_CallableSymbolicExpressionRing(QQ)
3608	False
3609	sage: var('x,y,z')
3610	(x, y, z)
3611	sage: is_CallableSymbolicExpressionRing(CallableSymbolicExpressionRing((x,y,z)))
3612	True
3613	"""
3614	return isinstance(x, CallableSymbolicExpressionRing_class)
3615
3616	class CallableSymbolicExpressionRing_class(CommutativeRing):
3617	def __init__(self, args):
3618	self._default_precision = 53 # default precision bits
3619	self._args = args
3620	ParentWithBase.__init__(self, RR)
3621
3622	def __call__(self, x):
3623	"""
3624
3625	TESTS:
3626	sage: f(x) = x+1; g(y) = y+1
3627	sage: f.parent()(g)
3628	x \|--> y + 1
3629	sage: g.parent()(f)
3630	y \|--> x + 1
3631	sage: f(x) = x+2y; g(y) = y+3x
3632	sage: f.parent()(g)
3633	x \|--> y + 3*x
3634	sage: g.parent()(f)
3635	y \|--> 2*y + x
3636	"""
3637	return self._coerce_impl(x)
3638
3639	def _coerce_impl(self, x):
3640	if isinstance(x, SymbolicExpression):
3641	if isinstance(x, CallableSymbolicExpression):
3642	x = x._expr
3643	return CallableSymbolicExpression(self, x)
3644	return self._coerce_try(x, [SR])
3645
3646	def _repr_(self):
3647	return "Callable function ring with arguments %s"%(self._args,)
3648
3649	def args(self):
3650	return self._args
3651
3652	def arguments(self):
3653	return self.args()
3654
3655	def zero_element(self):
3656	try:
3657	return self.__zero_element
3658	except AttributeError:
3659	z = CallableSymbolicExpression(self, SR.zero_element())
3660	self.__zero_element = z
3661	return z
3662
3663	def _an_element_impl(self):
3664	return CallableSymbolicExpression(self, SR._an_element())
3665
3666
3667	_cfr_cache = {}
3668	def CallableSymbolicExpressionRing(args, check=True):
3669	if check:
3670	if len(args) == 1 and isinstance(args[0], (list, tuple)):
3671	args = args[0]
3672	for arg in args:
3673	if not isinstance(arg, SymbolicVariable):
3674	raise TypeError, "Must construct a function with a tuple (or list) of" \
3675	+" SymbolicVariables."
3676	args = tuple(args)
3677	if _cfr_cache.has_key(args):
3678	R = _cfr_cache[args]()
3679	if not R is None:
3680	return R
3681	R = CallableSymbolicExpressionRing_class(args)
3682	_cfr_cache[args] = weakref.ref(R)
3683	return R
3684
3685	def is_CallableSymbolicExpression(x):
3686	"""
3687	Returns true if x is a callable symbolic expression.
3688
3689	EXAMPLES:
3690	sage: var('a x y z')
3691	(a, x, y, z)
3692	sage: f(x,y) = a + 2x + 3y + z
3693	sage: is_CallableSymbolicExpression(f)
3694	True
3695	sage: is_CallableSymbolicExpression(a+2*x)
3696	False
3697	sage: def foo(n): return n^2
3698	...
3699	sage: is_CallableSymbolicExpression(foo)
3700	False
3701	"""
3702	return isinstance(x, CallableSymbolicExpression)
3703
3704	class CallableSymbolicExpression(SymbolicExpression):
3705	r"""
3706	A callable symbolic expression that knows the ordered list of
3707	variables on which it depends.
3708
3709	EXAMPLES:
3710	sage: var('a, x, y, z')
3711	(a, x, y, z)
3712	sage: f(x,y) = a + 2x + 3y + z
3713	sage: f
3714	(x, y) \|--> z + 3y + 2x + a
3715	sage: f(1,2)
3716	z + a + 8
3717	"""
3718	def __init__(self, parent, expr):
3719	RingElement.__init__(self, parent)
3720	self._expr = expr
3721
3722	def variables(self):
3723	"""
3724	EXAMPLES:
3725	sage: a = var('a')
3726	sage: g(x) = sin(x) + a
3727	sage: g.variables()
3728	(a, x)
3729	sage: g.args()
3730	(x,)
3731	sage: g(y,x,z) = sin(x) + a - a
3732	sage: g
3733	(y, x, z) \|--> sin(x)
3734	sage: g.args()
3735	(y, x, z)
3736	"""
3737	return self._expr.variables()
3738
3739	def integral(self, x=None, a=None, b=None):
3740	"""
3741	Returns an integral of self.
3742	"""
3743	if a is None:
3744	return SymbolicExpression.integral(self, x, None, None)
3745	# if l. endpoint is None, compute an indefinite integral
3746	else:
3747	if x is None:
3748	x = self.default_variable()
3749	if not isinstance(x, SymbolicVariable):
3750	x = var(repr(x))
3751	# if we supplied an endpoint, then we want to return a number.
3752	return SR(self._maxima_().integrate(x, a, b))
3753
3754	integrate = integral
3755
3756	def expression(self):
3757	"""
3758	Return the underlying symbolic expression (i.e., forget the
3759	extra map structure).
3760	"""
3761	return self._expr
3762
3763	def args(self):
3764	return self.parent().args()
3765
3766	def arguments(self):
3767	return self.args()
3768
3769	def _maxima_init_(self):
3770	return self._expr._maxima_init_()
3771
3772	def __float__(self):
3773	return float(self._expr)
3774
3775	def __complex__(self):
3776	return complex(self._expr)
3777
3778	def _mpfr_(self, field):
3779	"""
3780	Coerce to a multiprecision real number.
3781
3782	EXAMPLES:
3783	sage: RealField(100)(SR(10))
3784	10.000000000000000000000000000
3785	"""
3786	return (self._expr)._mpfr_(field)
3787
3788	def _complex_mpfr_field_(self, field):
3789	return field(self._expr)
3790
3791	def _complex_double_(self, C):
3792	return C(self._expr)
3793
3794	def _real_double_(self, R):
3795	return R(self._expr)
3796
3797	def _real_rqdf_(self, R):
3798	return R(self._expr)
3799
3800	# TODO: should len(args) == len(vars)?
3801	def __call__(self, *args):
3802	vars = self.args()
3803	dct = dict( (vars[i], args[i]) for i in range(len(args)) )
3804	return self._expr.substitute(dct)
3805
3806	def _repr_(self, simplify=True):
3807	args = self.args()
3808	if len(args) == 1:
3809	return "%s \|--> %s" % (args[0], self._expr._repr_(simplify=simplify))
3810	else:
3811	args = ", ".join(map(str, args))
3812	return "(%s) \|--> %s" % (args, self._expr._repr_(simplify=simplify))
3813
3814	def _latex_(self):
3815	args = self.args()
3816	if len(args) == 1:
3817	return "%s \\ {\mapsto}\\ %s" % (args[0],
3818	self._expr._latex_())
3819	else:
3820	vars = ", ".join(map(latex, args))
3821	# the weird TeX is to workaround an apparent JsMath bug
3822	return "\\left(%s \\right)\\ {\\mapsto}\\ %s" % (args, self._expr._latex_())
3823
3824	def _neg_(self):
3825	return CallableSymbolicExpression(self.parent(), -self._expr)
3826
3827	def __add__(self, right):
3828	"""
3829	EXAMPLES:
3830	sage: var('x y z n m')
3831	(x, y, z, n, m)
3832	sage: f(x,n,y) = x^n + y^m; g(x,n,m,z) = x^n +z^m
3833	sage: f + g
3834	(x, n, m, y, z) \|--> z^m + y^m + 2*x^n
3835	sage: g + f
3836	(x, n, m, y, z) \|--> z^m + y^m + 2*x^n
3837	"""
3838	if isinstance(right, CallableSymbolicExpression):
3839	if self.parent() is right.parent():
3840	return self._add_(right)
3841	else:
3842	args = self._unify_args(right)
3843	R = CallableSymbolicExpressionRing(args)
3844	return R(self) + R(right)
3845	else:
3846	return RingElement.__add__(self, right)
3847
3848	def _add_(self, right):
3849	return CallableSymbolicExpression(self.parent(), self._expr + right._expr)
3850
3851	def __sub__(self, right):
3852	"""
3853	EXAMPLES:
3854	sage: var('x y z n m')
3855	(x, y, z, n, m)
3856	sage: f(x,n,y) = x^n + y^m; g(x,n,m,z) = x^n +z^m
3857	sage: f - g
3858	(x, n, m, y, z) \|--> y^m - z^m
3859	sage: g - f
3860	(x, n, m, y, z) \|--> z^m - y^m
3861	"""
3862	if isinstance(right, CallableSymbolicExpression):
3863	if self.parent() is right.parent():
3864	return self._sub_(right)
3865	else:
3866	args = self._unify_args(right)
3867	R = CallableSymbolicExpressionRing(args)
3868	return R(self) - R(right)
3869	else:
3870	return RingElement.__sub__(self, right)
3871
3872	def _sub_(self, right):
3873	return CallableSymbolicExpression(self.parent(), self._expr - right._expr)
3874
3875	def __mul__(self, right):
3876	"""
3877	EXAMPLES:
3878	sage: var('x y z a b c n m')
3879	(x, y, z, a, b, c, n, m)
3880
3881	sage: f(x) = x+2y; g(y) = y+3x
3882	sage: f*(2/3)
3883	x \|--> 2(2y + x)/3
3884	sage: f*g
3885	(x, y) \|--> (y + 3x)(2*y + x)
3886	sage: (2/3)*f
3887	x \|--> 2(2y + x)/3
3888
3889	sage: f(x,y,z,a,b) = x+y+z-a-b; f
3890	(x, y, z, a, b) \|--> z + y + x - b - a
3891	sage: f * (b*c)
3892	(x, y, z, a, b) \|--> bc(z + y + x - b - a)
3893	sage: g(x,y,w,t) = xyw*t
3894	sage: f*g
3895	(x, y, a, b, t, w, z) \|--> twxy(z + y + x - b - a)
3896	sage: (f*g)(2,3)
3897	6tw*(z - b - a + 5)
3898
3899	sage: f(x,n,y) = x^n + y^m; g(x,n,m,z) = x^n +z^m
3900	sage: f * g
3901	(x, n, m, y, z) \|--> (y^m + x^n)*(z^m + x^n)
3902	sage: g * f
3903	(x, n, m, y, z) \|--> (y^m + x^n)*(z^m + x^n)
3904	"""
3905	if isinstance(right, CallableSymbolicExpression):
3906	if self.parent() is right.parent():
3907	return self._mul_(right)
3908	else:
3909	args = self._unify_args(right)
3910	R = CallableSymbolicExpressionRing(args)
3911	return R(self)*R(right)
3912	else:
3913	return RingElement.__mul__(self, right)
3914
3915	def _mul_(self, right):
3916	return CallableSymbolicExpression(self.parent(), self._expr * right._expr)
3917
3918	def __div__(self, right):
3919	"""
3920	EXAMPLES:
3921	sage: var('x,y,z,m,n')
3922	(x, y, z, m, n)
3923	sage: f(x,n,y) = x^n + y^m; g(x,n,m,z) = x^n +z^m
3924	sage: f / g
3925	(x, n, m, y, z) \|--> (y^m + x^n)/(z^m + x^n)
3926	sage: g / f
3927	(x, n, m, y, z) \|--> (z^m + x^n)/(y^m + x^n)
3928	"""
3929	if isinstance(right, CallableSymbolicExpression):
3930	if self.parent() is right.parent():
3931	return self._div_(right)
3932	else:
3933	args = self._unify_args(right)
3934	R = CallableSymbolicExpressionRing(args)
3935	return R(self) / R(right)
3936	else:
3937	return RingElement.__div__(self, right)
3938
3939	def _div_(self, right):
3940	return CallableSymbolicExpression(self.parent(), self._expr / right._expr)
3941
3942	def __pow__(self, right):
3943	return CallableSymbolicExpression(self.parent(), self._expr ** right)
3944
3945	def _unify_args(self, x):
3946	r"""
3947	Takes the variable list from another CallableSymbolicExpression object and
3948	compares it with the current CallableSymbolicExpression object's variable list,
3949	combining them according to the following rules:
3950
3951	Let \code{a} be \code{self}'s variable list, let \code{b} be
3952	\code{y}'s variable list.
3953
3954	\begin{enumerate}
3955
3956	\item If \code{a == b}, then the variable lists are identical,
3957	so return that variable list.
3958
3959	\item If \code{a} $\neq$ \code{b}, then check if the first $n$
3960	items in \code{a} are the first $n$ items in \code{b},
3961	or vice-versa. If so, return a list with these $n$
3962	items, followed by the remaining items in \code{a} and
3963	\code{b} sorted together in alphabetical order.
3964
3965	\end{enumerate}
3966
3967	Note: When used for arithmetic between CallableSymbolicExpressions,
3968	these rules ensure that the set of CallableSymbolicExpressions will have
3969	certain properties. In particular, it ensures that the set is
3970	a \emph{commutative} ring, i.e., the order of the input
3971	variables is the same no matter in which order arithmetic is
3972	done.
3973
3974	INPUT:
3975	x -- A CallableSymbolicExpression
3976
3977	OUTPUT:
3978	A tuple of variables.
3979
3980	EXAMPLES:
3981	sage: f(x, y, z) = sin(x+y+z)
3982	sage: f
3983	(x, y, z) \|--> sin(z + y + x)
3984	sage: g(x, y) = y + 2*x
3985	sage: g
3986	(x, y) \|--> y + 2*x
3987	sage: f._unify_args(g)
3988	(x, y, z)
3989	sage: g._unify_args(f)
3990	(x, y, z)
3991
3992	sage: f(x, y, z) = sin(x+y+z)
3993	sage: g(w, t) = cos(w - t)
3994	sage: g
3995	(w, t) \|--> cos(w - t)
3996	sage: f._unify_args(g)
3997	(t, w, x, y, z)
3998
3999	sage: f(x, y, t) = y*(x^2-t)
4000	sage: f
4001	(x, y, t) \|--> (x^2 - t)*y
4002	sage: g(x, y, w) = x + y - cos(w)
4003	sage: f._unify_args(g)
4004	(x, y, t, w)
4005	sage: g._unify_args(f)
4006	(x, y, t, w)
4007	sage: f*g
4008	(x, y, t, w) \|--> (x^2 - t)y(y + x - cos(w))
4009
4010	sage: f(x,y, t) = x+y
4011	sage: g(x, y, w) = w + t
4012	sage: f._unify_args(g)
4013	(x, y, t, w)
4014	sage: g._unify_args(f)
4015	(x, y, t, w)
4016	sage: f + g
4017	(x, y, t, w) \|--> y + x + w + t
4018
4019	AUTHORS:
4020	-- Bobby Moretti, thanks to William Stein for the rules
4021
4022	"""
4023	a = self.args()
4024	b = x.args()
4025
4026	# Rule #1
4027	if [str(x) for x in a] == [str(x) for x in b]:
4028	return a
4029
4030	# Rule #2
4031	new_list = []
4032	done = False
4033	i = 0
4034	while not done and i < min(len(a), len(b)):
4035	if var_cmp(a[i], b[i]) == 0:
4036	new_list.append(a[i])
4037	i += 1
4038	else:
4039	done = True
4040
4041	temp = set([])
4042	# Sorting remaining variables.
4043	for j in range(i, len(a)):
4044	if not a[j] in temp:
4045	temp.add(a[j])
4046
4047	for j in range(i, len(b)):
4048	if not b[j] in temp:
4049	temp.add(b[j])
4050
4051	temp = list(temp)
4052	temp.sort(var_cmp)
4053	new_list.extend(temp)
4054	return tuple(new_list)
4055
4056
4057	#########################################################################################
4058	# End callable functions
4059	#########################################################################################
4060
4061	class SymbolicComposition(SymbolicOperation):
4062	r"""
4063	Represents the symbolic composition of $f \circ g$.
4064	"""
4065	def __init__(self, f, g):
4066	"""
4067	INPUT:
4068	f, g -- both must be in the symbolic expression ring.
4069	"""
4070	SymbolicOperation.__init__(self, [f,g])
4071
4072	def _recursive_sub(self, kwds):
4073	ops = self._operands
4074	return ops[0](ops[1]._recursive_sub(kwds))
4075
4076	def _recursive_sub_over_ring(self, kwds, ring):
4077	ops = self._operands
4078	return ring(ops[0](ops[1]._recursive_sub_over_ring(kwds, ring=ring)))
4079
4080	def _is_atomic(self):
4081	return True
4082
4083	def _repr_(self, simplify=True):
4084	if simplify:
4085	if hasattr(self, '_simp'):
4086	return self._simp._repr_(simplify=False)
4087	else:
4088	return self.simplify()._repr_(simplify=False)
4089	ops = self._operands
4090	try:
4091	return ops[0]._repr_evaled_(ops[1]._repr_(simplify=False))
4092	except AttributeError:
4093	return "%s(%s)"% (ops[0]._repr_(simplify=False), ops[1]._repr_(simplify=False))
4094
4095	def _maxima_init_(self):
4096	ops = self._operands
4097	try:
4098	return ops[0]._maxima_init_evaled_(ops[1]._maxima_init_())
4099	except AttributeError:
4100	return '%s(%s)' % (ops[0]._maxima_init_(), ops[1]._maxima_init_())
4101
4102	def _latex_(self):
4103	if not self.is_simplified():
4104	return self.simplify()._latex_()
4105	ops = self._operands
4106
4107	#Check to see if the function has a _latex_composition method
4108	if hasattr(ops[0], '_latex_composition'):
4109	return ops[0]._latex_composition(ops[1])
4110
4111	# certain functions (such as \sqrt) need braces in LaTeX
4112	if (ops[0]).tex_needs_braces():
4113	return r"%s{ %s }" % ( (ops[0])._latex_(), (ops[1])._latex_())
4114	# ... while others (such as \cos) don't
4115	return r"%s \left( %s \right)"%((ops[0])._latex_(),(ops[1])._latex_())
4116
4117	def _sys_init_(self, system):
4118	ops = self._operands
4119	return '%s(%s)' % (sys_init(ops[0],system), sys_init(ops[1],system))
4120
4121	def _mathematica_init_(self):
4122	system = 'mathematica'
4123	ops = self._operands
4124	return '%s[%s]' % (sys_init(ops[0],system).capitalize(), sys_init(ops[1],system))
4125
4126	def __float__(self):
4127	f = self._operands[0]
4128	g = self._operands[1]
4129	return float(f._approx_(float(g)))
4130
4131	def __complex__(self):
4132	f = self._operands[0]
4133	g = self._operands[1]
4134	return complex(f._approx_(float(g)))
4135
4136	def _mpfr_(self, field):
4137	"""
4138	Coerce to a multiprecision real number.
4139
4140	EXAMPLES:
4141	sage: RealField(100)(sin(2)+cos(2))
4142	0.49315059027853930839845163641
4143
4144	sage: RR(sin(pi))
4145	0.000000000000000122464679914735
4146
4147	sage: type(RR(sqrt(163)*pi))
4148	<type 'sage.rings.real_mpfr.RealNumber'>
4149
4150	sage: RR(coth(pi))
4151	1.00374187319732
4152	sage: RealField(100)(coth(pi))
4153	1.0037418731973212882015526912
4154	sage: RealField(200)(acos(1/10))
4155	1.4706289056333368228857985121870581235299087274579233690964
4156	"""
4157	f = self._operands[0]
4158	g = self._operands[1]
4159	x = f(g._mpfr_(field))
4160	if isinstance(x, SymbolicExpression):
4161	if field.prec() <= 53:
4162	return field(float(x))
4163	else:
4164	raise TypeError, "precision loss"
4165	else:
4166	return x
4167
4168
4169	def _complex_mpfr_field_(self, field):
4170	"""
4171	Coerce to a multiprecision complex number.
4172
4173	EXAMPLES:
4174	sage: ComplexField(100)(sin(2)+cos(2)+I)
4175	0.49315059027853930839845163641 + 1.0000000000000000000000000000*I
4176
4177	"""
4178	f = self._operands[0]
4179	g = self._operands[1]
4180	x = f(g._complex_mpfr_field_(field))
4181	if isinstance(x, SymbolicExpression):
4182	if field.prec() <= 53:
4183	return field(complex(x))
4184	else:
4185	raise TypeError, "precision loss"
4186	else:
4187	return x
4188
4189	def _complex_double_(self, field):
4190	"""
4191	Coerce to a complex double.
4192
4193	EXAMPLES:
4194	sage: CDF(sin(2)+cos(2)+I)
4195	0.493150590279 + 1.0*I
4196	sage: CDF(coth(pi))
4197	1.0037418732
4198	"""
4199	f = self._operands[0]
4200	g = self._operands[1]
4201	z = f(g._complex_double_(field))
4202	if isinstance(z, SymbolicExpression):
4203	return field(complex(z))
4204	return z
4205
4206	def _real_double_(self, field):
4207	"""
4208	Coerce to a real double.
4209
4210	EXAMPLES:
4211	sage: RDF(sin(2)+cos(2))
4212	0.493150590279
4213	"""
4214	if not self.is_simplified():
4215	return self.simplify()._real_double_(field)
4216	f = self._operands[0]
4217	g = self._operands[1]
4218	z = f(g._real_double_(field))
4219	if isinstance(z, SymbolicExpression):
4220	return field(float(z))
4221	return z
4222
4223	def _real_rqdf_(self, field):
4224	"""
4225	Coerce to a real qdrf.
4226
4227	EXAMPLES:
4228
4229	"""
4230	if not self.is_simplified():
4231	return self.simplify()._real_rqdf_(field)
4232	f = self._operands[0]
4233	g = self._operands[1]
4234	z = f(g._real_rqdf_(field))
4235	if isinstance(z, SymbolicExpression):
4236	raise TypeError, "precision loss"
4237	else:
4238	return z
4239
4240	def _algebraic_(self, field):
4241	"""
4242	Coerce to an algebraic number.
4243
4244	EXAMPLES:
4245	sage: QQbar(sqrt(2))
4246	[1.4142135623730949 .. 1.4142135623730952]
4247	sage: AA(abs(1+I))
4248	[1.4142135623730949 .. 1.4142135623730952]
4249	"""
4250	# We try to avoid simplifying, because maxima's simplify command
4251	# can change the value of a radical expression (by changing which
4252	# root is selected).
4253	f = self._operands[0]
4254	g = self._operands[1]
4255	try:
4256	return field(f(g._algebraic_(field)))
4257	except TypeError:
4258	if self.is_simplified():
4259	raise
4260	else:
4261	return self.simplify()._algebraic_(field)
4262
4263
4264	class PrimitiveFunction(SymbolicExpression):
4265	def __init__(self, needs_braces=False):
4266	SymbolicExpression.__init__(self)
4267	self._tex_needs_braces = needs_braces
4268
4269	def _recursive_sub(self, kwds):
4270	if kwds.has_key(self):
4271	return kwds[self]
4272	return self
4273
4274	def _recursive_sub_over_ring(self, kwds, ring):
4275	if kwds.has_key(self):
4276	return kwds[self]
4277	return self
4278
4279	def plot(self, args, *kwds):
4280	f = self(var('x'))
4281	return SymbolicExpression.plot(f, args, *kwds)
4282
4283	def _is_atomic(self):
4284	return True
4285
4286	def tex_needs_braces(self):
4287	return self._tex_needs_braces
4288
4289	def __call__(self, x, *args):
4290	if isinstance(x, float):
4291	return self._approx_(x)
4292	try:
4293	return getattr(x, self._repr_())(*args)
4294	except AttributeError:
4295	return SymbolicComposition(self, SR(x))
4296
4297	def _approx_(self, x): # must always be called with a float x as input.
4298	s = '%s(%s), numer'%(self._repr_(), float(x))
4299	return float(maxima.eval(s))
4300
4301	_syms = {}
4302
4303	class Function_erf(PrimitiveFunction):
4304	r"""
4305	The error function, defined as $\text{erf}(x) =
4306	\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2} dt$.
4307	"""
4308
4309	def _repr_(self, simplify=True):
4310	return "erf"
4311
4312	def _latex_(self):
4313	return "\\text{erf}"
4314
4315	def _approx_(self, x):
4316	return float(1 - pari(float(x)).erfc())
4317
4318
4319	erf = Function_erf()
4320	_syms['erf'] = erf
4321
4322	class Function_abs(PrimitiveFunction):
4323	"""
4324	The absolute value function.
4325
4326	EXAMPLES:
4327	sage: var('x y')
4328	(x, y)
4329	sage: abs(x)
4330	abs(x)
4331	sage: abs(x^2 + y^2)
4332	y^2 + x^2
4333	sage: abs(-2)
4334	2
4335	sage: sqrt(x^2)
4336	abs(x)
4337	sage: abs(sqrt(x))
4338	sqrt(x)
4339	"""
4340	def _repr_(self, simplify=True):
4341	return "abs"
4342
4343	def _latex_(self):
4344	return "\\mathrm{abs}"
4345
4346	def _latex_composition(self, x):
4347	"""
4348	sage: f = sage.calculus.calculus.Function_abs()
4349	sage: latex(f)
4350	\mathrm{abs}
4351	sage: latex(abs(x))
4352	\left\| x \right\|
4353	"""
4354	return "\\left\| " + latex(x) + " \\right\|"
4355
4356	def _approx_(self, x):
4357	return float(x.__abs__())
4358
4359	def __call__(self, x): # special case
4360	return SymbolicComposition(self, SR(x))
4361
4362	abs_symbolic = Function_abs()
4363	_syms['abs'] = abs_symbolic
4364
4365
4366
4367	class Function_ceil(PrimitiveFunction):
4368	"""
4369	The ceiling function.
4370
4371	EXAMPLES:
4372	sage: a = ceil(2/5 + x)
4373	sage: a
4374	ceil(x + 2/5)
4375	sage: a(4)
4376	5
4377	sage: a(4.0)
4378	5
4379	sage: ZZ(a(3))
4380	4
4381	sage: a = ceil(x^3 + x + 5/2)
4382	sage: a
4383	ceil(x^3 + x + 1/2) + 2
4384	sage: a(x=2)
4385	13
4386
4387	sage: ceil(5.4)
4388	6
4389	sage: type(ceil(5.4))
4390	<type 'sage.rings.integer.Integer'>
4391	"""
4392	def _repr_(self, simplify=True):
4393	return "ceil"
4394
4395	def _latex_(self):
4396	return "\\text{ceil}"
4397
4398	def _maxima_init_(self):
4399	return "ceiling"
4400
4401	_approx_ = math.ceil
4402
4403	def __call__(self, x):
4404	try:
4405	return x.ceil()
4406	except AttributeError:
4407	if isinstance(x, (float, int, long, complex)):
4408	return int(math.ceil(x))
4409	return SymbolicComposition(self, SR(x))
4410
4411	ceil = Function_ceil()
4412	_syms['ceiling'] = ceil # spelled ceiling in maxima
4413
4414
4415	class Function_floor(PrimitiveFunction):
4416	"""
4417	The floor function.
4418
4419	EXAMPLES:
4420	sage: floor(5.4)
4421	5
4422	sage: type(floor(5.4))
4423	<type 'sage.rings.integer.Integer'>
4424	sage: var('x')
4425	x
4426	sage: a = floor(5.4 + x); a
4427	floor(x + 0.4000000000000004) + 5
4428	sage: a(2)
4429	7
4430	"""
4431	def _repr_(self, simplify=True):
4432	return "floor"
4433
4434	def _latex_(self):
4435	return "\\text{floor}"
4436
4437	def _maxima_init_(self):
4438	return "floor"
4439
4440	_approx_ = math.floor
4441
4442	def __call__(self, x):
4443	try:
4444	return x.floor()
4445	except AttributeError:
4446	if isinstance(x, (float, int, long, complex)):
4447	return int(math.floor(x))
4448	return SymbolicComposition(self, SR(x))
4449
4450	floor = Function_floor()
4451	_syms['floor'] = floor # spelled ceiling in maxima
4452
4453
4454	class Function_sin(PrimitiveFunction):
4455	"""
4456	The sine function
4457	"""
4458	def _repr_(self, simplify=True):
4459	return "sin"
4460
4461	def _latex_(self):
4462	return "\\sin"
4463
4464	_approx_ = math.sin
4465
4466	def __call__(self, x):
4467	try:
4468	return x.sin()
4469	except AttributeError:
4470	if isinstance(x, float):
4471	return math.sin(x)
4472	return SymbolicComposition(self, SR(x))
4473
4474	sin = Function_sin()
4475	_syms['sin'] = sin
4476
4477	class Function_cos(PrimitiveFunction):
4478	"""
4479	The cosine function
4480	"""
4481	def _repr_(self, simplify=True):
4482	return "cos"
4483
4484	def _latex_(self):
4485	return "\\cos"
4486
4487	_approx_ = math.cos
4488
4489	def __call__(self, x):
4490	try:
4491	return x.cos()
4492	except AttributeError:
4493	if isinstance(x, float):
4494	return math.cos(x)
4495	return SymbolicComposition(self, SR(x))
4496
4497
4498	cos = Function_cos()
4499	_syms['cos'] = cos
4500
4501	class Function_sec(PrimitiveFunction):
4502	"""
4503	The secant function
4504
4505	EXAMPLES:
4506	sage: sec(pi/4)
4507	sqrt(2)
4508	sage: RR(sec(pi/4))
4509	1.41421356237309
4510	sage: n(sec(pi/4),100)
4511	1.4142135623730950488016887242
4512	sage: sec(1/2)
4513	sec(1/2)
4514	sage: sec(0.5)
4515	1.13949392732455
4516	"""
4517	def _repr_(self, simplify=True):
4518	return "sec"
4519
4520	def _latex_(self):
4521	return "\\sec"
4522
4523	def _approx_(self, x):
4524	return 1/math.cos(x)
4525
4526	sec = Function_sec()
4527	_syms['sec'] = sec
4528
4529	class Function_csc(PrimitiveFunction):
4530	"""
4531	The cosecant function.
4532
4533	EXAMPLES:
4534	sage: csc(pi/4)
4535	sqrt(2)
4536	sage: RR(csc(pi/4))
4537	1.41421356237310
4538	sage: n(csc(pi/4),100)
4539	1.4142135623730950488016887242
4540	sage: csc(1/2)
4541	csc(1/2)
4542	sage: csc(0.5)
4543	2.08582964293349
4544	"""
4545	def _repr_(self, simplify=True):
4546	return "csc"
4547
4548	def _latex_(self):
4549	return "\\csc"
4550
4551	def _approx_(self, x):
4552	return 1/math.sin(x)
4553
4554	csc = Function_csc()
4555	_syms['csc'] = csc
4556
4557	class Function_cot(PrimitiveFunction):
4558	"""
4559	The cotangent function.
4560
4561	EXAMPLES:
4562	sage: cot(pi/4)
4563	1
4564	sage: RR(cot(pi/4))
4565	1.00000000000000
4566	sage: n(cot(pi/4),100)
4567	1.0000000000000000000000000000
4568	sage: cot(1/2)
4569	cot(1/2)
4570	sage: cot(0.5)
4571	1.83048772171245
4572	"""
4573	def _repr_(self, simplify=True):
4574	return "cot"
4575
4576	def _latex_(self):
4577	return "\\cot"
4578
4579	def _approx_(self, x):
4580	return 1/math.tan(x)
4581
4582	cot = Function_cot()
4583	_syms['cot'] = cot
4584
4585	class Function_tan(PrimitiveFunction):
4586	"""
4587	The tangent function
4588
4589	EXAMPLES:
4590	sage: tan(pi)
4591	0
4592	sage: tan(3.1415)
4593	-0.0000926535900581913
4594	sage: tan(3.1415/4)
4595	0.999953674278156
4596	sage: tan(pi/4)
4597	1
4598	sage: tan(1/2)
4599	tan(1/2)
4600	sage: RR(tan(1/2))
4601	0.546302489843790
4602	"""
4603	def _repr_(self, simplify=True):
4604	return "tan"
4605
4606	def _latex_(self):
4607	return "\\tan"
4608
4609	def _approx_(self, x):
4610	return math.tan(x)
4611
4612	tan = Function_tan()
4613	_syms['tan'] = tan
4614
4615	def atan2(y, x):
4616	return atan(y/x)
4617
4618	_syms['atan2'] = atan2
4619
4620	class Function_asin(PrimitiveFunction):
4621	"""
4622	The arcsine function
4623
4624	EXAMPLES:
4625	sage: asin(0.5)
4626	0.523598775598299
4627	sage: asin(1/2)
4628	pi/6
4629	sage: asin(1 + I*1.0)
4630	1.061275061905036*I + 0.6662394324925153
4631	"""
4632	def _repr_(self, simplify=True):
4633	return "asin"
4634
4635	def _latex_(self):
4636	return "\\sin^{-1}"
4637
4638	def _approx_(self, x):
4639	return math.asin(x)
4640
4641	asin = Function_asin()
4642	_syms['asin'] = asin
4643
4644	class Function_asinh(PrimitiveFunction):
4645	"""
4646	The inverse of the hyperbolic sine function.
4647
4648	EXAMPLES:
4649	sage: asinh(0.5)
4650	0.481211825059603
4651	sage: asinh(1/2)
4652	asinh(1/2)
4653	sage: asinh(1 + I*1.0)
4654	0.6662394324925153*I + 1.061275061905036
4655	"""
4656	def _repr_(self, simplify=True):
4657	return "asinh"
4658
4659	def _latex_(self):
4660	return "\\sinh^{-1}"
4661
4662	def _approx_(self, x):
4663	return float(pari(float(x)).asinh())
4664
4665	asinh = Function_asinh()
4666	_syms['asinh'] = asinh
4667
4668	class Function_acosh(PrimitiveFunction):
4669	"""
4670	The inverse of the hyperbolic cose function.
4671
4672	EXAMPLES:
4673	sage: acosh(1/2)
4674	acosh(1/2)
4675	sage: acosh(1 + I*1.0)
4676	0.9045568943023813*I + 1.061275061905036
4677
4678	Warning: If the input is real the output will be real or NaN:
4679	sage: acosh(0.5)
4680	NaN
4681
4682	But evaluate where the input is in the complex field yields a complex output:
4683	sage: acosh(CC(0.5))
4684	1.04719755119660*I
4685
4686	"""
4687	def _repr_(self, simplify=True):
4688	return "acosh"
4689
4690	def _latex_(self):
4691	return "\\cosh^{-1}"
4692
4693	def _approx_(self, x):
4694	return float(pari(float(x)).acosh())
4695
4696	acosh = Function_acosh()
4697	_syms['acosh'] = acosh
4698
4699	class Function_atanh(PrimitiveFunction):
4700	"""
4701	The inverse of the hyperbolic tane function.
4702
4703	EXAMPLES:
4704	sage: atanh(0.5)
4705	0.549306144334055
4706	sage: atanh(1/2)
4707	atanh(1/2)
4708	sage: atanh(1 + I*1.0)
4709	1.017221967897851*I + 0.4023594781085251
4710	"""
4711	def _repr_(self, simplify=True):
4712	return "atanh"
4713
4714	def _latex_(self):
4715	return "\\tanh^{-1}"
4716
4717	def _approx_(self, x):
4718	return float(pari(float(x)).atanh())
4719
4720	atanh = Function_atanh()
4721	_syms['atanh'] = atanh
4722
4723	class Function_acos(PrimitiveFunction):
4724	"""
4725	The arccosine function
4726
4727	EXAMPLES:
4728	sage: acos(0.5)
4729	1.04719755119660
4730	sage: acos(1/2)
4731	pi/3
4732	sage: acos(1 + I*1.0)
4733	0.9045568943023813 - 1.061275061905036*I
4734	"""
4735	def _repr_(self, simplify=True):
4736	return "acos"
4737
4738	def _latex_(self):
4739	return "\\cos^{-1}"
4740
4741	def _approx_(self, x):
4742	return math.acos(x)
4743
4744	acos = Function_acos()
4745	_syms['acos'] = acos
4746
4747
4748	class Function_atan(PrimitiveFunction):
4749	"""
4750	The arctangent function.
4751
4752	EXAMPLES:
4753	sage: atan(1/2)
4754	atan(1/2)
4755	sage: RDF(atan(1/2))
4756	0.463647609001
4757	sage: atan(1 + I)
4758	atan(I + 1)
4759	"""
4760	def _repr_(self, simplify=True):
4761	return "atan"
4762
4763	def _latex_(self):
4764	return "\\tan^{-1}"
4765
4766	def _approx_(self, x):
4767	return math.atan(x)
4768
4769	atan = Function_atan()
4770	_syms['atan'] = atan
4771
4772
4773	#######
4774	# Hyperbolic functions
4775	#######
4776
4777	#tanh
4778	class Function_tanh(PrimitiveFunction):
4779	"""
4780	The hyperbolic tangent function.
4781
4782	EXAMPLES:
4783	sage: tanh(pi)
4784	tanh(pi)
4785	sage: tanh(3.1415)
4786	0.996271386633702
4787	sage: float(tanh(pi)) # random low-order bits
4788	0.99627207622074987
4789	sage: tan(3.1415/4)
4790	0.999953674278156
4791	sage: tanh(pi/4)
4792	tanh(pi/4)
4793	sage: RR(tanh(1/2))
4794	0.462117157260010
4795
4796	sage: CC(tanh(pi + I*e))
4797	0.997524731976164 - 0.00279068768100315*I
4798	sage: ComplexField(100)(tanh(pi + I*e))
4799	0.99752473197616361034204366446 - 0.0027906876810031453884245163923*I
4800	sage: CDF(tanh(pi + I*e))
4801	0.997524731976 - 0.002790687681*I
4802	"""
4803	def _repr_(self, simplify=True):
4804	return "tanh"
4805
4806	def _latex_(self):
4807	return "\\tanh"
4808
4809	def _approx_(self, x):
4810	return math.tanh(x)
4811
4812	tanh = Function_tanh()
4813	_syms['tanh'] = tanh
4814
4815	#sinh
4816	class Function_sinh(PrimitiveFunction):
4817	"""
4818	The hyperbolic sine function.
4819
4820	EXAMPLES:
4821	sage: sinh(pi)
4822	sinh(pi)
4823	sage: sinh(3.1415)
4824	11.5476653707437
4825	sage: float(sinh(pi)) # random low-order bits
4826	11.548739357257748
4827	sage: RR(sinh(pi))
4828	11.5487393572577
4829	"""
4830	def _repr_(self, simplify=True):
4831	return "sinh"
4832
4833	def _latex_(self):
4834	return "\\sinh"
4835
4836	def _approx_(self, x):
4837	return math.sinh(x)
4838
4839	sinh = Function_sinh()
4840	_syms['sinh'] = sinh
4841
4842	#cosh
4843	class Function_cosh(PrimitiveFunction):
4844	"""
4845	The hyperbolic cosine function.
4846
4847	EXAMPLES:
4848	sage: cosh(pi)
4849	cosh(pi)
4850	sage: cosh(3.1415)
4851	11.5908832931176
4852	sage: float(cosh(pi)) # random low order bits
4853	11.591953275521519
4854	sage: RR(cosh(1/2))
4855	1.12762596520638
4856	"""
4857	def _repr_(self, simplify=True):
4858	return "cosh"
4859
4860	def _latex_(self):
4861	return "\\cosh"
4862
4863	def _approx_(self, x):
4864	return math.cosh(x)
4865
4866	cosh = Function_cosh()
4867	_syms['cosh'] = cosh
4868
4869	#coth
4870	class Function_coth(PrimitiveFunction):
4871	"""
4872	The hyperbolic cotangent function.
4873
4874	EXAMPLES:
4875	sage: coth(pi)
4876	coth(pi)
4877	sage: coth(3.1415)
4878	1.00374256795520
4879	sage: float(coth(pi))
4880	1.0037418731973213
4881	sage: RR(coth(pi))
4882	1.00374187319732
4883	"""
4884	def _repr_(self, simplify=True):
4885	return "coth"
4886
4887	def _latex_(self):
4888	return "\\coth"
4889
4890	def _approx_(self, x):
4891	return 1/math.tanh(x)
4892
4893	coth = Function_coth()
4894	_syms['coth'] = coth
4895
4896	#sech
4897	class Function_sech(PrimitiveFunction):
4898	"""
4899	The hyperbolic secant function.
4900
4901	EXAMPLES:
4902	sage: sech(pi)
4903	sech(pi)
4904	sage: sech(3.1415)
4905	0.0862747018248192
4906	sage: float(sech(pi)) # random low order bits
4907	0.086266738334054432
4908	sage: RR(sech(pi))
4909	0.0862667383340544
4910	"""
4911	def _repr_(self, simplify=True):
4912	return "sech"
4913
4914	def _latex_(self):
4915	return "\\sech"
4916
4917	def _approx_(self, x):
4918	return 1/math.cosh(x)
4919
4920	sech = Function_sech()
4921	_syms['sech'] = sech
4922
4923
4924	#csch
4925	class Function_csch(PrimitiveFunction):
4926	"""
4927	The hyperbolic cosecant function.
4928
4929	EXAMPLES:
4930	sage: csch(pi)
4931	csch(pi)
4932	sage: csch(3.1415)
4933	0.0865975907592133
4934	sage: float(csch(pi)) # random low-order bits
4935	0.086589537530046945
4936	sage: RR(csch(pi))
4937	0.0865895375300470
4938	"""
4939	def _repr_(self, simplify=True):
4940	return "csch"
4941
4942	def _latex_(self):
4943	return "\\csch"
4944
4945	def _approx_(self, x):
4946	return 1/math.sinh(x)
4947
4948	csch = Function_csch()
4949	_syms['csch'] = csch
4950
4951	#############
4952	# log
4953	#############
4954
4955	class Function_log(PrimitiveFunction):
4956	"""
4957	The log function.
4958
4959	EXAMPLES:
4960	sage: log(e^2)
4961	2
4962	sage: log(1024, 2) # the following is ugly (for now)
4963	log(1024)/log(2)
4964	sage: log(10, 4)
4965	log(10)/log(4)
4966
4967	sage: RDF(log(10,2))
4968	3.32192809489
4969	sage: RDF(log(8, 2))
4970	3.0
4971	sage: log(RDF(10))
4972	2.30258509299
4973	sage: log(2.718)
4974	0.999896315728952
4975	"""
4976	def __init__(self):
4977	PrimitiveFunction.__init__(self)
4978
4979	def _repr_(self, simplify=True):
4980	return "log"
4981
4982	def _latex_(self):
4983	return "\\log"
4984
4985	_approx_ = math.log
4986
4987	function_log = Function_log()
4988	ln = function_log
4989
4990	def log(x, base=None):
4991	if base is None:
4992	try:
4993	return x.log()
4994	except AttributeError:
4995	return function_log(x)
4996	else:
4997	try:
4998	return x.log(base)
4999	except AttributeError:
5000	return function_log(x) / function_log(base)
5001
5002	#####################
5003	# The polylogarithm
5004	#####################
5005
5006
5007	class Function_polylog(PrimitiveFunction):
5008	r"""
5009	The polylog function $\text{Li}_n(z) = \sum_{k=1}^{\infty} z^k / k^n$.
5010
5011	INPUT:
5012	n -- object
5013	z -- object
5014
5015	EXAMPLES:
5016
5017	"""
5018	def __init__(self, n):
5019	PrimitiveFunction.__init__(self)
5020	self._n = n
5021
5022	def _repr_(self, simplify=True):
5023	return "polylog(%s)"%self._n
5024
5025	def _repr_evaled_(self, args):
5026	return 'polylog(%s, %s)'%(self._n, args)
5027
5028	def _maxima_init_(self):
5029	if self._n in [1,2,3]:
5030	return 'li[%s]'%self._n
5031	else:
5032	return 'polylog(%s)'%self._n
5033
5034	def _maxima_init_evaled_(self, args):
5035	if self._n in [1,2,3]:
5036	return 'li[%s](%s)'%(self._n, args)
5037	else:
5038	return 'polylog(%s, %s)'%(self._n, args)
5039
5040	def n(self):
5041	return self._n
5042
5043	def _latex_(self):
5044	return "\\text{Li}"
5045
5046	def _approx_(self, x):
5047	try:
5048	return float(pari(x).polylog(self._n))
5049	except PariError:
5050	raise TypeError, 'unable to coerce polylogarithm to float'
5051
5052
5053	def polylog(n, z):
5054	r"""
5055	The polylog function $\text{Li}_n(z) = \sum_{k=1}^{\infty} z^k / k^n$.
5056
5057	EXAMPLES:
5058	sage: polylog(2,1)
5059	pi^2/6
5060	sage: polylog(2,x^2+1)
5061	polylog(2, x^2 + 1)
5062	sage: polylog(4,0.5)
5063	polylog(4, 0.500000000000000)
5064	sage: float(polylog(4,0.5))
5065	0.51747906167389934
5066
5067	sage: var('z')
5068	z
5069	sage: polylog(2,z).taylor(z, 1/2, 3)
5070	-(6log(2)^2 - pi^2)/12 + 2log(2)(z - 1/2) + (-2log(2) + 2)(z - 1/2)^2 + (8log(2) - 4)*(z - 1/2)^3/3
5071	"""
5072	return Function_polylog(n)(z)
5073
5074	_syms['polylog2'] = Function_polylog(2)
5075	_syms['polylog3'] = Function_polylog(3)
5076
5077	##############################
5078	# square root
5079	##############################
5080
5081	class Function_sqrt(PrimitiveFunction):
5082	"""
5083	The square root function. This is a symbolic square root.
5084
5085	EXAMPLES:
5086	sage: sqrt(-1)
5087	I
5088	sage: sqrt(2)
5089	sqrt(2)
5090	sage: sqrt(x^2)
5091	abs(x)
5092	"""
5093	def __init__(self):
5094	PrimitiveFunction.__init__(self, needs_braces=True)
5095
5096	def _repr_(self, simplify=True):
5097	return "sqrt"
5098
5099	def _latex_(self):
5100	return "\\sqrt"
5101
5102	def _do_sqrt(self, x, prec=None, extend=True, all=False):
5103	if prec:
5104	return ComplexField(prec)(x).sqrt(all=all)
5105	z = SymbolicComposition(self, SR(x))
5106	if all:
5107	return [z, -z]
5108	return z
5109
5110	def __call__(self, x, args, *kwds):
5111	"""
5112
5113	POSSIBLE INPUTS INCLUDE:
5114	x -- a number
5115	prec -- integer (default: None): if None, returns an exact
5116	square root; otherwise returns a numerical square
5117	root if necessary, to the given bits of precision.
5118	extend -- bool (default: True); if True, return a square
5119	root in an extension ring, if necessary. Otherwise,
5120	raise a ValueError if the square is not in the base
5121	ring.
5122	all -- bool (default: False); if True, return all square
5123	roots of self, instead of just one.
5124	"""
5125	if isinstance(x, float):
5126	return math.sqrt(x)
5127	if not isinstance(x, (Integer, Rational)):
5128	try:
5129	return x.sqrt(args, *kwds)
5130	except AttributeError:
5131	pass
5132	return self._do_sqrt(x, args, *kwds)
5133
5134	def _approx_(self, x):
5135	return math.sqrt(x)
5136
5137	sqrt = Function_sqrt()
5138	_syms['sqrt'] = sqrt
5139
5140	class Function_exp(PrimitiveFunction):
5141	r"""
5142	The exponential function, $\exp(x) = e^x$.
5143
5144	EXAMPLES:
5145	sage: exp(-1)
5146	e^-1
5147	sage: exp(2)
5148	e^2
5149	sage: exp(x^2 + log(x))
5150	x*e^x^2
5151	sage: exp(2.5)
5152	12.1824939607035
5153	sage: exp(float(2.5)) # random low order bits
5154	12.182493960703473
5155	sage: exp(RDF('2.5'))
5156	12.1824939607
5157	"""
5158	def __init__(self):
5159	PrimitiveFunction.__init__(self, needs_braces=True)
5160
5161	def _repr_(self, simplify=True):
5162	return "exp"
5163
5164	def _latex_(self):
5165	return "\\exp"
5166
5167	def __call__(self, x):
5168	# if x is an integer or rational, never call the sqrt method
5169	if isinstance(x, float):
5170	return self._approx_(x)
5171	if not isinstance(x, (Integer, Rational)):
5172	try:
5173	return x.exp()
5174	except AttributeError:
5175	pass
5176	return SymbolicComposition(self, SR(x))
5177
5178
5179	def _approx_(self, x):
5180	return math.exp(x)
5181
5182	exp = Function_exp()
5183	_syms['exp'] = exp
5184
5185
5186	#######################################################
5187	# Symbolic functions
5188	#######################################################
5189
5190	class SymbolicFunction(PrimitiveFunction):
5191	"""
5192	A formal symbolic function.
5193
5194	EXAMPLES:
5195	sage: f = function('foo')
5196	sage: var('x,y,z')
5197	(x, y, z)
5198	sage: g = f(x,y,z)
5199	sage: g
5200	foo(x, y, z)
5201	sage: g(x=var('theta'))
5202	foo(theta, y, z)
5203	"""
5204	def __init__(self, name):
5205	PrimitiveFunction.__init__(self, needs_braces=True)
5206	self._name = str(name)
5207
5208	def __hash__(self):
5209	return hash(self._name)
5210
5211	def _repr_(self, simplify=True):
5212	return self._name
5213
5214
5215	def _is_atomic(self):
5216	return True
5217
5218	def _latex_(self):
5219	return "{\\rm %s}"%self._name
5220
5221	def _maxima_init_(self):
5222	return "'%s"%self._name
5223
5224	def _approx_(self, x):
5225	raise TypeError
5226
5227	def __call__(self, args, *kwds):
5228	return SymbolicFunctionEvaluation(self, [SR(x) for x in args])
5229
5230
5231
5232	class SymbolicFunction_delayed(SymbolicFunction):
5233	def simplify(self):
5234	return self
5235
5236	def is_simplified(self):
5237	return True
5238
5239	def _maxima_init_(self):
5240	return "%s"%self._name
5241
5242	def __call__(self, *args):
5243	return SymbolicFunctionEvaluation_delayed(self, [SR(x) for x in args])
5244
5245
5246
5247	class SymbolicFunctionEvaluation(SymbolicExpression):
5248	"""
5249	The result of evaluating a formal symbolic function.
5250
5251	EXAMPLES:
5252	sage: h = function('gfun', x); h
5253	gfun(x)
5254	sage: k = h.integral(x); k
5255	integrate(gfun(x), x)
5256	sage: k(gfun=sin)
5257	-cos(x)
5258	sage: k(gfun=cos)
5259	sin(x)
5260	sage: k.diff(x)
5261	gfun(x)
5262	"""
5263	def __init__(self, f, args=None, kwds=None):
5264	"""
5265	INPUT:
5266	f -- symbolic function
5267	args -- a tuple or list of symbolic expressions, at which
5268	f is formally evaluated.
5269	"""
5270	SymbolicExpression.__init__(self)
5271	self._f = f
5272	if not args is None:
5273	if not isinstance(args, tuple):
5274	args = tuple(args)
5275	self._args = args
5276	self._kwds = kwds
5277
5278	def __float__(self):
5279	return float(maxima(self))
5280
5281	def _is_atomic(self):
5282	return True
5283
5284	def arguments(self):
5285	return tuple(self._args)
5286
5287	def keyword_arguments(self):
5288	return self._kwds
5289
5290	def _repr_(self, simplify=True):
5291	if simplify:
5292	return self.simplify()._repr_(simplify=False)
5293	else:
5294	args = ', '.join([x._repr_(simplify=simplify) for x in
5295	self._args])
5296	if not self._kwds is None:
5297	kwds = ', '.join(["%s=%s" %(x, y) for x,y in self._kwds.iteritems()])
5298	return '%s(%s, %s)' % (self._f._name, args, kwds)
5299	else:
5300	return '%s(%s)' % (self._f._name, args)
5301
5302	def _latex_(self):
5303	return "{\\rm %s}(%s)"%(self._f._name, ', '.join([x._latex_() for
5304	x in self._args]))
5305
5306	def _maxima_init_(self):
5307	try:
5308	return self.__maxima_init
5309	except AttributeError:
5310	n = self._f._name
5311	if not (n in ['integrate', 'diff']):
5312	n = "'" + n
5313	s = "%s(%s)"%(n, ', '.join([x._maxima_init_()
5314	for x in self._args]))
5315	self.__maxima_init = s
5316	return s
5317
5318	def _recursive_sub(self, kwds):
5319	"""
5320	EXAMPLES:
5321	sage: y = var('y')
5322	sage: f = function('foo',x); f
5323	foo(x)
5324	sage: f(foo=sin)
5325	sin(x)
5326	sage: f(x+y)
5327	foo(y + x)
5328	sage: a = f(pi)
5329	sage: a.substitute(foo = sin)
5330	0
5331	sage: a = f(pi/2)
5332	sage: a.substitute(foo = sin)
5333	1
5334
5335	sage: b = f(pi/3) + x + y
5336	sage: b
5337	y + x + foo(pi/3)
5338	sage: b(foo = sin)
5339	y + x + sqrt(3)/2
5340	sage: b(foo = cos)
5341	y + x + 1/2
5342	sage: b(foo = cos, x=y)
5343	2*y + 1/2
5344	"""
5345	function_sub = False
5346	for x in kwds.keys():
5347	if repr(x) == self._f._name:
5348	g = kwds[x]
5349	function_sub = True
5350	break
5351	if function_sub:
5352	# Very important to make a copy, since we are mutating a dictionary
5353	# that will get used again by the calling function!
5354	kwds = dict(kwds)
5355	del kwds[x]
5356
5357	arg = tuple([SR(x._recursive_sub(kwds)) for x in self._args])
5358
5359	if function_sub:
5360	return g(*arg)
5361	else:
5362	return self.__class__(self._f, arg)
5363
5364	def _recursive_sub_over_ring(self, kwds, ring):
5365	raise TypeError, "no way to coerce the formal function to ring."
5366
5367	def variables(self):
5368	"""
5369	Return the variables appearing in the simplified form of self.
5370
5371	EXAMPLES:
5372	sage: foo = function('foo')
5373	sage: var('x,y,a,b,z,t')
5374	(x, y, a, b, z, t)
5375	sage: w = foo(x,y,a,b,z) + t
5376	sage: w
5377	foo(x, y, a, b, z) + t
5378	sage: w.variables()
5379	(a, b, t, x, y, z)
5380	"""
5381	try:
5382	return self.__variables
5383	except AttributeError:
5384	pass
5385	vars = sum([list(op.variables()) for op in self._args], [])
5386	vars.sort(var_cmp)
5387	vars = tuple(vars)
5388	self.__variables = vars
5389	return vars
5390
5391	class SymbolicFunctionEvaluation_delayed(SymbolicFunctionEvaluation):
5392	def simplify(self):
5393	return self
5394
5395	def is_simplified(self):
5396	return True
5397
5398	def __float__(self):
5399	return float(self._maxima_())
5400
5401	def __complex__(self):
5402	return complex(self._maxima_())
5403
5404	def _real_double_(self, R):
5405	return R(float(self))
5406
5407	def _real_rqdf_(self, R):
5408	raise TypeError
5409
5410	def _complex_double_(self, C):
5411	return C(float(self))
5412
5413	def _mpfr_(self, field):
5414	if field.prec() <= 53:
5415	return field(float(self))
5416	raise TypeError
5417
5418	def _complex_mpfr_field_(self, field):
5419	if field.prec() <= 53:
5420	return field(float(self))
5421	raise TypeError
5422
5423	def _maxima_init_(self):
5424	try:
5425	return self.__maxima_init
5426	except AttributeError:
5427	n = self._f._name
5428	s = "%s(%s)"%(n, ', '.join([x._maxima_init_()
5429	for x in self._args]))
5430	self.__maxima_init = s
5431	return s
5432
5433
5434
5435	_functions = {}
5436	def function(s, *args):
5437	"""
5438	Create a formal symbolic function with the name \emph{s}.
5439
5440	EXAMPLES:
5441	sage: var('a, b')
5442	(a, b)
5443	sage: f = function('cr', a)
5444	sage: g = f.diff(a).integral(b)
5445	sage: g
5446	diff(cr(a), a, 1)*b
5447	sage: g(cr=cos)
5448	-sin(a)*b
5449	sage: g(cr=sin(x) + cos(x))
5450	(cos(a) - sin(a))*b
5451
5452	Basic arithmetic:
5453	sage: x = var('x')
5454	sage: h = function('f',x)
5455	sage: 2*f
5456	2*f
5457	sage: 2*h
5458	2*f(x)
5459	"""
5460	if len(args) > 0:
5461	return function(s)(*args)
5462	if isinstance(s, SymbolicFunction):
5463	return s
5464	s = str(s)
5465	if ',' in s:
5466	return tuple([function(x.strip()) for x in s.split(',')])
5467	elif ' ' in s:
5468	return tuple([function(x.strip()) for x in s.split()])
5469	try:
5470	f = _functions[s]
5471	_syms[s] = f
5472	return f
5473	except KeyError:
5474	pass
5475	v = SymbolicFunction(s)
5476	_functions[s] = v
5477	_syms[s] = v
5478	return v
5479
5480	#######################################################
5481	# This is buggy as is:
5482	# sage: a = lim(exp(x^2)*(1-erf(x)), x=infinity)
5483	# sage: maxima(a)
5484	# 'limit(%e^x^2-%e^x^2*erf(x))
5485	# ??? -- what happened to the infinity in the above.
5486	# With the old version one gets
5487	# sage: maxima(a)
5488	# 'limit(%e^x^2-%e^x^2*erf(x),x,inf)
5489	# Which is right, since:
5490	# (%i3) limit(%e^x^2-%e^x^2*erf(x),x,inf);
5491	# (%o3) 'limit(%e^x^2-%e^x^2*erf(x),x,inf)
5492	#a
5493	def dummy_limit(*args):
5494	s = str(args[1])
5495	return SymbolicFunctionEvaluation(function('limit'), args=(args[0],), kwds={s: args[2]})
5496
5497	######################################i################
5498
5499
5500
5501
5502	#######################################################
5503
5504	symtable = {'%pi':'pi', '%e': 'e', '%i':'I', '%gamma':'euler_gamma',
5505	'li[2]':'polylog2', 'li[3]':'polylog3'}
5506
5507	from sage.rings.infinity import infinity, minus_infinity
5508
5509	_syms['inf'] = infinity
5510	_syms['minf'] = minus_infinity
5511
5512	from sage.misc.multireplace import multiple_replace
5513
5514	maxima_tick = re.compile("'[a-z\|A-Z\|0-9\|_]*")
5515
5516	maxima_qp = re.compile("\?\%[a-z\|A-Z\|0-9\|_]*") # e.g., ?%jacobi_cd
5517
5518	maxima_var = re.compile("\%[a-z\|A-Z\|0-9\|_]*") # e.g., ?%jacobi_cd
5519
5520	sci_not = re.compile("(-?(?:0\|[1-9]\d*))(\.\d+)?([eE][-+]\d+)")
5521
5522	def symbolic_expression_from_maxima_string(x, equals_sub=False, maxima=maxima):
5523	syms = dict(_syms)
5524
5525	if len(x) == 0:
5526	raise RuntimeError, "invalid symbolic expression -- ''"
5527	maxima.set('_tmp_',x)
5528
5529	# This is inefficient since it so rarely is needed:
5530	#r = maxima._eval_line('listofvars(_tmp_);')[1:-1]
5531
5532	s = maxima._eval_line('_tmp_;')
5533
5534	formal_functions = maxima_tick.findall(s)
5535	if len(formal_functions) > 0:
5536	for X in formal_functions:
5537	syms[X[1:]] = function(X[1:])
5538	# You might think there is a potential very subtle bug if 'foo is in a string literal --
5539	# but string literals should never ever be part of a symbolic expression.
5540	s = s.replace("'","")
5541
5542	delayed_functions = maxima_qp.findall(s)
5543	if len(delayed_functions) > 0:
5544	for X in delayed_functions:
5545	syms[X[2:]] = SymbolicFunction_delayed(X[2:])
5546	s = s.replace("?%","")
5547
5548	s = multiple_replace(symtable, s)
5549	s = s.replace("%","")
5550
5551	if equals_sub:
5552	s = s.replace('=','==')
5553
5554	#replace all instances of scientific notation
5555	#with regular notation
5556	search = sci_not.search(s)
5557	while not search is None:
5558	(start, end) = search.span()
5559	s = s.replace(s[start:end], str(RR(s[start:end])))
5560	search = sci_not.search(s)
5561
5562	# have to do this here, otherwise maxima_tick catches it
5563	syms['limit'] = dummy_limit
5564
5565	global is_simplified
5566	try:
5567	# use a global flag so all expressions obtained via
5568	# evaluation of maxima code are assumed pre-simplified
5569	is_simplified = True
5570	last_msg = ''
5571	while True:
5572	try:
5573	w = sage_eval(s, syms)
5574	except NameError, msg:
5575	if msg == last_msg:
5576	raise NameError, msg
5577	msg = str(msg)
5578	last_msg = msg
5579	i = msg.find("'")
5580	j = msg.rfind("'")
5581	nm = msg[i+1:j]
5582
5583	res = re.match('.' + nm + '\s\(.*\)', s)
5584	if res:
5585	syms[nm] = function(nm)
5586	else:
5587	syms[nm] = var(nm)
5588	else:
5589	break
5590	if isinstance(w, (list, tuple)):
5591	return w
5592	else:
5593	x = SR(w)
5594	return x
5595	except SyntaxError:
5596	raise TypeError, "unable to make sense of Maxima expression '%s' in SAGE"%s
5597	finally:
5598	is_simplified = False
5599
5600	def symbolic_expression_from_maxima_element(x):
5601	return symbolic_expression_from_maxima_string(x.name())
5602
5603	def evaled_symbolic_expression_from_maxima_string(x):
5604	return symbolic_expression_from_maxima_string(maxima.eval(x))
5605
5606
5607	# External access used by restore
5608	syms_cur = _syms
5609	syms_default = dict(syms_cur)

Note: See TracBrowser for help on using the repository browser.

Download in other formats: