| blob | 01b7449b6d2370fd1cbb5141fd764429af5f6313 |
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 definied 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 ;; z in {minf, inf, infinity}, use http://dlmf.nist.gov/8.11#i
434 ;; Revert a & z to the arguments of gamma_incomplete.
437 ;; Use gamma_incomplete(a,z) = z^(a-1)/exp(z) + ... for
438 ;; fixed a and |z| --> inf. Unfortunately, limit(x*exp(%i*x),x,inf) = und.
439 ;; And that makes limit(gamma_incomplete(2, -%i*x), x, inf) = und, not
440 ;; infinity (this is a test in rtest_limit) But this bug needs to be fixed
441 ;; elsewhere. When a isn't fixed, give up.
454 ;; z in {zerob, 0, zeroa}.
456 ;; Two cases: (i) a <= 0 & integer (ii) a not a negative integer.
457 ;; Revert a & z to the arguments of gamma_incomplete.
461 ;; gamma_incomplete(a,n) = - (-1)^(-a) log(z)/(-a)! + ...
467 ;; z is on the branch cut. We need to know if the imaginary part of
468 ;; (caddr exp) approaches zero from above or below. The incomplete
469 ;; gamma function is continuous from above its branch cut.
484 ;; Revert a & z to the arguments of gamma_incomplete. When z is
485 ;; off the negative real axis or a is not a negative integer,
486 ;; return ind.
499 ;; All other cases are handled by the simplifier of the function.
502 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
504 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
505 ;;;
506 ;;; Implementation of the Lower Incomplete gamma function:
507 ;;;
508 ;;; A&S 6.5.2
509 ;;;
510 ;;; z
511 ;;; /
512 ;;; [ a - 1 - t
513 ;;; gamma_incomplete_lower(a, z) = I t %e dt
514 ;;; ]
515 ;;; /
516 ;;; 0
517 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
518 ;;; distribute over bags (aggregates)
522 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
524 ;; (defprop %gamma_incomplete_lower ??? grad) WHAT TO PUT HERE ??
526 ;;
527 ;; Handles some special cases for the order a and simplifies it to an
528 ;; equivalent form, possibly involving erf and gamma_incomplete_lower
529 ;; to a lower order.
539 ($hypergeometric_representation
540 ;; We make use of the fact that gamma_incomplete_lower(a,z) +
541 ;; gamma_incomplete(a,z) = gamma(a). Thus,
542 ;; gamma_incomplete_lower(a,z) = gamma(a)-gamma_incomplete(a,z).
543 ;; And we know the hypergeometric form for gamma_incomplete.
549 ;; Try to express gamma_incomplete_lower(a,z) in simpler terms for
550 ;; special values of the order a. If we can't, return NIL to say so.
553 ;; gamma_incomplete_lower(0, x) = 1-exp(x)
556 ;; gamma_incomplete_lower(a,z) can be simplified if a is a
557 ;; positive integer. Since gamma_incomplete_lower(a,z) =
558 ;; gamma(a) - gamma_incomplete(a,z), use gamma_incomplete to
559 ;; get the result.
563 ;; A&S 6.5.12:
564 ;;
565 ;; gamma_incomplete_lower(1/2,x^2) = sqrt(%pi)*erf(x)
566 ;;
567 ;; gamma_incomplete_lower(1/2,z) = sqrt(%pi)*erf(sqrt(x))
568 ;;
572 ;; gamma_incomplete_lower(a+n,z), where n is an integer
577 ;; gamma_incomplete_lower(a+n,z). See DLMF 8.8.7:
578 ;;
579 ;; gamma_incomplete_lower(a+n,z)
580 ;; = pochhammer(a,n)*gamma_incomplete_lower(a,z)
581 ;; -z^a*exp(-z)*sum(gamma(a+n)gamma(a+k+1)*z^k,k,0,n-1)
600 ;; See DLMF 8.8.8. For simplicity let g(a,z) = gamma_incomplete_lower(a,z).
601 ;;
602 ;; g(a,z) = gamma(a)/gamma(a-n)*g(a-n,z)
603 ;; - z^(a-1)*exp(z)*sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
604 ;;
605 ;; Rewrite:
606 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
607 ;; + gamma(a-n)/gamma(a)*z^(a-1)*exp(-z)
608 ;; * sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
609 ;; Or
610 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
611 ;; + z^(a-1)*exp(-z)
612 ;; * sum(gamma(a-n)/gamma(a-k)*z^(-k),k,0,n-1)
633 ;; gamma_incomplete_lower of (numeric) rational order.
634 ;; Expand it out so that the resulting order is between 0 and
635 ;; 1.
638 ;; a = n + order where 0 <= order < 1.
642 ;; Nothing to do if the order is already between 0 and 1
645 ;; Use gamma_incomplete(a+n,z) above. and then substitute
646 ;; a=order. This works for n positive or negative.
651 ;; No expansion so return nil to indicate that
652 nil)))
654 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
656 ;;; Incomplete Gamma function is a simplifying function
663 ;; Check for specific values
666 ;; gamma_incomplete(v,0) is gamma(v) only if the realpart(v) >
667 ;; 0. If realpart(v) <= 0, gamma_incomplete is undefined. For
668 ;; all other cases, return the noun form.
674 a z))
683 ;; Check for numerical evaluation in Float or Bigfloat precision
688 ;; a and z are Maxima numbers, at least one has a float value
695 ;; a is zero or a negative float representing an integer.
696 ;; For these cases the numerical routines of gamma-incomplete
697 ;; do not work. Call the numerical routine for the Exponential
698 ;; Integral E(n,z). The routine is called with a positive integer!.
708 ;; a and z are Maxima numbers, at least one is a complex value and
709 ;; we have at least one float part
717 ;; Call expintegral-e. See comment above.
732 ;; Call bfloat-expintegral-e. See comment above.
753 ;; Call bfloat-expintegral-e. See comment above.
763 ;; Check for transformations and argument simplification
766 ;; Integer or a symbol declared to be an integer. Expand in a series.
788 ;; Simplify gamma, but do not expand to avoid division
789 ;; by zero.
817 ;; We have a half integral order and $gamma_expand is not NIL.
818 ;; We expand in a series with the Erfc function
868 ;; Handle gamma_incomplete(a+n, z), where n is a numerical integer
873 ;; See DLMF 8.8.9: https://dlmf.nist.gov/8.8.E9
874 ;;
875 ;; gamma_incomplete(a+n,z) = pochhammer(a,n)*gamma_incomplete(a,z)
876 ;; + z^a*exp(-z)*sum(gamma(a+n)/gamma(a+k+1)*z^k,k,0,n-1)
894 ;; See http://functions.wolfram.com/06.06.17.0004.01
895 ;;
896 ;; gamma_incomplate(a,z) =
897 ;; (-1)^n*pochhammer(1-a,n)
898 ;; *[gamma_incomplete(a-n,z)
899 ;; + z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)]
900 ;;
901 ;; Rerarrange this in terms of gamma_incomplete(a-n,z):
902 ;;
903 ;; gamma_incomplete(a-n,z) =
904 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
905 ;; -z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)
906 ;;
907 ;; Change the summation index to go from k = 0 to n-1:
908 ;;
909 ;; z^(a-n-1)*sum(z^k/pochhammer(a-n,k),k,1,n)
910 ;; = z^(a-n-1)*sum(z^(k+1)/pochhammer(a-n,k+1),k,0,n-1)
911 ;; = z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
912 ;;
913 ;; Thuus:
914 ;; gamma_incomplete(a-n,z) =
915 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
916 ;; -z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
918 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
927 ;; sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
937 ;; gamma_incomplete of (numeric) rational order. Expand it out
938 ;; so that the resulting order is between 0 and 1.
941 ;; a = n + order where 0 <= order < 1.
945 ;; Nothing to do if the order is already between 0 and 1
948 ;; Use gamma_incomplete(a+n,z) above. and then substitute
949 ;; a=order. This works for n positive or negative.
954 ($hypergeometric_representation
955 ;; See http://functions.wolfram.com/06.06.26.0002.01
956 ;;
957 ;; gamma_incomplete(a,z) = gamma(z) - z^a/a*%f[1,1]([a],[a+1},-z)
958 ;;
959 ;; However, hypergeomtric simplifies
960 ;; hypergeometric([a],[a+1],-z) to
961 ;; hypergeometric([1],[a+1],z)*exp(-z).
971 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
972 ;;; Numerical evaluation of the Incomplete Gamma function
973 ;;;
974 ;;; gamma-incomplete (a,z) - real and complex double float
975 ;;; bfloat-gamma-incomplete (a z) - bigfloat
976 ;;; complex-bfloat-gamma-incomplete (a z) - complex bigfloat
977 ;;;
978 ;;; Expansion in a power series for abs(x) < R, where R is
979 ;;; *gamma-radius* + real(a) if real(a) > 0 or *gamma-radius*
980 ;;; otherwise.
981 ;;;
982 ;;; (A&S 6.5.29):
983 ;;;
984 ;;; inf
985 ;;; ===
986 ;;; \ gamma(a)
987 ;;; gamma(a,z) = exp(-x)*z^a * > ------------ * z^n
988 ;;; / gamma(a+1+n)
989 ;;; ===
990 ;;; n=0
991 ;;;
992 ;;; This expansion does not work for an integer a<=0, because the Gamma function
993 ;;; in the denominator is not defined for a=0 and negative integers. For this
994 ;;; case we use the Exponential Integral E for numerically evaluation. The
995 ;;; Incomplete Gamma function and the Exponential integral are connected by
996 ;;;
997 ;;; gamma(a,z) = z^a * expintegral_e(1-a,z)
998 ;;;
999 ;;; When the series is not used, two forms of the continued fraction
1000 ;;; are used. When z is not near the negative real axis use the
1001 ;;; continued fractions (A&S 6.5.31):
1002 ;;;
1003 ;;; 1 1-a 1 2-a 2
1004 ;;; gamma(a,z) = exp(-z) z^a *( -- --- --- --- --- ... )
1005 ;;; z+ 1+ z+ 1+ z+
1006 ;;;
1007 ;;; The accuracy is controlled by *gamma-incomplete-eps* for double float
1008 ;;; precision. For bigfloat precision epsilon is 10^(-$fpprec). The expansions
1009 ;;; in a power series or continued fractions stops if *gamma-incomplete-maxit*
1010 ;;; is exceeded and an Maxima error is thrown.
1011 ;;;
1012 ;;; The above fraction does not converge on the negative real axis and
1013 ;;; converges very slowly near the axis. In this case, use the
1014 ;;; relationship
1015 ;;;
1016 ;;; gamma(a,z) = gamma(a) - gamma_lower(a,z)
1017 ;;;
1018 ;;; The continued fraction for gamma_incomplete_lower(a,z) is
1019 ;;; (http://functions.wolfram.com/06.06.10.0009.01):
1020 ;;;
1021 ;;; gamma_lower(a,z) = exp(-z) * z^a / cf(a,z)
1022 ;;;
1023 ;;; where
1024 ;;;
1025 ;;; -a*z z (a+1)*z 2*z (a+2)*z
1026 ;;; cf(a,z) = a + ---- ---- ------- ---- -------
1027 ;;; a+1+ a+2- a+3+ a+4- a+5+
1028 ;;;
1029 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1039 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1041 ;;; The numerical evaluation for CL float or complex values a and x
1042 ;;; When the flag regularized is T, the result is divided by gamma(a) and
1043 ;;; Maxima returns the numercial result for gamma_incomplete_regularized
1049 ;; Compute the factor needed to scale the series or continued
1050 ;; fraction. This is x^a*exp(-x) or x^a*exp(-x)/gamma(a)
1051 ;; depending on whether we want a non-regularized or
1052 ;; regularized form. We want to compute the factor carefully
1053 ;; to avoid unnecessary overflow if possible.
1057 ;; gammafloat is more accurate for real
1058 ;; values of a.
1065 ;; Easy way failed. Use logs to compute the
1066 ;; result. This loses some precision, especially
1067 ;; for large values of a and/or x because we use
1068 ;; many bits to hold the exponent part.
1070 x
1072 a
1078 ;; Easy way failed, so use the log form.
1080 x)))))))
1084 ;; %gamma-incomplete compute the lower incomplete gamma
1085 ;; function, so we need to subtract that from gamma(a),
1086 ;; more or less.
1094 ;; Continued fraction used. Just multiply by the factor
1095 ;; to get the final result.
1098 ;; Compute the key part of the gamma incomplete function using either
1099 ;; a series expression or a continued fraction expression. Two values
1100 ;; are returned: the value itself and a boolean, indicating what the
1101 ;; computed value is. If the boolean non-NIL, then the computed value
1102 ;; is the lower incomplete gamma function.
1110 ;; The series expansion is done for x within a circle of radius
1111 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1112 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1113 ;; continued fraction is used.
1119 ;; For x near the negative real axis, use the
1120 ;; relationship gamma_incomplete(a,z) = gamma(a) -
1121 ;; gamma_incomplete_lower(a,z), where
1122 ;; gamma_incomplete_lower(a,z) is the lower poart of the
1123 ;; incomplete gamma function. We can evaluate that
1124 ;; using a continued fraction from
1125 ;; http://functions.wolfram.com/06.06.10.0009.01. (Note
1126 ;; that the alternative fraction,
1127 ;; http://functions.wolfram.com/06.06.10.0007.01,
1128 ;; appears to be less accurate.)
1129 ;;
1130 ;; Also note that this appears to be valid for all
1131 ;; values of x (real or complex), but we don't want to
1132 ;; use this everywhere for gamma_incomplete. Consider
1133 ;; what happens for large real x. gamma_incomplete(a,x)
1134 ;; is small, but gamma_incomplete(a,x) = gamma(x) - cf
1135 ;; will have large roundoff errors.
1148 t)))
1150 ;; Expansion in continued fractions for gamma_incomplete.
1161 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1170 ;; Return nil to indicate we used the continued fraction.
1173 ;; Expansion in a series
1181 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1184 ;; Return T to indicate we used the series and the series
1185 ;; is for the integral from 0 to x, not x to inf.
1188 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1190 ;;; This function is called for a and x real
1198 ;; The series expansion is done for x within a circle of radius
1199 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1200 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1201 ;; continued fraction is used.
1208 ;; Expansion in continued fractions of the Incomplete Gamma function
1217 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1238 ;; Expand to multiply everything out.
1239 ($expand
1240 ;; Expansion in continued fraction for the lower incomplete gamma.
1242 ;; NOTE: We want (power x a) instead of bigfloat:expt
1243 ;; because this preserves how maxima computes x^a when
1244 ;; x is negative and a is rational. For, example
1245 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1246 ;; principal value.
1259 x))))))))))))))
1262 ;; Series expansion of the Incomplete Gamma function
1268 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1282 ($rectform
1287 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1299 ;; The series expansion is done for x within a circle of radius
1300 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1301 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1302 ;; continued fraction is used.
1314 ;; Expansion in continued fractions of the Incomplete Gamma function
1325 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1336 gm-eps))
1345 ;; Expand to multiply everything out.
1346 ($expand
1347 ;; Expansion in continued fraction for the lower incomplete gamma.
1349 ;; NOTE: We want (power x a) instead of bigfloat:expt
1350 ;; because this preserves how maxima computes x^a when
1351 ;; x is negative and a is rational. For, example
1352 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1353 ;; principal value.
1366 x))))))))))))))
1368 ;; Series expansion of the Incomplete Gamma function
1376 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1389 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1390 ;;;
1391 ;;; Implementation of the Generalized Incomplete Gamma function
1392 ;;;
1393 ;;; z2
1394 ;;; /
1395 ;;; [ a - 1 - t
1396 ;;; gamma_incomplete_generalized(a, z1, z2) = I t %e dt
1397 ;;; ]
1398 ;;; /
1399 ;;; z1
1400 ;;;
1401 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1403 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1404 ;;; on the negative real axis.
1405 ;;; We support a conjugate-function which test this case.
1408 conjugate-gamma-incomplete-generalized conjugate-function)
1413 ;; z1 and z2 definitely not on the negative real axis.
1414 ;; Mirror symmetry.
1422 ;; On the negative real axis or no information. Unsimplified.
1427 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1429 ;;; Generalized Incomplete Gamma distributes over bags
1433 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1435 ;;; Differentiation of Generalized Incomplete Gamma function
1439 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1440 ;; and the Generalized Incomplete Gamma function (functions.wolfram.com)
1462 ;; The derivative wrt z1
1466 ;; The derivative wrt z2
1470 grad)
1472 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1474 ;;; Generalized Incomplete Gamma function is a simplifying function
1480 ;; Check for specific values
1507 ;; Check for numerical evaluation in Float or Bigfloat precision
1508 ;; Use the numerical routines of the Incomplete Gamma function
1535 ;; Check for transformations and argument simplification
1538 ;; Expand gamma_incomplete_generalized(a+n,z1,z2) with n an integer
1588 ($hypergeometric_representation
1589 ;; Use the fact that gamma_incomplete_generalized(a,z1,z2) =
1590 ;; gamma_incomplete(a,z1) - gamma_incomplete(a,z2).
1595 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1596 ;;;
1597 ;;; Implementation of the Regularized Incomplete Gamma function
1598 ;;;
1599 ;;; A&S 6.5.1
1600 ;;;
1601 ;;; gamma_incomplete(a, z)
1602 ;;; gamma_incomplete_regularized(a, z) = ----------------------
1603 ;;; gamma(a)
1604 ;;;
1605 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1607 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1609 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1610 ;;; on the negative real axis.
1611 ;;; We support a conjugate-function which test this case.
1614 conjugate-gamma-incomplete-regularized conjugate-function)
1619 ;; z definitely not on the negative real axis. Mirror symmetry.
1626 ;; On the negative real axis or no information. Unsimplified.
1631 ;;; Regularized Incomplete Gamma distributes over bags
1635 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1637 ;;; Differentiation of Regularized Incomplete Gamma function
1641 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1642 ;; and the Regularized Generalized Incomplete Gamma function
1643 ;; (functions.wolfram.com)
1657 ;; The derivative wrt z
1659 -1
1663 grad)
1665 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1667 ;;; Regularized Incomplete Gamma function is a simplifying function
1675 ;; Check for specific values
1683 a z))
1690 ;; Check for numerical evaluation in Float or Bigfloat precision
1694 ;; gamma_incomplete returns a regularized result
1700 ;; gamma_incomplete returns a regularized result
1712 ($rectform
1717 ;; Check for transformations and argument simplification
1720 ;; An integer. Expand the expression.
1737 ;; We have a half integral order and $gamma_expand is not NIL.
1738 ;; We expand in a series with the Erfc function
1791 ;; We factor the second summand.
1792 ;; Some factors vanish and the result is more readable.
1793 ($factor
1808 ;; We factor the second summand.
1809 ($factor
1823 ;; gamma_incomplete_regularized of numeric rational order.
1824 ;; Expand it out so that the resulting order is between 0 and
1825 ;; 1. Use gamma_incomplete_regularized(a+n,z) to do the dirty
1826 ;; work.
1829 ;; a = n + order where 0 <= order < 1.
1833 ;; Nothing to do if the order is already between 0 and 1
1836 ;; Use gamma_incomplete_regularized(a+n,z) above. and
1837 ;; then substitute a=order. This works for n positive or
1838 ;; negative.
1843 ($hypergeometric_representation
1844 ;; gamma_incomplete_regularized(a,z)
1845 ;; = gamma_incomplete(a,z)/gamma(a)
1846 ;; = (gamma(a) - gamma_incomplete_lower(a,z))/gamma(a)
1847 ;; = 1 - gamma_incomplete_lower(a,z)/gamma(a)
1853 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1854 ;;;
1855 ;;; Implementation of the Logarithm of the Gamma function
1856 ;;;
1857 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1862 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1870 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1874 ;;; Logarithm of the Gamma function distributes over bags
1878 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1883 grad)
1885 ;; integrate(log_gamma(x),x) = psi[-2](x)
1891 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1898 ;; Check for specific values
1904 ;; We have zero, a negative integer or a float or bigfloat representation.
1910 ;; Check for numerical evaluation
1927 ;; Transform to Logarithm of Factorial for integer values
1928 ;; At this point the integer value is positive and not zero.
1935 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1936 ;;; The functions log-gamma-lanczos, bfloat-log-gamma and
1937 ;;; complex-bfloat-log-gamma are modified versions of the related functions
1938 ;;; gamma-lanczos, bffac and cbffac. The functions return the Logarithm of
1939 ;;; the Gamma function. If we have to calculate the quotient of Gamma functions,
1940 ;;; e. g. for the Beta function, it is much more appropriate to use the
1941 ;;; logarithmic versions to avoid overflow.
1942 ;;;
1943 ;;; Be careful log(gamma(z)) is only for realpart(z) positive equal to
1944 ;;; log_gamma(z). For a negative realpart(z) log_gamma differ by multiple of
1945 ;;; %pi from log(gamma(z)). But we always have exp(log_gamma(z))= gamma(z).
1946 ;;; The terms to get the transformation for log_gamma(-z) are taken from
1947 ;;; functions.wolfram.com.
1948 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1954 :initial-contents
1956 57.156235665862923517
1957 -59.597960355475491248
1958 14.136097974741747174
1959 -0.49191381609762019978
1960 .33994649984811888699e-4
1961 .46523628927048575665e-4
1962 -.98374475304879564677e-4
1963 .15808870322491248884e-3
1964 -.21026444172410488319e-3
1965 .21743961811521264320e-3
1966 -.16431810653676389022e-3
1967 .84418223983852743293e-4
1968 -.26190838401581408670e-4
1974 (+
1975 (*
1983 (*
1997 sum)
2003 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2024 bigfloat%pi
2027 bigfloat%pi '$%i
2045 y)))
2051 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2072 bigfloat%pi
2075 bigfloat%pi '$%i
2094 y))))
2095 ($rectform
2101 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2102 ;;;
2103 ;;; Implementation of the Error function Erf(z)
2104 ;;;
2105 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2107 ;;; erf has mirror symmetry
2111 ;;; erf is an odd function
2115 ;;; erf distributes over bags
2119 ;;; Derivative of the Error function erf
2126 grad)
2128 ;;; Integral of the Error function erf
2137 integral)
2139 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2142 ; See A&S 7.1.21 or http://functions.wolfram.com/06.25.26.0001.01
2144 z
2151 ;;; erf is a simplifying function
2156 ;; Check for specific values
2162 ;; Check for numerical evaluation
2175 ;; Argument simplification
2184 ;; Representation through more general functions
2186 ($hypergeometric_representation
2189 ;; Transformation to Erfc or Erfi
2194 (%erfc
2196 (%erfi
2204 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2207 ;; We use the slatec routine for float values.
2210 ;; Compute erf(z) using the relationship
2211 ;;
2212 ;; erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2,z^2)/sqrt(%pi))
2213 ;;
2214 ;; When z is real sqrt(z^2)/z is signum(z). For complex z,
2215 ;; sqrt(z^2)/z = 1 if -%pi/2 < arg(z) <= %pi/2 and -1 otherwise.
2216 ;;
2217 ;; This relationship has serious round-off issues when z is small
2218 ;; because gamma_incomplete(1/2,z^2)/sqrt(%pi) is near 1.
2219 ;;
2220 ;; complex-erf is for (lisp) complex numbers; bfloat-erf is for real
2221 ;; bfloats, and complex-bfloat-erf is for complex bfloats. Care is
2222 ;; taken to return real results for real arguments and imaginary
2223 ;; results for imaginary arguments
2226 (*
2228 -1
2234 ;; Pure real argument, the result is real
2237 ;; Pure imaginary argument, the result is pure imaginary
2240 result))))
2243 ;; Warning! This has round-off problems when abs(z) is very small.
2245 ;; The argument is real, the result is real too
2246 ($realpart
2255 ;; Warning! This has round-off problems when abs(z) is very small.
2269 ;; Pure real argument, the result is real
2272 ;; Pure imaginary argument, the result is pure imaginary
2275 ;; A general complex result
2276 result))))
2280 ;; Erf(z) for all z. Z must be a CL real or complex number or a
2281 ;; BIGFLOAT or COMPLEX-BIGFLOAT object. The result will be of the
2282 ;; same type as Z.
2285 ;; Use Slatec derf, which should be faster than the series.
2288 ;; Use the series A&S 7.1.5 for small x:
2289 ;;
2290 ;; erf(z) = 2*z/sqrt(%pi) * sum((-1)^n*z^(2*n)/n!/(2*n+1), n, 0, inf)
2291 ;;
2292 ;; The threshold is approximately erf(x) = 0.5. (Doesn't
2293 ;; have to be super accurate.) This gives max accuracy when
2294 ;; using the identity erf(x) = 1 - erfc(x).
2300 k
2304 ;; The general case.
2311 (bigfloat (maxima::$bfloat (maxima::$expand (maxima::complex-bfloat-erf (maxima::to z))))))))))
2314 ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
2315 ;; near 1. Wolfram says
2316 ;;
2317 ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
2318 ;;
2319 ;; For real(z) > 0, sqrt(z^2)/z is 1 so
2320 ;;
2321 ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2322 ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
2323 ;;
2324 ;; For real(z) < 0, sqrt(z^2)/z is -1 so
2325 ;;
2326 ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2327 ;; = 2 - 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
2349 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2350 ;;;
2351 ;;; Implementation of the Generalized Error function Erf(z1,z2)
2352 ;;;
2353 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2355 ;;; Generalized Erf has mirror symmetry
2359 ;;; Generalized Erf distributes over bags
2363 ;;; Generalized Erf is antisymmetric Erf(z1,z2) = - Erf(z2,z1)
2369 ;; Let's remove built-in symbols from list for user-defined properties.
2372 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2376 ;; derivative wrt z1
2380 ;; derviative wrt z2
2384 grad)
2386 ;;; ----------------------------------------------------------------------------
2391 ;; Look for the limit of the arguments.
2395 ;; Handle infinities at this place.
2409 ;; All other cases are handled by the simplifier of the function.
2412 ;;; ----------------------------------------------------------------------------
2417 ;; Check for specific values
2435 ;; Check for numerical evaluation. Use erf(z1,z2) = erf(z2)-erf(z1)
2458 ;; Argument simplification
2463 ;; Representation through more general functions
2465 ($hypergeometric_representation
2468 ;; Transformation to Erf
2470 ($erf_representation
2476 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2477 ;;;
2478 ;;; Implementation of the Complementary Error function Erfc(z)
2479 ;;;
2480 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2484 ;;; Complementary Error function distributes over bags
2488 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2495 grad)
2497 ;;; Integral of the Error function erfc
2506 integral)
2508 ;;; ----------------------------------------------------------------------------
2513 ;; Look for the limit of the arguments.
2516 ;; Handle infinities at this place.
2533 ;; All other cases are handled by the simplifier of the function.
2536 ;;; ----------------------------------------------------------------------------
2541 ;; Check for specific values
2547 ;; Check for numerical evaluation.
2552 ;; Argument simplification
2558 ;; Representation through more general functions
2560 ($hypergeometric_representation
2563 ;; Transformation to Erf or Erfi
2568 (%erf
2570 (%erfi
2578 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2579 ;;;
2580 ;;; Implementation of the Imaginary Error function Erfi(z)
2581 ;;;
2582 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2584 ;;; erfi has mirror symmetry
2588 ;;; erfi is an odd function
2592 ;;; erfi distributes over bags
2596 ;;; Derivative of the Error function erfi
2603 grad)
2605 ;;; Integral of the Error function erfi
2614 integral)
2616 ;;; ----------------------------------------------------------------------------
2621 ;; Look for the limit of the arguments.
2624 ;; Handle infinities at this place.
2628 ;; All other cases are handled by the simplifier of the function.
2631 ;;; ----------------------------------------------------------------------------
2636 ;; Use the relationship erfi(z) = -%i*erf(%i*z)
2640 ;; Take care to return real results when the argument is real.
2649 ;;; erfi is an simplifying function
2654 ;; Check for specific values
2660 ;; Check for numerical evaluation. Use erfi(z) = -%i*erf(%i*z).
2665 ;; Argument simplification
2673 ;; Representation through more general functions
2675 ($hypergeometric_representation
2678 ;; Transformation to Erf or Erfc
2683 (%erf
2685 (%erfc
2693 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2694 ;;;
2695 ;;; The implementation of the Inverse Error function
2696 ;;;
2697 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2699 ;;; The Inverse Error function distributes over bags
2703 ;;; inverse_erf is the inverse of the erf function
2708 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2710 ;;; Differentiation of the Inverse Error function
2718 grad)
2720 ;;; Integral of the Inverse Error function
2727 integral)
2729 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2731 ;;; We support a simplim%function. The function is looked up in simplimit and
2732 ;;; handles specific values of the function.
2737 ;; Look for the limit of the argument.
2740 ;; Handle an argument 1 at this place
2744 ;; All other cases are handled by the simplifier of the function.
2747 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2749 ;;; The Inverse Error function is a simplifying function
2764 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2765 ;;;
2766 ;;; The implementation of the Inverse Complementary Error function
2767 ;;;
2768 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2770 ;;; Inverse Complementary Error function distributes over bags
2774 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2776 ;;; inverse_erfc is the inverse of the erfc function
2782 ;;; Differentiation of the Inverse Complementary Error function
2790 grad)
2792 ;;; Integral of the Inverse Complementary Error function
2800 integral)
2802 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2804 ;;; We support a simplim%function. The function is looked up in simplimit and
2805 ;;; handles specific values of the function.
2810 ;; Look for the limit of the argument.
2813 ;; Handle an argument 0 at this place
2820 ;; All other cases are handled by the simplifier of the function.
2823 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2825 ;;; Inverse Complementary Error function is a simplifying function
2839 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2841 ;;; Implementation of the Newton algorithm to calculate numerical values
2842 ;;; of the Inverse Error functions in float or bigfloat precision.
2843 ;;; The algorithm is a modified version of the routine in newton1.mac.
2880 ;; Compute inverse_erf(z) for z a real or complex number, including
2881 ;; bigfloat objects. The value is computing using a Newton iteration
2882 ;; to solve erf(x) = z.
2885 z)
2889 z))
2891 ;; inverse_erf is odd because erf is.
2896 ;; Find an approximate solution for x = inverse_erf(z).
2899 ;; For small z, inverse_erf(z) = z*sqrt(%pi)/2
2900 ;; + O(z^3). Thus, x = z*sqrt(%pi)/2 is our
2901 ;; initial starting point.
2904 ;; For |z| > 1 and realpart(z) >= 0, we have
2905 ;; the asymptotic series z = erf(x) = 1 -
2906 ;; exp(-x^2)/x/sqrt(%pi).
2907 ;;
2908 ;; Then
2909 ;; x = sqrt(-log(x*sqrt(%pi)*(1-z))
2910 ;;
2911 ;; We can use this as a fixed-point iteration
2912 ;; to find x, and we can start the iteration at
2913 ;; x = 1. Just do one more iteration. I (RLT)
2914 ;; think that's close enough to get the Newton
2915 ;; algorithm to converge.
2922 result)))
2925 ;; Try to pick a reasonable epsilon value for the
2926 ;; Newton iteration.
2937 ;; Derivative of erf(x)
2940 #'diff
2941 z
2943 eps)))))))
2949 z))
2955 ;; Find an approximate solution for x =
2956 ;; inverse_erfc(z). We have inverse_erfc(z) =
2957 ;; inverse_erf(1-z), so that's a good starting point.
2958 ;; We don't need full precision, so a float value is
2959 ;; good enough. But if 1-z is 1, inverse_erf is
2960 ;; undefined, so we need to do something else.
2971 result)))
2977 ;; Derivative of erfc(x)
2980 #'diff
2981 z
2983 eps)))))))
2985 ;; Newton iteration for solving f(x) = z, given f and the derivative
2986 ;; of f.
2998 (intl:gettext "bf-newton: failed to converge after ~M iterations: delta = ~S, x = ~S eps = ~S")
2999 count delta x eps)
3000 x))
3008 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3009 ;;;
3010 ;;; Implementation of the Fresnel Integral S(z)
3011 ;;;
3012 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3014 ;;; fresnel_s distributes over bags
3018 ;;; fresnel_s has mirror symmetry
3022 ;;; fresnel_s is an odd function
3023 ;;;
3024 ;;; Maxima has two mechanism to define a reflection property
3025 ;;; 1. Declare the feature oddfun or evenfun
3026 ;;; 2. Put a reflection rule on the property list
3027 ;;;
3028 ;;; The second way is used for the trig functions. We use it here too.
3030 ;;; This would be the first way to give the property of an odd function.
3031 ;(eval-when
3032 ; (:load-toplevel :execute)
3033 ; (let (($context '$global) (context '$global))
3034 ; (meval '(($declare) %fresnel_s $oddfun))
3035 ; ;; Let's remove built-in symbols from list for user-defined properties.
3036 ; (setq $props (remove '%fresnel_s $props))))
3040 ;;; Differentiation of the Fresnel Integral S
3045 grad)
3047 ;;; Integration of the Fresnel Integral S
3056 integral)
3058 ;;; Limits of the Fresnel Integral S
3070 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3076 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3077 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3086 ;; For very small x, we have
3087 ;; fresnel_s(x) = %pi/6*z^3
3088 ;; fresnel_c(x) = x
3106 (maxima::merror (intl:gettext "fresnel: series expansion failed for (COMPLEX-BFLOAT-FRESNEL ~:M).") z))
3125 (* sqrt-pi
3126 z)))
3128 (* sqrt-pi
3129 z)))
3140 ;; A pure imaginary argument. Use fresnel_s(%i*x)=-%i*fresnel_s(x).
3146 fs))))
3150 ;; A pure imaginary argument. Use fresnel_c(%i*x)=%i*fresnel_c(x).
3156 fc))))
3160 ;;; fresnel_s is a simplifying function
3164 ;; Check for specific values
3170 ;; Check for numerical evaluation
3174 ;; Check for argument simplification
3181 ;; Representation through equivalent functions
3183 ($erf_representation
3198 ($hypergeometric_representation
3207 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3208 ;;;
3209 ;;; Implementation of the Fresnel Integral C(z)
3210 ;;;
3211 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3213 ;;; fresnel_c distributes over bags
3217 ;;; fresnel_c has mirror symmetry
3221 ;;; fresnel_c is an odd function
3222 ;;; Maxima has two mechanism to define a reflection property
3223 ;;; 1. Declare the feature oddfun or evenfun
3224 ;;; 2. Put a reflection rule on the property list
3225 ;;;
3226 ;;; The second way is used for the trig functions. We use it here too.
3230 ;;; Differentiation of the Fresnel Integral C
3235 grad)
3237 ;;; Integration of the Fresnel Integral C
3246 integral)
3248 ;;; Limits of the Fresnel Integral C
3260 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3262 ;;; fresnel_c is a simplifying function
3266 ;; Check for specific values
3272 ;; Check for numerical evaluation
3277 ;; Check for argument simplification
3284 ;; Representation through equivalent functions
3286 ($erf_representation
3301 ($hypergeometric_representation
3310 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3311 ;;;
3312 ;;; Implementation of the Beta function
3313 ;;;
3314 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3316 ;;; The code for the implementation of the beta function is in the files
3317 ;;; csimp2.lisp, simp.lisp and mactex.lisp.
3318 ;;; At this place we only implement the operator property SYMMETRIC.
3320 ;;; Beta is symmetric beta(a,b) = beta(b,a)
3326 ;; Let's remove built-in symbols from list for user-defined properties.
3329 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3330 ;;;
3331 ;;; Implementation of the Incomplete Beta function
3332 ;;;
3333 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3338 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3340 ;;; beta_incomplete distributes over bags
3344 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3348 ;; Derivative wrt a
3356 z)
3358 ;; Derivative wrt b
3375 ;; The derivative wrt z
3379 grad)
3381 ;;; Integral of the Incomplete Beta function
3385 nil
3386 nil
3390 integral)
3392 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3403 ;; Check for specific values
3411 a b z))
3413 z)
3418 ;; z=1 and realpart(b)>0. Simplify to a Beta function.
3419 ;; If we have to evaluate numerically preserve the type.
3434 ;; The argument a is zero or a is negative and the argument b is
3435 ;; not in a valid range. Beta incomplete is undefined.
3436 ;; It would be more correct to return Complex infinity.
3441 ;; Check for numerical evaluation in Float or Bigfloat precision
3449 ;; Negative integer a and b is in a valid range. Expand.
3450 ($rectform
3461 ;; Negative integer a and b is in a valid range. Expand.
3462 ($rectform
3466 ;; Argument simplifications and transformations
3472 ;; Expand for b a positive integer and a not a negative integer.
3486 ;; Expand for a a positive integer.
3502 ;; Expand for a a negative integer and b an integer with b <= -a.
3531 index))
3534 index)))
3560 index)))
3565 ($hypergeometric_representation
3566 ;; There are several cases to consider.
3567 ;;
3568 ;; For -a not an integer see
3569 ;; http://functions.wolfram.com/06.19.26.0005.01
3570 ;;
3571 ;; beta_incomplete(a,b,z) = z^a/a*%f[2,1]([a,1-b],[a+1],z)
3572 ;;
3573 ;; For -b not an integer see
3574 ;; http://functions.wolfram.com/06.19.26.0007.01
3575 ;;
3576 ;; beta_incomplete(a,b,z) = beta(a,b) -
3577 ;; (1-z)^b*z^a/b*%f[2,1]([1,a+b],[b+1],1-z)
3578 ;;
3579 ;; For a+b not a positive integer see
3580 ;; http://functions.wolfram.com/06.19.26.0008.01
3581 ;;
3582 ;; beta_incomplete(a,b,z) =
3583 ;; z^a*(-z)^(b-1)/(a+b-1)*%f[2,1]([1-b,-a-b+1],[-a-b+2],1/z)
3584 ;; + z^a*beta(1-a-b,a)*(-z)^(-a)
3585 ;;
3586 ;; We need to handle more cases here
3588 a)
3592 z)))
3608 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3615 ($rectform
3644 (merror (intl:gettext "beta_incomplete: continued fractions failed for beta_incomplete(~:M, ~:M, ~:M).") a b z))
3682 result))))))
3684 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3685 ;;;
3686 ;;; Implementation of the Generalized Incomplete Beta function
3687 ;;;
3688 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3690 ;;; beta_incomplete_generalized distributes over bags
3694 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3696 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
3697 ;;; but not on the negative real axis and for z1 or z2 real and > 1.
3698 ;;; We support a conjugate-function which test these cases.
3701 conjugate-beta-incomplete-generalized conjugate-function)
3725 a b z1 z2)))))))
3727 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3731 ;; Derivative wrt a
3743 z1)
3749 z2)
3752 ;; Derivative wrt b
3775 ;; The derivative wrt z1
3781 ;; The derivative wrt z2
3787 grad)
3789 ;;; Integral of the Incomplete Beta function
3793 nil
3794 nil
3795 ;; Integral wrt z1
3799 ;; Integral wrt z2
3804 integral)
3806 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3812 ;; Check for specific values
3820 a b z1 z2))
3832 a b z1 z2))
3854 ;; Check for numerical evaluation in Float or Bigfloat precision
3867 ;; Check for argument simplifications and transformations
3936 index))
3939 index)))
3968 index)))
3979 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3980 ;;;
3981 ;;; Implementation of the Regularized Incomplete Beta function
3982 ;;;
3983 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3985 ;;; beta_incomplete_regularized distributes over bags
3989 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3993 ;; Derivative wrt a
4009 ;; Derivative wrt b
4026 ;; The derivative wrt z
4031 grad)
4033 ;;; Integral of the Generalized Incomplete Beta function
4037 nil
4038 nil
4039 ;; Integral wrt z
4045 integral)
4047 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4053 ;; Check for specific values
4061 a b z))
4079 ;; Problem: for b a negative integer the Regularized Incomplete
4080 ;; Beta function is defined to be zero. BUT: When we calculate
4081 ;; e.g. beta_incomplete(1.0,-2.0,1/2)/beta(1.0,-2.0) we get the
4082 ;; result -3.0, because beta_incomplete and beta are defined for
4083 ;; for this case. How do we get a consistent behaviour?
4089 ;; TODO: The following line does not work for bigfloats.
4091 ;; Does $beta_incomplete or simpbeta underflow in this case?
4097 ;; Check for numerical evaluation in Float or Bigfloat precision
4101 beta ibeta )
4105 ;; In case of underflow (see bug #2999) or precision loss use bigfloats
4106 ;; and emporarily give some extra precision but avoid fpprec dependency.
4107 ;; Is this workaround correct for complex values?
4116 ($rectform
4120 ;; Check for argument simplifications and transformations
4144 index))
4147 index)))
4167 index)))
4175 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;