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

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

Fixed #847.

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