Introduction to Special Functions
=================================

Special function notation follows:

     bessel_j (index, expr)         Bessel function, 1st kind
     bessel_y (index, expr)         Bessel function, 2nd kind
     bessel_i (index, expr)         Modified Bessel function, 1st kind
     bessel_k (index, expr)         Modified Bessel function, 2nd kind

     hankel_1 (v,z)                 Hankel function of the 1st kind
     hankel_2 (v,z)                 Hankel function of the 2nd kind
     struve_h (v,z)                 Struve H function
     struve_l (v,z)                 Struve L function

     assoc_legendre_p[v,u] (z)      Legendre function of degree v and order u
     assoc_legendre_q[v,u] (z)      Legendre function, 2nd kind

     %f[p,q] ([], [], expr)         Generalized Hypergeometric function
     gamma (z)                      Gamma function
     gamma_incomplete_lower (a,z)   Lower incomplete gamma function
     gamma_incomplete (a,z)         Tail of incomplete gamma function
     hypergeometric (l1, l2, z)     Hypergeometric function
     %s[u,v] (z)                    Lommel "small" s function
     slommel[u,v] (z)               Lommel "big" S function
     %m[u,k] (z)                    Whittaker function, 1st kind
     %w[u,k] (z)                    Whittaker function, 2nd kind
     erfc (z)                       Complement of the erf function

     expintegral_e (v,z)            Exponential integral E
     expintegral_e1 (z)             Exponential integral E1
     expintegral_ei (z)             Exponential integral Ei
     expintegral_li (z)             Logarithmic integral Li
     expintegral_si (z)             Exponential integral Si
     expintegral_ci (z)             Exponential integral Ci
     expintegral_shi (z)            Exponential integral Shi
     expintegral_chi (z)            Exponential integral Chi

     kelliptic (z)                  Complete elliptic integral of the first
                                    kind (K)
     parabolic_cylinder_d (v,z)     Parabolic cylinder D function

Categories: Bessel functions Airy functions Special functions

Bessel Functions
================

 -- Function: bessel_j (V, Z)

     The Bessel function of the first kind of order v and argument z. See eqn 9.1.10, 360 and 10.2.E2.

     bessel_j is defined as

                          inf
                          ====       k  - v - 2 k  v + 2 k
                          \     (- 1)  2          z
                           >    --------------------------
                          /        k! gamma(v + k + 1)
                          ====
                          k = 0

     although the infinite series is not used for computations.

     When besselexpand is true, bessel_j is expanded in terms of elementary functions when the order v is half of an odd integer. See besselexpand

     Categories: Bessel functions Special functions

 -- Function: bessel_y (V, Z)

     The Bessel function of the second kind of order v and argument z. See eqn 9.1.2, 358 and 10.2.E3.

     bessel_y is defined as

                        cos(%pi v) bessel_j(v, z) - bessel_j(-v, z)
                        -------------------------------------------
                                       sin(%pi v)

     when v is not an integer.  When v is an integer n, the limit as v approaches n is taken.

     When besselexpand is true, bessel_y is expanded in terms of elementary functions when the order v is half of an odd integer. See besselexpand

     Categories: Bessel functions Special functions

 -- Function: bessel_i (V, Z)

     The modified Bessel function of the first kind of order v and argument z. See eqn 9.6.10, 375 and 10.25.E2.

     bessel_i is defined as

                              inf
                              ====   - v - 2 k  v + 2 k
                              \     2          z
                               >    -------------------
                              /     k! gamma(v + k + 1)
                              ====
                              k = 0

     although the infinite series is not used for computations.

     When besselexpand is true, bessel_i is expanded in terms of elementary functions when the order v is half of an odd integer. See besselexpand

     Categories: Bessel functions Special functions

 -- Function: bessel_k (V, Z)

     The modified Bessel function of the second kind of order v and argument z. See eqn 9.6.2, 375 and 10.27.E4.

     bessel_k is defined as

                     %pi csc(%pi v) (bessel_i(-v, z) - bessel_i(v, z))
                     -------------------------------------------------
                                            2

     when v is not an integer.  If v is an integer n, then the limit as v approaches n is taken.

     When besselexpand is true, bessel_k is expanded in terms of elementary functions when the order v is half of an odd integer. See besselexpand

     Categories: Bessel functions Special functions

 -- Function: hankel_1 (V, Z)

     The Hankel function of the first kind of order v and argument z. See eqn 9.1.3, 358 and 10.4.E3.

     hankel_1 is defined as

             bessel_j(v,z) + %i * bessel_y(v,z)

     Maxima evaluates hankel_1 numerically for a complex order v and complex argument z in float precision. The numerical evaluation in bigfloat precision is not supported.

     When besselexpand is true, hankel_1 is expanded in terms of elementary functions when the order v is half of an odd integer. See besselexpand

     Maxima knows the derivative of hankel_1 wrt the argument z.

     Examples:

     Numerical evaluation:

          (%i1) hankel_1(1,0.5);
          (%o1)        0.24226845767487 - 1.471472392670243 %i

          (%i2) hankel_1(1,0.5+%i);
          (%o2)       - 0.25582879948621 %i - 0.23957560188301

     Expansion of hankel_1 when besselexpand is true:

          (%i1) hankel_1(1/2,z),besselexpand:true;
                         sqrt(2) sin(z) - sqrt(2) %i cos(z)
          (%o1)          ----------------------------------
                                 sqrt(%pi) sqrt(z)

     Derivative of hankel_1 wrt the argument z. The derivative wrt the order v is not supported. Maxima returns a noun form:

          (%i1) diff(hankel_1(v,z),z);
                       hankel_1(v - 1, z) - hankel_1(v + 1, z)
          (%o1)        ---------------------------------------
                                          2

          (%i2) diff(hankel_1(v,z),v);
                                 d
          (%o2)                  -- (hankel_1(v, z))
                                 dv

     Categories: Bessel functions Special functions

 -- Function: hankel_2 (V, Z)

     The Hankel function of the second kind of order v and argument z. See eqn 9.1.4, 358 and 10.4.E3.

     hankel_2 is defined as

             bessel_j(v,z) - %i * bessel_y(v,z)

     Maxima evaluates hankel_2 numerically for a complex order v and complex argument z in float precision. The numerical evaluation in bigfloat precision is not supported.

     When besselexpand is true, hankel_2 is expanded in terms of elementary functions when the order v is half of an odd integer. See besselexpand

     Maxima knows the derivative of hankel_2 wrt the argument z.

     For examples see hankel_1.

     Categories: Bessel functions Special functions

 -- Option variable: besselexpand
     Default value: false

     Controls expansion of the Bessel, Hankel and Struve functions when the order is half of an odd integer.  In this case, the functions can be expanded in terms of other elementary functions.  When besselexpand is true, the Bessel function is expanded.

          (%i1) besselexpand: false$
          (%i2) bessel_j (3/2, z);
                                              3
          (%o2)                      bessel_j(-, z)
                                              2
          (%i3) besselexpand: true$
          (%i4) bessel_j (3/2, z);
                                                  sin(z)   cos(z)
                                 sqrt(2) sqrt(z) (------ - ------)
                                                     2       z
                                                    z
          (%o4)                  ---------------------------------
                                             sqrt(%pi)
          (%i5) bessel_y(3/2,z);
                                                  sin(z)    cos(z)
                              sqrt(2) sqrt(z) ((- ------) - ------)
                                                    z          2
                                                              z
          (%o5)               -------------------------------------
                                            sqrt(%pi)
          (%i6) bessel_i(3/2,z);
                                                cosh(z)   sinh(z)
                               sqrt(2) sqrt(z) (------- - -------)
                                                   z         2
                                                            z
          (%o6)                -----------------------------------
                                            sqrt(%pi)
          (%i7) bessel_k(3/2,z);
                                                1        - z
                                     sqrt(%pi) (- + 1) %e
                                                z
          (%o7)                      -----------------------
                                         sqrt(2) sqrt(z)

     Categories: Bessel functions Simplification flags and variables Special functions

 -- Function: scaled_bessel_i (V, Z)

     The scaled modified Bessel function of the first kind of order v and argument z.  That is,

     scaled_bessel_i(v,z) = exp(-abs(z)) * bessel_i(v,z).

     This function is particularly useful for calculating bessel_i for large z, which is large. However, maxima does not otherwise know much about this function.  For symbolic work, it is probably preferable to work with the expression exp(-abs(z))*bessel_i(v, z).

     Categories: Bessel functions

 -- Function: scaled_bessel_i0 (Z)

     Identical to scaled_bessel_i(0,z).

     Categories: Bessel functions Special functions

 -- Function: scaled_bessel_i1 (Z)

     Identical to scaled_bessel_i(1,z). Categories: Bessel functions Special functions

 -- Function: %s [U,V] (Z)
     Lommel's little s[u,v](z) function. (11.9.E3)(G&R 8.570.1).

     This Lommel function is the particular solution of the inhomogeneous Bessel differential equation:

     diff(w,z,2)+1/z*diff(w,z)+(1-v^2/z^2)*w = z^(u-1)

     This can be defined by the series

     %s[u,v](z) = z^(u+1)*sum((-1)^k*z^(2*k)/a[k+1](u,v)

     where

     a[k](u,v) = product(((u+2*m-1)^2-v^2),m,1,k) = 4^k*pochhammer((u-v+1)/2,k)*pochhammer((u+v+1)/2,k)

     Categories: Bessel functions Special functions

 -- Function: slommel [U,V] (Z)
     Lommel's big S[u,v](z) function. (11.9.E5)(G&R 8.570.2).

     Lommels big S function is another particular solution of the inhomogeneous Bessel differential equation (see %s) defined for all values of u and v , where

     slommel[u,v](z) = %s[u,v](z) + 2^(u-1)*gamma(u/2+v/2+1/2)*gamma(u/2-v/2+1/2) * (sin(1/2*(u-v)*%pi)*bessel_j(v,z) - cos(1/2*(u-v)*%pi)*bessel_y(v,z))

     When u+/-v ) is an odd negative integer, the limit must be used.

     Categories: Bessel functions Special functions

