Ticket #14801: trac_14801_piecewiese.patch

File trac_14801_piecewiese.patch, 141.0 KB (added by vbraun, 6 years ago)

Updated patch

  • sage/functions/all.py

    # HG changeset patch
    # User Volker Braun <vbraun@stp.dias.ie>
    # Date 1372346122 14400
    #      Thu Jun 27 11:15:22 2013 -0400
    # Node ID 530f6f96f708d3565c153752a7c6d04ed07c8bd0
    # Parent  3f607223351f015a7fef512ba512497ed782e0b9
    Piecewise functions using the symbolic expression framework
    
    diff --git a/sage/functions/all.py b/sage/functions/all.py
    a b  
    1 from piecewise import piecewise, Piecewise
     1from sage.misc.lazy_import import lazy_import
     2
     3lazy_import('sage.functions.piecewise', 'Piecewise')   # deprecated
     4lazy_import('sage.functions.piecewise', 'piecewise')
    25
    36from trig import ( sin, cos, sec, csc, cot, tan,
    47                   asin, acos, atan,
  • new file sage/functions/piecewise-old.py

    diff --git a/sage/functions/piecewise-old.py b/sage/functions/piecewise-old.py
    new file mode 100644
    - +  
     1r"""
     2Piecewise-defined Functions
     3
     4Sage implements a very simple class of piecewise-defined functions.
     5Functions may be any type of symbolic expression. Infinite
     6intervals are not supported. The endpoints of each interval must
     7line up.
     8
     9TODO:
     10
     11- Implement max/min location and values,
     12
     13- Need: parent object - ring of piecewise functions
     14
     15- This class should derive from an element-type class, and should
     16  define ``_add_``, ``_mul_``, etc. That will automatically take care
     17  of left multiplication and proper coercion. The coercion mentioned
     18  below for scalar mult on right is bad, since it only allows ints and
     19  rationals. The right way is to use an element class and only define
     20  ``_mul_``, and have a parent, so anything gets coerced properly.
     21
     22AUTHORS:
     23
     24- David Joyner (2006-04): initial version
     25
     26- David Joyner (2006-09): added __eq__, extend_by_zero_to, unextend,
     27  convolution, trapezoid, trapezoid_integral_approximation,
     28  riemann_sum, riemann_sum_integral_approximation, tangent_line fixed
     29  bugs in __mul__, __add__
     30
     31- David Joyner (2007-03): adding Hann filter for FS, added general FS
     32  filter methods for computing and plotting, added options to plotting
     33  of FS (eg, specifying rgb values are now allowed). Fixed bug in
     34  documentation reported by Pablo De Napoli.
     35
     36- David Joyner (2007-09): bug fixes due to behaviour of
     37  SymbolicArithmetic
     38
     39- David Joyner (2008-04): fixed docstring bugs reported by J Morrow; added
     40  support for Laplace transform of functions with infinite support.
     41
     42- David Joyner (2008-07): fixed a left multiplication bug reported by
     43  C. Boncelet (by defining __rmul__ = __mul__).
     44
     45- Paul Butler (2009-01): added indefinite integration and default_variable
     46
     47TESTS::
     48
     49    sage: R.<x> = QQ[]
     50    sage: f = Piecewise([[(0,1),1*x^0]])
     51    sage: 2*f
     52    Piecewise defined function with 1 parts, [[(0, 1), 2]]
     53"""
     54
     55#*****************************************************************************
     56#       Copyright (C) 2006 William Stein <wstein@gmail.com>
     57#                     2006 David Joyner <wdjoyner@gmail.com>
     58#
     59#  Distributed under the terms of the GNU General Public License (GPL)
     60#
     61#    This code is distributed in the hope that it will be useful,
     62#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     63#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     64#    General Public License for more details.
     65#
     66#  The full text of the GPL is available at:
     67#
     68#                  http://www.gnu.org/licenses/
     69#*****************************************************************************
     70
     71from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     72from sage.misc.sage_eval import sage_eval
     73from sage.rings.all import QQ, RR, Integer, Rational, infinity
     74from sage.calculus.functional import derivative
     75from sage.symbolic.expression import is_Expression
     76from sage.symbolic.assumptions import assume, forget
     77
     78from sage.calculus.calculus import SR, maxima
     79from sage.calculus.all import var
     80
     81def piecewise(list_of_pairs, var=None):
     82    """
     83    Returns a piecewise function from a list of (interval, function)
     84    pairs.
     85   
     86    ``list_of_pairs`` is a list of pairs (I, fcn), where
     87    fcn is a Sage function (such as a polynomial over RR, or functions
     88    using the lambda notation), and I is an interval such as I = (1,3).
     89    Two consecutive intervals must share a common endpoint.
     90   
     91    If the optional ``var`` is specified, then any symbolic expressions
     92    in the list will be converted to symbolic functions using
     93    ``fcn.function(var)``.  (This says which variable is considered to
     94    be "piecewise".)
     95
     96    We assume that these definitions are consistent (ie, no checking is
     97    done).
     98   
     99    EXAMPLES::
     100   
     101        sage: f1(x) = -1
     102        sage: f2(x) = 2
     103        sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
     104        sage: f(1)
     105        -1
     106        sage: f(3)
     107        2
     108        sage: f = Piecewise([[(0,1),x], [(1,2),x^2]], x); f
     109        Piecewise defined function with 2 parts, [[(0, 1), x |--> x], [(1, 2), x |--> x^2]]
     110        sage: f(0.9)
     111        0.900000000000000
     112        sage: f(1.1)
     113        1.21000000000000
     114    """
     115    return PiecewisePolynomial(list_of_pairs, var=var)
     116
     117Piecewise = piecewise
     118
     119class PiecewisePolynomial:
     120    """
     121    Returns a piecewise function from a list of (interval, function)
     122    pairs.
     123   
     124    EXAMPLES::
     125   
     126        sage: f1(x) = -1
     127        sage: f2(x) = 2
     128        sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
     129        sage: f(1)
     130        -1
     131        sage: f(3)
     132        2
     133    """
     134    def __init__(self, list_of_pairs, var=None):
     135        r"""
     136        ``list_of_pairs`` is a list of pairs (I, fcn), where
     137        fcn is a Sage function (such as a polynomial over RR, or functions
     138        using the lambda notation), and I is an interval such as I = (1,3).
     139        Two consecutive intervals must share a common endpoint.
     140       
     141        If the optional ``var`` is specified, then any symbolic
     142        expressions in the list will be converted to symbolic
     143        functions using ``fcn.function(var)``.  (This says which
     144        variable is considered to be "piecewise".)
     145
     146        We assume that these definitions are consistent (ie, no checking is
     147        done).
     148
     149        EXAMPLES::
     150
     151            sage: f1(x) = 1
     152            sage: f2(x) = 1 - x
     153            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     154            sage: f.list()
     155            [[(0, 1), x |--> 1], [(1, 2), x |--> -x + 1]]
     156            sage: f.length()
     157            2
     158        """
     159        self._length = len(list_of_pairs)
     160        self._intervals = [x[0] for x in list_of_pairs]
     161        functions = [x[1] for x in list_of_pairs]
     162        if var is not None:
     163            for i in range(len(functions)):
     164                if is_Expression(functions[i]):
     165                    functions[i] = functions[i].function(var)
     166        self._functions = functions
     167        # We regenerate self._list in case self._functions was modified
     168        # above.  This also protects us in case somebody mutates a list
     169        # after they use it as an argument to piecewise().
     170        self._list = [[self._intervals[i], self._functions[i]] for i in range(self._length)]
     171 
     172    def list(self):
     173        """
     174        Returns the pieces of this function as a list of functions.
     175       
     176        EXAMPLES::
     177
     178            sage: f1(x) = 1
     179            sage: f2(x) = 1 - x
     180            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     181            sage: f.list()
     182            [[(0, 1), x |--> 1], [(1, 2), x |--> -x + 1]]
     183        """
     184        return self._list
     185 
     186    def length(self):
     187        """
     188        Returns the number of pieces of this function.
     189       
     190        EXAMPLES::
     191       
     192            sage: f1(x) = 1
     193            sage: f2(x) = 1 - x
     194            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     195            sage: f.length()
     196            2
     197        """
     198        return self._length
     199 
     200    def __repr__(self):
     201        """
     202        EXAMPLES::
     203       
     204            sage: f1(x) = 1
     205            sage: f2(x) = 1 - x
     206            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]); f
     207            Piecewise defined function with 2 parts, [[(0, 1), x |--> 1], [(1, 2), x |--> -x + 1]]
     208        """
     209        return 'Piecewise defined function with %s parts, %s'%(
     210            self.length(),self.list())
     211 
     212    def _latex_(self):
     213        r"""
     214        EXAMPLES::
     215       
     216            sage: f1(x) = 1
     217            sage: f2(x) = 1 - x
     218            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     219            sage: latex(f)
     220            \begin{cases}
     221            x \ {\mapsto}\ 1 &\text{on $(0, 1)$}\cr
     222            x \ {\mapsto}\ -x + 1 &\text{on $(1, 2)$}\cr
     223            \end{cases}
     224       
     225        ::
     226       
     227            sage: f(x) = sin(x*pi/2)
     228            sage: g(x) = 1-(x-1)^2
     229            sage: h(x) = -x
     230            sage: P = Piecewise([[(0,1), f], [(1,3),g], [(3,5), h]])
     231            sage: latex(P)
     232            \begin{cases}
     233            x \ {\mapsto}\ \sin\left(\frac{1}{2} \, \pi x\right) &\text{on $(0, 1)$}\cr
     234            x \ {\mapsto}\ -{\left(x - 1\right)}^{2} + 1 &\text{on $(1, 3)$}\cr
     235            x \ {\mapsto}\ -x &\text{on $(3, 5)$}\cr
     236            \end{cases}
     237        """
     238        from sage.misc.latex import latex
     239        tex = ['\\begin{cases}\n']
     240        for (left, right), f in self.list():
     241            tex.append('%s &\\text{on $(%s, %s)$}\\cr\n' % (latex(f), left, right))
     242        tex.append(r'\end{cases}')
     243        return ''.join(tex)
     244
     245    def intervals(self):
     246        """
     247        A piecewise non-polynomial example.
     248       
     249        EXAMPLES::
     250       
     251            sage: f1(x) = 1
     252            sage: f2(x) = 1-x
     253            sage: f3(x) = exp(x)
     254            sage: f4(x) = sin(2*x)
     255            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     256            sage: f.intervals()
     257            [(0, 1), (1, 2), (2, 3), (3, 10)]
     258        """
     259        return self._intervals
     260 
     261    def domain(self):
     262        """
     263        Returns the domain of the function.
     264       
     265        EXAMPLES::
     266       
     267            sage: f1(x) = 1
     268            sage: f2(x) = 1-x
     269            sage: f3(x) = exp(x)
     270            sage: f4(x) = sin(2*x)
     271            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     272            sage: f.domain()
     273            (0, 10)
     274        """
     275        endpoints = sum(self.intervals(), ())
     276        return (min(endpoints), max(endpoints))
     277
     278    def functions(self):
     279        """
     280        Returns the list of functions (the "pieces").
     281       
     282        EXAMPLES::
     283       
     284            sage: f1(x) = 1
     285            sage: f2(x) = 1-x
     286            sage: f3(x) = exp(x)
     287            sage: f4(x) = sin(2*x)
     288            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     289            sage: f.functions()
     290            [x |--> 1, x |--> -x + 1, x |--> e^x, x |--> sin(2*x)]
     291        """
     292        return self._functions
     293       
     294    def extend_by_zero_to(self,xmin=-1000,xmax=1000):
     295        """
     296        This function simply returns the piecewise defined function which
     297        is extended by 0 so it is defined on all of (xmin,xmax). This is
     298        needed to add two piecewise functions in a reasonable way.
     299       
     300        EXAMPLES::
     301       
     302            sage: f1(x) = 1
     303            sage: f2(x) = 1 - x
     304            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     305            sage: f.extend_by_zero_to(-1, 3)
     306            Piecewise defined function with 4 parts, [[(-1, 0), 0], [(0, 1), x |--> 1], [(1, 2), x |--> -x + 1], [(2, 3), 0]]
     307        """
     308        zero = QQ['x'](0)
     309        list_of_pairs = self.list()
     310        a, b = self.domain()
     311        if xmin < a:
     312            list_of_pairs = [[(xmin, a), zero]] + list_of_pairs
     313        if xmax > b:
     314            list_of_pairs = list_of_pairs + [[(b, xmax), zero]]
     315        return Piecewise(list_of_pairs)
     316
     317    def unextend(self):
     318        """
     319        This removes any parts in the front or back of the function which
     320        is zero (the inverse to extend_by_zero_to).
     321       
     322        EXAMPLES::
     323       
     324            sage: R.<x> = QQ[]
     325            sage: f = Piecewise([[(-3,-1),1+2+x],[(-1,1),1-x^2]])
     326            sage: e = f.extend_by_zero_to(-10,10); e
     327            Piecewise defined function with 4 parts, [[(-10, -3), 0], [(-3, -1), x + 3], [(-1, 1), -x^2 + 1], [(1, 10), 0]]
     328            sage: d = e.unextend(); d
     329            Piecewise defined function with 2 parts, [[(-3, -1), x + 3], [(-1, 1), -x^2 + 1]]
     330            sage: d==f
     331            True
     332        """
     333        list_of_pairs = self.list()
     334        funcs = self.functions()
     335        if funcs[0] == 0:
     336            list_of_pairs = list_of_pairs[1:]
     337        if funcs[-1] == 0:
     338            list_of_pairs = list_of_pairs[:-1]
     339        return Piecewise(list_of_pairs)
     340
     341    def _riemann_sum_helper(self, N, func, initial=0):
     342        """
     343        A helper function for computing Riemann sums.
     344       
     345        INPUT:
     346       
     347       
     348        -  ``N`` - the number of subdivisions
     349       
     350        -  ``func`` - a function to apply to the endpoints of
     351           each subdivision
     352       
     353        -  ``initial`` - the starting value
     354       
     355       
     356        EXAMPLES::
     357       
     358            sage: f1(x) = x^2                   ## example 1
     359            sage: f2(x) = 5-x^2
     360            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     361            sage: f._riemann_sum_helper(6, lambda x0, x1: (x1-x0)*f(x1))
     362            19/6
     363        """
     364        a,b = self.domain()
     365        rsum = initial
     366        h = (b-a)/N
     367        for i in range(N):
     368            x0 = a+i*h
     369            x1 = a+(i+1)*h
     370            rsum += func(x0, x1)
     371        return rsum
     372
     373    def riemann_sum_integral_approximation(self,N,mode=None):
     374        """
     375        Returns the piecewise line function defined by the Riemann sums in
     376        numerical integration based on a subdivision into N subintervals.
     377       
     378        Set mode="midpoint" for the height of the rectangles to be
     379        determined by the midpoint of the subinterval; set mode="right" for
     380        the height of the rectangles to be determined by the right-hand
     381        endpoint of the subinterval; the default is mode="left" (the height
     382        of the rectangles to be determined by the left-hand endpoint of
     383        the subinterval).
     384       
     385        EXAMPLES::
     386       
     387            sage: f1(x) = x^2                   ## example 1
     388            sage: f2(x) = 5-x^2
     389            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     390            sage: f.riemann_sum_integral_approximation(6)
     391            17/6
     392            sage: f.riemann_sum_integral_approximation(6,mode="right")
     393            19/6
     394            sage: f.riemann_sum_integral_approximation(6,mode="midpoint")
     395            3
     396            sage: f.integral(definite=True)
     397            3
     398        """
     399        if mode is None:
     400            return self._riemann_sum_helper(N, lambda x0, x1: (x1-x0)*self(x0))
     401        elif mode == "right":
     402            return self._riemann_sum_helper(N, lambda x0, x1: (x1-x0)*self(x1))
     403        elif mode == "midpoint":
     404            return self._riemann_sum_helper(N, lambda x0, x1: (x1-x0)*self((x0+x1)/2))
     405        else:
     406            raise ValueError, "invalid mode"
     407
     408    def riemann_sum(self,N,mode=None):
     409        """
     410        Returns the piecewise line function defined by the Riemann sums in
     411        numerical integration based on a subdivision into N subintervals.
     412        Set mode="midpoint" for the height of the rectangles to be
     413        determined by the midpoint of the subinterval; set mode="right" for
     414        the height of the rectangles to be determined by the right-hand
     415        endpoint of the subinterval; the default is mode="left" (the height
     416        of the rectangles to be determined by the left-hand endpoint of
     417        the subinterval).
     418       
     419        EXAMPLES::
     420       
     421            sage: f1(x) = x^2
     422            sage: f2(x) = 5-x^2
     423            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     424            sage: f.riemann_sum(6,mode="midpoint")
     425            Piecewise defined function with 6 parts, [[(0, 1/3), 1/36], [(1/3, 2/3), 1/4], [(2/3, 1), 25/36], [(1, 4/3), 131/36], [(4/3, 5/3), 11/4], [(5/3, 2), 59/36]]
     426       
     427        ::
     428       
     429            sage: f = Piecewise([[(-1,1),(1-x^2).function(x)]])
     430            sage: rsf = f.riemann_sum(7)
     431            sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
     432            sage: Q = rsf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
     433            sage: L = add([line([[a,0],[a,f(x=a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in rsf.list()])
     434            sage: P + Q + L
     435       
     436        ::
     437       
     438            sage: f = Piecewise([[(-1,1),(1/2+x-x^3)]], x) ## example 3
     439            sage: rsf = f.riemann_sum(8)
     440            sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
     441            sage: Q = rsf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
     442            sage: L = add([line([[a,0],[a,f(x=a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in rsf.list()])
     443            sage: P + Q + L
     444        """
     445        if mode is None:
     446            rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self(x0))]],
     447                                            initial=[])
     448        elif mode == "right":
     449            rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self(x1))]],
     450                                            initial=[])
     451        elif mode == "midpoint":
     452            rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self((x0+x1)/2))]],
     453                                            initial=[])
     454        else:
     455            raise ValueError, "invalid mode"
     456        return Piecewise(rsum)
     457
     458    def trapezoid(self,N):
     459        """
     460        Returns the piecewise line function defined by the trapezoid rule
     461        for numerical integration based on a subdivision into N
     462        subintervals.
     463       
     464        EXAMPLES::
     465       
     466            sage: R.<x> = QQ[]
     467            sage: f1 = x^2
     468            sage: f2 = 5-x^2
     469            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     470            sage: f.trapezoid(4)
     471            Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*x], [(1/2, 1), 9/2*x - 2], [(1, 3/2), 1/2*x + 2], [(3/2, 2), -7/2*x + 8]]
     472       
     473        ::
     474       
     475            sage: R.<x> = QQ[]
     476            sage: f = Piecewise([[(-1,1),1-x^2]])
     477            sage: tf = f.trapezoid(4)
     478            sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
     479            sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
     480            sage: L = add([line([[a,0],[a,f(a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in tf.list()])
     481            sage: P+Q+L
     482       
     483        ::
     484       
     485            sage: R.<x> = QQ[]
     486            sage: f = Piecewise([[(-1,1),1/2+x-x^3]]) ## example 3
     487            sage: tf = f.trapezoid(6)
     488            sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
     489            sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
     490            sage: L = add([line([[a,0],[a,f(a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in tf.list()])
     491            sage: P+Q+L
     492
     493        TESTS:
     494
     495        Use variables other than x (:trac:`13836`)::
     496
     497            sage: R.<y> = QQ[]
     498            sage: f1 = y^2
     499            sage: f2 = 5-y^2
     500            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     501            sage: f.trapezoid(4)
     502            Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*y], [(1/2, 1), 9/2*y - 2], [(1, 3/2), 1/2*y + 2], [(3/2, 2), -7/2*y + 8]]
     503
     504        """
     505        x = QQ[self.default_variable()].gen()
     506        def f(x0, x1):
     507            f0, f1 = self(x0), self(x1)
     508            return [[(x0,x1),f0+(f1-f0)*(x1-x0)**(-1)*(x-x0)]]
     509        rsum = self._riemann_sum_helper(N, f, initial=[])
     510        return Piecewise(rsum)
     511
     512    def trapezoid_integral_approximation(self,N):
     513        """
     514        Returns the approximation given by the trapezoid rule for numerical
     515        integration based on a subdivision into N subintervals.
     516       
     517        EXAMPLES::
     518       
     519            sage: f1(x) = x^2                      ## example 1
     520            sage: f2(x) = 1-(1-x)^2
     521            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     522            sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
     523            sage: tf = f.trapezoid(6)
     524            sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
     525            sage: ta = f.trapezoid_integral_approximation(6)
     526            sage: t = text('trapezoid approximation = %s'%ta, (1.5, 0.25))
     527            sage: a = f.integral(definite=True)
     528            sage: tt = text('area under curve = %s'%a, (1.5, -0.5))
     529            sage: P + Q + t + tt
     530       
     531        ::
     532       
     533            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])  ## example 2
     534            sage: tf = f.trapezoid(4)
     535            sage: ta = f.trapezoid_integral_approximation(4)
     536            sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
     537            sage: t = text('trapezoid approximation = %s'%ta, (1.5, 0.25))
     538            sage: a = f.integral(definite=True)
     539            sage: tt = text('area under curve = %s'%a, (1.5, -0.5))
     540            sage: P+Q+t+tt
     541        """
     542        def f(x0, x1):
     543            f0, f1 = self(x0), self(x1)
     544            return ((f1+f0)/2)*(x1-x0)
     545        return self._riemann_sum_helper(N, f)
     546
     547    def critical_points(self):
     548        """
     549        Return the critical points of this piecewise function.
     550       
     551        .. warning::
     552
     553           Uses maxima, which prints the warning to use results with
     554           caution. Only works for piecewise functions whose parts are
     555           polynomials with real critical not occurring on the
     556           interval endpoints.
     557       
     558        EXAMPLES::
     559       
     560            sage: R.<x> = QQ[]
     561            sage: f1 = x^0
     562            sage: f2 = 10*x - x^2
     563            sage: f3 = 3*x^4 - 156*x^3 + 3036*x^2 - 26208*x
     564            sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]])
     565            sage: expected = [5, 12, 13, 14]
     566            sage: all(abs(e-a) < 0.001 for e,a in zip(expected, f.critical_points()))
     567            True
     568
     569        TESTS:
     570
     571        Use variables other than x (:trac:`13836`)::
     572
     573            sage: R.<y> = QQ[]
     574            sage: f1 = y^0
     575            sage: f2 = 10*y - y^2
     576            sage: f3 = 3*y^4 - 156*y^3 + 3036*y^2 - 26208*y
     577            sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]])
     578            sage: expected = [5, 12, 13, 14]
     579            sage: all(abs(e-a) < 0.001 for e,a in zip(expected, f.critical_points()))
     580            True
     581        """
     582        from sage.calculus.calculus import maxima
     583        x = QQ[self.default_variable()].gen()
     584        crit_pts = []
     585        for (a,b), f in self.list():
     586            for root in maxima.allroots(SR(f).diff(x)==0):
     587                root = float(root.rhs())
     588                if a < root < b:
     589                    crit_pts.append(root)
     590        return crit_pts
     591       
     592    def base_ring(self):
     593        """
     594        Returns the base ring of the function pieces.   This
     595        is useful when this class is extended.
     596
     597        EXAMPLES::
     598       
     599            sage: f1(x) = 1
     600            sage: f2(x) = 1-x
     601            sage: f3(x) = x^2-5
     602            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3]])
     603            sage: base_ring(f)
     604            Symbolic Ring
     605
     606        ::
     607
     608            sage: R.<x> = QQ[]
     609            sage: f1 = x^0
     610            sage: f2 = 10*x - x^2
     611            sage: f3 = 3*x^4 - 156*x^3 + 3036*x^2 - 26208*x
     612            sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]])
     613            sage: f.base_ring()
     614            Rational Field
     615        """
     616        return (self.functions()[0]).base_ring()
     617
     618    def end_points(self):
     619        """
     620        Returns a list of all interval endpoints for this function.
     621       
     622        EXAMPLES::
     623       
     624            sage: f1(x) = 1
     625            sage: f2(x) = 1-x
     626            sage: f3(x) = x^2-5
     627            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3]])
     628            sage: f.end_points()
     629            [0, 1, 2, 3]
     630        """
     631        intervals = self.intervals()
     632        return [ intervals[0][0] ] + [b for a,b in intervals]
     633
     634    def __call__(self,x0):
     635        """
     636        Evaluates self at x0. Returns the average value of the jump if x0
     637        is an interior endpoint of one of the intervals of self and the
     638        usual value otherwise.
     639       
     640        EXAMPLES::
     641       
     642            sage: f1(x) = 1
     643            sage: f2(x) = 1-x
     644            sage: f3(x) = exp(x)
     645            sage: f4(x) = sin(2*x)
     646            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     647            sage: f(0.5)
     648            1
     649            sage: f(5/2)
     650            e^(5/2)
     651            sage: f(5/2).n()
     652            12.1824939607035
     653            sage: f(1)
     654            1/2
     655        """
     656        #x0 = QQ(x0) ## does not allow for evaluation at pi
     657        n = self.length()
     658        endpts = self.end_points()
     659        for i in range(1,n):
     660            if x0 == endpts[i]:
     661                return (self.functions()[i-1](x0) + self.functions()[i](x0))/2
     662        if x0 == endpts[0]:
     663            return self.functions()[0](x0)
     664        if x0 == endpts[n]:
     665            return self.functions()[n-1](x0)
     666        for i in range(n):
     667            if endpts[i] < x0 < endpts[i+1]:
     668                return self.functions()[i](x0)
     669        raise ValueError,"Value not defined outside of domain."
     670
     671    def which_function(self,x0):
     672        """
     673        Returns the function piece used to evaluate self at x0.
     674       
     675        EXAMPLES::
     676       
     677            sage: f1(z) = z
     678            sage: f2(x) = 1-x
     679            sage: f3(y) = exp(y)
     680            sage: f4(t) = sin(2*t)
     681            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     682            sage: f.which_function(3/2)
     683            x |--> -x + 1
     684        """
     685        for (a,b), f in self.list():
     686            if a <= x0 <= b:
     687                return f
     688        raise ValueError,"Function not defined outside of domain."
     689
     690    def default_variable(self):
     691        r"""
     692        Return the default variable. The default variable is defined as the
     693        first variable in the first piece that has a variable. If no pieces have
     694        a variable (each piece is a constant value), `x` is returned.
     695
     696        The result is cached.
     697
     698        AUTHOR: Paul Butler
     699
     700        EXAMPLES::
     701       
     702            sage: f1(x) = 1
     703            sage: f2(x) = 5*x
     704            sage: p = Piecewise([[(0,1),f1],[(1,4),f2]])
     705            sage: p.default_variable()
     706            x
     707
     708            sage: f1 = 3*var('y')
     709            sage: p = Piecewise([[(0,1),4],[(1,4),f1]])
     710            sage: p.default_variable()
     711            y
     712
     713        """
     714        try:
     715            return self.__default_variable
     716        except AttributeError:
     717            pass
     718        for _, fun in self._list:
     719            try:
     720                fun = SR(fun)
     721                if fun.variables():
     722                    v = fun.variables()[0]
     723                    self.__default_variable = v
     724                    return v
     725            except TypeError:
     726                # pass if fun is lambda function
     727                pass
     728        # default to x
     729        v = var('x')
     730        self.__default_value = v
     731        return v
     732
     733    def integral(self, x=None, a=None, b=None, definite=False):
     734        r"""
     735        By default, returns the indefinite integral of the function.
     736        If definite=True is given, returns the definite integral.
     737
     738        AUTHOR:
     739
     740        - Paul Butler
     741
     742        EXAMPLES::
     743
     744            sage: f1(x) = 1-x
     745            sage: f = Piecewise([[(0,1),1],[(1,2),f1]])
     746            sage: f.integral(definite=True)
     747            1/2
     748       
     749        ::
     750       
     751            sage: f1(x) = -1
     752            sage: f2(x) = 2
     753            sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
     754            sage: f.integral(definite=True)
     755            1/2*pi
     756           
     757            sage: f1(x) = 2
     758            sage: f2(x) = 3 - x
     759            sage: f = Piecewise([[(-2, 0), f1], [(0, 3), f2]])
     760            sage: f.integral()
     761            Piecewise defined function with 2 parts, [[(-2, 0), x |--> 2*x + 4], [(0, 3), x |--> -1/2*x^2 + 3*x + 4]]
     762
     763            sage: f1(y) = -1
     764            sage: f2(y) = y + 3
     765            sage: f3(y) = -y - 1
     766            sage: f4(y) = y^2 - 1
     767            sage: f5(y) = 3
     768            sage: f = Piecewise([[(-4,-3),f1],[(-3,-2),f2],[(-2,0),f3],[(0,2),f4],[(2,3),f5]])
     769            sage: F = f.integral(y)
     770            sage: F
     771            Piecewise defined function with 5 parts, [[(-4, -3), y |--> -y - 4], [(-3, -2), y |--> 1/2*y^2 + 3*y + 7/2], [(-2, 0), y |--> -1/2*y^2 - y - 1/2], [(0, 2), y |--> 1/3*y^3 - y - 1/2], [(2, 3), y |--> 3*y - 35/6]]
     772           
     773        Ensure results are consistent with FTC::
     774
     775            sage: F(-3) - F(-4)
     776            -1
     777            sage: F(-1) - F(-3)
     778            1
     779            sage: F(2) - F(0)
     780            2/3
     781            sage: f.integral(y, 0, 2)
     782            2/3
     783            sage: F(3) - F(-4)
     784            19/6
     785            sage: f.integral(y, -4, 3)
     786            19/6
     787            sage: f.integral(definite=True)
     788            19/6
     789           
     790        ::
     791
     792            sage: f1(y) = (y+3)^2
     793            sage: f2(y) = y+3
     794            sage: f3(y) = 3
     795            sage: f = Piecewise([[(-infinity, -3), f1], [(-3, 0), f2], [(0, infinity), f3]])
     796            sage: f.integral()
     797            Piecewise defined function with 3 parts, [[(-Infinity, -3), y |--> 1/3*y^3 + 3*y^2 + 9*y + 9], [(-3, 0), y |--> 1/2*y^2 + 3*y + 9/2], [(0, +Infinity), y |--> 3*y + 9/2]]
     798
     799        ::
     800
     801            sage: f1(x) = e^(-abs(x))
     802            sage: f = Piecewise([[(-infinity, infinity), f1]])
     803            sage: f.integral(definite=True)
     804            2
     805            sage: f.integral()
     806            Piecewise defined function with 1 parts, [[(-Infinity, +Infinity), x |--> -1/2*((sgn(x) - 1)*e^(2*x) - 2*e^x*sgn(x) + sgn(x) + 1)*e^(-x) - 1]]
     807
     808        ::
     809
     810            sage: f = Piecewise([((0, 5), cos(x))])
     811            sage: f.integral()
     812            Piecewise defined function with 1 parts, [[(0, 5), x |--> sin(x)]]
     813
     814
     815        TESTS:
     816
     817        Verify that piecewise integrals of zero work (trac #10841)::
     818
     819            sage: f0(x) = 0
     820            sage: f = Piecewise([[(0,1),f0]])
     821            sage: f.integral(x,0,1)
     822            0
     823            sage: f = Piecewise([[(0,1), 0]])
     824            sage: f.integral(x,0,1)
     825            0
     826            sage: f = Piecewise([[(0,1), SR(0)]])
     827            sage: f.integral(x,0,1)
     828            0
     829
     830        """
     831        if a != None and b != None:
     832            F = self.integral(x)
     833            return F(b) - F(a)
     834
     835        if a != None or b != None:
     836            raise TypeError, 'only one endpoint given'
     837
     838        area = 0 # cumulative definite integral of parts to the left of the current interval
     839        integrand_pieces = self.list()
     840        integrand_pieces.sort()
     841        new_pieces = []
     842
     843        if x == None:
     844            x = self.default_variable()
     845       
     846        # The integral is computed by iterating over the pieces in order.
     847        # The definite integral for each piece is calculated and accumulated in `area`.
     848        # Thus at any time, `area` represents the definite integral of all the pieces
     849        # encountered so far. The indefinite integral of each piece is also calculated,
     850        # and the `area` before each piece is added to the piece.
     851        #
     852        # If a definite integral is requested, `area` is returned.
     853        # Otherwise, a piecewise function is constructed from the indefinite integrals
     854        # and returned.
     855        #
     856        # An exception is made if integral is called on a piecewise function
     857        # that starts at -infinity. In this case, we do not try to calculate the
     858        # definite integral of the first piece, and the value of `area` remains 0
     859        # after the first piece.
     860
     861        for (start, end), fun in integrand_pieces:
     862            fun = SR(fun)
     863            if start == -infinity and not definite:
     864                fun_integrated = fun.integral(x, end, x)
     865            else:
     866                try:
     867                    assume(start < x)
     868                except ValueError: # Assumption is redundant
     869                    pass
     870                fun_integrated = fun.integral(x, start, x) + area
     871                forget(start < x)
     872                if definite or end != infinity:
     873                    area += fun.integral(x, start, end)
     874            new_pieces.append([(start, end), SR(fun_integrated).function(x)])
     875
     876        if definite:
     877            return SR(area)
     878        else:
     879            return Piecewise(new_pieces)
     880
     881    def convolution(self, other):
     882        """
     883        Returns the convolution function,
     884        `f*g(t)=\int_{-\infty}^\infty f(u)g(t-u)du`, for compactly
     885        supported `f,g`.
     886       
     887        EXAMPLES::
     888       
     889            sage: x = PolynomialRing(QQ,'x').gen()
     890            sage: f = Piecewise([[(0,1),1*x^0]])  ## example 0
     891            sage: g = f.convolution(f)
     892            sage: h = f.convolution(g)
     893            sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1));
     894            sage: # Type show(P+Q+R) to view
     895            sage: f = Piecewise([[(0,1),1*x^0],[(1,2),2*x^0],[(2,3),1*x^0]])  ## example 1
     896            sage: g = f.convolution(f)
     897            sage: h = f.convolution(g)
     898            sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1));
     899            sage: # Type show(P+Q+R) to view
     900            sage: f = Piecewise([[(-1,1),1]])                             ## example 2
     901            sage: g = Piecewise([[(0,3),x]])
     902            sage: f.convolution(g)
     903            Piecewise defined function with 3 parts, [[(-1, 1), 0], [(1, 2), -3/2*x], [(2, 4), -3/2*x]]
     904            sage: g = Piecewise([[(0,3),1*x^0],[(3,4),2*x^0]])
     905            sage: f.convolution(g)
     906            Piecewise defined function with 5 parts, [[(-1, 1), x + 1], [(1, 2), 3], [(2, 3), x], [(3, 4), -x + 8], [(4, 5), -2*x + 10]]
     907        """
     908        f = self
     909        g = other
     910        M = min(min(f.end_points()),min(g.end_points()))
     911        N = max(max(f.end_points()),max(g.end_points()))
     912        R2 = PolynomialRing(QQ,2,names=["tt","uu"])
     913        tt,uu = R2.gens()
     914        conv = 0
     915        f0 = f.functions()[0]
     916        g0 = g.functions()[0]
     917        R1 = f0.parent()
     918        xx = R1.gen()
     919        var = repr(xx)
     920        if len(f.intervals())==1 and len(g.intervals())==1:
     921            f = f.unextend()
     922            g = g.unextend()
     923            a1 = f.intervals()[0][0]
     924            a2 = f.intervals()[0][1]
     925            b1 = g.intervals()[0][0]
     926            b2 = g.intervals()[0][1]
     927            i1 = repr(f0).replace(var,repr(uu))
     928            i2 = repr(g0).replace(var,"("+repr(tt-uu)+")")
     929            cmd1 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, a1,    tt-b1)    ## if a1+b1 < tt < a2+b1
     930            cmd2 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, tt-b2, tt-b1)    ## if a1+b2 < tt < a2+b1
     931            cmd3 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, tt-b2, a2)       ## if a1+b2 < tt < a2+b2
     932            cmd4 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, a1, a2)          ## if a2+b1 < tt < a1+b2
     933            conv1 = maxima.eval(cmd1)
     934            conv2 = maxima.eval(cmd2)
     935            conv3 = maxima.eval(cmd3)
     936            conv4 = maxima.eval(cmd4)
     937            # this is a very, very, very ugly hack
     938            x = PolynomialRing(QQ,'x').gen()
     939            fg1 = sage_eval(conv1.replace("tt",var), {'x':x}) ## should be = R2(conv1)
     940            fg2 = sage_eval(conv2.replace("tt",var), {'x':x}) ## should be = R2(conv2)
     941            fg3 = sage_eval(conv3.replace("tt",var), {'x':x}) ## should be = R2(conv3)
     942            fg4 = sage_eval(conv4.replace("tt",var), {'x':x}) ## should be = R2(conv4)
     943            if a1-b1<a2-b2:
     944                if a2+b1!=a1+b2:
     945                    h = Piecewise([[(a1+b1,a1+b2),fg1],[(a1+b2,a2+b1),fg4],[(a2+b1,a2+b2),fg3]])
     946                else:
     947                    h = Piecewise([[(a1+b1,a1+b2),fg1],[(a1+b2,a2+b2),fg3]])
     948            else:
     949                if a1+b2!=a2+b1:
     950                    h = Piecewise([[(a1+b1,a2+b1),fg1],[(a2+b1,a1+b2),fg2],[(a1+b2,a2+b2),fg3]])
     951                else:
     952                    h = Piecewise([[(a1+b1,a2+b1),fg1],[(a2+b1,a2+b2),fg3]])
     953            return h
     954       
     955        if len(f.intervals())>1 or len(g.intervals())>1:
     956            z = Piecewise([[(-3*abs(N-M),3*abs(N-M)),0*xx**0]])
     957            ff = f.functions()
     958            gg = g.functions()
     959            intvlsf = f.intervals()
     960            intvlsg = g.intervals()
     961            for i in range(len(ff)):
     962                for j in range(len(gg)):
     963                    f0 = Piecewise([[intvlsf[i],ff[i]]])
     964                    g0 = Piecewise([[intvlsg[j],gg[j]]])
     965                    h = g0.convolution(f0)
     966                    z = z + h
     967            return z.unextend()
     968   
     969    def derivative(self):
     970        r"""
     971        Returns the derivative (as computed by maxima)
     972        Piecewise(I,`(d/dx)(self|_I)`), as I runs over the
     973        intervals belonging to self. self must be piecewise polynomial.
     974       
     975        EXAMPLES::
     976       
     977            sage: f1(x) = 1
     978            sage: f2(x) = 1-x
     979            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     980            sage: f.derivative()
     981            Piecewise defined function with 2 parts, [[(0, 1), x |--> 0], [(1, 2), x |--> -1]]
     982            sage: f1(x) = -1
     983            sage: f2(x) = 2
     984            sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
     985            sage: f.derivative()
     986            Piecewise defined function with 2 parts, [[(0, 1/2*pi), x |--> 0], [(1/2*pi, pi), x |--> 0]]
     987           
     988        ::
     989       
     990            sage: f = Piecewise([[(0,1), (x * 2)]], x)
     991            sage: f.derivative()
     992            Piecewise defined function with 1 parts, [[(0, 1), x |--> 2]]
     993        """
     994        x = self.default_variable()
     995        dlist = [[(a, b), derivative(f(x), x).function(x)] for (a,b),f in self.list()]
     996        return Piecewise(dlist)
     997 
     998    def tangent_line(self, pt):
     999        """
     1000        Computes the linear function defining the tangent line of the
     1001        piecewise function self.
     1002       
     1003        EXAMPLES::
     1004       
     1005            sage: f1(x) = x^2
     1006            sage: f2(x) = 5-x^3+x
     1007            sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
     1008            sage: tf = f.tangent_line(0.9) ## tangent line at x=0.9
     1009            sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
     1010            sage: Q = tf.plot(rgbcolor=(0.7,0.2,0.2), plot_points=40)
     1011            sage: P + Q
     1012        """
     1013        pt = QQ(pt)
     1014        R = QQ[self.default_variable()]
     1015        x = R.gen()
     1016        der = self.derivative()
     1017        tanline = (x-pt)*der(pt)+self(pt)
     1018        dlist = [[(a, b), tanline] for (a,b),f in self.list()]
     1019        return Piecewise(dlist)
     1020       
     1021    def plot(self, *args, **kwds):
     1022        """
     1023        Returns the plot of self.
     1024       
     1025        Keyword arguments are passed onto the plot command for each piece
     1026        of the function. E.g., the plot_points keyword affects each
     1027        segment of the plot.
     1028       
     1029        EXAMPLES::
     1030       
     1031            sage: f1(x) = 1
     1032            sage: f2(x) = 1-x
     1033            sage: f3(x) = exp(x)
     1034            sage: f4(x) = sin(2*x)
     1035            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     1036            sage: P = f.plot(rgbcolor=(0.7,0.1,0), plot_points=40)
     1037            sage: P
     1038       
     1039        Remember: to view this, type show(P) or P.save("path/myplot.png")
     1040        and then open it in a graphics viewer such as GIMP.
     1041
     1042        TESTS:
     1043
     1044        We should not add each piece to the legend individually, since
     1045        this creates duplicates (:trac:`12651`). This tests that only
     1046        one of the graphics objects in the plot has a non-``None``
     1047        ``legend_label``::
     1048
     1049            sage: f1 = sin(x)
     1050            sage: f2 = cos(x)
     1051            sage: f = piecewise([[(-1,0), f1],[(0,1), f2]])
     1052            sage: p = f.plot(legend_label='$f(x)$')
     1053            sage: lines = [
     1054            ...     line
     1055            ...     for line in p._objects
     1056            ...     if line.options()['legend_label'] is not None ]
     1057            sage: len(lines)
     1058            1
     1059        """
     1060        from sage.plot.all import plot, Graphics
     1061
     1062        g = Graphics()
     1063
     1064        for i, ((a,b), f) in enumerate(self.list()):
     1065            # If it's the first piece, pass all arguments. Otherwise,
     1066            # filter out 'legend_label' so that we don't add each
     1067            # piece to the legend separately (trac #12651).
     1068            if i != 0 and 'legend_label' in kwds:
     1069                del kwds['legend_label']
     1070
     1071            g += plot(f, a, b, *args, **kwds)
     1072
     1073        return g
     1074
     1075    def fourier_series_cosine_coefficient(self,n,L):
     1076        r"""
     1077        Returns the n-th Fourier series coefficient of
     1078        `\cos(n\pi x/L)`, `a_n`.
     1079       
     1080        INPUT:
     1081       
     1082       
     1083        -  ``self`` - the function f(x), defined over -L x L
     1084       
     1085        -  ``n`` - an integer n=0
     1086       
     1087        -  ``L`` - (the period)/2
     1088       
     1089       
     1090        OUTPUT:
     1091        `a_n = \frac{1}{L}\int_{-L}^L f(x)\cos(n\pi x/L)dx`
     1092       
     1093        EXAMPLES::
     1094       
     1095            sage: f(x) = x^2
     1096            sage: f = Piecewise([[(-1,1),f]])
     1097            sage: f.fourier_series_cosine_coefficient(2,1)
     1098            pi^(-2)
     1099            sage: f(x) = x^2
     1100            sage: f = Piecewise([[(-pi,pi),f]])
     1101            sage: f.fourier_series_cosine_coefficient(2,pi)
     1102            1
     1103            sage: f1(x) = -1
     1104            sage: f2(x) = 2
     1105            sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
     1106            sage: f.fourier_series_cosine_coefficient(5,pi)
     1107            -3/5/pi
     1108        """
     1109        from sage.all import cos, pi
     1110        x = var('x')
     1111        result = sum([(f(x)*cos(pi*x*n/L)/L).integrate(x, a, b)
     1112                      for (a,b), f in self.list()])
     1113        if is_Expression(result):
     1114            return result.simplify_trig()
     1115        return result
     1116   
     1117    def fourier_series_sine_coefficient(self,n,L):
     1118        r"""
     1119        Returns the n-th Fourier series coefficient of
     1120        `\sin(n\pi x/L)`, `b_n`.
     1121       
     1122        INPUT:
     1123       
     1124       
     1125        -  ``self`` - the function f(x), defined over -L x L
     1126       
     1127        -  ``n`` - an integer n0
     1128       
     1129        -  ``L`` - (the period)/2
     1130       
     1131       
     1132        OUTPUT:
     1133        `b_n = \frac{1}{L}\int_{-L}^L f(x)\sin(n\pi x/L)dx`
     1134       
     1135        EXAMPLES::
     1136       
     1137            sage: f(x) = x^2
     1138            sage: f = Piecewise([[(-1,1),f]])
     1139            sage: f.fourier_series_sine_coefficient(2,1)  # L=1, n=2
     1140            0
     1141        """
     1142        from sage.all import sin, pi
     1143        x = var('x')
     1144        result = sum([(f(x)*sin(pi*x*n/L)/L).integrate(x, a, b)
     1145                      for (a,b), f in self.list()])
     1146        if is_Expression(result):
     1147            return result.simplify_trig()
     1148        return result
     1149
     1150    def _fourier_series_helper(self, N, L, scale_function):
     1151        r"""
     1152        A helper function for the construction of Fourier series. The
     1153        argument scale_function is a function which takes in n,
     1154        representing the `n^{th}` coefficient, and return an
     1155        expression to scale the sine and cosine coefficients by.
     1156       
     1157        EXAMPLES::
     1158       
     1159            sage: f(x) = x^2
     1160            sage: f = Piecewise([[(-1,1),f]])
     1161            sage: f._fourier_series_helper(3, 1, lambda n: 1)
     1162            cos(2*pi*x)/pi^2 - 4*cos(pi*x)/pi^2 + 1/3
     1163        """
     1164        from sage.all import pi, sin, cos, srange
     1165        x = self.default_variable()
     1166        a0 = self.fourier_series_cosine_coefficient(0,L)
     1167        result = a0/2 + sum([(self.fourier_series_cosine_coefficient(n,L)*cos(n*pi*x/L) +
     1168                              self.fourier_series_sine_coefficient(n,L)*sin(n*pi*x/L))*
     1169                             scale_function(n)
     1170                             for n in srange(1,N)])
     1171        return result.expand()
     1172       
     1173
     1174    def fourier_series_partial_sum(self,N,L):
     1175        r"""
     1176        Returns the partial sum
     1177       
     1178        .. math::
     1179       
     1180           f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
     1181       
     1182        as a string.
     1183       
     1184        EXAMPLE::
     1185       
     1186            sage: f(x) = x^2
     1187            sage: f = Piecewise([[(-1,1),f]])
     1188            sage: f.fourier_series_partial_sum(3,1)
     1189            cos(2*pi*x)/pi^2 - 4*cos(pi*x)/pi^2 + 1/3
     1190            sage: f1(x) = -1
     1191            sage: f2(x) = 2
     1192            sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
     1193            sage: f.fourier_series_partial_sum(3,pi)
     1194            -3*cos(x)/pi - 3*sin(2*x)/pi + 3*sin(x)/pi - 1/4
     1195        """
     1196        return self._fourier_series_helper(N, L, lambda n: 1)
     1197 
     1198    def fourier_series_partial_sum_cesaro(self,N,L):
     1199        r"""
     1200        Returns the Cesaro partial sum
     1201       
     1202        .. math::
     1203       
     1204           f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N (1-n/N)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
     1205       
     1206       
     1207        as a string. This is a "smoother" partial sum - the Gibbs
     1208        phenomenon is mollified.
     1209       
     1210        EXAMPLE::
     1211       
     1212            sage: f(x) = x^2
     1213            sage: f = Piecewise([[(-1,1),f]])
     1214            sage: f.fourier_series_partial_sum_cesaro(3,1)
     1215            1/3*cos(2*pi*x)/pi^2 - 8/3*cos(pi*x)/pi^2 + 1/3
     1216            sage: f1(x) = -1
     1217            sage: f2(x) = 2
     1218            sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
     1219            sage: f.fourier_series_partial_sum_cesaro(3,pi)
     1220            -2*cos(x)/pi - sin(2*x)/pi + 2*sin(x)/pi - 1/4
     1221        """
     1222        return self._fourier_series_helper(N, L, lambda n: 1-n/N)
     1223
     1224    def fourier_series_partial_sum_hann(self,N,L):
     1225        r"""
     1226        Returns the Hann-filtered partial sum (named after von Hann, not
     1227        Hamming)
     1228       
     1229        .. math::
     1230       
     1231           f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N H_N(n)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
     1232       
     1233        as a string, where `H_N(x) = (1+\cos(\pi x/N))/2`. This is
     1234        a "smoother" partial sum - the Gibbs phenomenon is mollified.
     1235       
     1236        EXAMPLE::
     1237       
     1238            sage: f(x) = x^2
     1239            sage: f = Piecewise([[(-1,1),f]])
     1240            sage: f.fourier_series_partial_sum_hann(3,1)
     1241            1/4*cos(2*pi*x)/pi^2 - 3*cos(pi*x)/pi^2 + 1/3
     1242            sage: f1(x) = -1
     1243            sage: f2(x) = 2
     1244            sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
     1245            sage: f.fourier_series_partial_sum_hann(3,pi)
     1246            -9/4*cos(x)/pi - 3/4*sin(2*x)/pi + 9/4*sin(x)/pi - 1/4
     1247        """
     1248        from sage.all import cos, pi
     1249        return self._fourier_series_helper(N, L, lambda n: (1+cos(pi*n/N))/2)
     1250
     1251    def fourier_series_partial_sum_filtered(self,N,L,F):
     1252        r"""
     1253        Returns the "filtered" partial sum
     1254       
     1255        .. math::
     1256       
     1257           f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N F_n*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
     1258       
     1259        as a string, where `F = [F_1,F_2, ..., F_{N}]` is a list
     1260        of length `N` consisting of real numbers. This can be used
     1261        to plot FS solutions to the heat and wave PDEs.
     1262       
     1263        EXAMPLE::
     1264       
     1265            sage: f(x) = x^2
     1266            sage: f = Piecewise([[(-1,1),f]])
     1267            sage: f.fourier_series_partial_sum_filtered(3,1,[1,1,1])
     1268            cos(2*pi*x)/pi^2 - 4*cos(pi*x)/pi^2 + 1/3
     1269            sage: f1(x) = -1
     1270            sage: f2(x) = 2
     1271            sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
     1272            sage: f.fourier_series_partial_sum_filtered(3,pi,[1,1,1])
     1273            -3*cos(x)/pi - 3*sin(2*x)/pi + 3*sin(x)/pi - 1/4
     1274        """
     1275        return self._fourier_series_helper(N, L, lambda n: F[n])
     1276       
     1277    def plot_fourier_series_partial_sum(self,N,L,xmin,xmax, **kwds):
     1278        r"""
     1279        Plots the partial sum
     1280       
     1281        .. math::
     1282       
     1283           f(x) \sim \frac{a_0}{2} +  sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
     1284       
     1285        over xmin x xmin.
     1286       
     1287        EXAMPLE::
     1288       
     1289            sage: f1(x) = -2
     1290            sage: f2(x) = 1
     1291            sage: f3(x) = -1
     1292            sage: f4(x) = 2
     1293            sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]])
     1294            sage: P = f.plot_fourier_series_partial_sum(3,pi,-5,5)    # long time
     1295            sage: f1(x) = -1
     1296            sage: f2(x) = 2
     1297            sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
     1298            sage: P = f.plot_fourier_series_partial_sum(15,pi,-5,5)   # long time
     1299       
     1300        Remember, to view this type show(P) or P.save("path/myplot.png")
     1301        and then open it in a graphics viewer such as GIMP.
     1302        """
     1303        from sage.plot.all import plot
     1304        return plot(self.fourier_series_partial_sum(N,L), xmin, xmax, **kwds)
     1305
     1306    def plot_fourier_series_partial_sum_cesaro(self,N,L,xmin,xmax, **kwds):
     1307        r"""
     1308        Plots the partial sum
     1309       
     1310        .. math::
     1311       
     1312                     f(x) \sim \frac{a_0}{2} +                     \sum_{n=1}^N (1-n/N)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
     1313       
     1314       
     1315        over xmin x xmin. This is a "smoother" partial sum - the Gibbs
     1316        phenomenon is mollified.
     1317       
     1318        EXAMPLE::
     1319       
     1320            sage: f1(x) = -2
     1321            sage: f2(x) = 1
     1322            sage: f3(x) = -1
     1323            sage: f4(x) = 2
     1324            sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]])
     1325            sage: P = f.plot_fourier_series_partial_sum_cesaro(3,pi,-5,5)    # long time
     1326            sage: f1(x) = -1
     1327            sage: f2(x) = 2
     1328            sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
     1329            sage: P = f.plot_fourier_series_partial_sum_cesaro(15,pi,-5,5)   # long time
     1330       
     1331        Remember, to view this type show(P) or P.save("path/myplot.png")
     1332        and then open it in a graphics viewer such as GIMP.
     1333        """     
     1334        from sage.plot.all import plot
     1335        return plot(self.fourier_series_partial_sum_cesaro(N,L), xmin, xmax, **kwds)
     1336   
     1337    def plot_fourier_series_partial_sum_hann(self,N,L,xmin,xmax, **kwds):
     1338        r"""
     1339        Plots the partial sum
     1340       
     1341        .. math::
     1342       
     1343           f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N H_N(n)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
     1344       
     1345       
     1346        over xmin x xmin, where H_N(x) = (0.5)+(0.5)\*cos(x\*pi/N) is the
     1347        N-th Hann filter.
     1348       
     1349        EXAMPLE::
     1350       
     1351            sage: f1(x) = -2
     1352            sage: f2(x) = 1
     1353            sage: f3(x) = -1
     1354            sage: f4(x) = 2
     1355            sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]])
     1356            sage: P = f.plot_fourier_series_partial_sum_hann(3,pi,-5,5)    # long time
     1357            sage: f1(x) = -1
     1358            sage: f2(x) = 2
     1359            sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
     1360            sage: P = f.plot_fourier_series_partial_sum_hann(15,pi,-5,5)   # long time
     1361       
     1362        Remember, to view this type show(P) or P.save("path/myplot.png")
     1363        and then open it in a graphics viewer such as GIMP.
     1364        """     
     1365        from sage.plot.all import plot
     1366        return plot(self.fourier_series_partial_sum_hann(N,L), xmin, xmax, **kwds)
     1367       
     1368    def plot_fourier_series_partial_sum_filtered(self,N,L,F,xmin,xmax, **kwds):
     1369        r"""
     1370        Plots the partial sum
     1371       
     1372        .. math::
     1373       
     1374                     f(x) \sim \frac{a_0}{2} +                     \sum_{n=1}^N F_n*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
     1375       
     1376       
     1377        over xmin x xmin, where `F = [F_1,F_2, ..., F_{N}]` is a
     1378        list of length `N` consisting of real numbers. This can be
     1379        used to plot FS solutions to the heat and wave PDEs.
     1380       
     1381        EXAMPLE::
     1382       
     1383            sage: f1(x) = -2
     1384            sage: f2(x) = 1
     1385            sage: f3(x) = -1
     1386            sage: f4(x) = 2
     1387            sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]])
     1388            sage: P = f.plot_fourier_series_partial_sum_filtered(3,pi,[1]*3,-5,5)    # long time
     1389            sage: f1(x) = -1
     1390            sage: f2(x) = 2
     1391            sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f1],[(pi/2,pi),f2]])
     1392            sage: P = f.plot_fourier_series_partial_sum_filtered(15,pi,[1]*15,-5,5)   # long time
     1393       
     1394        Remember, to view this type show(P) or P.save("path/myplot.png")
     1395        and then open it in a graphics viewer such as GIMP.
     1396        """
     1397        from sage.plot.all import plot
     1398        return plot(self.fourier_series_partial_sum_filtered(N,L,F), xmin, xmax, **kwds)
     1399               
     1400    def fourier_series_value(self,x,L):
     1401        r"""
     1402        Returns the value of the Fourier series coefficient of self at
     1403        `x`,
     1404       
     1405       
     1406        .. math::
     1407       
     1408                     f(x) \sim \frac{a_0}{2} +                     \sum_{n=1}^\infty [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         \ \ \ -L<x<L.         
     1409       
     1410       
     1411        This method applies to piecewise non-polynomial functions as well.
     1412       
     1413        INPUT:
     1414       
     1415       
     1416        -  ``self`` - the function f(x), defined over -L x L
     1417       
     1418        -  ``x`` - a real number
     1419       
     1420        -  ``L`` - (the period)/2
     1421       
     1422       
     1423        OUTPUT: `(f^*(x+)+f^*(x-)/2`, where `f^*` denotes
     1424        the function `f` extended to `\RR` with period
     1425        `2L` (Dirichlet's Theorem for Fourier series).
     1426       
     1427        EXAMPLES::
     1428       
     1429            sage: f1(x) = 1
     1430            sage: f2(x) = 1-x
     1431            sage: f3(x) = exp(x)
     1432            sage: f4(x) = sin(2*x)
     1433            sage: f = Piecewise([[(-10,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     1434            sage: f.fourier_series_value(101,10) 
     1435            1/2
     1436            sage: f.fourier_series_value(100,10)
     1437            1
     1438            sage: f.fourier_series_value(10,10)
     1439            1/2*sin(20) + 1/2
     1440            sage: f.fourier_series_value(20,10)
     1441            1
     1442            sage: f.fourier_series_value(30,10)
     1443            1/2*sin(20) + 1/2
     1444            sage: f1(x) = -1
     1445            sage: f2(x) = 2
     1446            sage: f = Piecewise([[(-pi,0),lambda x:0],[(0,pi/2),f1],[(pi/2,pi),f2]])
     1447            sage: f.fourier_series_value(-1,pi)
     1448            0
     1449            sage: f.fourier_series_value(20,pi)
     1450            -1
     1451            sage: f.fourier_series_value(pi/2,pi)
     1452            1/2
     1453        """
     1454        xnew = x - int(RR(x/(2*L)))*2*L
     1455        endpts = self.end_points()
     1456        if xnew == endpts[0] or xnew == endpts[-1]:
     1457            return (self.functions()[0](endpts[0]) + self.functions()[-1](endpts[-1]))/2
     1458        else:
     1459            return self(xnew)
     1460
     1461    def cosine_series_coefficient(self,n,L):
     1462        r"""
     1463        Returns the n-th cosine series coefficient of
     1464        `\cos(n\pi x/L)`, `a_n`.
     1465       
     1466        INPUT:
     1467       
     1468       
     1469        -  ``self`` - the function f(x), defined over 0 x L (no
     1470           checking is done to insure this)
     1471       
     1472        -  ``n`` - an integer n=0
     1473       
     1474        -  ``L`` - (the period)/2
     1475       
     1476       
     1477        OUTPUT:
     1478        `a_n = \frac{2}{L}\int_{-L}^L f(x)\cos(n\pi x/L)dx` such
     1479        that
     1480       
     1481        .. math::
     1482       
     1483                     f(x) \sim \frac{a_0}{2} +                     \sum_{n=1}^\infty a_n\cos(\frac{n\pi x}{L}),\ \ 0<x<L.         
     1484       
     1485       
     1486       
     1487        EXAMPLES::
     1488       
     1489            sage: f(x) = x
     1490            sage: f = Piecewise([[(0,1),f]])
     1491            sage: f.cosine_series_coefficient(2,1) 
     1492            0
     1493            sage: f.cosine_series_coefficient(3,1)
     1494            -4/9/pi^2
     1495            sage: f1(x) = -1
     1496            sage: f2(x) = 2
     1497            sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
     1498            sage: f.cosine_series_coefficient(2,pi)
     1499            0
     1500            sage: f.cosine_series_coefficient(3,pi)
     1501            2/pi
     1502            sage: f.cosine_series_coefficient(111,pi)
     1503            2/37/pi
     1504            sage: f1 = lambda x: x*(pi-x)
     1505            sage: f = Piecewise([[(0,pi),f1]])
     1506            sage: f.cosine_series_coefficient(0,pi)
     1507            1/3*pi^2
     1508
     1509        """
     1510        from sage.all import cos, pi
     1511        x = var('x')
     1512        result = sum([(2*f(x)*cos(pi*x*n/L)/L).integrate(x, a, b)
     1513                      for (a,b), f in self.list()])
     1514        if is_Expression(result):
     1515            return result.simplify_trig()
     1516        return result
     1517
     1518
     1519    def sine_series_coefficient(self,n,L):
     1520        r"""
     1521        Returns the n-th sine series coefficient of
     1522        `\sin(n\pi x/L)`, `b_n`.
     1523       
     1524        INPUT:
     1525       
     1526        -  ``self`` - the function f(x), defined over 0 x L (no
     1527           checking is done to insure this)
     1528       
     1529        -  ``n`` - an integer n0
     1530       
     1531        -  ``L`` - (the period)/2
     1532       
     1533        OUTPUT:
     1534
     1535        `b_n = \frac{2}{L}\int_{-L}^L f(x)\sin(n\pi x/L)dx` such
     1536        that
     1537       
     1538        .. math::
     1539       
     1540           f(x) \sim \sum_{n=1}^\infty b_n\sin(\frac{n\pi x}{L}),\ \ 0<x<L.         
     1541       
     1542        EXAMPLES::
     1543       
     1544            sage: f(x) = 1
     1545            sage: f = Piecewise([[(0,1),f]])
     1546            sage: f.sine_series_coefficient(2,1) 
     1547            0
     1548            sage: f.sine_series_coefficient(3,1)
     1549            4/3/pi
     1550        """
     1551        from sage.all import sin, pi
     1552        x = var('x')
     1553        result = sum([(2*f(x)*sin(pi*x*n/L)/L).integrate(x, a, b)
     1554                      for (a,b), f in self.list()])
     1555        if is_Expression(result):
     1556            return result.simplify_trig()
     1557        return result
     1558
     1559    def laplace(self, x='x', s='t'):
     1560        r"""
     1561        Returns the Laplace transform of self with respect to the variable
     1562        var.
     1563       
     1564        INPUT:
     1565       
     1566       
     1567        -  ``x`` - variable of self
     1568       
     1569        -  ``s`` - variable of Laplace transform.
     1570       
     1571       
     1572        We assume that a piecewise function is 0 outside of its domain and
     1573        that the left-most endpoint of the domain is 0.
     1574       
     1575        EXAMPLES::
     1576       
     1577            sage: x, s, w = var('x, s, w')
     1578            sage: f = Piecewise([[(0,1),1],[(1,2), 1-x]])
     1579            sage: f.laplace(x, s)
     1580            -e^(-s)/s + (s + 1)*e^(-2*s)/s^2 + 1/s - e^(-s)/s^2
     1581            sage: f.laplace(x, w)
     1582            -e^(-w)/w + (w + 1)*e^(-2*w)/w^2 + 1/w - e^(-w)/w^2
     1583           
     1584        ::
     1585       
     1586            sage: y, t = var('y, t')
     1587            sage: f = Piecewise([[(1,2), 1-y]])
     1588            sage: f.laplace(y, t)
     1589            (t + 1)*e^(-2*t)/t^2 - e^(-t)/t^2
     1590           
     1591        ::
     1592       
     1593            sage: s = var('s')
     1594            sage: t = var('t')
     1595            sage: f1(t) = -t
     1596            sage: f2(t) = 2
     1597            sage: f = Piecewise([[(0,1),f1],[(1,infinity),f2]])
     1598            sage: f.laplace(t,s)
     1599            (s + 1)*e^(-s)/s^2 + 2*e^(-s)/s - 1/s^2
     1600        """
     1601        from sage.all import assume, exp, forget
     1602        x = var(x)
     1603        s = var(s)
     1604        assume(s>0)
     1605        result =  sum([(SR(f)*exp(-s*x)).integral(x,a,b)
     1606                       for (a,b),f in self.list()])
     1607        forget(s>0)
     1608        return result
     1609
     1610    def _make_compatible(self, other):
     1611        """
     1612        Returns self and other extended to be defined on the same domain as
     1613        well as a refinement of their intervals. This is used for adding
     1614        and multiplying piecewise functions.
     1615       
     1616        EXAMPLES::
     1617       
     1618            sage: R.<x> = QQ[]
     1619            sage: f1 = Piecewise([[(0, 2), x]])
     1620            sage: f2 = Piecewise([[(1, 3), x^2]])
     1621            sage: f1._make_compatible(f2)
     1622            (Piecewise defined function with 2 parts, [[(0, 2), x], [(2, 3), 0]],
     1623            Piecewise defined function with 2 parts, [[(0, 1), 0], [(1, 3), x^2]],
     1624            [(0, 1), (1, 2), (2, 3)])
     1625        """
     1626        a1, b1 = self.domain()
     1627        a2, b2 = other.domain()
     1628        a = min(a1, a2)
     1629        b = max(b1, b2)
     1630        F = self.extend_by_zero_to(a,b)
     1631        G = other.extend_by_zero_to(a,b)
     1632        endpts = list(set(F.end_points()).union(set(G.end_points())))
     1633        endpts.sort()
     1634        return F, G, zip(endpts, endpts[1:])
     1635   
     1636    def __add__(self,other):
     1637        """
     1638        Returns the piecewise defined function which is the sum of self and
     1639        other. Does not require both domains be the same.
     1640       
     1641        EXAMPLES::
     1642       
     1643            sage: x = PolynomialRing(QQ,'x').gen()
     1644            sage: f1 = x^0
     1645            sage: f2 = 1-x
     1646            sage: f3 = 2*x
     1647            sage: f4 = 10-x
     1648            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     1649            sage: g1 = x-2
     1650            sage: g2 = x-5
     1651            sage: g = Piecewise([[(0,5),g1],[(5,10),g2]])
     1652            sage: h = f+g
     1653            sage: h
     1654            Piecewise defined function with 5 parts, [[(0, 1), x - 1], [(1, 2), -1], [(2, 3), 3*x - 2], [(3, 5), 8], [(5, 10), 5]]
     1655       
     1656        Note that in this case the functions must be defined using
     1657        polynomial expressions *not* using the lambda notation.
     1658        """
     1659        F, G, intervals = self._make_compatible(other)
     1660        fcn = []
     1661        for a,b in intervals:
     1662            fcn.append([(a,b), F.which_function(b)+G.which_function(b)])       
     1663        return Piecewise(fcn)
     1664       
     1665    def __mul__(self,other):
     1666        r"""
     1667        Returns the piecewise defined function which is the product of one
     1668        piecewise function (self) with another one (other).
     1669       
     1670        EXAMPLES::
     1671       
     1672            sage: x = PolynomialRing(QQ,'x').gen()
     1673            sage: f1 = x^0
     1674            sage: f2 = 1-x
     1675            sage: f3 = 2*x
     1676            sage: f4 = 10-x
     1677            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     1678            sage: g1 = x-2
     1679            sage: g2 = x-5
     1680            sage: g = Piecewise([[(0,5),g1],[(5,10),g2]])
     1681            sage: h = f*g
     1682            sage: h
     1683            Piecewise defined function with 5 parts, [[(0, 1), x - 2], [(1, 2), -x^2 + 3*x - 2], [(2, 3), 2*x^2 - 4*x], [(3, 5), -x^2 + 12*x - 20], [(5, 10), -x^2 + 15*x - 50]]
     1684            sage: g*(11/2)
     1685            Piecewise defined function with 2 parts, [[(0, 5), 11/2*x - 11], [(5, 10), 11/2*x - 55/2]]
     1686           
     1687        Note that in this method the functions must be defined using
     1688        polynomial expressions *not* using the lambda notation.
     1689        """
     1690        ## needed for scalar multiplication
     1691        if isinstance(other,Rational) or isinstance(other,Integer):
     1692            return Piecewise([[(a,b), other*f] for (a,b),f in self.list()])
     1693        else:
     1694            F, G, intervals = self._make_compatible(other)
     1695            fcn = []
     1696            for a,b in intervals:
     1697                fcn.append([(a,b),F.which_function(b)*G.which_function(b)])     
     1698            return Piecewise(fcn)
     1699
     1700    __rmul__ = __mul__
     1701
     1702    def __eq__(self,other):
     1703        """
     1704        Implements Boolean == operator.
     1705
     1706        EXAMPLES::
     1707       
     1708            sage: f1 = x^0
     1709            sage: f2 = 1-x
     1710            sage: f3 = 2*x
     1711            sage: f4 = 10-x
     1712            sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     1713            sage: g = Piecewise([[(0,1),1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
     1714            sage: f==g
     1715            True
     1716        """
     1717        return self.list()==other.list()
  • sage/functions/piecewise.py

    diff --git a/sage/functions/piecewise.py b/sage/functions/piecewise.py
    a b  
    11r"""
    22Piecewise-defined Functions
    33
    4 Sage implements a very simple class of piecewise-defined functions.
    5 Functions may be any type of symbolic expression. Infinite
    6 intervals are not supported. The endpoints of each interval must
    7 line up.
     4This module implement piecewise functions in a single variable. See
     5:mod:`sage.sets.real_set` for more information about how to construct
     6subsets of the real line for the domains.
     7
     8EXAMPLES::
     9
     10    sage: f = piecewise([((0,1), x^3), ([-1,0], -x^2)]);  f
     11    piecewise(x|-->x^3 on (0, 1), x|-->-x^2 on [-1, 0]; x)
     12    sage: 2*f
     13    2*piecewise(x|-->x^3 on (0, 1), x|-->-x^2 on [-1, 0]; x)
     14    sage: f(x=1/2)
     15    1/8
     16    sage: plot(f)    # not tested
    817
    918TODO:
    1019
    1120- Implement max/min location and values,
    1221
    13 - Need: parent object - ring of piecewise functions
    14 
    15 - This class should derive from an element-type class, and should
    16   define ``_add_``, ``_mul_``, etc. That will automatically take care
    17   of left multiplication and proper coercion. The coercion mentioned
    18   below for scalar mult on right is bad, since it only allows ints and
    19   rationals. The right way is to use an element class and only define
    20   ``_mul_``, and have a parent, so anything gets coerced properly.
    21 
    2222AUTHORS:
    2323
    2424- David Joyner (2006-04): initial version
     
    4444
    4545- Paul Butler (2009-01): added indefinite integration and default_variable
    4646
     47- Volker Braun (2013): Complete rewrite
     48
    4749TESTS::
    4850
    49     sage: R.<x> = QQ[]
    50     sage: f = Piecewise([[(0,1),1*x^0]])
    51     sage: 2*f
    52     Piecewise defined function with 1 parts, [[(0, 1), 2]]
     51    sage: fast_callable(f, vars=[x])(0.5)
     52    0.125000000000...
    5353"""
    5454
    5555#*****************************************************************************
    5656#       Copyright (C) 2006 William Stein <wstein@gmail.com>
    5757#                     2006 David Joyner <wdjoyner@gmail.com>
     58#                     2013 Volker Braun <vbraun.name@gmail.com>
    5859#
    59 #  Distributed under the terms of the GNU General Public License (GPL)
    60 #
    61 #    This code is distributed in the hope that it will be useful,
    62 #    but WITHOUT ANY WARRANTY; without even the implied warranty of
    63 #    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    64 #    General Public License for more details.
    65 #
    66 #  The full text of the GPL is available at:
    67 #
     60#  Distributed under the terms of the GNU General Public License (GPL),
     61#  version 2 or any later version.  The full text of the GPL is available at:
    6862#                  http://www.gnu.org/licenses/
    6963#*****************************************************************************
    7064
    71 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    72 from sage.misc.sage_eval import sage_eval
    73 from sage.rings.all import QQ, RR, Integer, Rational, infinity
    74 from sage.calculus.functional import derivative
    75 from sage.symbolic.expression import is_Expression
    76 from sage.symbolic.assumptions import assume, forget
     65from sage.symbolic.function import BuiltinFunction
     66from sage.sets.real_set import RealSet
     67from sage.symbolic.ring import SR
    7768
    78 from sage.calculus.calculus import SR, maxima
    79 from sage.calculus.all import var
    8069
    81 def piecewise(list_of_pairs, var=None):
    82     """
    83     Returns a piecewise function from a list of (interval, function)
    84     pairs.
    85    
    86     ``list_of_pairs`` is a list of pairs (I, fcn), where
    87     fcn is a Sage function (such as a polynomial over RR, or functions
    88     using the lambda notation), and I is an interval such as I = (1,3).
    89     Two consecutive intervals must share a common endpoint.
    90    
    91     If the optional ``var`` is specified, then any symbolic expressions
    92     in the list will be converted to symbolic functions using
    93     ``fcn.function(var)``.  (This says which variable is considered to
    94     be "piecewise".)
    9570
    96     We assume that these definitions are consistent (ie, no checking is
    97     done).
    98    
    99     EXAMPLES::
    100    
    101         sage: f1(x) = -1
    102         sage: f2(x) = 2
    103         sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
    104         sage: f(1)
    105         -1
    106         sage: f(3)
    107         2
    108         sage: f = Piecewise([[(0,1),x], [(1,2),x^2]], x); f
    109         Piecewise defined function with 2 parts, [[(0, 1), x |--> x], [(1, 2), x |--> x^2]]
    110         sage: f(0.9)
    111         0.900000000000000
    112         sage: f(1.1)
    113         1.21000000000000
    114     """
    115     return PiecewisePolynomial(list_of_pairs, var=var)
    11671
    117 Piecewise = piecewise
     72class PiecewiseFunction(BuiltinFunction):
    11873
    119 class PiecewisePolynomial:
    120     """
    121     Returns a piecewise function from a list of (interval, function)
    122     pairs.
    123    
    124     EXAMPLES::
    125    
    126         sage: f1(x) = -1
    127         sage: f2(x) = 2
    128         sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
    129         sage: f(1)
    130         -1
    131         sage: f(3)
    132         2
    133     """
    134     def __init__(self, list_of_pairs, var=None):
    135         r"""
    136         ``list_of_pairs`` is a list of pairs (I, fcn), where
    137         fcn is a Sage function (such as a polynomial over RR, or functions
    138         using the lambda notation), and I is an interval such as I = (1,3).
    139         Two consecutive intervals must share a common endpoint.
    140        
    141         If the optional ``var`` is specified, then any symbolic
    142         expressions in the list will be converted to symbolic
    143         functions using ``fcn.function(var)``.  (This says which
    144         variable is considered to be "piecewise".)
    145 
    146         We assume that these definitions are consistent (ie, no checking is
    147         done).
     74    def __init__(self):
     75        """
     76        Piecewise function
    14877
    14978        EXAMPLES::
    15079
    151             sage: f1(x) = 1
    152             sage: f2(x) = 1 - x
    153             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    154             sage: f.list()
    155             [[(0, 1), x |--> 1], [(1, 2), x |--> -x + 1]]
    156             sage: f.length()
    157             2
     80            sage: var('x, y')
     81            (x, y)
     82            sage: f = piecewise([((0,1), x^2*y), ([-1,0], -x*y^2)], var=x);  f
     83            piecewise(x|-->x^2*y on (0, 1), x|-->-x*y^2 on [-1, 0]; x)
     84            sage: f(1/2)
     85            1/4*y
     86            sage: f(-1/2)
     87            1/2*y^2
    15888        """
    159         self._length = len(list_of_pairs)
    160         self._intervals = [x[0] for x in list_of_pairs]
    161         functions = [x[1] for x in list_of_pairs]
    162         if var is not None:
    163             for i in range(len(functions)):
    164                 if is_Expression(functions[i]):
    165                     functions[i] = functions[i].function(var)
    166         self._functions = functions
    167         # We regenerate self._list in case self._functions was modified
    168         # above.  This also protects us in case somebody mutates a list
    169         # after they use it as an argument to piecewise().
    170         self._list = [[self._intervals[i], self._functions[i]] for i in range(self._length)]
     89        BuiltinFunction.__init__(self, "piecewise",
     90                                 latex_name="piecewise",
     91                                 conversions=dict(), nargs=2)
     92
     93    def __call__(self, function_pieces, **kwds):
     94        r"""
     95        Piecewise functions
     96
     97        INPUT:
     98   
     99        - ``function_pieces`` -- a list of pairs consisting of a
     100          domain and a symbolic function.
    171101 
    172     def list(self):
     102        - ``var=x`` -- a symbolic variable or ``None`` (default). The
     103        real variable in which the function is piecewise in.
     104
     105        OUTPUT:
     106
     107        A piecewise-defined function. A ``ValueError`` will be raised
     108        if the domains of the pieces are not pairwise disjoint.
     109   
     110        EXAMPLES::
     111       
     112            sage: my_abs = piecewise([((-1, 0), -x), ([0, 1], x)], var=x);  my_abs
     113            piecewise(x|-->-x on (-1, 0), x|-->x on [0, 1]; x)
     114            sage: [ my_abs(i/5) for i in range(-4, 5)]
     115            [4/5, 3/5, 2/5, 1/5, 0, 1/5, 2/5, 3/5, 4/5]
     116
     117        TESTS::
     118
     119            sage: piecewise([([-1, 0], -x), ([0, 1], x)], var=x)
     120            Traceback (most recent call last):
     121            ...
     122            ValueError: domains must be pairwise disjoint
     123
     124            sage: step = piecewise([((-1, 0), -1), ([0, 0], 0), ((0, 1), 1)], var=x);  step
     125            piecewise(x|-->-1 on (-1, 0), x|-->0 on {0}, x|-->1 on (0, 1); x)
     126            sage: step(-1/2), step(0), step(1/2)
     127            (-1, 0, 1)
    173128        """
    174         Returns the pieces of this function as a list of functions.
     129        #print 'pf_call', function_pieces, kwds
     130        var = kwds.pop('var', None)
     131        parameters = []
     132        domain_list = []
     133        for piece in function_pieces:
     134            domain, function = piece
     135            if not isinstance(domain, RealSet):
     136                domain = RealSet(domain)
     137            if domain.is_empty():
     138                continue
     139            function = SR(function)
     140            if var is None and len(function.variables()) > 0:
     141                var = function.variables()[0]
     142            parameters.append((domain, function))
     143            domain_list.append(domain)
     144        if not RealSet.are_pairwise_disjoint(*domain_list):
     145            raise ValueError('domains must be pairwise disjoint')
     146        if var is None:
     147            var = self.default_variable()
     148        parameters = SR._force_pyobject(tuple(parameters), recursive=False)
     149        return BuiltinFunction.__call__(self, parameters, var, **kwds)
     150
     151    def _print_(self, parameters, variable):
     152        """
     153        Return a string representation
     154       
     155        OUTPUT:
     156       
     157        String.
    175158       
    176159        EXAMPLES::
    177160
    178             sage: f1(x) = 1
    179             sage: f2(x) = 1 - x
    180             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    181             sage: f.list()
    182             [[(0, 1), x |--> 1], [(1, 2), x |--> -x + 1]]
     161       
    183162        """
    184         return self._list
    185  
    186     def length(self):
     163        s = 'piecewise('
     164        args = []
     165        for domain, func in parameters:
     166            args.append('{0}|-->{1} on {2}'.format(str(variable), str(func), str(domain)))
     167        s += ', '.join(args) + '; {0})'.format(str(variable))
     168        return s
     169
     170    def _subs_(self, subs_map, options, parameters, x):
    187171        """
    188         Returns the number of pieces of this function.
     172        Callback from Pynac `subs()`
    189173       
    190         EXAMPLES::
     174        EXAMPLES:
     175
     176        If the substitution changes the piecewise variable, it must
     177        evaluate to a number so that we know which component we are
     178        on::
     179
     180            sage: p = piecewise([((-2, 0), -x), ([0, 2], x)], var=x)
     181            sage: p.subs(x=-1)
     182            1
     183            sage: (10+p).subs(x=-1)
     184            11
     185   
     186        Auxiliary variables can be substituted arbitrarily::
     187
     188            sage: var('x,y')
     189            (x, y)
     190            sage: p = piecewise([((-2, 0), -x^y), ([0, 2], x-y)], var=x);  p
     191            piecewise(x|-->-x^y on (-2, 0), x|-->x - y on [0, 2]; x)
     192            sage: p.subs(y=sin(y))
     193            piecewise(x|-->-x^sin(y) on (-2, 0), x|-->x - sin(y) on [0, 2]; x)
     194        """
     195        point = subs_map.apply_to(x, 0)
     196        # print 'point =', point
     197        if point == x:
     198            # substitution only in auxiliary variables
     199            new_params = []
     200            for domain, func in parameters:
     201                new_params.append((domain, subs_map.apply_to(func, 0)))
     202            return piecewise(new_params, var=x)
     203        if not (point.is_numeric() and point.is_real()):
     204            raise ValueError('substituting the piecewise variable must result in real number')
    191205       
    192             sage: f1(x) = 1
    193             sage: f2(x) = 1 - x
    194             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    195             sage: f.length()
    196             2
     206        for domain, func in parameters:
     207            if domain.contains(point):
     208                return subs_map.apply_to(func, 0)
     209        raise ValueError('point is not in the domain')
     210
     211    @staticmethod
     212    def in_operands(ex):
    197213        """
    198         return self._length
    199  
    200     def __repr__(self):
    201         """
    202         EXAMPLES::
    203        
    204             sage: f1(x) = 1
    205             sage: f2(x) = 1 - x
    206             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]]); f
    207             Piecewise defined function with 2 parts, [[(0, 1), x |--> 1], [(1, 2), x |--> -x + 1]]
    208         """
    209         return 'Piecewise defined function with %s parts, %s'%(
    210             self.length(),self.list())
    211  
    212     def _latex_(self):
    213         r"""
    214         EXAMPLES::
    215        
    216             sage: f1(x) = 1
    217             sage: f2(x) = 1 - x
    218             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    219             sage: latex(f)
    220             \begin{cases}
    221             x \ {\mapsto}\ 1 &\text{on $(0, 1)$}\cr
    222             x \ {\mapsto}\ -x + 1 &\text{on $(1, 2)$}\cr
    223             \end{cases}
    224        
    225         ::
    226        
    227             sage: f(x) = sin(x*pi/2)
    228             sage: g(x) = 1-(x-1)^2
    229             sage: h(x) = -x
    230             sage: P = Piecewise([[(0,1), f], [(1,3),g], [(3,5), h]])
    231             sage: latex(P)
    232             \begin{cases}
    233             x \ {\mapsto}\ \sin\left(\frac{1}{2} \, \pi x\right) &\text{on $(0, 1)$}\cr
    234             x \ {\mapsto}\ -{\left(x - 1\right)}^{2} + 1 &\text{on $(1, 3)$}\cr
    235             x \ {\mapsto}\ -x &\text{on $(3, 5)$}\cr
    236             \end{cases}
    237         """
    238         from sage.misc.latex import latex
    239         tex = ['\\begin{cases}\n']
    240         for (left, right), f in self.list():
    241             tex.append('%s &\\text{on $(%s, %s)$}\\cr\n' % (latex(f), left, right))
    242         tex.append(r'\end{cases}')
    243         return ''.join(tex)
     214        Return whether a symbolic expression contains a piecewise
     215        function as operand
    244216
    245     def intervals(self):
    246         """
    247         A piecewise non-polynomial example.
    248        
    249         EXAMPLES::
    250        
    251             sage: f1(x) = 1
    252             sage: f2(x) = 1-x
    253             sage: f3(x) = exp(x)
    254             sage: f4(x) = sin(2*x)
    255             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    256             sage: f.intervals()
    257             [(0, 1), (1, 2), (2, 3), (3, 10)]
    258         """
    259         return self._intervals
    260  
    261     def domain(self):
    262         """
    263         Returns the domain of the function.
    264        
    265         EXAMPLES::
    266        
    267             sage: f1(x) = 1
    268             sage: f2(x) = 1-x
    269             sage: f3(x) = exp(x)
    270             sage: f4(x) = sin(2*x)
    271             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    272             sage: f.domain()
    273             (0, 10)
    274         """
    275         endpoints = sum(self.intervals(), ())
    276         return (min(endpoints), max(endpoints))
    277 
    278     def functions(self):
    279         """
    280         Returns the list of functions (the "pieces").
    281        
    282         EXAMPLES::
    283        
    284             sage: f1(x) = 1
    285             sage: f2(x) = 1-x
    286             sage: f3(x) = exp(x)
    287             sage: f4(x) = sin(2*x)
    288             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    289             sage: f.functions()
    290             [x |--> 1, x |--> -x + 1, x |--> e^x, x |--> sin(2*x)]
    291         """
    292         return self._functions
    293        
    294     def extend_by_zero_to(self,xmin=-1000,xmax=1000):
    295         """
    296         This function simply returns the piecewise defined function which
    297         is extended by 0 so it is defined on all of (xmin,xmax). This is
    298         needed to add two piecewise functions in a reasonable way.
    299        
    300         EXAMPLES::
    301        
    302             sage: f1(x) = 1
    303             sage: f2(x) = 1 - x
    304             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    305             sage: f.extend_by_zero_to(-1, 3)
    306             Piecewise defined function with 4 parts, [[(-1, 0), 0], [(0, 1), x |--> 1], [(1, 2), x |--> -x + 1], [(2, 3), 0]]
    307         """
    308         zero = QQ['x'](0)
    309         list_of_pairs = self.list()
    310         a, b = self.domain()
    311         if xmin < a:
    312             list_of_pairs = [[(xmin, a), zero]] + list_of_pairs
    313         if xmax > b:
    314             list_of_pairs = list_of_pairs + [[(b, xmax), zero]]
    315         return Piecewise(list_of_pairs)
    316 
    317     def unextend(self):
    318         """
    319         This removes any parts in the front or back of the function which
    320         is zero (the inverse to extend_by_zero_to).
    321        
    322         EXAMPLES::
    323        
    324             sage: R.<x> = QQ[]
    325             sage: f = Piecewise([[(-3,-1),1+2+x],[(-1,1),1-x^2]])
    326             sage: e = f.extend_by_zero_to(-10,10); e
    327             Piecewise defined function with 4 parts, [[(-10, -3), 0], [(-3, -1), x + 3], [(-1, 1), -x^2 + 1], [(1, 10), 0]]
    328             sage: d = e.unextend(); d
    329             Piecewise defined function with 2 parts, [[(-3, -1), x + 3], [(-1, 1), -x^2 + 1]]
    330             sage: d==f
    331             True
    332         """
    333         list_of_pairs = self.list()
    334         funcs = self.functions()
    335         if funcs[0] == 0:
    336             list_of_pairs = list_of_pairs[1:]
    337         if funcs[-1] == 0:
    338             list_of_pairs = list_of_pairs[:-1]
    339         return Piecewise(list_of_pairs)
    340 
    341     def _riemann_sum_helper(self, N, func, initial=0):
    342         """
    343         A helper function for computing Riemann sums.
    344        
    345217        INPUT:
    346218       
    347        
    348         -  ``N`` - the number of subdivisions
    349        
    350         -  ``func`` - a function to apply to the endpoints of
    351            each subdivision
    352        
    353         -  ``initial`` - the starting value
    354        
    355        
    356         EXAMPLES::
    357        
    358             sage: f1(x) = x^2                   ## example 1
    359             sage: f2(x) = 5-x^2
    360             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    361             sage: f._riemann_sum_helper(6, lambda x0, x1: (x1-x0)*f(x1))
    362             19/6
    363         """
    364         a,b = self.domain()
    365         rsum = initial
    366         h = (b-a)/N
    367         for i in range(N):
    368             x0 = a+i*h
    369             x1 = a+(i+1)*h
    370             rsum += func(x0, x1)
    371         return rsum
     219        - ``ex`` -- a symbolic expression.
    372220
    373     def riemann_sum_integral_approximation(self,N,mode=None):
    374         """
    375         Returns the piecewise line function defined by the Riemann sums in
    376         numerical integration based on a subdivision into N subintervals.
    377        
    378         Set mode="midpoint" for the height of the rectangles to be
    379         determined by the midpoint of the subinterval; set mode="right" for
    380         the height of the rectangles to be determined by the right-hand
    381         endpoint of the subinterval; the default is mode="left" (the height
    382         of the rectangles to be determined by the left-hand endpoint of
    383         the subinterval).
    384        
    385         EXAMPLES::
    386        
    387             sage: f1(x) = x^2                   ## example 1
    388             sage: f2(x) = 5-x^2
    389             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    390             sage: f.riemann_sum_integral_approximation(6)
    391             17/6
    392             sage: f.riemann_sum_integral_approximation(6,mode="right")
    393             19/6
    394             sage: f.riemann_sum_integral_approximation(6,mode="midpoint")
    395             3
    396             sage: f.integral(definite=True)
    397             3
    398         """
    399         if mode is None:
    400             return self._riemann_sum_helper(N, lambda x0, x1: (x1-x0)*self(x0))
    401         elif mode == "right":
    402             return self._riemann_sum_helper(N, lambda x0, x1: (x1-x0)*self(x1))
    403         elif mode == "midpoint":
    404             return self._riemann_sum_helper(N, lambda x0, x1: (x1-x0)*self((x0+x1)/2))
    405         else:
    406             raise ValueError, "invalid mode"
     221        OUTPUT:
    407222
    408     def riemann_sum(self,N,mode=None):
    409         """
    410         Returns the piecewise line function defined by the Riemann sums in
    411         numerical integration based on a subdivision into N subintervals.
    412         Set mode="midpoint" for the height of the rectangles to be
    413         determined by the midpoint of the subinterval; set mode="right" for
    414         the height of the rectangles to be determined by the right-hand
    415         endpoint of the subinterval; the default is mode="left" (the height
    416         of the rectangles to be determined by the left-hand endpoint of
    417         the subinterval).
    418        
    419         EXAMPLES::
    420        
    421             sage: f1(x) = x^2
    422             sage: f2(x) = 5-x^2
    423             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    424             sage: f.riemann_sum(6,mode="midpoint")
    425             Piecewise defined function with 6 parts, [[(0, 1/3), 1/36], [(1/3, 2/3), 1/4], [(2/3, 1), 25/36], [(1, 4/3), 131/36], [(4/3, 5/3), 11/4], [(5/3, 2), 59/36]]
    426        
    427         ::
    428        
    429             sage: f = Piecewise([[(-1,1),(1-x^2).function(x)]])
    430             sage: rsf = f.riemann_sum(7)
    431             sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
    432             sage: Q = rsf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
    433             sage: L = add([line([[a,0],[a,f(x=a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in rsf.list()])
    434             sage: P + Q + L
    435        
    436         ::
    437        
    438             sage: f = Piecewise([[(-1,1),(1/2+x-x^3)]], x) ## example 3
    439             sage: rsf = f.riemann_sum(8)
    440             sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
    441             sage: Q = rsf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
    442             sage: L = add([line([[a,0],[a,f(x=a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in rsf.list()])
    443             sage: P + Q + L
    444         """
    445         if mode is None:
    446             rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self(x0))]],
    447                                             initial=[])
    448         elif mode == "right":
    449             rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self(x1))]],
    450                                             initial=[])
    451         elif mode == "midpoint":
    452             rsum = self._riemann_sum_helper(N, lambda x0,x1: [[(x0,x1),SR(self((x0+x1)/2))]],
    453                                             initial=[])
    454         else:
    455             raise ValueError, "invalid mode"
    456         return Piecewise(rsum)
    457 
    458     def trapezoid(self,N):
    459         """
    460         Returns the piecewise line function defined by the trapezoid rule
    461         for numerical integration based on a subdivision into N
    462         subintervals.
    463        
    464         EXAMPLES::
    465        
    466             sage: R.<x> = QQ[]
    467             sage: f1 = x^2
    468             sage: f2 = 5-x^2
    469             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    470             sage: f.trapezoid(4)
    471             Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*x], [(1/2, 1), 9/2*x - 2], [(1, 3/2), 1/2*x + 2], [(3/2, 2), -7/2*x + 8]]
    472        
    473         ::
    474        
    475             sage: R.<x> = QQ[]
    476             sage: f = Piecewise([[(-1,1),1-x^2]])
    477             sage: tf = f.trapezoid(4)
    478             sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
    479             sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
    480             sage: L = add([line([[a,0],[a,f(a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in tf.list()])
    481             sage: P+Q+L
    482        
    483         ::
    484        
    485             sage: R.<x> = QQ[]
    486             sage: f = Piecewise([[(-1,1),1/2+x-x^3]]) ## example 3
    487             sage: tf = f.trapezoid(6)
    488             sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
    489             sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
    490             sage: L = add([line([[a,0],[a,f(a)]],rgbcolor=(0.7,0.6,0.6)) for (a,b),f in tf.list()])
    491             sage: P+Q+L
    492 
    493         TESTS:
    494 
    495         Use variables other than x (:trac:`13836`)::
    496 
    497             sage: R.<y> = QQ[]
    498             sage: f1 = y^2
    499             sage: f2 = 5-y^2
    500             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    501             sage: f.trapezoid(4)
    502             Piecewise defined function with 4 parts, [[(0, 1/2), 1/2*y], [(1/2, 1), 9/2*y - 2], [(1, 3/2), 1/2*y + 2], [(3/2, 2), -7/2*y + 8]]
    503 
    504         """
    505         x = QQ[self.default_variable()].gen()
    506         def f(x0, x1):
    507             f0, f1 = self(x0), self(x1)
    508             return [[(x0,x1),f0+(f1-f0)*(x1-x0)**(-1)*(x-x0)]]
    509         rsum = self._riemann_sum_helper(N, f, initial=[])
    510         return Piecewise(rsum)
    511 
    512     def trapezoid_integral_approximation(self,N):
    513         """
    514         Returns the approximation given by the trapezoid rule for numerical
    515         integration based on a subdivision into N subintervals.
    516        
    517         EXAMPLES::
    518        
    519             sage: f1(x) = x^2                      ## example 1
    520             sage: f2(x) = 1-(1-x)^2
    521             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    522             sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
    523             sage: tf = f.trapezoid(6)
    524             sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
    525             sage: ta = f.trapezoid_integral_approximation(6)
    526             sage: t = text('trapezoid approximation = %s'%ta, (1.5, 0.25))
    527             sage: a = f.integral(definite=True)
    528             sage: tt = text('area under curve = %s'%a, (1.5, -0.5))
    529             sage: P + Q + t + tt
    530        
    531         ::
    532        
    533             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])  ## example 2
    534             sage: tf = f.trapezoid(4)
    535             sage: ta = f.trapezoid_integral_approximation(4)
    536             sage: Q = tf.plot(rgbcolor=(0.7,0.6,0.6), plot_points=40)
    537             sage: t = text('trapezoid approximation = %s'%ta, (1.5, 0.25))
    538             sage: a = f.integral(definite=True)
    539             sage: tt = text('area under curve = %s'%a, (1.5, -0.5))
    540             sage: P+Q+t+tt
    541         """
    542         def f(x0, x1):
    543             f0, f1 = self(x0), self(x1)
    544             return ((f1+f0)/2)*(x1-x0)
    545         return self._riemann_sum_helper(N, f)
    546 
    547     def critical_points(self):
    548         """
    549         Return the critical points of this piecewise function.
    550        
    551         .. warning::
    552 
    553            Uses maxima, which prints the warning to use results with
    554            caution. Only works for piecewise functions whose parts are
    555            polynomials with real critical not occurring on the
    556            interval endpoints.
    557        
    558         EXAMPLES::
    559        
    560             sage: R.<x> = QQ[]
    561             sage: f1 = x^0
    562             sage: f2 = 10*x - x^2
    563             sage: f3 = 3*x^4 - 156*x^3 + 3036*x^2 - 26208*x
    564             sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]])
    565             sage: expected = [5, 12, 13, 14]
    566             sage: all(abs(e-a) < 0.001 for e,a in zip(expected, f.critical_points()))
    567             True
    568 
    569         TESTS:
    570 
    571         Use variables other than x (:trac:`13836`)::
    572 
    573             sage: R.<y> = QQ[]
    574             sage: f1 = y^0
    575             sage: f2 = 10*y - y^2
    576             sage: f3 = 3*y^4 - 156*y^3 + 3036*y^2 - 26208*y
    577             sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]])
    578             sage: expected = [5, 12, 13, 14]
    579             sage: all(abs(e-a) < 0.001 for e,a in zip(expected, f.critical_points()))
    580             True
    581         """
    582         from sage.calculus.calculus import maxima
    583         x = QQ[self.default_variable()].gen()
    584         crit_pts = []
    585         for (a,b), f in self.list():
    586             for root in maxima.allroots(SR(f).diff(x)==0):
    587                 root = float(root.rhs())
    588                 if a < root < b:
    589                     crit_pts.append(root)
    590         return crit_pts
    591        
    592     def base_ring(self):
    593         """
    594         Returns the base ring of the function pieces.   This
    595         is useful when this class is extended.
     223        Boolean
    596224
    597225        EXAMPLES::
    598226       
    599             sage: f1(x) = 1
    600             sage: f2(x) = 1-x
    601             sage: f3(x) = x^2-5
    602             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3]])
    603             sage: base_ring(f)
    604             Symbolic Ring
     227            sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]);  f
     228            piecewise(x|-->sin(x) on {0}, x|-->cos(x) on (0, 2); x)
     229            sage: piecewise.in_operands(f)
     230            True
     231            sage: piecewise.in_operands(1+sin(f))
     232            True
     233            sage: piecewise.in_operands(1+sin(0*f))
     234            False
     235        """
     236        def is_piecewise(ex):
     237            result = ex.operator() is piecewise
     238            for op in ex.operands():
     239                result = result or is_piecewise(op)
     240            return result
     241        return is_piecewise(ex)
    605242
    606         ::
    607243
    608             sage: R.<x> = QQ[]
    609             sage: f1 = x^0
    610             sage: f2 = 10*x - x^2
    611             sage: f3 = 3*x^4 - 156*x^3 + 3036*x^2 - 26208*x
    612             sage: f = Piecewise([[(0,3),f1],[(3,10),f2],[(10,20),f3]])
    613             sage: f.base_ring()
    614             Rational Field
     244    @staticmethod
     245    def simplify():
    615246        """
    616         return (self.functions()[0]).base_ring()
    617 
    618     def end_points(self):
    619         """
    620         Returns a list of all interval endpoints for this function.
    621        
    622         EXAMPLES::
    623        
    624             sage: f1(x) = 1
    625             sage: f2(x) = 1-x
    626             sage: f3(x) = x^2-5
    627             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3]])
    628             sage: f.end_points()
    629             [0, 1, 2, 3]
    630         """
    631         intervals = self.intervals()
    632         return [ intervals[0][0] ] + [b for a,b in intervals]
    633 
    634     def __call__(self,x0):
    635         """
    636         Evaluates self at x0. Returns the average value of the jump if x0
    637         is an interior endpoint of one of the intervals of self and the
    638         usual value otherwise.
    639        
    640         EXAMPLES::
    641        
    642             sage: f1(x) = 1
    643             sage: f2(x) = 1-x
    644             sage: f3(x) = exp(x)
    645             sage: f4(x) = sin(2*x)
    646             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    647             sage: f(0.5)
    648             1
    649             sage: f(5/2)
    650             e^(5/2)
    651             sage: f(5/2).n()
    652             12.1824939607035
    653             sage: f(1)
    654             1/2
    655         """
    656         #x0 = QQ(x0) ## does not allow for evaluation at pi
    657         n = self.length()
    658         endpts = self.end_points()
    659         for i in range(1,n):
    660             if x0 == endpts[i]:
    661                 return (self.functions()[i-1](x0) + self.functions()[i](x0))/2
    662         if x0 == endpts[0]:
    663             return self.functions()[0](x0)
    664         if x0 == endpts[n]:
    665             return self.functions()[n-1](x0)
    666         for i in range(n):
    667             if endpts[i] < x0 < endpts[i+1]:
    668                 return self.functions()[i](x0)
    669         raise ValueError,"Value not defined outside of domain."
    670 
    671     def which_function(self,x0):
    672         """
    673         Returns the function piece used to evaluate self at x0.
    674        
    675         EXAMPLES::
    676        
    677             sage: f1(z) = z
    678             sage: f2(x) = 1-x
    679             sage: f3(y) = exp(y)
    680             sage: f4(t) = sin(2*t)
    681             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    682             sage: f.which_function(3/2)
    683             x |--> -x + 1
    684         """
    685         for (a,b), f in self.list():
    686             if a <= x0 <= b:
    687                 return f
    688         raise ValueError,"Function not defined outside of domain."
    689 
    690     def default_variable(self):
    691         r"""
    692         Return the default variable. The default variable is defined as the
    693         first variable in the first piece that has a variable. If no pieces have
    694         a variable (each piece is a constant value), `x` is returned.
    695 
    696         The result is cached.
    697 
    698         AUTHOR: Paul Butler
    699 
    700         EXAMPLES::
    701        
    702             sage: f1(x) = 1
    703             sage: f2(x) = 5*x
    704             sage: p = Piecewise([[(0,1),f1],[(1,4),f2]])
    705             sage: p.default_variable()
    706             x
    707 
    708             sage: f1 = 3*var('y')
    709             sage: p = Piecewise([[(0,1),4],[(1,4),f1]])
    710             sage: p.default_variable()
    711             y
    712 
    713         """
    714         try:
    715             return self.__default_variable
    716         except AttributeError:
    717             pass
    718         for _, fun in self._list:
    719             try:
    720                 fun = SR(fun)
    721                 if fun.variables():
    722                     v = fun.variables()[0]
    723                     self.__default_variable = v
    724                     return v
    725             except TypeError:
    726                 # pass if fun is lambda function
    727                 pass
    728         # default to x
    729         v = var('x')
    730         self.__default_value = v
    731         return v
    732 
    733     def integral(self, x=None, a=None, b=None, definite=False):
    734         r"""
    735         By default, returns the indefinite integral of the function.
    736         If definite=True is given, returns the definite integral.
    737 
    738         AUTHOR:
    739 
    740         - Paul Butler
    741 
    742         EXAMPLES::
    743 
    744             sage: f1(x) = 1-x
    745             sage: f = Piecewise([[(0,1),1],[(1,2),f1]])
    746             sage: f.integral(definite=True)
    747             1/2
    748        
    749         ::
    750        
    751             sage: f1(x) = -1
    752             sage: f2(x) = 2
    753             sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
    754             sage: f.integral(definite=True)
    755             1/2*pi
    756            
    757             sage: f1(x) = 2
    758             sage: f2(x) = 3 - x
    759             sage: f = Piecewise([[(-2, 0), f1], [(0, 3), f2]])
    760             sage: f.integral()
    761             Piecewise defined function with 2 parts, [[(-2, 0), x |--> 2*x + 4], [(0, 3), x |--> -1/2*x^2 + 3*x + 4]]
    762 
    763             sage: f1(y) = -1
    764             sage: f2(y) = y + 3
    765             sage: f3(y) = -y - 1
    766             sage: f4(y) = y^2 - 1
    767             sage: f5(y) = 3
    768             sage: f = Piecewise([[(-4,-3),f1],[(-3,-2),f2],[(-2,0),f3],[(0,2),f4],[(2,3),f5]])
    769             sage: F = f.integral(y)
    770             sage: F
    771             Piecewise defined function with 5 parts, [[(-4, -3), y |--> -y - 4], [(-3, -2), y |--> 1/2*y^2 + 3*y + 7/2], [(-2, 0), y |--> -1/2*y^2 - y - 1/2], [(0, 2), y |--> 1/3*y^3 - y - 1/2], [(2, 3), y |--> 3*y - 35/6]]
    772            
    773         Ensure results are consistent with FTC::
    774 
    775             sage: F(-3) - F(-4)
    776             -1
    777             sage: F(-1) - F(-3)
    778             1
    779             sage: F(2) - F(0)
    780             2/3
    781             sage: f.integral(y, 0, 2)
    782             2/3
    783             sage: F(3) - F(-4)
    784             19/6
    785             sage: f.integral(y, -4, 3)
    786             19/6
    787             sage: f.integral(definite=True)
    788             19/6
    789            
    790         ::
    791 
    792             sage: f1(y) = (y+3)^2
    793             sage: f2(y) = y+3
    794             sage: f3(y) = 3
    795             sage: f = Piecewise([[(-infinity, -3), f1], [(-3, 0), f2], [(0, infinity), f3]])
    796             sage: f.integral()
    797             Piecewise defined function with 3 parts, [[(-Infinity, -3), y |--> 1/3*y^3 + 3*y^2 + 9*y + 9], [(-3, 0), y |--> 1/2*y^2 + 3*y + 9/2], [(0, +Infinity), y |--> 3*y + 9/2]]
    798 
    799         ::
    800 
    801             sage: f1(x) = e^(-abs(x))
    802             sage: f = Piecewise([[(-infinity, infinity), f1]])
    803             sage: f.integral(definite=True)
    804             2
    805             sage: f.integral()
    806             Piecewise defined function with 1 parts, [[(-Infinity, +Infinity), x |--> -1/2*((sgn(x) - 1)*e^(2*x) - 2*e^x*sgn(x) + sgn(x) + 1)*e^(-x) - 1]]
    807 
    808         ::
    809 
    810             sage: f = Piecewise([((0, 5), cos(x))])
    811             sage: f.integral()
    812             Piecewise defined function with 1 parts, [[(0, 5), x |--> sin(x)]]
    813 
    814 
    815         TESTS:
    816 
    817         Verify that piecewise integrals of zero work (trac #10841)::
    818 
    819             sage: f0(x) = 0
    820             sage: f = Piecewise([[(0,1),f0]])
    821             sage: f.integral(x,0,1)
    822             0
    823             sage: f = Piecewise([[(0,1), 0]])
    824             sage: f.integral(x,0,1)
    825             0
    826             sage: f = Piecewise([[(0,1), SR(0)]])
    827             sage: f.integral(x,0,1)
    828             0
    829 
    830         """
    831         if a != None and b != None:
    832             F = self.integral(x)
    833             return F(b) - F(a)
    834 
    835         if a != None or b != None:
    836             raise TypeError, 'only one endpoint given'
    837 
    838         area = 0 # cumulative definite integral of parts to the left of the current interval
    839         integrand_pieces = self.list()
    840         integrand_pieces.sort()
    841         new_pieces = []
    842 
    843         if x == None:
    844             x = self.default_variable()
    845        
    846         # The integral is computed by iterating over the pieces in order.
    847         # The definite integral for each piece is calculated and accumulated in `area`.
    848         # Thus at any time, `area` represents the definite integral of all the pieces
    849         # encountered so far. The indefinite integral of each piece is also calculated,
    850         # and the `area` before each piece is added to the piece.
    851         #
    852         # If a definite integral is requested, `area` is returned.
    853         # Otherwise, a piecewise function is constructed from the indefinite integrals
    854         # and returned.
    855         #
    856         # An exception is made if integral is called on a piecewise function
    857         # that starts at -infinity. In this case, we do not try to calculate the
    858         # definite integral of the first piece, and the value of `area` remains 0
    859         # after the first piece.
    860 
    861         for (start, end), fun in integrand_pieces:
    862             fun = SR(fun)
    863             if start == -infinity and not definite:
    864                 fun_integrated = fun.integral(x, end, x)
    865             else:
    866                 try:
    867                     assume(start < x)
    868                 except ValueError: # Assumption is redundant
    869                     pass
    870                 fun_integrated = fun.integral(x, start, x) + area
    871                 forget(start < x)
    872                 if definite or end != infinity:
    873                     area += fun.integral(x, start, end)
    874             new_pieces.append([(start, end), SR(fun_integrated).function(x)])
    875 
    876         if definite:
    877             return SR(area)
    878         else:
    879             return Piecewise(new_pieces)
    880 
    881     def convolution(self, other):
    882         """
    883         Returns the convolution function,
    884         `f*g(t)=\int_{-\infty}^\infty f(u)g(t-u)du`, for compactly
    885         supported `f,g`.
    886        
    887         EXAMPLES::
    888        
    889             sage: x = PolynomialRing(QQ,'x').gen()
    890             sage: f = Piecewise([[(0,1),1*x^0]])  ## example 0
    891             sage: g = f.convolution(f)
    892             sage: h = f.convolution(g)
    893             sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1));
    894             sage: # Type show(P+Q+R) to view
    895             sage: f = Piecewise([[(0,1),1*x^0],[(1,2),2*x^0],[(2,3),1*x^0]])  ## example 1
    896             sage: g = f.convolution(f)
    897             sage: h = f.convolution(g)
    898             sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1));
    899             sage: # Type show(P+Q+R) to view
    900             sage: f = Piecewise([[(-1,1),1]])                             ## example 2
    901             sage: g = Piecewise([[(0,3),x]])
    902             sage: f.convolution(g)
    903             Piecewise defined function with 3 parts, [[(-1, 1), 0], [(1, 2), -3/2*x], [(2, 4), -3/2*x]]
    904             sage: g = Piecewise([[(0,3),1*x^0],[(3,4),2*x^0]])
    905             sage: f.convolution(g)
    906             Piecewise defined function with 5 parts, [[(-1, 1), x + 1], [(1, 2), 3], [(2, 3), x], [(3, 4), -x + 8], [(4, 5), -2*x + 10]]
    907         """
    908         f = self
    909         g = other
    910         M = min(min(f.end_points()),min(g.end_points()))
    911         N = max(max(f.end_points()),max(g.end_points()))
    912         R2 = PolynomialRing(QQ,2,names=["tt","uu"])
    913         tt,uu = R2.gens()
    914         conv = 0
    915         f0 = f.functions()[0]
    916         g0 = g.functions()[0]
    917         R1 = f0.parent()
    918         xx = R1.gen()
    919         var = repr(xx)
    920         if len(f.intervals())==1 and len(g.intervals())==1:
    921             f = f.unextend()
    922             g = g.unextend()
    923             a1 = f.intervals()[0][0]
    924             a2 = f.intervals()[0][1]
    925             b1 = g.intervals()[0][0]
    926             b2 = g.intervals()[0][1]
    927             i1 = repr(f0).replace(var,repr(uu))
    928             i2 = repr(g0).replace(var,"("+repr(tt-uu)+")")
    929             cmd1 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, a1,    tt-b1)    ## if a1+b1 < tt < a2+b1
    930             cmd2 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, tt-b2, tt-b1)    ## if a1+b2 < tt < a2+b1
    931             cmd3 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, tt-b2, a2)       ## if a1+b2 < tt < a2+b2
    932             cmd4 = "integrate((%s)*(%s),%s,%s,%s)"%(i1,i2, uu, a1, a2)          ## if a2+b1 < tt < a1+b2
    933             conv1 = maxima.eval(cmd1)
    934             conv2 = maxima.eval(cmd2)
    935             conv3 = maxima.eval(cmd3)
    936             conv4 = maxima.eval(cmd4)
    937             # this is a very, very, very ugly hack
    938             x = PolynomialRing(QQ,'x').gen()
    939             fg1 = sage_eval(conv1.replace("tt",var), {'x':x}) ## should be = R2(conv1)
    940             fg2 = sage_eval(conv2.replace("tt",var), {'x':x}) ## should be = R2(conv2)
    941             fg3 = sage_eval(conv3.replace("tt",var), {'x':x}) ## should be = R2(conv3)
    942             fg4 = sage_eval(conv4.replace("tt",var), {'x':x}) ## should be = R2(conv4)
    943             if a1-b1<a2-b2:
    944                 if a2+b1!=a1+b2:
    945                     h = Piecewise([[(a1+b1,a1+b2),fg1],[(a1+b2,a2+b1),fg4],[(a2+b1,a2+b2),fg3]])
    946                 else:
    947                     h = Piecewise([[(a1+b1,a1+b2),fg1],[(a1+b2,a2+b2),fg3]])
    948             else:
    949                 if a1+b2!=a2+b1:
    950                     h = Piecewise([[(a1+b1,a2+b1),fg1],[(a2+b1,a1+b2),fg2],[(a1+b2,a2+b2),fg3]])
    951                 else:
    952                     h = Piecewise([[(a1+b1,a2+b1),fg1],[(a2+b1,a2+b2),fg3]])
    953             return h
    954        
    955         if len(f.intervals())>1 or len(g.intervals())>1:
    956             z = Piecewise([[(-3*abs(N-M),3*abs(N-M)),0*xx**0]])
    957             ff = f.functions()
    958             gg = g.functions()
    959             intvlsf = f.intervals()
    960             intvlsg = g.intervals()
    961             for i in range(len(ff)):
    962                 for j in range(len(gg)):
    963                     f0 = Piecewise([[intvlsf[i],ff[i]]])
    964                     g0 = Piecewise([[intvlsg[j],gg[j]]])
    965                     h = g0.convolution(f0)
    966                     z = z + h
    967             return z.unextend()
    968    
    969     def derivative(self):
    970         r"""
    971         Returns the derivative (as computed by maxima)
    972         Piecewise(I,`(d/dx)(self|_I)`), as I runs over the
    973         intervals belonging to self. self must be piecewise polynomial.
    974        
    975         EXAMPLES::
    976        
    977             sage: f1(x) = 1
    978             sage: f2(x) = 1-x
    979             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    980             sage: f.derivative()
    981             Piecewise defined function with 2 parts, [[(0, 1), x |--> 0], [(1, 2), x |--> -1]]
    982             sage: f1(x) = -1
    983             sage: f2(x) = 2
    984             sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
    985             sage: f.derivative()
    986             Piecewise defined function with 2 parts, [[(0, 1/2*pi), x |--> 0], [(1/2*pi, pi), x |--> 0]]
    987            
    988         ::
    989        
    990             sage: f = Piecewise([[(0,1), (x * 2)]], x)
    991             sage: f.derivative()
    992             Piecewise defined function with 1 parts, [[(0, 1), x |--> 2]]
    993         """
    994         x = self.default_variable()
    995         dlist = [[(a, b), derivative(f(x), x).function(x)] for (a,b),f in self.list()]
    996         return Piecewise(dlist)
    997  
    998     def tangent_line(self, pt):
    999         """
    1000         Computes the linear function defining the tangent line of the
    1001         piecewise function self.
    1002        
    1003         EXAMPLES::
    1004        
    1005             sage: f1(x) = x^2
    1006             sage: f2(x) = 5-x^3+x
    1007             sage: f = Piecewise([[(0,1),f1],[(1,2),f2]])
    1008             sage: tf = f.tangent_line(0.9) ## tangent line at x=0.9
    1009             sage: P = f.plot(rgbcolor=(0.7,0.1,0.5), plot_points=40)
    1010             sage: Q = tf.plot(rgbcolor=(0.7,0.2,0.2), plot_points=40)
    1011             sage: P + Q
    1012         """
    1013         pt = QQ(pt)
    1014         R = QQ[self.default_variable()]
    1015         x = R.gen()
    1016         der = self.derivative()
    1017         tanline = (x-pt)*der(pt)+self(pt)
    1018         dlist = [[(a, b), tanline] for (a,b),f in self.list()]
    1019         return Piecewise(dlist)
    1020        
    1021     def plot(self, *args, **kwds):
    1022         """
    1023         Returns the plot of self.
    1024        
    1025         Keyword arguments are passed onto the plot command for each piece
    1026         of the function. E.g., the plot_points keyword affects each
    1027         segment of the plot.
    1028        
    1029         EXAMPLES::
    1030        
    1031             sage: f1(x) = 1
    1032             sage: f2(x) = 1-x
    1033             sage: f3(x) = exp(x)
    1034             sage: f4(x) = sin(2*x)
    1035             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    1036             sage: P = f.plot(rgbcolor=(0.7,0.1,0), plot_points=40)
    1037             sage: P
    1038        
    1039         Remember: to view this, type show(P) or P.save("path/myplot.png")
    1040         and then open it in a graphics viewer such as GIMP.
    1041 
    1042         TESTS:
    1043 
    1044         We should not add each piece to the legend individually, since
    1045         this creates duplicates (:trac:`12651`). This tests that only
    1046         one of the graphics objects in the plot has a non-``None``
    1047         ``legend_label``::
    1048 
    1049             sage: f1 = sin(x)
    1050             sage: f2 = cos(x)
    1051             sage: f = piecewise([[(-1,0), f1],[(0,1), f2]])
    1052             sage: p = f.plot(legend_label='$f(x)$')
    1053             sage: lines = [
    1054             ...     line
    1055             ...     for line in p._objects
    1056             ...     if line.options()['legend_label'] is not None ]
    1057             sage: len(lines)
    1058             1
    1059         """
    1060         from sage.plot.all import plot, Graphics
    1061 
    1062         g = Graphics()
    1063 
    1064         for i, ((a,b), f) in enumerate(self.list()):
    1065             # If it's the first piece, pass all arguments. Otherwise,
    1066             # filter out 'legend_label' so that we don't add each
    1067             # piece to the legend separately (trac #12651).
    1068             if i != 0 and 'legend_label' in kwds:
    1069                 del kwds['legend_label']
    1070 
    1071             g += plot(f, a, b, *args, **kwds)
    1072 
    1073         return g
    1074 
    1075     def fourier_series_cosine_coefficient(self,n,L):
    1076         r"""
    1077         Returns the n-th Fourier series coefficient of
    1078         `\cos(n\pi x/L)`, `a_n`.
    1079        
    1080         INPUT:
    1081        
    1082        
    1083         -  ``self`` - the function f(x), defined over -L x L
    1084        
    1085         -  ``n`` - an integer n=0
    1086        
    1087         -  ``L`` - (the period)/2
    1088        
    1089        
    1090         OUTPUT:
    1091         `a_n = \frac{1}{L}\int_{-L}^L f(x)\cos(n\pi x/L)dx`
    1092        
    1093         EXAMPLES::
    1094        
    1095             sage: f(x) = x^2
    1096             sage: f = Piecewise([[(-1,1),f]])
    1097             sage: f.fourier_series_cosine_coefficient(2,1)
    1098             pi^(-2)
    1099             sage: f(x) = x^2
    1100             sage: f = Piecewise([[(-pi,pi),f]])
    1101             sage: f.fourier_series_cosine_coefficient(2,pi)
    1102             1
    1103             sage: f1(x) = -1
    1104             sage: f2(x) = 2
    1105             sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
    1106             sage: f.fourier_series_cosine_coefficient(5,pi)
    1107             -3/5/pi
    1108         """
    1109         from sage.all import cos, pi
    1110         x = var('x')
    1111         result = sum([(f(x)*cos(pi*x*n/L)/L).integrate(x, a, b)
    1112                       for (a,b), f in self.list()])
    1113         if is_Expression(result):
    1114             return result.simplify_trig()
    1115         return result
    1116    
    1117     def fourier_series_sine_coefficient(self,n,L):
    1118         r"""
    1119         Returns the n-th Fourier series coefficient of
    1120         `\sin(n\pi x/L)`, `b_n`.
    1121        
    1122         INPUT:
    1123        
    1124        
    1125         -  ``self`` - the function f(x), defined over -L x L
    1126        
    1127         -  ``n`` - an integer n0
    1128        
    1129         -  ``L`` - (the period)/2
    1130        
    1131        
    1132         OUTPUT:
    1133         `b_n = \frac{1}{L}\int_{-L}^L f(x)\sin(n\pi x/L)dx`
    1134        
    1135         EXAMPLES::
    1136        
    1137             sage: f(x) = x^2
    1138             sage: f = Piecewise([[(-1,1),f]])
    1139             sage: f.fourier_series_sine_coefficient(2,1)  # L=1, n=2
    1140             0
    1141         """
    1142         from sage.all import sin, pi
    1143         x = var('x')
    1144         result = sum([(f(x)*sin(pi*x*n/L)/L).integrate(x, a, b)
    1145                       for (a,b), f in self.list()])
    1146         if is_Expression(result):
    1147             return result.simplify_trig()
    1148         return result
    1149 
    1150     def _fourier_series_helper(self, N, L, scale_function):
    1151         r"""
    1152         A helper function for the construction of Fourier series. The
    1153         argument scale_function is a function which takes in n,
    1154         representing the `n^{th}` coefficient, and return an
    1155         expression to scale the sine and cosine coefficients by.
    1156        
    1157         EXAMPLES::
    1158        
    1159             sage: f(x) = x^2
    1160             sage: f = Piecewise([[(-1,1),f]])
    1161             sage: f._fourier_series_helper(3, 1, lambda n: 1)
    1162             cos(2*pi*x)/pi^2 - 4*cos(pi*x)/pi^2 + 1/3
    1163         """
    1164         from sage.all import pi, sin, cos, srange
    1165         x = self.default_variable()
    1166         a0 = self.fourier_series_cosine_coefficient(0,L)
    1167         result = a0/2 + sum([(self.fourier_series_cosine_coefficient(n,L)*cos(n*pi*x/L) +
    1168                               self.fourier_series_sine_coefficient(n,L)*sin(n*pi*x/L))*
    1169                              scale_function(n)
    1170                              for n in srange(1,N)])
    1171         return result.expand()
    1172        
    1173 
    1174     def fourier_series_partial_sum(self,N,L):
    1175         r"""
    1176         Returns the partial sum
    1177        
    1178         .. math::
    1179        
    1180            f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
    1181        
    1182         as a string.
    1183        
    1184         EXAMPLE::
    1185        
    1186             sage: f(x) = x^2
    1187             sage: f = Piecewise([[(-1,1),f]])
    1188             sage: f.fourier_series_partial_sum(3,1)
    1189             cos(2*pi*x)/pi^2 - 4*cos(pi*x)/pi^2 + 1/3
    1190             sage: f1(x) = -1
    1191             sage: f2(x) = 2
    1192             sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
    1193             sage: f.fourier_series_partial_sum(3,pi)
    1194             -3*cos(x)/pi - 3*sin(2*x)/pi + 3*sin(x)/pi - 1/4
    1195         """
    1196         return self._fourier_series_helper(N, L, lambda n: 1)
    1197  
    1198     def fourier_series_partial_sum_cesaro(self,N,L):
    1199         r"""
    1200         Returns the Cesaro partial sum
    1201        
    1202         .. math::
    1203        
    1204            f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N (1-n/N)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
    1205        
    1206        
    1207         as a string. This is a "smoother" partial sum - the Gibbs
    1208         phenomenon is mollified.
    1209        
    1210         EXAMPLE::
    1211        
    1212             sage: f(x) = x^2
    1213             sage: f = Piecewise([[(-1,1),f]])
    1214             sage: f.fourier_series_partial_sum_cesaro(3,1)
    1215             1/3*cos(2*pi*x)/pi^2 - 8/3*cos(pi*x)/pi^2 + 1/3
    1216             sage: f1(x) = -1
    1217             sage: f2(x) = 2
    1218             sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
    1219             sage: f.fourier_series_partial_sum_cesaro(3,pi)
    1220             -2*cos(x)/pi - sin(2*x)/pi + 2*sin(x)/pi - 1/4
    1221         """
    1222         return self._fourier_series_helper(N, L, lambda n: 1-n/N)
    1223 
    1224     def fourier_series_partial_sum_hann(self,N,L):
    1225         r"""
    1226         Returns the Hann-filtered partial sum (named after von Hann, not
    1227         Hamming)
    1228        
    1229         .. math::
    1230        
    1231            f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N H_N(n)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
    1232        
    1233         as a string, where `H_N(x) = (1+\cos(\pi x/N))/2`. This is
    1234         a "smoother" partial sum - the Gibbs phenomenon is mollified.
    1235        
    1236         EXAMPLE::
    1237        
    1238             sage: f(x) = x^2
    1239             sage: f = Piecewise([[(-1,1),f]])
    1240             sage: f.fourier_series_partial_sum_hann(3,1)
    1241             1/4*cos(2*pi*x)/pi^2 - 3*cos(pi*x)/pi^2 + 1/3
    1242             sage: f1(x) = -1
    1243             sage: f2(x) = 2
    1244             sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
    1245             sage: f.fourier_series_partial_sum_hann(3,pi)
    1246             -9/4*cos(x)/pi - 3/4*sin(2*x)/pi + 9/4*sin(x)/pi - 1/4
    1247         """
    1248         from sage.all import cos, pi
    1249         return self._fourier_series_helper(N, L, lambda n: (1+cos(pi*n/N))/2)
    1250 
    1251     def fourier_series_partial_sum_filtered(self,N,L,F):
    1252         r"""
    1253         Returns the "filtered" partial sum
    1254        
    1255         .. math::
    1256        
    1257            f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N F_n*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
    1258        
    1259         as a string, where `F = [F_1,F_2, ..., F_{N}]` is a list
    1260         of length `N` consisting of real numbers. This can be used
    1261         to plot FS solutions to the heat and wave PDEs.
    1262        
    1263         EXAMPLE::
    1264        
    1265             sage: f(x) = x^2
    1266             sage: f = Piecewise([[(-1,1),f]])
    1267             sage: f.fourier_series_partial_sum_filtered(3,1,[1,1,1])
    1268             cos(2*pi*x)/pi^2 - 4*cos(pi*x)/pi^2 + 1/3
    1269             sage: f1(x) = -1
    1270             sage: f2(x) = 2
    1271             sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
    1272             sage: f.fourier_series_partial_sum_filtered(3,pi,[1,1,1])
    1273             -3*cos(x)/pi - 3*sin(2*x)/pi + 3*sin(x)/pi - 1/4
    1274         """
    1275         return self._fourier_series_helper(N, L, lambda n: F[n])
    1276        
    1277     def plot_fourier_series_partial_sum(self,N,L,xmin,xmax, **kwds):
    1278         r"""
    1279         Plots the partial sum
    1280        
    1281         .. math::
    1282        
    1283            f(x) \sim \frac{a_0}{2} +  sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
    1284        
    1285         over xmin x xmin.
    1286        
    1287         EXAMPLE::
    1288        
    1289             sage: f1(x) = -2
    1290             sage: f2(x) = 1
    1291             sage: f3(x) = -1
    1292             sage: f4(x) = 2
    1293             sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]])
    1294             sage: P = f.plot_fourier_series_partial_sum(3,pi,-5,5)    # long time
    1295             sage: f1(x) = -1
    1296             sage: f2(x) = 2
    1297             sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
    1298             sage: P = f.plot_fourier_series_partial_sum(15,pi,-5,5)   # long time
    1299        
    1300         Remember, to view this type show(P) or P.save("path/myplot.png")
    1301         and then open it in a graphics viewer such as GIMP.
    1302         """
    1303         from sage.plot.all import plot
    1304         return plot(self.fourier_series_partial_sum(N,L), xmin, xmax, **kwds)
    1305 
    1306     def plot_fourier_series_partial_sum_cesaro(self,N,L,xmin,xmax, **kwds):
    1307         r"""
    1308         Plots the partial sum
    1309        
    1310         .. math::
    1311        
    1312                      f(x) \sim \frac{a_0}{2} +                     \sum_{n=1}^N (1-n/N)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
    1313        
    1314        
    1315         over xmin x xmin. This is a "smoother" partial sum - the Gibbs
    1316         phenomenon is mollified.
    1317        
    1318         EXAMPLE::
    1319        
    1320             sage: f1(x) = -2
    1321             sage: f2(x) = 1
    1322             sage: f3(x) = -1
    1323             sage: f4(x) = 2
    1324             sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]])
    1325             sage: P = f.plot_fourier_series_partial_sum_cesaro(3,pi,-5,5)    # long time
    1326             sage: f1(x) = -1
    1327             sage: f2(x) = 2
    1328             sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
    1329             sage: P = f.plot_fourier_series_partial_sum_cesaro(15,pi,-5,5)   # long time
    1330        
    1331         Remember, to view this type show(P) or P.save("path/myplot.png")
    1332         and then open it in a graphics viewer such as GIMP.
    1333         """     
    1334         from sage.plot.all import plot
    1335         return plot(self.fourier_series_partial_sum_cesaro(N,L), xmin, xmax, **kwds)
    1336    
    1337     def plot_fourier_series_partial_sum_hann(self,N,L,xmin,xmax, **kwds):
    1338         r"""
    1339         Plots the partial sum
    1340        
    1341         .. math::
    1342        
    1343            f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N H_N(n)*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
    1344        
    1345        
    1346         over xmin x xmin, where H_N(x) = (0.5)+(0.5)\*cos(x\*pi/N) is the
    1347         N-th Hann filter.
    1348        
    1349         EXAMPLE::
    1350        
    1351             sage: f1(x) = -2
    1352             sage: f2(x) = 1
    1353             sage: f3(x) = -1
    1354             sage: f4(x) = 2
    1355             sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]])
    1356             sage: P = f.plot_fourier_series_partial_sum_hann(3,pi,-5,5)    # long time
    1357             sage: f1(x) = -1
    1358             sage: f2(x) = 2
    1359             sage: f = Piecewise([[(-pi,pi/2),f1],[(pi/2,pi),f2]])
    1360             sage: P = f.plot_fourier_series_partial_sum_hann(15,pi,-5,5)   # long time
    1361        
    1362         Remember, to view this type show(P) or P.save("path/myplot.png")
    1363         and then open it in a graphics viewer such as GIMP.
    1364         """     
    1365         from sage.plot.all import plot
    1366         return plot(self.fourier_series_partial_sum_hann(N,L), xmin, xmax, **kwds)
    1367        
    1368     def plot_fourier_series_partial_sum_filtered(self,N,L,F,xmin,xmax, **kwds):
    1369         r"""
    1370         Plots the partial sum
    1371        
    1372         .. math::
    1373        
    1374                      f(x) \sim \frac{a_0}{2} +                     \sum_{n=1}^N F_n*[a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         
    1375        
    1376        
    1377         over xmin x xmin, where `F = [F_1,F_2, ..., F_{N}]` is a
    1378         list of length `N` consisting of real numbers. This can be
    1379         used to plot FS solutions to the heat and wave PDEs.
    1380        
    1381         EXAMPLE::
    1382        
    1383             sage: f1(x) = -2
    1384             sage: f2(x) = 1
    1385             sage: f3(x) = -1
    1386             sage: f4(x) = 2
    1387             sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f3],[(pi/2,pi),f4]])
    1388             sage: P = f.plot_fourier_series_partial_sum_filtered(3,pi,[1]*3,-5,5)    # long time
    1389             sage: f1(x) = -1
    1390             sage: f2(x) = 2
    1391             sage: f = Piecewise([[(-pi,-pi/2),f1],[(-pi/2,0),f2],[(0,pi/2),f1],[(pi/2,pi),f2]])
    1392             sage: P = f.plot_fourier_series_partial_sum_filtered(15,pi,[1]*15,-5,5)   # long time
    1393        
    1394         Remember, to view this type show(P) or P.save("path/myplot.png")
    1395         and then open it in a graphics viewer such as GIMP.
    1396         """
    1397         from sage.plot.all import plot
    1398         return plot(self.fourier_series_partial_sum_filtered(N,L,F), xmin, xmax, **kwds)
    1399                
    1400     def fourier_series_value(self,x,L):
    1401         r"""
    1402         Returns the value of the Fourier series coefficient of self at
    1403         `x`,
    1404        
    1405        
    1406         .. math::
    1407        
    1408                      f(x) \sim \frac{a_0}{2} +                     \sum_{n=1}^\infty [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})],         \ \ \ -L<x<L.         
    1409        
    1410        
    1411         This method applies to piecewise non-polynomial functions as well.
    1412        
    1413         INPUT:
    1414        
    1415        
    1416         -  ``self`` - the function f(x), defined over -L x L
    1417        
    1418         -  ``x`` - a real number
    1419        
    1420         -  ``L`` - (the period)/2
    1421        
    1422        
    1423         OUTPUT: `(f^*(x+)+f^*(x-)/2`, where `f^*` denotes
    1424         the function `f` extended to `\RR` with period
    1425         `2L` (Dirichlet's Theorem for Fourier series).
    1426        
    1427         EXAMPLES::
    1428        
    1429             sage: f1(x) = 1
    1430             sage: f2(x) = 1-x
    1431             sage: f3(x) = exp(x)
    1432             sage: f4(x) = sin(2*x)
    1433             sage: f = Piecewise([[(-10,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    1434             sage: f.fourier_series_value(101,10) 
    1435             1/2
    1436             sage: f.fourier_series_value(100,10)
    1437             1
    1438             sage: f.fourier_series_value(10,10)
    1439             1/2*sin(20) + 1/2
    1440             sage: f.fourier_series_value(20,10)
    1441             1
    1442             sage: f.fourier_series_value(30,10)
    1443             1/2*sin(20) + 1/2
    1444             sage: f1(x) = -1
    1445             sage: f2(x) = 2
    1446             sage: f = Piecewise([[(-pi,0),lambda x:0],[(0,pi/2),f1],[(pi/2,pi),f2]])
    1447             sage: f.fourier_series_value(-1,pi)
    1448             0
    1449             sage: f.fourier_series_value(20,pi)
    1450             -1
    1451             sage: f.fourier_series_value(pi/2,pi)
    1452             1/2
    1453         """
    1454         xnew = x - int(RR(x/(2*L)))*2*L
    1455         endpts = self.end_points()
    1456         if xnew == endpts[0] or xnew == endpts[-1]:
    1457             return (self.functions()[0](endpts[0]) + self.functions()[-1](endpts[-1]))/2
    1458         else:
    1459             return self(xnew)
    1460 
    1461     def cosine_series_coefficient(self,n,L):
    1462         r"""
    1463         Returns the n-th cosine series coefficient of
    1464         `\cos(n\pi x/L)`, `a_n`.
    1465        
    1466         INPUT:
    1467        
    1468        
    1469         -  ``self`` - the function f(x), defined over 0 x L (no
    1470            checking is done to insure this)
    1471        
    1472         -  ``n`` - an integer n=0
    1473        
    1474         -  ``L`` - (the period)/2
    1475        
    1476        
    1477         OUTPUT:
    1478         `a_n = \frac{2}{L}\int_{-L}^L f(x)\cos(n\pi x/L)dx` such
    1479         that
    1480        
    1481         .. math::
    1482        
    1483                      f(x) \sim \frac{a_0}{2} +                     \sum_{n=1}^\infty a_n\cos(\frac{n\pi x}{L}),\ \ 0<x<L.         
    1484        
    1485        
    1486        
    1487         EXAMPLES::
    1488        
    1489             sage: f(x) = x
    1490             sage: f = Piecewise([[(0,1),f]])
    1491             sage: f.cosine_series_coefficient(2,1) 
    1492             0
    1493             sage: f.cosine_series_coefficient(3,1)
    1494             -4/9/pi^2
    1495             sage: f1(x) = -1
    1496             sage: f2(x) = 2
    1497             sage: f = Piecewise([[(0,pi/2),f1],[(pi/2,pi),f2]])
    1498             sage: f.cosine_series_coefficient(2,pi)
    1499             0
    1500             sage: f.cosine_series_coefficient(3,pi)
    1501             2/pi
    1502             sage: f.cosine_series_coefficient(111,pi)
    1503             2/37/pi
    1504             sage: f1 = lambda x: x*(pi-x)
    1505             sage: f = Piecewise([[(0,pi),f1]])
    1506             sage: f.cosine_series_coefficient(0,pi)
    1507             1/3*pi^2
    1508 
    1509         """
    1510         from sage.all import cos, pi
    1511         x = var('x')
    1512         result = sum([(2*f(x)*cos(pi*x*n/L)/L).integrate(x, a, b)
    1513                       for (a,b), f in self.list()])
    1514         if is_Expression(result):
    1515             return result.simplify_trig()
    1516         return result
    1517 
    1518 
    1519     def sine_series_coefficient(self,n,L):
    1520         r"""
    1521         Returns the n-th sine series coefficient of
    1522         `\sin(n\pi x/L)`, `b_n`.
    1523        
    1524         INPUT:
    1525        
    1526         -  ``self`` - the function f(x), defined over 0 x L (no
    1527            checking is done to insure this)
    1528        
    1529         -  ``n`` - an integer n0
    1530        
    1531         -  ``L`` - (the period)/2
     247        Combine piecewise operands into single piecewise function
    1532248       
    1533249        OUTPUT:
    1534250
    1535         `b_n = \frac{2}{L}\int_{-L}^L f(x)\sin(n\pi x/L)dx` such
    1536         that
     251        A piecewise function whose operands are not piecewiese if
     252        possible, that is, as long as the piecewise variable is the same.
     253
     254        """
     255        raise NotImplementedError
     256
     257
     258    class EvaluationMethods:
     259
     260        def expression_at(cls, self, parameters, variable, point):
     261            """
     262            Return the expression defining the piecewise function at
     263            ``value``
     264
     265            INPUT:
     266           
     267            - ``point`` -- a real number.
     268
     269            OUTPUT:
     270
     271            The symbolic expression defining the function value at the
     272            given ``point``.
     273
     274            EXAMPLES::
     275
     276                sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]);  f
     277                piecewise(x|-->sin(x) on {0}, x|-->cos(x) on (0, 2); x)
     278                sage: f.expression_at(0)
     279                sin(x)
     280                sage: f.expression_at(1)
     281                cos(x)
     282                sage: f.expression_at(2)
     283                Traceback (most recent call last):
     284                ...
     285                ValueError: point is not in the domain
     286            """
     287            for domain, func in parameters:
     288                if domain.contains(point):
     289                    return func
     290            raise ValueError('point is not in the domain')
     291
     292        which_function = expression_at
     293
     294        def domains(cls, self, parameters, variable):
     295            """
     296            Return the individual domains
     297           
     298            See also :meth:`~expressions`.
     299
     300            OUTPUT:
     301
     302            The collection of domains of the component functions as a
     303            tuple of :class:`~sage.sets.real_set.RealSet`.
     304           
     305            EXAMPLES::
     306
     307                sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]);  f
     308                piecewise(x|-->sin(x) on {0}, x|-->cos(x) on (0, 2); x)
     309                sage: f.domains()
     310                ({0}, (0, 2))
     311            """
     312            return tuple(dom for dom, fun in parameters)
     313
     314        def domain(cls, self, parameters, variable):
     315            """
     316            Return the domain
     317           
     318            OUTPUT:
     319
     320            The union of the domains of the individual pieces as a
     321            :class:`~sage.sets.real_set.RealSet`.
     322           
     323            EXAMPLES::
     324
     325                sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]);  f
     326                piecewise(x|-->sin(x) on {0}, x|-->cos(x) on (0, 2); x)
     327                sage: f.domain()
     328                [0, 2)
     329            """
     330            intervals = []
     331            for domain, func in parameters:
     332                intervals += list(domain)
     333            return RealSet(*intervals)
     334
     335        def __len__(cls, self, parameters, variable):
     336            """
     337            Return the number of "pieces"
     338
     339            OUTPUT:
     340
     341            Integer.
     342
     343            EXAMPLES::
     344
     345                sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]);  f
     346                piecewise(x|-->sin(x) on {0}, x|-->cos(x) on (0, 2); x)
     347                sage: len(f)
     348                2
     349            """
     350            return len(parameters)
     351
     352        def expressions(cls, self, parameters, variable):
     353            """
     354            Return the individual domains
     355           
     356            See also :meth:`~domains`.
     357
     358            OUTPUT:
     359
     360            The collection of expressions of the component functions.
     361           
     362            EXAMPLES::
     363
     364                sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]);  f
     365                piecewise(x|-->sin(x) on {0}, x|-->cos(x) on (0, 2); x)
     366                sage: f.expressions()
     367                (sin(x), cos(x))
     368            """
     369            return tuple(fun for dom, fun in parameters)
     370
     371        def iteritems(cls, self, parameters, variable):
     372            for pair in parameters:
     373                yield pair
     374
     375        def __call__(cls, self, parameters, variable, value=None, **kwds):
     376            """
     377            Call the piecewise function
     378
     379            EXAMPLES::
     380
     381                sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))]);  f
     382                piecewise(x|-->sin(x) on {0}, x|-->cos(x) on (0, 2); x)
     383                sage: f(0)
     384                0
     385                sage: f(1)
     386                cos(1)
     387                sage: f(2)
     388                Traceback (most recent call last):
     389                ...
     390                ValueError: point is not in the domain
     391            """
     392            self = piecewise(parameters, var=variable)
     393            substitution = dict()
     394            for k, v in kwds.iteritems():
     395                substitution[SR.var(k)] = v
     396            if value is not None:
     397                substitution[variable] = value
     398            return self.subs(substitution)
     399
     400        def _fast_float_(cls, self, *args):
     401            """
     402            Do not support the old ``fast_float``
     403
     404            OUTPUT:
     405
     406            This method raises ``NotImplementedError`` so that
     407            plotting uses the newer `fast_callable` implementation.
     408
     409            EXAMPLES::
     410           
     411                sage: f = piecewise([([0,0], sin(x)), ((0,2), cos(x))])
     412                sage: f._fast_float_()
     413                Traceback (most recent call last):
     414                ...
     415                NotImplementedError
     416            """
     417            raise NotImplementedError
     418
     419        def _fast_callable_(cls, self, parameters, variable, etb):
     420            """
     421            Override the ``fast_callable``
     422
     423            OUTPUT:
     424
     425            A :class:`~sage.ext.fast_callable.ExpressionCall`
     426            representing the piecewise function in the expression
     427            tree.
     428           
     429            EXAMPLES::
     430
     431                sage: p = piecewise([((-1, 0), -x), ([0, 1], x)], var=x)
     432                sage: from sage.ext.fast_callable import ExpressionTreeBuilder
     433                sage: etb = ExpressionTreeBuilder(vars=['x'])
     434                sage: p._fast_callable_(etb)
     435                {CommutativeRings.element_class}(v_0)
     436            """
     437            # print 'ev_fast_cal', parameters, variable, etb
     438            self = piecewise(parameters, var=variable)
     439            return etb.call(self, variable)
     440
     441        def restriction(cls, self, parameters, variable, restricted_domain):
     442            """
     443            Restrict the domain
     444
     445            INPUT:
     446           
     447            - ``restricted_domain`` -- a
     448              :class:`~sage.sets.real_set.RealSet` or something that
     449              defines one.
     450
     451            OUTPUT:
     452
     453            A new piecewise function obtained by restricting the domain.
     454
     455            EXAMPLES::
     456
     457                sage: f = piecewise([((-oo, oo), x)]);  f
     458                piecewise(x|-->x on (-oo, +oo); x)
     459                sage: f.restriction([[-1,1], [3,3]])
     460                piecewise(x|-->x on [-1, 1] + {3}; x)
     461            """
     462            restricted_domain = RealSet(*restricted_domain)
     463            new_param = []
     464            for domain, func in parameters:
     465                domain = domain.intersection(restricted_domain)
     466                new_param.append((domain, func))
     467            return piecewise(new_param, var=variable)
     468
     469        def extension(cls, self, parameters, variable, extension, extension_domain=None):
     470            """
     471            Extend the function
     472
     473            INPUT:
     474
     475            - ``extension`` -- a symbolic expression
     476
     477            - ``extension_domain`` -- a
     478              :class:`~sage.sets.real_set.RealSet` or ``None``
     479              (default). The domain of the extension. By default, the
     480              entire complement of the current domain.
     481
     482            EXAMPLES::
     483
     484                sage: f = piecewise([((-1,1), x)]);  f
     485                piecewise(x|-->x on (-1, 1); x)
     486                sage: f(3)
     487                Traceback (most recent call last):
     488                ...
     489                ValueError: point is not in the domain
     490
     491                sage: g = f.extension(0);  g
     492                piecewise(x|-->x on (-1, 1), x|-->0 on (-oo, -1] + [1, +oo); x)
     493                sage: g(3)
     494                0
     495
     496                sage: h = f.extension(1, RealSet.unbounded_above_closed(1));  h
     497                piecewise(x|-->x on (-1, 1), x|-->1 on [1, +oo); x)
     498                sage: h(3)
     499                1
     500            """
     501            self = piecewise(parameters, var=variable)
     502            if extension_domain is None:
     503                extension_domain = self.domain().complement()
     504            ext = ((extension_domain, SR(extension)),)
     505            return piecewise(parameters + ext, var=variable)
    1537506       
    1538         .. math::
    1539        
    1540            f(x) \sim \sum_{n=1}^\infty b_n\sin(\frac{n\pi x}{L}),\ \ 0<x<L.         
    1541        
    1542         EXAMPLES::
    1543        
    1544             sage: f(x) = 1
    1545             sage: f = Piecewise([[(0,1),f]])
    1546             sage: f.sine_series_coefficient(2,1) 
    1547             0
    1548             sage: f.sine_series_coefficient(3,1)
    1549             4/3/pi
    1550         """
    1551         from sage.all import sin, pi
    1552         x = var('x')
    1553         result = sum([(2*f(x)*sin(pi*x*n/L)/L).integrate(x, a, b)
    1554                       for (a,b), f in self.list()])
    1555         if is_Expression(result):
    1556             return result.simplify_trig()
    1557         return result
     507        def pieces(cls, self, parameters, variable):
     508            """
     509            Return the "pieces"
    1558510
    1559     def laplace(self, x='x', s='t'):
    1560         r"""
    1561         Returns the Laplace transform of self with respect to the variable
    1562         var.
    1563        
    1564         INPUT:
    1565        
    1566        
    1567         -  ``x`` - variable of self
    1568        
    1569         -  ``s`` - variable of Laplace transform.
    1570        
    1571        
    1572         We assume that a piecewise function is 0 outside of its domain and
    1573         that the left-most endpoint of the domain is 0.
    1574        
    1575         EXAMPLES::
    1576        
    1577             sage: x, s, w = var('x, s, w')
    1578             sage: f = Piecewise([[(0,1),1],[(1,2), 1-x]])
    1579             sage: f.laplace(x, s)
    1580             -e^(-s)/s + (s + 1)*e^(-2*s)/s^2 + 1/s - e^(-s)/s^2
    1581             sage: f.laplace(x, w)
    1582             -e^(-w)/w + (w + 1)*e^(-2*w)/w^2 + 1/w - e^(-w)/w^2
    1583            
    1584         ::
    1585        
    1586             sage: y, t = var('y, t')
    1587             sage: f = Piecewise([[(1,2), 1-y]])
    1588             sage: f.laplace(y, t)
    1589             (t + 1)*e^(-2*t)/t^2 - e^(-t)/t^2
    1590            
    1591         ::
    1592        
    1593             sage: s = var('s')
    1594             sage: t = var('t')
    1595             sage: f1(t) = -t
    1596             sage: f2(t) = 2
    1597             sage: f = Piecewise([[(0,1),f1],[(1,infinity),f2]])
    1598             sage: f.laplace(t,s)
    1599             (s + 1)*e^(-s)/s^2 + 2*e^(-s)/s - 1/s^2
    1600         """
    1601         from sage.all import assume, exp, forget
    1602         x = var(x)
    1603         s = var(s)
    1604         assume(s>0)
    1605         result =  sum([(SR(f)*exp(-s*x)).integral(x,a,b)
    1606                        for (a,b),f in self.list()])
    1607         forget(s>0)
    1608         return result
     511            OUTPUT:
    1609512
    1610     def _make_compatible(self, other):
    1611         """
    1612         Returns self and other extended to be defined on the same domain as
    1613         well as a refinement of their intervals. This is used for adding
    1614         and multiplying piecewise functions.
    1615        
    1616         EXAMPLES::
    1617        
    1618             sage: R.<x> = QQ[]
    1619             sage: f1 = Piecewise([[(0, 2), x]])
    1620             sage: f2 = Piecewise([[(1, 3), x^2]])
    1621             sage: f1._make_compatible(f2)
    1622             (Piecewise defined function with 2 parts, [[(0, 2), x], [(2, 3), 0]],
    1623             Piecewise defined function with 2 parts, [[(0, 1), 0], [(1, 3), x^2]],
    1624             [(0, 1), (1, 2), (2, 3)])
    1625         """
    1626         a1, b1 = self.domain()
    1627         a2, b2 = other.domain()
    1628         a = min(a1, a2)
    1629         b = max(b1, b2)
    1630         F = self.extend_by_zero_to(a,b)
    1631         G = other.extend_by_zero_to(a,b)
    1632         endpts = list(set(F.end_points()).union(set(G.end_points())))
    1633         endpts.sort()
    1634         return F, G, zip(endpts, endpts[1:])
    1635    
    1636     def __add__(self,other):
    1637         """
    1638         Returns the piecewise defined function which is the sum of self and
    1639         other. Does not require both domains be the same.
    1640        
    1641         EXAMPLES::
    1642        
    1643             sage: x = PolynomialRing(QQ,'x').gen()
    1644             sage: f1 = x^0
    1645             sage: f2 = 1-x
    1646             sage: f3 = 2*x
    1647             sage: f4 = 10-x
    1648             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    1649             sage: g1 = x-2
    1650             sage: g2 = x-5
    1651             sage: g = Piecewise([[(0,5),g1],[(5,10),g2]])
    1652             sage: h = f+g
    1653             sage: h
    1654             Piecewise defined function with 5 parts, [[(0, 1), x - 1], [(1, 2), -1], [(2, 3), 3*x - 2], [(3, 5), 8], [(5, 10), 5]]
    1655        
    1656         Note that in this case the functions must be defined using
    1657         polynomial expressions *not* using the lambda notation.
    1658         """
    1659         F, G, intervals = self._make_compatible(other)
    1660         fcn = []
    1661         for a,b in intervals:
    1662             fcn.append([(a,b), F.which_function(b)+G.which_function(b)])       
    1663         return Piecewise(fcn)
    1664        
    1665     def __mul__(self,other):
    1666         r"""
    1667         Returns the piecewise defined function which is the product of one
    1668         piecewise function (self) with another one (other).
    1669        
    1670         EXAMPLES::
    1671        
    1672             sage: x = PolynomialRing(QQ,'x').gen()
    1673             sage: f1 = x^0
    1674             sage: f2 = 1-x
    1675             sage: f3 = 2*x
    1676             sage: f4 = 10-x
    1677             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    1678             sage: g1 = x-2
    1679             sage: g2 = x-5
    1680             sage: g = Piecewise([[(0,5),g1],[(5,10),g2]])
    1681             sage: h = f*g
    1682             sage: h
    1683             Piecewise defined function with 5 parts, [[(0, 1), x - 2], [(1, 2), -x^2 + 3*x - 2], [(2, 3), 2*x^2 - 4*x], [(3, 5), -x^2 + 12*x - 20], [(5, 10), -x^2 + 15*x - 50]]
    1684             sage: g*(11/2)
    1685             Piecewise defined function with 2 parts, [[(0, 5), 11/2*x - 11], [(5, 10), 11/2*x - 55/2]]
    1686            
    1687         Note that in this method the functions must be defined using
    1688         polynomial expressions *not* using the lambda notation.
    1689         """
    1690         ## needed for scalar multiplication
    1691         if isinstance(other,Rational) or isinstance(other,Integer):
    1692             return Piecewise([[(a,b), other*f] for (a,b),f in self.list()])
    1693         else:
    1694             F, G, intervals = self._make_compatible(other)
    1695             fcn = []
    1696             for a,b in intervals:
    1697                 fcn.append([(a,b),F.which_function(b)*G.which_function(b)])     
    1698             return Piecewise(fcn)
     513            A tuple of piecewise functions, each having only a single
     514            expression.
    1699515
    1700     __rmul__ = __mul__
     516            EXAMPLES::
    1701517
    1702     def __eq__(self,other):
    1703         """
    1704         Implements Boolean == operator.
     518                sage: p = piecewise([((-1, 0), -x), ([0, 1], x)], var=x)
     519                sage: p.pieces()
     520                (piecewise(x|-->-x on (-1, 0); x),
     521                 piecewise(x|-->x on [0, 1]; x))
     522            """
     523            result = []
     524            for domain, func in parameters:
     525                result.append(piecewise([(domain, func)], var=variable))
     526            return tuple(result)
    1705527
    1706         EXAMPLES::
    1707        
    1708             sage: f1 = x^0
    1709             sage: f2 = 1-x
    1710             sage: f3 = 2*x
    1711             sage: f4 = 10-x
    1712             sage: f = Piecewise([[(0,1),f1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    1713             sage: g = Piecewise([[(0,1),1],[(1,2),f2],[(2,3),f3],[(3,10),f4]])
    1714             sage: f==g
    1715             True
    1716         """
    1717         return self.list()==other.list()
     528
     529
     530
     531piecewise = PiecewiseFunction()
     532
     533def Piecewise(*args, **kwds):
     534    """
     535    Deprecated spelling of ``piecewise``
     536
     537    EXAMPLES::
     538
     539        sage: Piecewise([((-1, 0), -x), ([0, 1], x)], var=x)
     540        doctest:...: DeprecationWarning: use lower-case piecewise instead
     541        See http://trac.sagemath.org/14801 for details.
     542        piecewise(x|-->-x on (-1, 0), x|-->x on [0, 1]; x)
     543    """
     544    from sage.misc.superseded import deprecation
     545    deprecation(14801, 'use lower-case piecewise instead')
     546    return piecewise(*args, **kwds)
     547
     548
     549
  • sage/symbolic/expression_conversions.py

    diff --git a/sage/symbolic/expression_conversions.py b/sage/symbolic/expression_conversions.py
    a b  
    2020from sage.symbolic.pynac import I
    2121from sage.functions.all import exp
    2222from sage.symbolic.operators import arithmetic_operators, relation_operators, FDerivativeOperator
     23from sage.functions.piecewise import piecewise
    2324from sage.rings.number_field.number_field_element_quadratic import NumberFieldElement_quadratic
    2425GaussianField = I.pyobject().parent()
    2526
     
    11811182        try:
    11821183            return self.ff.fast_float_constant(float(ex))
    11831184        except TypeError:
    1184             raise ValueError, "free variable: %s" % repr(ex)
     1185            raise NotImplementedError, "free variable: %s" % repr(ex)
    11851186
    11861187    def arithmetic(self, ex, operator):
    11871188        """