# source:sage/calculus/calculus.py@6377:cb38f5253f8b

Revision 6377:cb38f5253f8b, 148.0 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

work in progress on algebraic number theory.

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