Airy Functions
==============

The Airy functions Ai(x) and Bi(x) are defined in Abramowitz and Stegun, Handbook of Mathematical Functions, Section 10.4, 446 and 9.

The Airy differential equation is:

diff (y(x), x, 2) - x y(x) = 0

The numerically satisfactory pair of solutions (9.2#T1) on the real line are y = Ai(x) and y = Bi(x).

These two solutions are oscillatory for x < 0. Ai(x) is the solution subject to the condition that y->0 as x -> +inf, and Bi(x) is the second solution with the same amplitude as Ai(x) as x->minf which differs in phase by %pi/2. Also, Bi(x) is unbounded as x->+inf.

If the argument x is a real or complex floating point number, the numerical value of the function is returned.

 -- Function: airy_ai (X)
     The Airy function Ai(x). See eqn 10.4.2, 446 and 9.

     See also airy_bi airy_dai and airy_dbi

     Categories: Airy functions Special functions

 -- Function: airy_dai (X)
     The derivative of the Airy function Ai(x) :

     airy_dai(x) = diff(airy_ai(x),x)

     See airy_ai

     Categories: Airy functions Special functions

 -- Function: airy_bi (X)
     The Airy function Bi(x) .  See eqn 10.4.3, 446 and  9.

     See airy_ai and airy_dbi

     Categories: Airy functions Special functions

 -- Function: airy_dbi (X)
     The derivative of the Airy function Bi(x)>>> :

     airy_dbi(x) = diff(airy_bi(x),x)

     See airy_ai and airy_bi

     Categories: Airy functions Special functions

Gamma and Factorial Functions
=============================

The gamma function and the related beta, psi and incomplete gamma functions are defined in Abramowitz and Stegun, Handbook of Mathematical Functions, Chapter 6.

 -- Function: bffac (EXPR, N)

     Bigfloat version of the factorial (shifted gamma) function.  The second argument is how many digits to retain and return, it's a good idea to request a couple of extra.

          (%i1) bffac(1/2,16);
          (%o1)                        8.862269254527584b-1
          (%i2) (1/2)!,numer;
          (%o2)                          0.886226925452758
          (%i3) bffac(1/2,32);
          (%o3)                 8.862269254527580136490837416707b-1

     Categories: Gamma and factorial functions Numerical evaluation

 -- Function: bfpsi (N, Z, FPPREC)
 -- Function: bfpsi0 (Z, FPPREC)

     bfpsi is the polygamma function of real argument Z and integer order N.  See psi for further information.  bfpsi0 is the digamma function. bfpsi0(Z, FPPREC) is equivalent to bfpsi(0, Z, FPPREC).

     These functions return bigfloat values. FPPREC is the bigfloat precision of the return value.

          (%i1) bfpsi0(1/3, 15);
          (%o1)                        - 3.13203378002081b0
          (%i2) bfpsi0(1/3, 32);
          (%o2)                - 3.1320337800208063229964190742873b0
          (%i3) bfpsi(0,1/3,32);
          (%o3)                - 3.1320337800208063229964190742873b0
          (%i4) psi[0](1/3);
                                    3 log(3)       %pi
          (%o4)                  (- --------) - --------- - %gamma
                                       2        2 sqrt(3)
          (%i5) float(%);
          (%o5)                         - 3.132033780020806

     Categories: Gamma and factorial functions Numerical evaluation

 -- Function: cbffac (Z, FPPREC)
     Complex bigfloat factorial.

     load ("bffac") loads this function.

          (%i1) cbffac(1+%i,16);
          (%o1)           3.430658398165453b-1 %i + 6.529654964201666b-1
          (%i2) (1+%i)!,numer;
          (%o2)             0.3430658398165453 %i + 0.6529654964201667

     Categories: Gamma and factorial functions Complex variables Numerical evaluation

 -- Function: gamma (Z)

     The basic definition of the gamma function (5.2.E1 and eqn 6.1.1, 255) is

                                   inf
                                  /
                                  [     z - 1   - t
                       gamma(z) = I    t      %e    dt
                                  ]
                                  /
                                   0

     Maxima simplifies gamma for positive integer and positive and negative rational numbers. For half integral values the result is a rational number times sqrt(%pi) . The simplification for integer values is controlled by factlim. For integers greater than factlim the numerical result of the factorial function, which is used to calculate gamma, will overflow. The simplification for rational numbers is controlled by gammalim to avoid internal overflow. See factlim and gammalim.

     For negative integers gamma is not defined.

     Maxima can evaluate gamma numerically for real and complex values in float and bigfloat precision.

     gamma has mirror symmetry.

     When gamma_expand is true, Maxima expands gamma for arguments z+n and z-n where n is an integer.

     Maxima knows the derivative of gamma.

     Examples:

     Simplification for integer, half integral, and rational numbers:

          (%i1) map('gamma,[1,2,3,4,5,6,7,8,9]);
          (%o1)        [1, 1, 2, 6, 24, 120, 720, 5040, 40320]
          (%i2) map('gamma,[1/2,3/2,5/2,7/2]);
                              sqrt(%pi)  3 sqrt(%pi)  15 sqrt(%pi)
          (%o2)   [sqrt(%pi), ---------, -----------, ------------]
                                  2           4            8
          (%i3) map('gamma,[2/3,5/3,7/3]);
                                            2           1
                                    2 gamma(-)  4 gamma(-)
                                2           3           3
          (%o3)          [gamma(-), ----------, ----------]
                                3       3           9

     Numerical evaluation for real and complex values:

          (%i4) map('gamma,[2.5,2.5b0]);
          (%o4)     [1.329340388179137, 1.3293403881791370205b0]
          (%i5) map('gamma,[1.0+%i,1.0b0+%i]);
          (%o5) [0.498015668118356 - .1549498283018107 %i,
                    4.9801566811835604272b-1 - 1.5494982830181068513b-1 %i]

     gamma has mirror symmetry:

          (%i6) declare(z,complex)$
          (%i7) conjugate(gamma(z));
          (%o7)                  gamma(conjugate(z))

     Maxima expands gamma(z+n) and gamma(z-n), when gamma_expand is true:

          (%i8) gamma_expand:true$

          (%i9) [gamma(z+1),gamma(z-1),gamma(z+2)/gamma(z+1)];
                                         gamma(z)
          (%o9)             [z gamma(z), --------, z + 1]
                                          z - 1

     The derivative of gamma:

          (%i10) diff(gamma(z),z);
          (%o10)                  psi (z) gamma(z)
                                     0

     See also makegamma

     The Euler-Mascheroni constant is %gamma.

     Categories: Gamma and factorial functions Special functions

 -- Function: log_gamma (Z)

     The natural logarithm of the gamma function.

          (%i1) gamma(6);
          (%o1)                                 120
          (%i2) log_gamma(6);
          (%o2)                              log(120)
          (%i3) log_gamma(0.5);
          (%o3)                         0.5723649429247004

     Categories: Gamma and factorial functions Special functions

 -- Function: gamma_incomplete_lower (A, Z)

     The lower incomplete gamma function (8.2.E1 and eqn 6.5.2, 260):

                                              z
                                             /
                                             [  a - 1   - t
              gamma_incomplete_lower(a, z) = I t      %e    dt
                                             ]
                                             /
                                              0

     See also gamma_incomplete (upper incomplete gamma function).

     Categories: Gamma and factorial functions Special functions

 -- Function: gamma_incomplete (A, Z)

     The incomplete upper gamma function (8.2.E2 and eqn 6.5.3, 260):

                                        inf
                                       /
                                       [     a - 1   - t
              gamma_incomplete(a, z) = I    t      %e    dt
                                       ]
                                       /
                                        z

     See also gamma_expand for controlling how gamma_incomplete is expressed in terms of elementary functions and erfc.

     Also see the related functions gamma_incomplete_regularized and gamma_incomplete_generalized.

     Categories: Gamma and factorial functions Special functions

 -- Function: gamma_incomplete_regularized (A, Z)

     The regularized incomplete upper gamma function (8.2.E4):

          gamma_incomplete_regularized(a, z) =
                                                  gamma_incomplete(a, z)
                                                  ----------------------
                                                         gamma(a)

     See also gamma_expand for controlling how gamma_incomplete is expressed in terms of elementary functions and erfc

     Also see gamma_incomplete.

     Categories: Gamma and factorial functions Special functions

 -- Function: gamma_incomplete_generalized (A, Z1, Z1)

     The generalized incomplete gamma function.

          gamma_incomplete_generalized(a, z1, z2) =
                                                         z2
                                                        /
                                                        [    a - 1   - t
                                                        I   t      %e    dt
                                                        ]
                                                        /
                                                         z1

     Also see gamma_incomplete and gamma_incomplete_regularized.

     Categories: Gamma and factorial functions Special functions

 -- Option variable: gamma_expand
     Default value: false

     gamma_expand controls expansion of gamma_incomplete. When gamma_expand is true, gamma_incomplete(v,z) is expanded in terms of z, exp(z), and gamma_incomplete or erfc when possible.

          (%i1) gamma_incomplete(2,z);
          (%o1)                       gamma_incomplete(2, z)
          (%i2) gamma_expand:true;
          (%o2)                                true
          (%i3) gamma_incomplete(2,z);
                                                     - z
          (%o3)                            (z + 1) %e

          (%i4) gamma_incomplete(3/2,z);
                                        - z   sqrt(%pi) erfc(sqrt(z))
          (%o4)               sqrt(z) %e    + -----------------------
                                                         2

          (%i5) gamma_incomplete(4/3,z);
                                                              1
                                             gamma_incomplete(-, z)
                                 1/3   - z                    3
          (%o5)                 z    %e    + ----------------------
                                                       3

          (%i6) gamma_incomplete(a+2,z);
                       a               - z
          (%o6)       z  (z + a + 1) %e    + a (a + 1) gamma_incomplete(a, z)
          (%i7) gamma_incomplete(a-2, z);
                  gamma_incomplete(a, z)    a - 2         z            1      - z
          (%o7)   ---------------------- - z      (--------------- + -----) %e
                     (1 - a) (2 - a)               (a - 2) (a - 1)   a - 2

     Categories: Gamma and factorial functions Simplification flags and variables

 -- Option variable: gammalim
     Default value: 10000

     gammalim controls simplification of the gamma function for integral and rational number arguments.  If the absolute value of the argument is not greater than gammalim, then simplification will occur.  Note that the factlim switch controls simplification of the result of gamma of an integer argument as well.

     Categories: Gamma and factorial functions Simplification flags and variables

 -- Function: makegamma (EXPR)
     Transforms instances of binomial, factorial, and beta functions in EXPR into gamma functions.

     See also makefact

          (%i1) makegamma(binomial(n,k));
                                           gamma(n + 1)
          (%o1)                    -----------------------------
                                   gamma(k + 1) gamma(n - k + 1)
          (%i2) makegamma(x!);
          (%o2)                            gamma(x + 1)
          (%i3) makegamma(beta(a,b));
                                         gamma(a) gamma(b)
          (%o3)                          -----------------
                                           gamma(b + a)

     Categories: Gamma and factorial functions

 -- Function: beta (A, B)
     The beta function is defined as gamma(a) gamma(b)/gamma(a+b)

     (5.12.E1 and eqn 6.2.1, 258).

     Maxima simplifies the beta function for positive integers and rational numbers, which sum to an integer. When beta_args_sum_to_integer is true, Maxima simplifies also general expressions which sum to an integer.

     For A or B equal to zero the beta function is not defined.

     In general the beta function is not defined for negative integers as an argument. The exception is for A=-N, N a positive integer and B a positive integer with b<=n, it is possible to define an analytic continuation. Maxima gives for this case a result.

     When beta_expand is true, expressions like beta(a+n,b) and beta(a-n,b) or beta(a,b+n) and beta(a,b-n) with n an integer are simplified.

     Maxima can evaluate the beta function for real and complex values in float and bigfloat precision. For numerical evaluation Maxima uses log_gamma:

                     - log_gamma(b + a) + log_gamma(b) + log_gamma(a)
                   %e

     Maxima knows that the beta function is symmetric and has mirror symmetry.

     Maxima knows the derivatives of the beta function with respect to A or B.

     To express the beta function as a ratio of gamma functions see makegamma.

     Examples:

     Simplification, when one of the arguments is an integer:

          (%i1) [beta(2,3),beta(2,1/3),beta(2,a)];
                                         1   9      1
          (%o1)                         [--, -, ---------]
                                         12  4  a (a + 1)

     Simplification for two rational numbers as arguments which sum to an integer:

          (%i2) [beta(1/2,5/2),beta(1/3,2/3),beta(1/4,3/4)];
                                    3 %pi   2 %pi
          (%o2)                    [-----, -------, sqrt(2) %pi]
                                      8    sqrt(3)

     When setting beta_args_sum_to_integer to true more general expression are simplified, when the sum of the arguments is an integer:

          (%i3) beta_args_sum_to_integer:true$
          (%i4) beta(a+1,-a+2);
                                          %pi (a - 1) a
          (%o4)                         ------------------
                                        2 sin(%pi (2 - a))

     The possible results, when one of the arguments is a negative integer:

          (%i5) [beta(-3,1),beta(-3,2),beta(-3,3)];
                                              1  1    1
          (%o5)                            [- -, -, - -]
                                              3  6    3

     beta(a+n,b) or beta(a-n,b) with n an integer simplifies when beta_expand is true:

          (%i6) beta_expand:true$
          (%i7) [beta(a+1,b),beta(a-1,b),beta(a+1,b)/beta(a,b+1)];
                              a beta(a, b)  beta(a, b) (b + a - 1)  a
          (%o7)              [------------, ----------------------, -]
                                 b + a              a - 1           b

     Beta is not defined, when one of the arguments is zero:

          (%i7) beta(0,b);
          beta: expected nonzero arguments; found 0, b
           -- an error.  To debug this try debugmode(true);

     Numerical evaluation for real and complex arguments in float or bigfloat precision:

          (%i8) beta(2.5,2.3);
          (%o8) .08694748611299981

          (%i9) beta(2.5,1.4+%i);
          (%o9) 0.0640144950796695 - .1502078053286415 %i

          (%i10) beta(2.5b0,2.3b0);
          (%o10) 8.694748611299969b-2

          (%i11) beta(2.5b0,1.4b0+%i);
          (%o11) 6.401449507966944b-2 - 1.502078053286415b-1 %i

     Beta is symmetric and has mirror symmetry:

          (%i14) beta(a,b)-beta(b,a);
          (%o14)                                 0
          (%i15) declare(a,complex,b,complex)$
          (%i16) conjugate(beta(a,b));
          (%o16)                 beta(conjugate(a), conjugate(b))

     The derivative of the beta function wrt a:

          (%i17) diff(beta(a,b),a);
          (%o17)               - beta(a, b) (psi (b + a) - psi (a))
                                                0             0

     Categories: Gamma and factorial functions

 -- Function: beta_incomplete (A, B, Z)

     The basic definition of the incomplete beta function (8.17.E1 and eqn 6.6.1, 263) is

                  z
                 /
                 [         b - 1  a - 1
                 I  (1 - t)      t      dt
                 ]
                 /
                  0

     This definition is possible for realpart(a)>0 and realpart(b)>0 and abs(z)<1 . For other values the incomplete beta function can be defined through a generalized hypergeometric function:

             gamma(a) hypergeometric_generalized([a, 1 - b], [a + 1], z) z

     (See https://functions.wolfram.com/GammaBetaErf/Beta3/ for a complete definition of the incomplete beta function.)

     For negative integers a = -n and positive integers b=m with m<=n the incomplete beta function is defined through

                                      m - 1           k
                                      ====  (1 - m)  z
                                n - 1 \            k
                               z       >    -----------
                                      /     k! (n - k)
                                      ====
                                      k = 0

     Maxima uses this definition to simplify beta_incomplete for A a negative integer.

     For A a positive integer, beta_incomplete simplifies for any argument B and Z and for B a positive integer for any argument A and Z, with the exception of A a negative integer.

     For z=0 and realpart(a)>0 , beta_incomplete has the specific value zero. For z=1 and realpart(b)>0 , beta_incomplete simplifies to the beta function beta(a,b).

     Maxima evaluates beta_incomplete numerically for real and complex values in float or bigfloat precision. For the numerical evaluation an expansion of the incomplete beta function in continued fractions is used.

     When the option variable beta_expand is true, Maxima expands expressions like beta_incomplete(a+n,b,z) and beta_incomplete(a-n,b,z) where n is a positive integer.

     Maxima knows the derivatives of beta_incomplete with respect to the variables A, B and Z and the integral with respect to the variable Z.

     Examples:

     Simplification for A a positive integer:

          (%i1) beta_incomplete(2,b,z);
                                                 b
                                      1 - (1 - z)  (b z + 1)
          (%o1)                       ----------------------
                                            b (b + 1)

     Simplification for B a positive integer:

          (%i2) beta_incomplete(a,2,z);
                                                         a
                                        (a (1 - z) + 1) z
          (%o2)                         ------------------
                                            a (a + 1)

     Simplification for A and B a positive integer:

          (%i3) beta_incomplete(3,2,z);

                                                         3
                                        (3 (1 - z) + 1) z
          (%o3)                         ------------------
                                                12

     A is a negative integer and b<=(-a), Maxima simplifies:

          (%i4) beta_incomplete(-3,1,z);
                                                 1
          (%o4)                              - ----
                                                  3
                                               3 z

     For the specific values z=0 and z=1, Maxima simplifies:

          (%i5) assume(a>0,b>0)$
          (%i6) beta_incomplete(a,b,0);
          (%o6)                                 0
          (%i7) beta_incomplete(a,b,1);
          (%o7)                            beta(a, b)

     Numerical evaluation in float or bigfloat precision:

          (%i8) beta_incomplete(0.25,0.50,0.9);
          (%o8)                          4.594959440269333
          (%i9)  fpprec:25$
          (%i10) beta_incomplete(0.25,0.50,0.9b0);
          (%o10)                    4.594959440269324086971203b0

     For abs(z)>1 beta_incomplete returns a complex result:

          (%i11) beta_incomplete(0.25,0.50,1.7);
          (%o11)              5.244115108584249 - 1.45518047787844 %i

     Results for more general complex arguments:

          (%i14) beta_incomplete(0.25+%i,1.0+%i,1.7+%i);
          (%o14)             2.726960675662536 - .3831175704269199 %i
          (%i15) beta_incomplete(1/2,5/4*%i,2.8+%i);
          (%o15)             13.04649635168716 %i - 5.802067956270001
          (%i16)

     Expansion, when beta_expand is true:

          (%i23) beta_incomplete(a+1,b,z),beta_expand:true;
                                                                 b  a
                             a beta_incomplete(a, b, z)   (1 - z)  z
          (%o23)             -------------------------- - -----------
                                       b + a                 b + a

          (%i24) beta_incomplete(a-1,b,z),beta_expand:true;
                                                                 b  a - 1
                 beta_incomplete(a, b, z) (- b - a + 1)   (1 - z)  z
          (%o24) -------------------------------------- - ---------------
                                 1 - a                         1 - a

     Derivative and integral for beta_incomplete:

          (%i34) diff(beta_incomplete(a, b, z), z);

                                        b - 1  a - 1
          (%o34)                 (1 - z)      z

          (%i35) integrate(beta_incomplete(a, b, z), z);
                        b  a
                 (1 - z)  z
          (%o35) ----------- + beta_incomplete(a, b, z) z
                    b + a
                                                 a beta_incomplete(a, b, z)
                                               - --------------------------
                                                           b + a
          (%i36) factor(diff(%, z));
          (%o36)              beta_incomplete(a, b, z)

     Categories: Gamma and factorial functions

 -- Function: beta_incomplete_regularized (A, B, Z)

     The regularized incomplete beta function (8.17.E2 and eqn 6.6.2, 263), defined as

          beta_incomplete_regularized(a, b, z) =
                                                beta_incomplete(a, b, z)
                                                ------------------------
                                                       beta(a, b)

     As for beta_incomplete this definition is not complete. See https://functions.wolfram.com/GammaBetaErf/BetaRegularized/ for a complete definition of beta_incomplete_regularized.

     beta_incomplete_regularized simplifies A or B a positive integer.

     For z=0 and realpart(a) > 0 , beta_incomplete_regularized has the specific value 0. For z=1 and realpart(b)>0 , beta_incomplete_regularized simplifies to 1.

     Maxima can evaluate beta_incomplete_regularized for real and complex arguments in float and bigfloat precision.

     When beta_expand is true, Maxima expands beta_incomplete_regularized for arguments a+n or a-n, where n is an integer.

     Maxima knows the derivatives of beta_incomplete_regularized with respect to the variables A, B, and Z and the integral with respect to the variable Z.

     Examples:

     Simplification for A or B a positive integer:

          (%i1) beta_incomplete_regularized(2,b,z);
                                                 b
          (%o1)                       1 - (1 - z)  (b z + 1)

          (%i2) beta_incomplete_regularized(a,2,z);
                                                         a
          (%o2)                         (a (1 - z) + 1) z

          (%i3) beta_incomplete_regularized(3,2,z);
                                                         3
          (%o3)                         (3 (1 - z) + 1) z

     For the specific values z=0 and z=1, Maxima simplifies:

          (%i4) assume(a>0,b>0)$
          (%i5) beta_incomplete_regularized(a,b,0);
          (%o5)                                 0
          (%i6) beta_incomplete_regularized(a,b,1);
          (%o6)                                 1

     Numerical evaluation for real and complex arguments in float and bigfloat precision:

          (%i7) beta_incomplete_regularized(0.12,0.43,0.9);
          (%o7)                         .9114011367359802
          (%i8) fpprec:32$
          (%i9) beta_incomplete_regularized(0.12,0.43,0.9b0);
          (%o9)               9.1140113673598075519946998779975b-1
          (%i10) beta_incomplete_regularized(1+%i,3/3,1.5*%i);
          (%o10)             .2865367499935403 %i - 0.122995963334684
          (%i11) fpprec:20$
          (%i12) beta_incomplete_regularized(1+%i,3/3,1.5b0*%i);
          (%o12)      2.8653674999354036142b-1 %i - 1.2299596333468400163b-1

     Expansion, when beta_expand is true:

          (%i13) beta_incomplete_regularized(a+1,b,z);
                                                               b  a
                                                        (1 - z)  z
          (%o13) beta_incomplete_regularized(a, b, z) - ------------
                                                        a beta(a, b)
          (%i14) beta_incomplete_regularized(a-1,b,z);
          (%o14) beta_incomplete_regularized(a, b, z)
                                                               b  a - 1
                                                        (1 - z)  z
                                                   - ----------------------
                                                     beta(a, b) (b + a - 1)

     The derivative and the integral wrt Z:

          (%i15) diff(beta_incomplete_regularized(a,b,z),z);
                                        b - 1  a - 1
                                 (1 - z)      z
          (%o15)                 -------------------
                                     beta(a, b)
          (%i16) integrate(beta_incomplete_regularized(a,b,z),z);
          (%o16) beta_incomplete_regularized(a, b, z) z
                                                                     b  a
                                                              (1 - z)  z
                    a (beta_incomplete_regularized(a, b, z) - ------------)
                                                              a beta(a, b)
                  - -------------------------------------------------------
                                             b + a

     Categories: Gamma and factorial functions

 -- Function: beta_incomplete_generalized (A, B, Z1, Z2)

     The basic definition of the generalized incomplete beta function is

                       z2
                     /
                     [          b - 1  a - 1
                     I   (1 - t)      t      dt
                     ]
                     /
                      z1

     Maxima simplifies beta_incomplete_regularized for A and B a positive integer.

     For realpart(a)>0 and z1=0 or z2=0 , Maxima simplifies beta_incomplete_generalized to beta_incomplete. For realpart(b)>0 and z1=1 or z2=1 , Maxima simplifies to an expression with beta and beta_incomplete.

     Maxima evaluates beta_incomplete_regularized for real and complex values in float and bigfloat precision.

     When beta_expand is true, Maxima expands beta_incomplete_generalized for a+n and a-n, N a positive integer.

     Maxima knows the derivative of beta_incomplete_generalized with respect to the variables A, B, Z1, and Z2 and the integrals with respect to the variables Z1 and Z2.

     Examples:

     Maxima simplifies beta_incomplete_generalized for A and B a positive integer:

          (%i1) beta_incomplete_generalized(2,b,z1,z2);
                             b                      b
                     (1 - z1)  (b z1 + 1) - (1 - z2)  (b z2 + 1)
          (%o1)      -------------------------------------------
                                      b (b + 1)
          (%i2) beta_incomplete_generalized(a,2,z1,z2);

                                        a                      a
                     (a (1 - z2) + 1) z2  - (a (1 - z1) + 1) z1
          (%o2)      -------------------------------------------
                                      a (a + 1)

          (%i3) beta_incomplete_generalized(3,2,z1,z2);
                        2      2                       2      2
                (1 - z1)  (3 z1  + 2 z1 + 1) - (1 - z2)  (3 z2  + 2 z2 + 1)
          (%o3) -----------------------------------------------------------
                                            12

     Simplification for specific values z1=0, z2=0, z1=1, or z2=1:

          (%i4) assume(a > 0, b > 0)$
          (%i5) beta_incomplete_generalized(a,b,z1,0);
          (%o5)                    - beta_incomplete(a, b, z1)

          (%i6) beta_incomplete_generalized(a,b,0,z2);
          (%o6)                    - beta_incomplete(a, b, z2)

          (%i7) beta_incomplete_generalized(a,b,z1,1);
          (%o7)              beta(a, b) - beta_incomplete(a, b, z1)

          (%i8) beta_incomplete_generalized(a,b,1,z2);
          (%o8)              beta_incomplete(a, b, z2) - beta(a, b)

     Numerical evaluation for real arguments in float or bigfloat precision:

          (%i9) beta_incomplete_generalized(1/2,3/2,0.25,0.31);
          (%o9)                        .09638178086368676

          (%i10) fpprec:32$
          (%i10) beta_incomplete_generalized(1/2,3/2,0.25,0.31b0);
          (%o10)               9.6381780863686935309170054689964b-2

     Numerical evaluation for complex arguments in float or bigfloat precision:

          (%i11) beta_incomplete_generalized(1/2+%i,3/2+%i,0.25,0.31);
          (%o11)           - .09625463003205376 %i - .003323847735353769
          (%i12) fpprec:20$
          (%i13) beta_incomplete_generalized(1/2+%i,3/2+%i,0.25,0.31b0);
          (%o13)     - 9.6254630032054178691b-2 %i - 3.3238477353543591914b-3

     Expansion for a+n or a-n, N a positive integer, when beta_expand is true:

          (%i14) beta_expand:true$

          (%i15) beta_incomplete_generalized(a+1,b,z1,z2);

                         b   a           b   a
                 (1 - z1)  z1  - (1 - z2)  z2
          (%o15) -----------------------------
                             b + a
                                a beta_incomplete_generalized(a, b, z1, z2)
                              + -------------------------------------------
                                                   b + a
          (%i16) beta_incomplete_generalized(a-1,b,z1,z2);

                 beta_incomplete_generalized(a, b, z1, z2) (- b - a + 1)
          (%o16) -------------------------------------------------------
                                          1 - a
                                              b   a - 1           b   a - 1
                                      (1 - z2)  z2      - (1 - z1)  z1
                                    - -------------------------------------
                                                      1 - a

     Derivative wrt the variable Z1 and integrals wrt Z1 and Z2:

          (%i17) diff(beta_incomplete_generalized(a,b,z1,z2),z1);
                                         b - 1   a - 1
          (%o17)               - (1 - z1)      z1
          (%i18) integrate(beta_incomplete_generalized(a,b,z1,z2),z1);
          (%o18) beta_incomplete_generalized(a, b, z1, z2) z1
                                            + beta_incomplete(a + 1, b, z1)
          (%i19) integrate(beta_incomplete_generalized(a,b,z1,z2),z2);
          (%o19) beta_incomplete_generalized(a, b, z1, z2) z2
                                            - beta_incomplete(a + 1, b, z2)

     Categories: Gamma and factorial functions

 -- Option variable: beta_expand
     Default value: false

     When beta_expand is true, beta(a,b) and related functions are expanded for arguments like a+n or a-n, where n is an integer.

     See beta for examples.

     Categories: Gamma and factorial functions Simplification flags and variables

 -- Option variable: beta_args_sum_to_integer
     Default value: false

     When beta_args_sum_to_integer is true, Maxima simplifies beta(a,b), when the arguments A and B sum to an integer.

     See beta for examples.

     Categories: Gamma and factorial functions Simplification flags and variables

 -- Function: psi [N](X)

     psi[n](x) is the polygamma function (5.2E2, 5.15, eqn 6.3.1, 258 and eqn 6.4.1, 260) defined by diff(log(gamma(x)), x, n+1)

     Thus, psi[0](X) is the first derivative, psi[1](X) is the second derivative, etc.

     Maxima can compute some exact values for rational args as well for float and bfloat args.  Several variables control what range of rational args psi[n](x) ) will return an exact value, if possible.  See maxpsiposint, maxpsinegint, maxpsifracnum, and maxpsifracdenom That is, x must lie between maxpsinegint and maxpsiposint.  If the absolute value of the fractional part of x is rational and has a numerator less than maxpsifracnum and has a denominator less than maxpsifracdenom, psi[0](x) will return an exact value.

     The function bfpsi in the bffac package can compute numerical values.

          (%i1) psi[0](.25);
          (%o1)                        - 4.227453533376265
          (%i2) psi[0](1/4);
                                                  %pi
          (%o2)                    (- 3 log(2)) - --- - %gamma
                                                   2
          (%i3) float(%);
          (%o3)                        - 4.227453533376265
          (%i4) psi[2](0.75);
          (%o4)                        - 5.30263321633764
          (%i5) psi[2](3/4);
                                             1         3
          (%o5)                         psi (-) + 4 %pi
                                           2 4
          (%i6) float(%);
          (%o6)                        - 5.30263321633764

     Categories: Gamma and factorial functions

 -- Option variable: maxpsiposint
     Default value: 20

     maxpsiposint is the largest positive integer value for which psi[n](x) gives an exact value for rational x.

          (%i1) psi[0](20);
                                       275295799
          (%o1)                        --------- - %gamma
                                       77597520
          (%i2) psi[0](21);
          (%o2)                             psi (21)
                                               0
          (%i3) psi[2](20);
                                1683118856778495358491487
          (%o3)              2 (------------------------- - zeta(3))
                                1401731326612193601024000
          (%i4) psi[2](21);
          (%o4)                            psi (21)
                                                2

     Categories: Gamma and factorial functions

 -- Option variable: maxpsinegint
     Default value: -10

     maxpsinegint is the most negative value for which psi[0](x) will try to compute an exact value for rational x.  That is if x is less than maxpsinegint, psi[n](x) will not return simplified answer, even if it could.

          (%i1) psi[0](-100/9);
                                                  100
          (%o1)                            psi (- ---)
                                              0    9
          (%i2) psi[0](-100/11);
                                   100 %pi         1     5231385863539
          (%o2)            %pi cot(-------) + psi (--) + -------------
                                     11          0 11    381905105400

          (%i3) psi[2](-100/9);
                                                  100
          (%o3)                            psi (- ---)
                                              2    9
          (%i4) psi[2](-100/11);
                     3     100 %pi     2 100 %pi         1
          (%o4) 2 %pi  cot(-------) csc (-------) + psi (--)
                             11            11          2 11
                                                   74191313259470963498957651385614962459
                                                 + --------------------------------------
                                                    27850718060013605318710152732000000

     Categories: Gamma and factorial functions

 -- Option variable: maxpsifracnum
     Default value: 6

     Let x be a rational number of the form p/q. If p is greater than maxpsifracnum, then psi[0](x) will not try to return a simplified value.

          (%i1) psi[0](3/4);
                                                  %pi
          (%o1)                    (- 3 log(2)) + --- - %gamma
                                                   2
          (%i2) psi[2](3/4);
                                             1         3
          (%o2)                         psi (-) + 4 %pi
                                           2 4
          (%i3) maxpsifracnum:2;
          (%o3)                                 2
          (%i4) psi[0](3/4);
                                                  3
          (%o4)                              psi (-)
                                                0 4
          (%i5) psi[2](3/4);
                                             1         3
          (%o5)                         psi (-) + 4 %pi
                                           2 4

     Categories: Gamma and factorial functions

 -- Option variable: maxpsifracdenom
     Default value: 6

     Let x be a rational number of the form p/q. If q is greater than maxpsifracdenom, then psi[0](X) will not try to return a simplified value.

          (%i1) psi[0](3/4);
                                                  %pi
          (%o1)                    (- 3 log(2)) + --- - %gamma
                                                   2
          (%i2) psi[2](3/4);
                                             1         3
          (%o2)                         psi (-) + 4 %pi
                                           2 4
          (%i3) maxpsifracdenom:2;
          (%o3)                                 2
          (%i4) psi[0](3/4);
                                                  3
          (%o4)                              psi (-)
                                                0 4
          (%i5) psi[2](3/4);
                                             1         3
          (%o5)                         psi (-) + 4 %pi
                                           2 4

     Categories: Gamma and factorial functions

 -- Function: makefact (EXPR)
     Transforms instances of binomial, gamma, and beta functions in EXPR into factorials.

     See also makegamma

          (%i1) makefact(binomial(n,k));
                                                n!
          (%o1)                             -----------
                                            k! (n - k)!
          (%i2) makefact(gamma(x));
          (%o2)                              (x - 1)!
          (%i3) makefact(beta(a,b));
                                         (a - 1)! (b - 1)!
          (%o3)                          -----------------
                                           (b + a - 1)!

     Categories: Gamma and factorial functions

 -- Function: numfactor (EXPR)
     Returns the numerical factor multiplying the expression EXPR, which should be a single term.

     content returns the greatest common divisor (gcd) of all terms in a sum.

          (%i1) gamma (7/2);
                                    15 sqrt(%pi)
          (%o1)                     ------------
                                         8
          (%i2) numfactor (%);
                                         15
          (%o2)                          --
                                         8

     Categories: Expressions

Exponential Integrals
=====================

The Exponential Integral and related functions are defined in Abramowitz and Stegun, Handbook of Mathematical Functions, Chapter 5, 228.

 -- Function: expintegral_e1 (Z)
     The Exponential Integral E1(z) defined as

     integrate(exp(-t)/t, t, z, inf)

     with abs(arg z) < %pi . (eqn 5.1.1, 228) and (6.2E2)

     This can be written in terms of other functions.  See expintrep for examples. Categories: Exponential Integrals Special functions

 -- Function: expintegral_ei (X)
     The Exponential Integral Ei(x) defined as

     -integrate(%e^(-t)/t, t, -x, inf) = integrate(%e^(t)/t, t, -inf, x)

     with x real and x > 0. (eqn 5.1.2, 228) and (6.2E5)

     This can be written in terms of other functions.  See expintrep for examples.

     Categories: Exponential Integrals Special functions

 -- Function: expintegral_li (X)
     The Exponential Integral li(x) defined as

     integrate(1/ln(t), t, 0, x)

     with x real and x > 1. (eqn 5.1.3, 228) and (6.2E8)

     This can be written in terms of other functions.  See expintrep for examples.

     Categories: Exponential Integrals Special functions

 -- Function: expintegral_e (N,Z)
     The Exponential Integral En(z) (eqn 5.1.4, 228) defined as

     integrate(exp(-z*t)/t^n, t, 1, inf)

     with realpart(z) > 1 and n a non-negative integer.

     For half-integral orders, this can be written in terms of erfc or erf.  See expintexpand for examples.

     Categories: Exponential Integrals Special functions

 -- Function: expintegral_si (Z)
     The Exponential Integral Si(z) (eqn 5.2.1, 231 and 6.2#E9) defined as

     integrate(sin(t)/t, t, 0, z)

     This can be written in terms of other functions.  See expintrep for examples.

     Categories: Exponential Integrals Special functions

 -- Function: expintegral_ci (Z)
     The Exponential Integral Ci(z) (eqn 5.2.2, 231 and 6.2#E13) defined as

     %gamma + log(z) + integrate((cos(t) - 1)/t, t, 0, z)

     with abs(arg z) < %pi .

     This can be written in terms of other functions.  See expintrep for examples.

     Categories: Exponential Integrals Special functions

 -- Function: expintegral_shi (Z)
     The Exponential Integral Shi(z) (eqn 5.2.3, 231 and 6.2#E15) defined as

     integrate(sinh(t)/t, t, 0, z)

     This can be written in terms of other functions.  See expintrep for examples.

     Categories: Exponential Integrals Special functions

 -- Function: expintegral_chi (Z)
     The Exponential Integral Chi(z) (eqn 5.2.4, 231 and 6.2#E16) defined as

     %gamma + log(z) + integrate((cosh(t) - 1)/t, t, 0, z)

     with abs(arg z) < %pi .

     This can be written in terms of other functions.  See expintrep for examples.

     Categories: Exponential Integrals Special functions

 -- Option variable: expintrep
     Default value: false

     Change the representation of one of the exponential integrals, expintegral_e(m, z), expintegral_e1, or expintegral_ei to an equivalent form if possible.

     Possible values for expintrep are false, gamma_incomplete, expintegral_e1, expintegral_ei, expintegral_li, expintegral_trig, or expintegral_hyp.

     false means that the representation is not changed.  Other values indicate the representation is to be changed to use the function specified where expintegral_trig means expintegral_si, expintegral_ci; and expintegral_hyp means expintegral_shi or expintegral_chi.

     Here are some examples for expintrep set to gamma_incomplete:

          (%i1) expintrep:'gamma_incomplete;
          (%o1)                          gamma_incomplete
          (%i2) expintegral_e1(z);
          (%o2)                       gamma_incomplete(0, z)
          (%i3) expintegral_ei(z);
          (%o3)            log(z) - log(- z) - gamma_incomplete(0, - z)
          (%i4) expintegral_li(z);
          (%o4)     log(log(z)) - log(- log(z)) - gamma_incomplete(0, - log(z))
          (%i5) expintegral_e(n,z);
                                                             n - 1
          (%o5)                  gamma_incomplete(1 - n, z) z
          (%i6) expintegral_si(z);
          (%o6) (%i ((- log(%i z)) + log(- %i z) - gamma_incomplete(0, %i z)
                                                        + gamma_incomplete(0, - %i z)))/2
          (%i7) expintegral_ci(z);
          (%o7) log(z) - (log(%i z) + log(- %i z) + gamma_incomplete(0, %i z)
                                                         + gamma_incomplete(0, - %i z))/2
          (%i8) expintegral_shi(z);
                log(z) - log(- z) + gamma_incomplete(0, z) - gamma_incomplete(0, - z)
          (%o8) ---------------------------------------------------------------------
                                                  2
          (%i9) expintegral_chi(z);
          (%o9)
                (- log(z)) + log(- z) + gamma_incomplete(0, z) + gamma_incomplete(0, - z)
              - -------------------------------------------------------------------------
                                                    2

     For expintrep set to expintegral_e1:

          (%i1) expintrep:'expintegral_e1;
          (%o1)                           expintegral_e1
          (%i2) expintegral_ei(z);
          (%o2)               log(z) - log(- z) - expintegral_e1(- z)
          (%i3) expintegral_li(z);
          (%o3)       log(log(z)) - log(- log(z)) - expintegral_e1(- log(z))
          (%i4) expintegral_e(n,z);
          (%o4)                         expintegral_e(n, z)
          (%i5) expintegral_si(z);
          (%o5) (%i ((- log(%i z)) - expintegral_e1(%i z) + log(- %i z)
                                                             + expintegral_e1(- %i z)))/2
          (%i6) expintegral_ci(z);
          (%o6) log(z)
                    log(- %i z) (expintegral_e1(%i z) + expintegral_e1(- %i z)) log(%i z)
                  - ---------------------------------------------------------------------
                                                      2
          (%i7) expintegral_shi(z);
                    log(z) + expintegral_e1(z) - log(- z) - expintegral_e1(- z)
          (%o7)     -----------------------------------------------------------
                                                 2
          (%i8) expintegral_chi(z);
                   (- log(z)) + expintegral_e1(z) + log(- z) + expintegral_e1(- z)
          (%o8)  - ---------------------------------------------------------------
                                                  2

     For expintrep set to expintegral_ei:

          (%i1) expintrep:'expintegral_ei;
          (%o1)                           expintegral_ei
          (%i2) expintegral_e1(z);
                                                    1
                                   log(- z) - log(- -)
                                                    z
          (%o2)       (- log(z)) + ------------------- - expintegral_ei(- z)
                                            2
          (%i3) expintegral_ei(z);
          (%o3)                          expintegral_ei(z)
          (%i4) expintegral_li(z);
          (%o4)                       expintegral_ei(log(z))
          (%i5) expintegral_e(n,z);
          (%o5)                         expintegral_e(n, z)
          (%i6) expintegral_si(z);
          (%o6) (%i (log(%i z) + 2 (expintegral_ei(- %i z) - expintegral_ei(%i z))
                                                                      %i          %i
                                                  - log(- %i z) + log(--) - log(- --)))/4
                                                                      z           z
          (%i7) expintegral_ci(z);
          (%o7) ((- log(%i z)) + 2 (expintegral_ei(%i z) + expintegral_ei(- %i z))
                                                              %i          %i
                                          - log(- %i z) + log(--) + log(- --))/4 + log(z)
                                                              z           z
          (%i8) expintegral_shi(z);
          (%o8) ((- 2 log(z)) + 2 (expintegral_ei(z) - expintegral_ei(- z)) + log(- z)
                                                                                    1
                                                                            - log(- -))/4
                                                                                    z
          (%i9) expintegral_chi(z);
          (%o9)
                                                                                       1
             2 log(z) + 2 (expintegral_ei(z) + expintegral_ei(- z)) - log(- z) + log(- -)
                                                                                       z
             ----------------------------------------------------------------------------
                                                  4

     For expintrep set to expintegral_li:

          (%i1) expintrep:'expintegral_li;
          (%o1)                           expintegral_li
          (%i2) expintegral_e1(z);
                                                                           1
                                                          log(- z) - log(- -)
                                         - z                               z
          (%o2)      (- expintegral_li(%e   )) - log(z) + -------------------
                                                                   2
          (%i3) expintegral_ei(z);
                                                         z
          (%o3)                         expintegral_li(%e )
          (%i4) expintegral_li(z);
          (%o4)                          expintegral_li(z)
          (%i5) expintegral_e(n,z);
          (%o5)                         expintegral_e(n, z)
          (%i6) expintegral_si(z);
                                       %i z                     - %e z    %pi signum(z)
                  %i (expintegral_li(%e    ) - expintegral_li(%e      ) - -------------)
                                                                                2
          (%o6) - ----------------------------------------------------------------------
                                                    2
          (%i7) expintegral_ci(z);
                                  %i z                     - %i z
                 expintegral_li(%e    ) + expintegral_li(%e      )
          (%o7)  ------------------------------------------------- - signum(z) + 1
                                         2
          (%i8) expintegral_shi(z);
                                             z                     - z
                            expintegral_li(%e ) - expintegral_li(%e   )
          (%o8)             -------------------------------------------
                                                 2
          (%i9) expintegral_chi(z);
                                             z                     - z
                            expintegral_li(%e ) + expintegral_li(%e   )
          (%o9)             -------------------------------------------
                                                 2

     For expintrep set to expintegral_trig:

          (%i1) expintrep:'expintegral_trig;
          (%o1)                          expintegral_trig
          (%i2) expintegral_e1(z);
          (%o2) log(%i z) - %i expintegral_si(%i z) - expintegral_ci(%i z) - log(z)
          (%i3) expintegral_ei(z);
          (%o3) (- log(%i z)) - %i expintegral_si(%i z) + expintegral_ci(%i z) + log(z)
          (%i4) expintegral_li(z);
          (%o4) (- log(%i log(z))) - %i expintegral_si(%i log(z))
                                                + expintegral_ci(%i log(z)) + log(log(z))
          (%i5) expintegral_e(n,z);
          (%o5)                         expintegral_e(n, z)
          (%i6) expintegral_si(z);
          (%o6)                          expintegral_si(z)
          (%i7) expintegral_ci(z);
          (%o7)                          expintegral_ci(z)
          (%i8) expintegral_shi(z);
          (%o8)                      - %i expintegral_si(%i z)
          (%i9) expintegral_chi(z);
          (%o9)            (- log(%i z)) + expintegral_ci(%i z) + log(z)

     For expintrep set to expintegral_hyp:

          (%i1) expintrep:'expintegral_hyp;
          (%o1)                           expintegral_hyp
          (%i2) expintegral_e1(z);
          (%o2)               expintegral_shi(z) - expintegral_chi(z)
          (%i3) expintegral_ei(z);
          (%o3)               expintegral_shi(z) + expintegral_chi(z)
          (%i4) expintegral_li(z);
          (%o4)          expintegral_shi(log(z)) + expintegral_chi(log(z))
          (%i5) expintegral_e(n,z);
          (%o5)                         expintegral_e(n, z)
          (%i6) expintegral_si(z);
          (%o6)                     - %i expintegral_shi(%i z)
          (%i7) expintegral_ci(z);
          (%o7)           (- log(%i z)) + expintegral_chi(%i z) + log(z)
          (%i8) expintegral_shi(z);
          (%o8)                         expintegral_shi(z)
          (%i9) expintegral_chi(z);
          (%o9)                         expintegral_chi(z)

     Categories: Exponential Integrals

 -- Option variable: expintexpand
     Default value: false

     Expand expintegral_e(n,z) for half integral values in terms of erfc or erf and for positive integers in terms of expintegral_ei.

          (%i1) expintegral_e(1/2,z);
                                                      1
          (%o1)                         expintegral_e(-, z)
                                                      2
          (%i2) expintegral_e(1,z);
          (%o2)                         expintegral_e(1, z)
          (%i3) expintexpand:true;
          (%o3)                                true
          (%i4) expintegral_e(1/2,z);
                                      sqrt(%pi) erfc(sqrt(z))
          (%o4)                       -----------------------
                                              sqrt(z)
          (%i5) expintegral_e(1,z);
                                         1
                                   log(- -) - log(- z)
                                         z
          (%o5)       (- log(z)) - ------------------- - expintegral_ei(- z)
                                            2

     Categories: Exponential Integrals

Error Function
==============

The Error function and related functions are defined in Abramowitz and Stegun, Handbook of Mathematical Functions, Chapter 7, 297 and (7)

 -- Function: erf (Z)

     The Error Function erf(z): erf(z) = 2/sqrt(%pi)*integrate(exp(-t^2), t, 0, z)

     (eqn 7.1.1, 297) and (7.2.E1).

     See also flag erfflag  This can also be expressed in terms of a hypergeometric function.  See hypergeometric_representation. Categories: Special functions

 -- Function: erfc (Z)
     The Complementary Error Function erfc(z): erfc(z) = 1-erf(z)

     (eqn 7.1.2, 297) and (7.2.E2).

     This can also be expressed in terms of a hypergeometric function.  See hypergeometric_representation. Categories: Special functions

 -- Function: erfi (Z)
     The Imaginary Error Function. erfi(z) = -%i*erf(%i*z)

     Categories: Special functions

 -- Function: erf_generalized (Z1,Z2)
     Generalized Error function Erf(z1,z2): erf(z) = 2/sqrt(%pi)*integrate(exp(-t^2), t, z1, z2)

     This can also be expressed in terms of a hypergeometric function.  See hypergeometric_representation. Categories: Special functions

 -- Function: fresnel_c (Z)
     The Fresnel Integral

     C(z) = integrate(cos((%pi/2)*t^2),t,0,z)

     (eqn 7.3.1, 300) and (7.2.E7).

     The simplification fresnel_c(-x) = -fresnel_c(x) is applied when flag trigsign is true.

     The simplification fresnel_c(%i*x) =  %i*fresnel_c(x) is applied when flag %iargs is true.

     See flags erf_representation and hypergeometric_representation. Categories: Special functions

 -- Function: fresnel_s (Z)
     The Fresnel Integral S(z) = integrate(sin((%pi/2)*t^2),t,0,z)

     (eqn 7.3.2, 300) and (7.2.E8).

     The simplification fresnel_s(-x) = -fresnel_s(x) is applied when flag trigsign is true.

     The simplification fresnel_s(%i*x) =  -%i*fresnel_s(x) is applied when flag %iargs is true.

     See flags erf_representation and hypergeometric_representation. Categories: Special functions

 -- Option variable: erf_representation
     Default value: false

     erf_representation controls how the error functions are represented.  It must be set to one of false, erf, erfc, or erfi.  When set to false, the error functions are not modified.  When set to erf, all error functions (erfc, erfi, erf_generalized, fresnel_s and fresnel_c) are converted to erf functions.  Similary, erfc converts error functions to erfc.  Finally erfi converts the functions to erfi.

     Converting to erf:

          (%i1) erf_representation:erf;
          (%o1)                                true
          (%i2) erfc(z);
          (%o2)                               erfc(z)
          (%i3) erfi(z);
          (%o3)                               erfi(z)
          (%i4) erf_generalized(z1,z2);
          (%o4)                          erf(z2) - erf(z1)
          (%i5) fresnel_c(z);
                              sqrt(%pi) (%i + 1) z           sqrt(%pi) (1 - %i) z
                (1 - %i) (erf(--------------------) + %i erf(--------------------))
                                       2                              2
          (%o5) -------------------------------------------------------------------
                                                 4
          (%i6) fresnel_s(z);
                              sqrt(%pi) (%i + 1) z           sqrt(%pi) (1 - %i) z
                (%i + 1) (erf(--------------------) - %i erf(--------------------))
                                       2                              2
          (%o6) -------------------------------------------------------------------
                                                 4

     Converting to erfc:

          (%i1) erf_representation:erfc;
          (%o1)                                erfc
          (%i2) erf(z);
          (%o2)                             1 - erfc(z)
          (%i3) erfc(z);
          (%o3)                               erfc(z)
          (%i4) erf_generalized(z1,z2);
          (%o4)                         erfc(z1) - erfc(z2)
          (%i5) fresnel_s(c);
                                   sqrt(%pi) (%i + 1) c
          (%o5) ((%i + 1) ((- erfc(--------------------))
                                            2
                                                           sqrt(%pi) (1 - %i) c
                                            - %i (1 - erfc(--------------------)) + 1))/4
                                                                    2
          (%i6) fresnel_c(c);
                                   sqrt(%pi) (%i + 1) c
          (%o6) ((1 - %i) ((- erfc(--------------------))
                                            2
                                                           sqrt(%pi) (1 - %i) c
                                            + %i (1 - erfc(--------------------)) + 1))/4
                                                                    2

     Converting to erfc:

          (%i1) erf_representation:erfi;
          (%o1)                                erfi
          (%i2) erf(z);
          (%o2)                           - %i erfi(%i z)
          (%i3) erfc(z);
          (%o3)                          %i erfi(%i z) + 1
          (%i4) erfi(z);
          (%o4)                               erfi(z)
          (%i5) erf_generalized(z1,z2);
          (%o5)                   %i erfi(%i z1) - %i erfi(%i z2)
          (%i6) fresnel_s(z);
                                      sqrt(%pi) %i (%i + 1) z
          (%o6) ((%i + 1) ((- %i erfi(-----------------------))
                                                 2
                                                             sqrt(%pi) (1 - %i) %i z
                                                      - erfi(-----------------------)))/4
                                                                        2
          (%i7) fresnel_c(z);
          (%o7)
                             sqrt(%pi) (1 - %i) %i z            sqrt(%pi) %i (%i + 1) z
              (1 - %i) (erfi(-----------------------) - %i erfi(-----------------------))
                                        2                                  2
              ---------------------------------------------------------------------------
                                                   4

 -- Option variable: hypergeometric_representation
     Default value: false

     Enables transformation to a Hypergeometric representation for fresnel_s and fresnel_c and other error functions.

          (%i1) hypergeometric_representation:true;
          (%o1)                                true
          (%i2) fresnel_s(z);
                                                                2  4
                                               3    3  7     %pi  z    3
                           %pi hypergeometric([-], [-, -], - -------) z
                                               4    2  4       16
          (%o2)            ---------------------------------------------
                                                 6
          (%i3) fresnel_c(z);
                                                              2  4
                                             1    1  5     %pi  z
          (%o3)              hypergeometric([-], [-, -], - -------) z
                                             4    2  4       16
          (%i4) erf(z);
                                                  1    3      2
                                2 hypergeometric([-], [-], - z ) z
                                                  2    2
          (%o4)                 ----------------------------------
                                            sqrt(%pi)
          (%i5) erfi(z);
                                                   1    3    2
                                 2 hypergeometric([-], [-], z ) z
                                                   2    2
          (%o5)                  --------------------------------
                                            sqrt(%pi)
          (%i6) erfc(z);
                                                    1    3      2
                                  2 hypergeometric([-], [-], - z ) z
                                                    2    2
          (%o6)               1 - ----------------------------------
                                              sqrt(%pi)
          (%i7) erf_generalized(z1,z2);
                                  1    3       2
                2 hypergeometric([-], [-], - z2 ) z2
                                  2    2
          (%o7) ------------------------------------
                             sqrt(%pi)
                                                                       1    3       2
                                                     2 hypergeometric([-], [-], - z1 ) z1
                                                                       2    2
                                                   - ------------------------------------
                                                                  sqrt(%pi)

Struve Functions
================

The Struve functions are defined in Abramowitz and Stegun, Handbook of Mathematical Functions, Chapter 12, 496 and (11). The Struve Function H[v](z) is a particular solution of the differential equation: z^2*diff(w,z,2) + z*diff(w,z)+(z^2-v^2)*w = 4*(z/2)^(v+1)/(sqrt(%pi)*gamma(v+1/2))

which has the general soution w = a*bessel_j(v,z) + b*bessel_y(v,z) + struve_h(v,z)

 -- Function: struve_h (V, Z)
     The Struve Function H of order v and argument z:

     struve_h(v,z) = (z/2)^(v+1)*sum((-1)^k*(z/2)^(2*k)/(gamma(k+3/2)*gamma(k+v+3/2)), k, 0, inf)

     (eqn 12.1.3, 496) and (11.2.E1).

     When besselexpand is true, struve_h is expanded in terms of elementary functions when the order v is half of an odd integer. See besselexpand

     Categories: Special functions

 -- Function: struve_l (V, Z)
     The Modified Struve Function L of order v and argument z: struve_l(v,z) = %i*exp(-%i*v*%pi/2)*struve_h(v,z)

     (eqn 12.2.1, 498) and (11.2.E2).

     When besselexpand is true, struve_l is expanded in terms of elementary functions when the order v is half of an odd integer. See besselexpand

     Categories: Special functions

Hypergeometric Functions
========================

The Hypergeometric Functions are defined in Abramowitz and Stegun, Handbook of Mathematical Functions, Chapters 13, 504 and 15, 555.

Maxima has very limited knowledge of these functions.  They can be returned from function hgfred.

 -- Function: %m [K,U] (Z)
     Whittaker M function (eqn 13.1.32, 505):

     %m[k,u](z) = exp(-z/2)*z^(1/2+u)*M(1/2+u-k,1+2*u,z)

     where M(a,b,z) is Kummer's solution of the confluent hypergeometric equation.

     This can also be expressed by the series (13.14.E6): %m[k,u](z) = %e^(-z/2)*z^(1/2+u)*sum(pochhammer(1/2+u-k,s)/(pochhammer(1+2*u,s)*s!)*z^s, s, 0, inf)

     Categories: Special functions

 -- Function: %w [K,U] (Z)
     Whittaker W function (eqn 13.1.33, 505): %w[k,v](z) = exp(-z/2)*z^(1/2+u)*U(1/2+u-k, 1+2*u, z)

     where U(a,b,z) is Kummer's second solution of the confluent hypergeometric equation.

     Categories: Special functions

 -- Function: %f [P,Q] ([A],[B],Z)
     The pFq(a_1,a_2,...,a_p;b_1,b_2,...,b_q;z) hypergeometric function, where A a list of length P and B a list of length Q. Categories: Bessel functions Special functions

 -- Function: hypergeometric ([A1, ..., AP],[B1, ... ,BQ], x)
     The hypergeometric function. Unlike Maxima's %f hypergeometric function, the function hypergeometric is a simplifying function; also, hypergeometric supports complex double and big floating point evaluation. For the Gauss hypergeometric function, that is p = 2 and q = 1, floating point evaluation outside the unit circle is supported, but in general, it is not supported.

     When the option variable expand_hypergeometric is true (default is false) and one of the arguments a1 through ap is a negative integer (a polynomial case), hypergeometric returns an expanded polynomial.

     Examples:

          (%i1)  hypergeometric([],[],x);
          (%o1) %e^x

     Polynomial cases automatically expand when expand_hypergeometric is true:

          (%i2) hypergeometric([-3],[7],x);
          (%o2) hypergeometric([-3],[7],x)

          (%i3) hypergeometric([-3],[7],x), expand_hypergeometric : true;
          (%o3) -x^3/504+3*x^2/56-3*x/7+1

     Both double float and big float evaluation is supported:

          (%i4) hypergeometric([5.1],[7.1 + %i],0.42);
          (%o4)       1.346250786375334 - 0.0559061414208204 %i
          (%i5) hypergeometric([5,6],[8], 5.7 - %i);
          (%o5)     .007375824009774946 - .001049813688578674 %i
          (%i6) hypergeometric([5,6],[8], 5.7b0 - %i), fpprec : 30;
          (%o6) 7.37582400977494674506442010824b-3
                                    - 1.04981368857867315858055393376b-3 %i

 -- Function: hypergeometric_simp (E)

     hypergeometric_simp simplifies hypergeometric functions by applying hgfred to the arguments of any hypergeometric functions in the expression E.

     Only instances of hypergeometric are affected; any %f, %w, and %m in the expression E are not affected. Any unsimplified hypergeometric functions are returned unchanged (instead of changing to %f as hgfred would).

     load("hypergeometric"); loads this function.

     See also hgfred

     Examples:

          (%i1) load ("hypergeometric") $
          (%i2) foo : [hypergeometric([1,1], [2], z), hypergeometric([1/2], [1], z)];
          (%o2) [hypergeometric([1, 1], [2], z),
                                                               1
                                               hypergeometric([-], [1], z)]
                                                               2
          (%i3) hypergeometric_simp (foo);
                           log(1 - z)              z    z/2
          (%o3)         [- ----------, bessel_i(0, -) %e   ]
                               z                   2
          (%i4) bar : hypergeometric([n], [m], z + 1);
          (%o4)            hypergeometric([n], [m], z + 1)
          (%i5) hypergeometric_simp (bar);
          (%o5)            hypergeometric([n], [m], z + 1)

     Categories: Hypergeometric functions Simplification functions Special functions

 -- Function: hgfred (A, B, T)

     Simplify the generalized hypergeometric function in terms of other, simpler, forms.  A is a list of numerator parameters and B is a list of the denominator parameters.

     If hgfred cannot simplify the hypergeometric function, it returns an expression of the form %f[p,q]([a], [b], x) where P is the number of elements in A, and Q is the number of elements in B.  This is the usual pFq generalized hypergeometric function.

          (%i1) assume(not(equal(z,0)));
          (%o1)                          [notequal(z, 0)]
          (%i2) hgfred([v+1/2],[2*v+1],2*%i*z);

                               v/2                               %i z
                              4    bessel_j(v, z) gamma(v + 1) %e
          (%o2)               ---------------------------------------
                                                 v
                                                z
          (%i3) hgfred([1,1],[2],z);

                                             log(1 - z)
          (%o3)                            - ----------
                                                 z
          (%i4) hgfred([a,a+1/2],[3/2],z^2);

                                         1 - 2 a          1 - 2 a
                                  (z + 1)        - (1 - z)
          (%o4)                   -------------------------------
                                           2 (1 - 2 a) z

     It can be beneficial to load orthopoly too as the following example shows.  Note that L is the generalized Laguerre polynomial.

          (%i5) load("orthopoly")$
          (%i6) hgfred([-2],[a],z);

                                              (a - 1)
                                           2 L       (z)
                                              2
          (%o6)                            -------------
                                             a (a + 1)

          (%i7) ev(%);

                                            2
                                           z        2 z
          (%o7)                         --------- - --- + 1
                                        a (a + 1)    a

Parabolic Cylinder Functions
============================

The Parabolic Cylinder Functions are defined in Abramowitz and Stegun, Handbook of Mathematical Functions, Chapter 19, 686.

Maxima has very limited knowledge of these functions.  They can be returned from function specint.

 -- Function: parabolic_cylinder_d (V, Z)
     The parabolic cylinder function parabolic_cylinder_d(v,z). (eqn 19.3.1, 687).

     The solution of the Weber differential equation diff(y(z), z, 2) + (v+1/2-z^2/4)*y(z) = 0

     has two independent solutions, one of which is parabolic_cylinder_d(v,z) , the parabolic cylinder d function.

     Function specint can return expressions containing parabolic_cylinder_d(v,z) if the option variable prefer_d is true.

     Categories: Special functions

Functions and Variables for Special Functions
=============================================

 -- Function: lambert_w (Z)
     The principal branch of Lambert's W function W(z) (4.13), the solution of z = W(z) * exp(W(z))

     Categories: Special functions

 -- Function: generalized_lambert_w (K, Z)
     The K-th branch of Lambert's W function W(z) (4.13), the solution of z = W(z) * exp(W(z)) .

     The principal branch, denoted Wp(z) in DLMF, is lambert_w(z) = generalized_lambert_w(0,z).

     The other branch with real values, denoted Wm(z) in DLMF, is generalized_lambert_w(-1,z). Categories: Special functions

 -- Function: kbateman [V] (X)
     The Bateman k function

     kbateman[v](x) = (2/%pi) integrate(cos(x*tan(t)-v*t),t,0,%pi/2)

     It is a special case of the confluent hypergeometric function. Maxima can calculate the Laplace transform of kbateman using laplace or specint but has no other knowledge of this function.

     Categories: Special functions

 -- Function: nzeta (Z)
     The Plasma Dispersion Function nzeta(z) = %i*sqrt(%pi)*exp(-z^2)*(1-erf(-%i*z))

     Categories: Special functions

 -- Function: nzetar (Z)
     Returns realpart(nzeta(z)). Categories: Special functions

 -- Function: nzetai (Z)
     Returns imagpart(nzeta(z)). Categories: Special functions