| blob | f9fb3ab7209a0db609c30b34fcc897edf94614e6 |
1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2 ;;;
3 ;;; Double Factorial, Incomplete Gamma function, ...
4 ;;;
5 ;;; This file will be extended with further functions related to the
6 ;;; Factorial and Gamma functions.
7 ;;;
8 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
9 ;;;
10 ;;; This file contains the following Maxima User functions:
11 ;;;
12 ;;; double_factorial(z)
13 ;;;
14 ;;; gamma_incomplete(a,z)
15 ;;; gamma_incomplete_generalized(a,z1,z2)
16 ;;; gamma_incomplete_regularized(a,z)
17 ;;;
18 ;;; log_gamma(z)
19 ;;;
20 ;;; erf(z)
21 ;;; erfc(z)
22 ;;; erfi(z)
23 ;;; erf_generalized(z1,z2)
24 ;;;
25 ;;; inverse_erf(z)
26 ;;; inverse_erfc(z)
27 ;;;
28 ;;; fresnel_s(z)
29 ;;; fresnel_c(z)
30 ;;;
31 ;;; beta_incomplete(a,b,z)
32 ;;; beta_incomplete_generalized(a,b,z1,z2)
33 ;;; beta_incomplete_regularized(a,b,z)
34 ;;;
35 ;;; Maxima User variable:
36 ;;;
37 ;;; $factorial_expand - Allows argument simplificaton for expressions like
38 ;;; double_factorial(n-1) and double_factorial(2*k+n)
39 ;;; $beta_expand - Switch on further expansions for the Beta functions
40 ;;;
41 ;;; $erf_representation - When T erfc, erfi, erf_generalized, fresnel_s
42 ;;; and fresnel_c are transformed to erf.
43 ;;; $erf_%iargs - Enable simplification of Erf and Erfi for
44 ;;; imaginary arguments
45 ;;; $hypergeometric_representation
46 ;;; - Enables transformation to a Hypergeometric
47 ;;; representation for fresnel_s and fresnel_c
48 ;;;
49 ;;; Maxima User variable (not defined in this file):
50 ;;;
51 ;;; $factlim - biggest integer for numerically evaluation
52 ;;; of the Double factorial
53 ;;; $gamma_expand - Expansion of the Gamma und Incomplete Gamma
54 ;;; function for some special cases
55 ;;;
56 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
57 ;;; This library is free software; you can redistribute it and/or modify it
58 ;;; under the terms of the GNU General Public License as published by the
59 ;;; Free Software Foundation; either version 2 of the License, or (at
60 ;;; your option) any later version.
61 ;;;
62 ;;; This library is distributed in the hope that it will be useful, but
63 ;;; WITHOUT ANY WARRANTY; without even the implied warranty of
64 ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
65 ;;; Library General Public License for more details.
66 ;;;
67 ;;; You should have received a copy of the GNU General Public License along
68 ;;; with this library; if not, write to the Free Software
69 ;;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
70 ;;;
71 ;;; Copyright (C) 2008 Dieter Kaiser
72 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
76 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
79 "When not false, error functions are transformed to erf, erfc, or erfi
80 depending on the value of $erf_representation."
81 ;; Need to use the noun form of these.
87 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
88 ;;; The following functions test if numerical evaluation has to be done.
89 ;;; The functions should help to test for numerical evaluation more consistent
90 ;;; and without complicated conditional tests including more than one or two
91 ;;; arguments.
92 ;;;
93 ;;; The functions take a list of arguments. All arguments have to be a CL or
94 ;;; Maxima number. If all arguments are numbers we have two cases:
95 ;;; 1. $numer is T we return T. The function has to be evaluated numerically.
96 ;;; 2. One of the args is a float or a bigfloat. Evaluate numerically.
97 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
99 ;;; Test for numerically evaluation in float precision
109 ;;; Test for numerically evaluation in complex float precision
112 "Determine if ARGS consists of numerical values by determining if
113 the real and imaginary parts of each arg are nuemrical (but not
114 bigfloats). A non-NIL result is returned if at least one of args is
115 a floating-point value or if numer is true. If the result is
116 non-NIL, it is a list of the arguments reduced via COMPLEX-NUMBER-P"
123 ;; Always save the result from complex-number-p. But for backward
124 ;; compatibility, only set the flag if any item is a float.
128 ;; Return the values in the same order as the args!
131 ;;; Test for numerically evaluation in bigfloat precision
142 ;;; Test for numerically evaluation in complex bigfloat precision
145 "Determine if ARGS consists of numerical values by determining if
146 the real and imaginary parts of each arg are nuemrical (including
147 bigfloats). A non-NIL result is returned if at least one of args is
148 a floating-point value or if numer is true. If the result is
149 non-NIL, it is a list of the arguments reduced via COMPLEX-NUMBER-P."
157 ;; Always save the result from complex-number-p. But for backward
158 ;; compatibility, only set the flag if any item is a bfloat.
163 ;; Return the values in the same order as the args!
166 ;;; Test for numerical evaluation in any precision, real or complex.
173 ;;; Check for an integer or a float or bigfloat representation. When we
174 ;;; have a float or bigfloat representation return the integer value.
183 val)
186 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
188 ;;; The changes to the parser to connect the operator !! to double_factorial(z)
190 ;(def-mheader |$!!| (%double_factorial))
192 ;(def-led (|$!!| 160.) (op left)
193 ; (list '$expr
194 ; (mheader '$!!)
195 ; (convert left '$expr)))
197 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
198 ;;;
199 ;;; The implementation of the function Double factorial
200 ;;;
201 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
203 ;;; Double factorial distributes over bags
207 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
209 ;;; Double factorial has mirror symmetry
213 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
215 ;;; Differentiation of Double factorial
231 grad)
233 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
235 ;;; Double factorial is a simplifying function
240 ;; Positive Integer less then $factlim or $factlim is -1. Call gfact.
246 ;; Even negative integer or real representation. Not defined.
257 ;; Odd negative integer, float or complex float.
264 ;; bigfloat or complex bigfloat.
268 ;; double_factorial(inf) -> inf
278 ;; Special case double_factorial(n-1)
279 ;; Not sure if this simplification is useful.
283 ;; Special case double_factorial(2*k+n), k integer
296 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
297 ;;; Double factorial for a complex float argument. The result is a CL complex.
298 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
302 (*
309 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
310 ;;; Double factorial for a bigfloat or complex bigfloat argument
311 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
328 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
329 ;;;
330 ;;; The implementation of the Incomplete Gamma function
331 ;;;
332 ;;; A&S 6.5.3:
333 ;;;
334 ;;; inf
335 ;;; /
336 ;;; [ a - 1 - t
337 ;;; gamma_incomplete(a, z) = I t %e dt
338 ;;; ]
339 ;;; /
340 ;;; z
341 ;;;
342 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
346 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
349 ;;; Incomplete Gamma distributes over bags
353 ;;; Incomplete Gamma function has not mirror symmetry for z on the negative
354 ;;; real axis. We support a conjugate-function which test this case.
361 ;; Definitely not on the negative real axis for z. Mirror symmetry.
368 ;; On the negative real axis or no information. Unsimplified.
373 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
375 ;;; Derivative of the Incomplete Gamma function
381 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2
382 ;; function and the Generalized Incomplete Gamma function
383 ;; (functions.wolfram.com), only for a>0.
399 ;; No derivative. Maxima generates a noun form.
400 nil)))
401 ;; The derivative wrt z
407 ;;; Integral of the Incomplete Gamma function
411 nil
415 integral)
417 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
419 ;;; We support a simplim%function. The function is looked up in simplimit and
420 ;;; handles specific values of the function.
425 ;; Set preserve-direction to true and find the limit of each argument.
430 ;; when either a or z is und, return und
432 ;; z in {minf, inf, infinity}, use http://dlmf.nist.gov/8.11#i
436 ;; Revert a & z to the arguments of gamma_incomplete.
439 ;; Use gamma_incomplete(a,z) = z^(a-1)/exp(z) + ... for
440 ;; fixed a and |z| --> inf. Unfortunately, limit(x*exp(%i*x),x,inf) = und.
441 ;; And that makes limit(gamma_incomplete(2, -%i*x), x, inf) = und, not
442 ;; infinity (this is a test in rtest_limit) But this bug needs to be fixed
443 ;; elsewhere. When a isn't fixed, give up.
456 ;; z in {zerob, 0, zeroa}.
458 ;; Two cases: (i) a <= 0 & integer (ii) a not a negative integer.
459 ;; Revert a & z to the arguments of gamma_incomplete.
463 ;; gamma_incomplete(a,n) = - (-1)^(-a) log(z)/(-a)! + ...
469 ;; z is on the branch cut. We need to know if the imaginary part of
470 ;; (caddr exp) approaches zero from above or below. The incomplete
471 ;; gamma function is continuous from above its branch cut. The check for
472 ;; $ind is needed to avoid calling sign on $ind.
486 ;; Revert a & z to the arguments of gamma_incomplete. When z is
487 ;; off the negative real axis or a is not a negative integer,
488 ;; return ind.
501 ;; All other cases are handled by the simplifier of the function.
504 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
506 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
507 ;;;
508 ;;; Implementation of the Lower Incomplete gamma function:
509 ;;;
510 ;;; A&S 6.5.2
511 ;;;
512 ;;; z
513 ;;; /
514 ;;; [ a - 1 - t
515 ;;; gamma_incomplete_lower(a, z) = I t %e dt
516 ;;; ]
517 ;;; /
518 ;;; 0
519 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
520 ;;; distribute over bags (aggregates)
524 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
526 ;; (defprop %gamma_incomplete_lower ??? grad) WHAT TO PUT HERE ??
528 ;;
529 ;; Handles some special cases for the order a and simplifies it to an
530 ;; equivalent form, possibly involving erf and gamma_incomplete_lower
531 ;; to a lower order.
541 ($hypergeometric_representation
542 ;; We make use of the fact that gamma_incomplete_lower(a,z) +
543 ;; gamma_incomplete(a,z) = gamma(a). Thus,
544 ;; gamma_incomplete_lower(a,z) = gamma(a)-gamma_incomplete(a,z).
545 ;; And we know the hypergeometric form for gamma_incomplete.
551 ;; Try to express gamma_incomplete_lower(a,z) in simpler terms for
552 ;; special values of the order a. If we can't, return NIL to say so.
555 ;; gamma_incomplete_lower(0, x) = 1-exp(x)
558 ;; gamma_incomplete_lower(a,z) can be simplified if a is a
559 ;; positive integer. Since gamma_incomplete_lower(a,z) =
560 ;; gamma(a) - gamma_incomplete(a,z), use gamma_incomplete to
561 ;; get the result.
565 ;; A&S 6.5.12:
566 ;;
567 ;; gamma_incomplete_lower(1/2,x^2) = sqrt(%pi)*erf(x)
568 ;;
569 ;; gamma_incomplete_lower(1/2,z) = sqrt(%pi)*erf(sqrt(x))
570 ;;
574 ;; gamma_incomplete_lower(a+n,z), where n is an integer
579 ;; gamma_incomplete_lower(a+n,z). See DLMF 8.8.7:
580 ;;
581 ;; gamma_incomplete_lower(a+n,z)
582 ;; = pochhammer(a,n)*gamma_incomplete_lower(a,z)
583 ;; -z^a*exp(-z)*sum(gamma(a+n)gamma(a+k+1)*z^k,k,0,n-1)
602 ;; See DLMF 8.8.8. For simplicity let g(a,z) = gamma_incomplete_lower(a,z).
603 ;;
604 ;; g(a,z) = gamma(a)/gamma(a-n)*g(a-n,z)
605 ;; - z^(a-1)*exp(z)*sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
606 ;;
607 ;; Rewrite:
608 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
609 ;; + gamma(a-n)/gamma(a)*z^(a-1)*exp(-z)
610 ;; * sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
611 ;; Or
612 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
613 ;; + z^(a-1)*exp(-z)
614 ;; * sum(gamma(a-n)/gamma(a-k)*z^(-k),k,0,n-1)
635 ;; gamma_incomplete_lower of (numeric) rational order.
636 ;; Expand it out so that the resulting order is between 0 and
637 ;; 1.
640 ;; a = n + order where 0 <= order < 1.
644 ;; Nothing to do if the order is already between 0 and 1
647 ;; Use gamma_incomplete(a+n,z) above. and then substitute
648 ;; a=order. This works for n positive or negative.
653 ;; No expansion so return nil to indicate that
654 nil)))
656 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
658 ;;; Incomplete Gamma function is a simplifying function
665 ;; Check for specific values
668 ;; gamma_incomplete(v,0) is gamma(v) only if the realpart(v) >
669 ;; 0. If realpart(v) <= 0, gamma_incomplete is undefined. For
670 ;; all other cases, return the noun form.
676 a z))
685 ;; Check for numerical evaluation in Float or Bigfloat precision
690 ;; a and z are Maxima numbers, at least one has a float value
697 ;; a is zero or a negative float representing an integer.
698 ;; For these cases the numerical routines of gamma-incomplete
699 ;; do not work. Call the numerical routine for the Exponential
700 ;; Integral E(n,z). The routine is called with a positive integer!.
710 ;; a and z are Maxima numbers, at least one is a complex value and
711 ;; we have at least one float part
719 ;; Call expintegral-e. See comment above.
734 ;; Call bfloat-expintegral-e. See comment above.
755 ;; Call bfloat-expintegral-e. See comment above.
765 ;; Check for transformations and argument simplification
768 ;; Integer or a symbol declared to be an integer. Expand in a series.
790 ;; Simplify gamma, but do not expand to avoid division
791 ;; by zero.
819 ;; We have a half integral order and $gamma_expand is not NIL.
820 ;; We expand in a series with the Erfc function
870 ;; Handle gamma_incomplete(a+n, z), where n is a numerical integer
875 ;; See DLMF 8.8.9: https://dlmf.nist.gov/8.8.E9
876 ;;
877 ;; gamma_incomplete(a+n,z) = pochhammer(a,n)*gamma_incomplete(a,z)
878 ;; + z^a*exp(-z)*sum(gamma(a+n)/gamma(a+k+1)*z^k,k,0,n-1)
896 ;; See http://functions.wolfram.com/06.06.17.0004.01
897 ;;
898 ;; gamma_incomplate(a,z) =
899 ;; (-1)^n*pochhammer(1-a,n)
900 ;; *[gamma_incomplete(a-n,z)
901 ;; + z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)]
902 ;;
903 ;; Rerarrange this in terms of gamma_incomplete(a-n,z):
904 ;;
905 ;; gamma_incomplete(a-n,z) =
906 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
907 ;; -z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)
908 ;;
909 ;; Change the summation index to go from k = 0 to n-1:
910 ;;
911 ;; z^(a-n-1)*sum(z^k/pochhammer(a-n,k),k,1,n)
912 ;; = z^(a-n-1)*sum(z^(k+1)/pochhammer(a-n,k+1),k,0,n-1)
913 ;; = z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
914 ;;
915 ;; Thuus:
916 ;; gamma_incomplete(a-n,z) =
917 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
918 ;; -z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
920 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
929 ;; sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
939 ;; gamma_incomplete of (numeric) rational order. Expand it out
940 ;; so that the resulting order is between 0 and 1.
943 ;; a = n + order where 0 <= order < 1.
947 ;; Nothing to do if the order is already between 0 and 1
950 ;; Use gamma_incomplete(a+n,z) above. and then substitute
951 ;; a=order. This works for n positive or negative.
956 ($hypergeometric_representation
957 ;; See http://functions.wolfram.com/06.06.26.0002.01
958 ;;
959 ;; gamma_incomplete(a,z) = gamma(z) - z^a/a*%f[1,1]([a],[a+1},-z)
960 ;;
961 ;; However, hypergeomtric simplifies
962 ;; hypergeometric([a],[a+1],-z) to
963 ;; hypergeometric([1],[a+1],z)*exp(-z).
973 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
974 ;;; Numerical evaluation of the Incomplete Gamma function
975 ;;;
976 ;;; gamma-incomplete (a,z) - real and complex double float
977 ;;; bfloat-gamma-incomplete (a z) - bigfloat
978 ;;; complex-bfloat-gamma-incomplete (a z) - complex bigfloat
979 ;;;
980 ;;; Expansion in a power series for abs(x) < R, where R is
981 ;;; *gamma-radius* + real(a) if real(a) > 0 or *gamma-radius*
982 ;;; otherwise.
983 ;;;
984 ;;; (A&S 6.5.29):
985 ;;;
986 ;;; inf
987 ;;; ===
988 ;;; \ gamma(a)
989 ;;; gamma(a,z) = exp(-x)*z^a * > ------------ * z^n
990 ;;; / gamma(a+1+n)
991 ;;; ===
992 ;;; n=0
993 ;;;
994 ;;; This expansion does not work for an integer a<=0, because the Gamma function
995 ;;; in the denominator is not defined for a=0 and negative integers. For this
996 ;;; case we use the Exponential Integral E for numerically evaluation. The
997 ;;; Incomplete Gamma function and the Exponential integral are connected by
998 ;;;
999 ;;; gamma(a,z) = z^a * expintegral_e(1-a,z)
1000 ;;;
1001 ;;; When the series is not used, two forms of the continued fraction
1002 ;;; are used. When z is not near the negative real axis use the
1003 ;;; continued fractions (A&S 6.5.31):
1004 ;;;
1005 ;;; 1 1-a 1 2-a 2
1006 ;;; gamma(a,z) = exp(-z) z^a *( -- --- --- --- --- ... )
1007 ;;; z+ 1+ z+ 1+ z+
1008 ;;;
1009 ;;; The accuracy is controlled by *gamma-incomplete-eps* for double float
1010 ;;; precision. For bigfloat precision epsilon is 10^(-$fpprec). The expansions
1011 ;;; in a power series or continued fractions stops if *gamma-incomplete-maxit*
1012 ;;; is exceeded and an Maxima error is thrown.
1013 ;;;
1014 ;;; The above fraction does not converge on the negative real axis and
1015 ;;; converges very slowly near the axis. In this case, use the
1016 ;;; relationship
1017 ;;;
1018 ;;; gamma(a,z) = gamma(a) - gamma_lower(a,z)
1019 ;;;
1020 ;;; The continued fraction for gamma_incomplete_lower(a,z) is
1021 ;;; (http://functions.wolfram.com/06.06.10.0009.01):
1022 ;;;
1023 ;;; gamma_lower(a,z) = exp(-z) * z^a / cf(a,z)
1024 ;;;
1025 ;;; where
1026 ;;;
1027 ;;; -a*z z (a+1)*z 2*z (a+2)*z
1028 ;;; cf(a,z) = a + ---- ---- ------- ---- -------
1029 ;;; a+1+ a+2- a+3+ a+4- a+5+
1030 ;;;
1031 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1041 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1043 ;;; The numerical evaluation for CL float or complex values a and x
1044 ;;; When the flag regularized is T, the result is divided by gamma(a) and
1045 ;;; Maxima returns the numerical result for gamma_incomplete_regularized
1051 ;; Compute the factor needed to scale the series or continued
1052 ;; fraction. This is x^a*exp(-x) or x^a*exp(-x)/gamma(a)
1053 ;; depending on whether we want a non-regularized or
1054 ;; regularized form. We want to compute the factor carefully
1055 ;; to avoid unnecessary overflow if possible.
1059 ;; gammafloat is more accurate for real
1060 ;; values of a.
1067 ;; Easy way failed. Use logs to compute the
1068 ;; result. This loses some precision, especially
1069 ;; for large values of a and/or x because we use
1070 ;; many bits to hold the exponent part.
1072 x
1074 a
1080 ;; Easy way failed, so use the log form.
1082 x)))))))
1086 ;; %gamma-incomplete compute the lower incomplete gamma
1087 ;; function, so we need to subtract that from gamma(a),
1088 ;; more or less.
1096 ;; Continued fraction used. Just multiply by the factor
1097 ;; to get the final result.
1100 ;; Compute the key part of the gamma incomplete function using either
1101 ;; a series expression or a continued fraction expression. Two values
1102 ;; are returned: the value itself and a boolean, indicating what the
1103 ;; computed value is. If the boolean non-NIL, then the computed value
1104 ;; is the lower incomplete gamma function.
1112 ;; The series expansion is done for x within a circle of radius
1113 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1114 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1115 ;; continued fraction is used.
1121 ;; For x near the negative real axis, use the
1122 ;; relationship gamma_incomplete(a,z) = gamma(a) -
1123 ;; gamma_incomplete_lower(a,z), where
1124 ;; gamma_incomplete_lower(a,z) is the lower poart of the
1125 ;; incomplete gamma function. We can evaluate that
1126 ;; using a continued fraction from
1127 ;; http://functions.wolfram.com/06.06.10.0009.01. (Note
1128 ;; that the alternative fraction,
1129 ;; http://functions.wolfram.com/06.06.10.0007.01,
1130 ;; appears to be less accurate.)
1131 ;;
1132 ;; Also note that this appears to be valid for all
1133 ;; values of x (real or complex), but we don't want to
1134 ;; use this everywhere for gamma_incomplete. Consider
1135 ;; what happens for large real x. gamma_incomplete(a,x)
1136 ;; is small, but gamma_incomplete(a,x) = gamma(x) - cf
1137 ;; will have large roundoff errors.
1150 t)))
1152 ;; Expansion in continued fractions for gamma_incomplete.
1163 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1172 ;; Return nil to indicate we used the continued fraction.
1175 ;; Expansion in a series
1183 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1186 ;; Return T to indicate we used the series and the series
1187 ;; is for the integral from 0 to x, not x to inf.
1190 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1192 ;;; This function is called for a and x real
1200 ;; The series expansion is done for x within a circle of radius
1201 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1202 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1203 ;; continued fraction is used.
1210 ;; Expansion in continued fractions of the Incomplete Gamma function
1219 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1240 ;; Expand to multiply everything out.
1241 ($expand
1242 ;; Expansion in continued fraction for the lower incomplete gamma.
1244 ;; NOTE: We want (power x a) instead of bigfloat:expt
1245 ;; because this preserves how maxima computes x^a when
1246 ;; x is negative and a is rational. For, example
1247 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1248 ;; principal value.
1261 x))))))))))))))
1264 ;; Series expansion of the Incomplete Gamma function
1270 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1284 ($rectform
1289 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1301 ;; The series expansion is done for x within a circle of radius
1302 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1303 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1304 ;; continued fraction is used.
1316 ;; Expansion in continued fractions of the Incomplete Gamma function
1327 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1338 gm-eps))
1347 ;; Expand to multiply everything out.
1348 ($expand
1349 ;; Expansion in continued fraction for the lower incomplete gamma.
1351 ;; NOTE: We want (power x a) instead of bigfloat:expt
1352 ;; because this preserves how maxima computes x^a when
1353 ;; x is negative and a is rational. For, example
1354 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1355 ;; principal value.
1368 x))))))))))))))
1370 ;; Series expansion of the Incomplete Gamma function
1378 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1391 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1392 ;;;
1393 ;;; Implementation of the Generalized Incomplete Gamma function
1394 ;;;
1395 ;;; z2
1396 ;;; /
1397 ;;; [ a - 1 - t
1398 ;;; gamma_incomplete_generalized(a, z1, z2) = I t %e dt
1399 ;;; ]
1400 ;;; /
1401 ;;; z1
1402 ;;;
1403 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1405 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1406 ;;; on the negative real axis.
1407 ;;; We support a conjugate-function which test this case.
1410 conjugate-gamma-incomplete-generalized conjugate-function)
1415 ;; z1 and z2 definitely not on the negative real axis.
1416 ;; Mirror symmetry.
1424 ;; On the negative real axis or no information. Unsimplified.
1429 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1431 ;;; Generalized Incomplete Gamma distributes over bags
1435 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1437 ;;; Differentiation of Generalized Incomplete Gamma function
1441 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1442 ;; and the Generalized Incomplete Gamma function (functions.wolfram.com)
1464 ;; The derivative wrt z1
1468 ;; The derivative wrt z2
1472 grad)
1474 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1476 ;;; Generalized Incomplete Gamma function is a simplifying function
1482 ;; Check for specific values
1509 ;; Check for numerical evaluation in Float or Bigfloat precision
1510 ;; Use the numerical routines of the Incomplete Gamma function
1537 ;; Check for transformations and argument simplification
1540 ;; Expand gamma_incomplete_generalized(a+n,z1,z2) with n an integer
1590 ($hypergeometric_representation
1591 ;; Use the fact that gamma_incomplete_generalized(a,z1,z2) =
1592 ;; gamma_incomplete(a,z1) - gamma_incomplete(a,z2).
1597 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1598 ;;;
1599 ;;; Implementation of the Regularized Incomplete Gamma function
1600 ;;;
1601 ;;; A&S 6.5.1
1602 ;;;
1603 ;;; gamma_incomplete(a, z)
1604 ;;; gamma_incomplete_regularized(a, z) = ----------------------
1605 ;;; gamma(a)
1606 ;;;
1607 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1609 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1611 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1612 ;;; on the negative real axis.
1613 ;;; We support a conjugate-function which test this case.
1616 conjugate-gamma-incomplete-regularized conjugate-function)
1621 ;; z definitely not on the negative real axis. Mirror symmetry.
1628 ;; On the negative real axis or no information. Unsimplified.
1633 ;;; Regularized Incomplete Gamma distributes over bags
1637 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1639 ;;; Differentiation of Regularized Incomplete Gamma function
1643 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1644 ;; and the Regularized Generalized Incomplete Gamma function
1645 ;; (functions.wolfram.com)
1659 ;; The derivative wrt z
1661 -1
1665 grad)
1667 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1669 ;;; Regularized Incomplete Gamma function is a simplifying function
1677 ;; Check for specific values
1685 a z))
1692 ;; Check for numerical evaluation in Float or Bigfloat precision
1696 ;; gamma_incomplete returns a regularized result
1702 ;; gamma_incomplete returns a regularized result
1714 ($rectform
1719 ;; Check for transformations and argument simplification
1722 ;; An integer. Expand the expression.
1739 ;; We have a half integral order and $gamma_expand is not NIL.
1740 ;; We expand in a series with the Erfc function
1793 ;; We factor the second summand.
1794 ;; Some factors vanish and the result is more readable.
1795 ($factor
1810 ;; We factor the second summand.
1811 ($factor
1825 ;; gamma_incomplete_regularized of numeric rational order.
1826 ;; Expand it out so that the resulting order is between 0 and
1827 ;; 1. Use gamma_incomplete_regularized(a+n,z) to do the dirty
1828 ;; work.
1831 ;; a = n + order where 0 <= order < 1.
1835 ;; Nothing to do if the order is already between 0 and 1
1838 ;; Use gamma_incomplete_regularized(a+n,z) above. and
1839 ;; then substitute a=order. This works for n positive or
1840 ;; negative.
1845 ($hypergeometric_representation
1846 ;; gamma_incomplete_regularized(a,z)
1847 ;; = gamma_incomplete(a,z)/gamma(a)
1848 ;; = (gamma(a) - gamma_incomplete_lower(a,z))/gamma(a)
1849 ;; = 1 - gamma_incomplete_lower(a,z)/gamma(a)
1855 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1856 ;;;
1857 ;;; Implementation of the Logarithm of the Gamma function
1858 ;;;
1859 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1864 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1872 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1876 ;;; Logarithm of the Gamma function distributes over bags
1880 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1885 grad)
1887 ;; integrate(log_gamma(x),x) = psi[-2](x)
1893 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1900 ;; Check for specific values
1906 ;; We have zero, a negative integer or a float or bigfloat representation.
1912 ;; Check for numerical evaluation
1929 ;; Transform to Logarithm of Factorial for integer values
1930 ;; At this point the integer value is positive and not zero.
1937 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1938 ;;; The functions log-gamma-lanczos, bfloat-log-gamma and
1939 ;;; complex-bfloat-log-gamma are modified versions of the related functions
1940 ;;; gamma-lanczos, bffac and cbffac. The functions return the Logarithm of
1941 ;;; the Gamma function. If we have to calculate the quotient of Gamma functions,
1942 ;;; e. g. for the Beta function, it is much more appropriate to use the
1943 ;;; logarithmic versions to avoid overflow.
1944 ;;;
1945 ;;; Be careful log(gamma(z)) is only for realpart(z) positive equal to
1946 ;;; log_gamma(z). For a negative realpart(z) log_gamma differ by multiple of
1947 ;;; %pi from log(gamma(z)). But we always have exp(log_gamma(z))= gamma(z).
1948 ;;; The terms to get the transformation for log_gamma(-z) are taken from
1949 ;;; functions.wolfram.com.
1950 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1956 :initial-contents
1958 57.156235665862923517
1959 -59.597960355475491248
1960 14.136097974741747174
1961 -0.49191381609762019978
1962 .33994649984811888699e-4
1963 .46523628927048575665e-4
1964 -.98374475304879564677e-4
1965 .15808870322491248884e-3
1966 -.21026444172410488319e-3
1967 .21743961811521264320e-3
1968 -.16431810653676389022e-3
1969 .84418223983852743293e-4
1970 -.26190838401581408670e-4
1976 (+
1977 (*
1985 (*
1999 sum)
2005 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2026 bigfloat%pi
2029 bigfloat%pi '$%i
2047 y)))
2053 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2074 bigfloat%pi
2077 bigfloat%pi '$%i
2096 y))))
2097 ($rectform
2103 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2104 ;;;
2105 ;;; Implementation of the Error function Erf(z)
2106 ;;;
2107 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2109 ;;; erf has mirror symmetry
2113 ;;; erf is an odd function
2117 ;;; erf distributes over bags
2121 ;;; Derivative of the Error function erf
2128 grad)
2130 ;;; Integral of the Error function erf
2139 integral)
2141 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2144 ; See A&S 7.1.21 or http://functions.wolfram.com/06.25.26.0001.01
2146 z
2153 ;;; erf is a simplifying function
2158 ;; Check for specific values
2164 ;; Check for numerical evaluation
2177 ;; Argument simplification
2186 ;; Representation through more general functions
2188 ($hypergeometric_representation
2191 ;; Transformation to Erfc or Erfi
2196 (%erfc
2198 (%erfi
2206 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2209 ;; We use the slatec routine for float values.
2212 ;; Compute erf(z) using the relationship
2213 ;;
2214 ;; erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2,z^2)/sqrt(%pi))
2215 ;;
2216 ;; When z is real sqrt(z^2)/z is signum(z). For complex z,
2217 ;; sqrt(z^2)/z = 1 if -%pi/2 < arg(z) <= %pi/2 and -1 otherwise.
2218 ;;
2219 ;; This relationship has serious round-off issues when z is small
2220 ;; because gamma_incomplete(1/2,z^2)/sqrt(%pi) is near 1.
2221 ;;
2222 ;; complex-erf is for (lisp) complex numbers; bfloat-erf is for real
2223 ;; bfloats, and complex-bfloat-erf is for complex bfloats. Care is
2224 ;; taken to return real results for real arguments and imaginary
2225 ;; results for imaginary arguments
2228 (*
2230 -1
2236 ;; Pure real argument, the result is real
2239 ;; Pure imaginary argument, the result is pure imaginary
2242 result))))
2245 ;; Warning! This has round-off problems when abs(z) is very small.
2247 ;; The argument is real, the result is real too
2248 ($realpart
2257 ;; Warning! This has round-off problems when abs(z) is very small.
2271 ;; Pure real argument, the result is real
2274 ;; Pure imaginary argument, the result is pure imaginary
2277 ;; A general complex result
2278 result))))
2282 ;; Erf(z) for all z. Z must be a CL real or complex number or a
2283 ;; BIGFLOAT or COMPLEX-BIGFLOAT object. The result will be of the
2284 ;; same type as Z.
2287 ;; Use Slatec derf, which should be faster than the series.
2290 ;; Use the series A&S 7.1.5 for small x:
2291 ;;
2292 ;; erf(z) = 2*z/sqrt(%pi) * sum((-1)^n*z^(2*n)/n!/(2*n+1), n, 0, inf)
2293 ;;
2294 ;; The threshold is approximately erf(x) = 0.5. (Doesn't
2295 ;; have to be super accurate.) This gives max accuracy when
2296 ;; using the identity erf(x) = 1 - erfc(x).
2302 k
2306 ;; The general case.
2313 (bigfloat (maxima::$bfloat (maxima::$expand (maxima::complex-bfloat-erf (maxima::to z))))))))))
2316 ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
2317 ;; near 1. Wolfram says
2318 ;;
2319 ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
2320 ;;
2321 ;; For real(z) > 0, sqrt(z^2)/z is 1 so
2322 ;;
2323 ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2324 ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
2325 ;;
2326 ;; For real(z) < 0, sqrt(z^2)/z is -1 so
2327 ;;
2328 ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2329 ;; = 2 - 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
2351 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2352 ;;;
2353 ;;; Implementation of the Generalized Error function Erf(z1,z2)
2354 ;;;
2355 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2357 ;;; Generalized Erf has mirror symmetry
2361 ;;; Generalized Erf distributes over bags
2365 ;;; Generalized Erf is antisymmetric Erf(z1,z2) = - Erf(z2,z1)
2371 ;; Let's remove built-in symbols from list for user-defined properties.
2374 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2378 ;; derivative wrt z1
2382 ;; derivative wrt z2
2386 grad)
2388 ;;; ----------------------------------------------------------------------------
2393 ;; Look for the limit of the arguments.
2397 ;; Handle infinities at this place.
2411 ;; All other cases are handled by the simplifier of the function.
2414 ;;; ----------------------------------------------------------------------------
2419 ;; Check for specific values
2437 ;; Check for numerical evaluation. Use erf(z1,z2) = erf(z2)-erf(z1)
2460 ;; Argument simplification
2465 ;; Representation through more general functions
2467 ($hypergeometric_representation
2470 ;; Transformation to Erf
2472 ($erf_representation
2478 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2479 ;;;
2480 ;;; Implementation of the Complementary Error function Erfc(z)
2481 ;;;
2482 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2486 ;;; Complementary Error function distributes over bags
2490 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2497 grad)
2499 ;;; Integral of the Error function erfc
2508 integral)
2510 ;;; ----------------------------------------------------------------------------
2515 ;; Look for the limit of the arguments.
2518 ;; Handle infinities at this place.
2535 ;; All other cases are handled by the simplifier of the function.
2538 ;;; ----------------------------------------------------------------------------
2543 ;; Check for specific values
2549 ;; Check for numerical evaluation.
2554 ;; Argument simplification
2560 ;; Representation through more general functions
2562 ($hypergeometric_representation
2565 ;; Transformation to Erf or Erfi
2570 (%erf
2572 (%erfi
2580 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2581 ;;;
2582 ;;; Implementation of the Imaginary Error function Erfi(z)
2583 ;;;
2584 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2586 ;;; erfi has mirror symmetry
2590 ;;; erfi is an odd function
2594 ;;; erfi distributes over bags
2598 ;;; Derivative of the Error function erfi
2605 grad)
2607 ;;; Integral of the Error function erfi
2616 integral)
2618 ;;; ----------------------------------------------------------------------------
2623 ;; Look for the limit of the arguments.
2626 ;; Handle extended reals
2634 ;; All other cases are handled by the simplifier of the function.
2637 ;;; ----------------------------------------------------------------------------
2642 ;; Use the relationship erfi(z) = -%i*erf(%i*z)
2646 ;; Take care to return real results when the argument is real.
2655 ;;; erfi is an simplifying function
2660 ;; Check for specific values
2666 ;; Check for numerical evaluation. Use erfi(z) = -%i*erf(%i*z).
2671 ;; Argument simplification
2679 ;; Representation through more general functions
2681 ($hypergeometric_representation
2684 ;; Transformation to Erf or Erfc
2689 (%erf
2691 (%erfc
2699 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2700 ;;;
2701 ;;; The implementation of the Inverse Error function
2702 ;;;
2703 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2705 ;;; The Inverse Error function distributes over bags
2709 ;;; inverse_erf is the inverse of the erf function
2714 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2716 ;;; Differentiation of the Inverse Error function
2724 grad)
2726 ;;; Integral of the Inverse Error function
2733 integral)
2735 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2737 ;;; We support a simplim%function. The function is looked up in simplimit and
2738 ;;; handles specific values of the function.
2743 ;; Look for the limit of the argument.
2746 ;; Handle an argument 1 at this place
2750 ;; All other cases are handled by the simplifier of the function.
2753 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2755 ;;; The Inverse Error function is a simplifying function
2770 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2771 ;;;
2772 ;;; The implementation of the Inverse Complementary Error function
2773 ;;;
2774 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2776 ;;; Inverse Complementary Error function distributes over bags
2780 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2782 ;;; inverse_erfc is the inverse of the erfc function
2788 ;;; Differentiation of the Inverse Complementary Error function
2796 grad)
2798 ;;; Integral of the Inverse Complementary Error function
2806 integral)
2808 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2810 ;;; We support a simplim%function. The function is looked up in simplimit and
2811 ;;; handles specific values of the function.
2816 ;; Look for the limit of the argument.
2819 ;; Handle an argument 0 at this place
2826 ;; All other cases are handled by the simplifier of the function.
2829 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2831 ;;; Inverse Complementary Error function is a simplifying function
2845 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2847 ;;; Implementation of the Newton algorithm to calculate numerical values
2848 ;;; of the Inverse Error functions in float or bigfloat precision.
2849 ;;; The algorithm is a modified version of the routine in newton1.mac.
2886 ;; Compute inverse_erf(z) for z a real or complex number, including
2887 ;; bigfloat objects. The value is computing using a Newton iteration
2888 ;; to solve erf(x) = z.
2891 z)
2895 z))
2897 ;; inverse_erf is odd because erf is.
2902 ;; Find an approximate solution for x = inverse_erf(z).
2905 ;; For small z, inverse_erf(z) = z*sqrt(%pi)/2
2906 ;; + O(z^3). Thus, x = z*sqrt(%pi)/2 is our
2907 ;; initial starting point.
2910 ;; For |z| > 1 and realpart(z) >= 0, we have
2911 ;; the asymptotic series z = erf(x) = 1 -
2912 ;; exp(-x^2)/x/sqrt(%pi).
2913 ;;
2914 ;; Then
2915 ;; x = sqrt(-log(x*sqrt(%pi)*(1-z))
2916 ;;
2917 ;; We can use this as a fixed-point iteration
2918 ;; to find x, and we can start the iteration at
2919 ;; x = 1. Just do one more iteration. I (RLT)
2920 ;; think that's close enough to get the Newton
2921 ;; algorithm to converge.
2928 result)))
2931 ;; Try to pick a reasonable epsilon value for the
2932 ;; Newton iteration.
2943 ;; Derivative of erf(x)
2946 #'diff
2947 z
2949 eps)))))))
2955 z))
2961 ;; Find an approximate solution for x =
2962 ;; inverse_erfc(z). We have inverse_erfc(z) =
2963 ;; inverse_erf(1-z), so that's a good starting point.
2964 ;; We don't need full precision, so a float value is
2965 ;; good enough. But if 1-z is 1, inverse_erf is
2966 ;; undefined, so we need to do something else.
2977 result)))
2983 ;; Derivative of erfc(x)
2986 #'diff
2987 z
2989 eps)))))))
2991 ;; Newton iteration for solving f(x) = z, given f and the derivative
2992 ;; of f.
3004 (intl:gettext "bf-newton: failed to converge after ~M iterations: delta = ~S, x = ~S eps = ~S")
3005 count delta x eps)
3006 x))
3014 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3015 ;;;
3016 ;;; Implementation of the Fresnel Integral S(z)
3017 ;;;
3018 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3020 ;;; fresnel_s distributes over bags
3024 ;;; fresnel_s has mirror symmetry
3028 ;;; fresnel_s is an odd function
3029 ;;;
3030 ;;; Maxima has two mechanism to define a reflection property
3031 ;;; 1. Declare the feature oddfun or evenfun
3032 ;;; 2. Put a reflection rule on the property list
3033 ;;;
3034 ;;; The second way is used for the trig functions. We use it here too.
3036 ;;; This would be the first way to give the property of an odd function.
3037 ;(eval-when
3038 ; (:load-toplevel :execute)
3039 ; (let (($context '$global) (context '$global))
3040 ; (meval '(($declare) %fresnel_s $oddfun))
3041 ; ;; Let's remove built-in symbols from list for user-defined properties.
3042 ; (setq $props (remove '%fresnel_s $props))))
3046 ;;; Differentiation of the Fresnel Integral S
3051 grad)
3053 ;;; Integration of the Fresnel Integral S
3062 integral)
3064 ;;; Limits of the Fresnel Integral S
3076 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3082 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3083 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3092 ;; For very small x, we have
3093 ;; fresnel_s(x) = %pi/6*z^3
3094 ;; fresnel_c(x) = x
3112 (maxima::merror (intl:gettext "fresnel: series expansion failed for (COMPLEX-BFLOAT-FRESNEL ~:M).") z))
3131 (* sqrt-pi
3132 z)))
3134 (* sqrt-pi
3135 z)))
3146 ;; A pure imaginary argument. Use fresnel_s(%i*x)=-%i*fresnel_s(x).
3152 fs))))
3156 ;; A pure imaginary argument. Use fresnel_c(%i*x)=%i*fresnel_c(x).
3162 fc))))
3166 ;;; fresnel_s is a simplifying function
3170 ;; Check for specific values
3176 ;; Check for numerical evaluation
3180 ;; Check for argument simplification
3187 ;; Representation through equivalent functions
3189 ($erf_representation
3204 ($hypergeometric_representation
3213 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3214 ;;;
3215 ;;; Implementation of the Fresnel Integral C(z)
3216 ;;;
3217 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3219 ;;; fresnel_c distributes over bags
3223 ;;; fresnel_c has mirror symmetry
3227 ;;; fresnel_c is an odd function
3228 ;;; Maxima has two mechanism to define a reflection property
3229 ;;; 1. Declare the feature oddfun or evenfun
3230 ;;; 2. Put a reflection rule on the property list
3231 ;;;
3232 ;;; The second way is used for the trig functions. We use it here too.
3236 ;;; Differentiation of the Fresnel Integral C
3241 grad)
3243 ;;; Integration of the Fresnel Integral C
3252 integral)
3254 ;;; Limits of the Fresnel Integral C
3266 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3268 ;;; fresnel_c is a simplifying function
3272 ;; Check for specific values
3278 ;; Check for numerical evaluation
3283 ;; Check for argument simplification
3290 ;; Representation through equivalent functions
3292 ($erf_representation
3307 ($hypergeometric_representation
3316 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3317 ;;;
3318 ;;; Implementation of the Beta function
3319 ;;;
3320 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3322 ;;; The code for the implementation of the beta function is in the files
3323 ;;; csimp2.lisp, simp.lisp and mactex.lisp.
3324 ;;; At this place we only implement the operator property SYMMETRIC.
3326 ;;; Beta is symmetric beta(a,b) = beta(b,a)
3332 ;; Let's remove built-in symbols from list for user-defined properties.
3335 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3336 ;;;
3337 ;;; Implementation of the Incomplete Beta function
3338 ;;;
3339 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3344 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3346 ;;; beta_incomplete distributes over bags
3350 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3354 ;; Derivative wrt a
3362 z)
3364 ;; Derivative wrt b
3381 ;; The derivative wrt z
3385 grad)
3387 ;;; Integral of the Incomplete Beta function
3391 nil
3392 nil
3396 integral)
3398 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3409 ;; Check for specific values
3417 a b z))
3419 z)
3424 ;; z=1 and realpart(b)>0. Simplify to a Beta function.
3425 ;; If we have to evaluate numerically preserve the type.
3440 ;; The argument a is zero or a is negative and the argument b is
3441 ;; not in a valid range. Beta incomplete is undefined.
3442 ;; It would be more correct to return Complex infinity.
3447 ;; Check for numerical evaluation in Float or Bigfloat precision
3455 ;; Negative integer a and b is in a valid range. Expand.
3456 ($rectform
3467 ;; Negative integer a and b is in a valid range. Expand.
3468 ($rectform
3472 ;; Argument simplifications and transformations
3478 ;; Expand for b a positive integer and a not a negative integer.
3492 ;; Expand for a a positive integer.
3508 ;; Expand for a a negative integer and b an integer with b <= -a.
3537 index))
3540 index)))
3566 index)))
3571 ($hypergeometric_representation
3572 ;; There are several cases to consider.
3573 ;;
3574 ;; For -a not an integer see
3575 ;; http://functions.wolfram.com/06.19.26.0005.01
3576 ;;
3577 ;; beta_incomplete(a,b,z) = z^a/a*%f[2,1]([a,1-b],[a+1],z)
3578 ;;
3579 ;; For -b not an integer see
3580 ;; http://functions.wolfram.com/06.19.26.0007.01
3581 ;;
3582 ;; beta_incomplete(a,b,z) = beta(a,b) -
3583 ;; (1-z)^b*z^a/b*%f[2,1]([1,a+b],[b+1],1-z)
3584 ;;
3585 ;; For a+b not a positive integer see
3586 ;; http://functions.wolfram.com/06.19.26.0008.01
3587 ;;
3588 ;; beta_incomplete(a,b,z) =
3589 ;; z^a*(-z)^(b-1)/(a+b-1)*%f[2,1]([1-b,-a-b+1],[-a-b+2],1/z)
3590 ;; + z^a*beta(1-a-b,a)*(-z)^(-a)
3591 ;;
3592 ;; We need to handle more cases here
3594 a)
3598 z)))
3614 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3621 ($rectform
3650 (merror (intl:gettext "beta_incomplete: continued fractions failed for beta_incomplete(~:M, ~:M, ~:M).") a b z))
3688 result))))))
3690 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3691 ;;;
3692 ;;; Implementation of the Generalized Incomplete Beta function
3693 ;;;
3694 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3696 ;;; beta_incomplete_generalized distributes over bags
3700 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3702 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
3703 ;;; but not on the negative real axis and for z1 or z2 real and > 1.
3704 ;;; We support a conjugate-function which test these cases.
3707 conjugate-beta-incomplete-generalized conjugate-function)
3731 a b z1 z2)))))))
3733 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3737 ;; Derivative wrt a
3749 z1)
3755 z2)
3758 ;; Derivative wrt b
3781 ;; The derivative wrt z1
3787 ;; The derivative wrt z2
3793 grad)
3795 ;;; Integral of the Incomplete Beta function
3799 nil
3800 nil
3801 ;; Integral wrt z1
3805 ;; Integral wrt z2
3810 integral)
3812 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3818 ;; Check for specific values
3826 a b z1 z2))
3838 a b z1 z2))
3860 ;; Check for numerical evaluation in Float or Bigfloat precision
3873 ;; Check for argument simplifications and transformations
3942 index))
3945 index)))
3974 index)))
3985 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3986 ;;;
3987 ;;; Implementation of the Regularized Incomplete Beta function
3988 ;;;
3989 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3991 ;;; beta_incomplete_regularized distributes over bags
3995 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3999 ;; Derivative wrt a
4015 ;; Derivative wrt b
4032 ;; The derivative wrt z
4037 grad)
4039 ;;; Integral of the Generalized Incomplete Beta function
4043 nil
4044 nil
4045 ;; Integral wrt z
4051 integral)
4053 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4059 ;; Check for specific values
4067 a b z))
4085 ;; Problem: for b a negative integer the Regularized Incomplete
4086 ;; Beta function is defined to be zero. BUT: When we calculate
4087 ;; e.g. beta_incomplete(1.0,-2.0,1/2)/beta(1.0,-2.0) we get the
4088 ;; result -3.0, because beta_incomplete and beta are defined for
4089 ;; for this case. How do we get a consistent behaviour?
4095 ;; TODO: The following line does not work for bigfloats.
4097 ;; Does $beta_incomplete or simpbeta underflow in this case?
4103 ;; Check for numerical evaluation in Float or Bigfloat precision
4107 beta ibeta )
4111 ;; In case of underflow (see bug #2999) or precision loss use bigfloats
4112 ;; and emporarily give some extra precision but avoid fpprec dependency.
4113 ;; Is this workaround correct for complex values?
4122 ($rectform
4126 ;; Check for argument simplifications and transformations
4150 index))
4153 index)))
4173 index)))
4181 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;