| blob | 8f5f30ddc08e360ec2f61820cd383facb6a2b70b |
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., 675 Mass Ave, Cambridge, MA 02139, USA.
70 ;;;
71 ;;; Copyright (C) 2008 Dieter Kaiser
72 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
78 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
92 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
93 ;;; The following functions test if numerical evaluation has to be done.
94 ;;; The functions should help to test for numerical evaluation more consitent
95 ;;; and without complicated conditional tests including more than one or two
96 ;;; arguments.
97 ;;;
98 ;;; The functions take a list of arguments. All arguments have to be a CL or
99 ;;; Maxima number. If all arguments are numbers we have two cases:
100 ;;; 1. $numer is T we return T. The function has to be evaluated numerically.
101 ;;; 2. One of the args is a float or a bigfloat. Evaluate numerically.
102 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
104 ;;; Test for numerically evaluation in float precision
114 ;;; Test for numerically evaluation in complex float precision
117 "Determine if ARGS consists of numerical values by determining if
118 the real and imaginary parts of each arg are nuemrical (but not
119 bigfloats). A non-NIL result is returned if at least one of args is
120 a floating-point value or if numer is true. If the result is
121 non-NIL, it is a list of the arguments reduced via COMPLEX-NUMBER-P"
128 ;; Always save the result from complex-number-p. But for backward
129 ;; compatibility, only set the flag if any item is a float.
133 ;; Return the values in the same order as the args!
136 ;;; Test for numerically evaluation in bigfloat precision
147 ;;; Test for numerically evaluation in complex bigfloat precision
150 "Determine if ARGS consists of numerical values by determining if
151 the real and imaginary parts of each arg are nuemrical (including
152 bigfloats). A non-NIL result is returned if at least one of args is
153 a floating-point value or if numer is true. If the result is
154 non-NIL, it is a list of the arguments reduced via COMPLEX-NUMBER-P."
162 ;; Always save the result from complex-number-p. But for backward
163 ;; compatibility, only set the flag if any item is a bfloat.
168 ;; Return the values in the same order as the args!
171 ;;; Test for numerical evaluation in any precision, real or complex.
178 ;;; Check for an integer or a float or bigfloat representation. When we
179 ;;; have a float or bigfloat representation return the integer value.
188 val)
191 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
193 ;;; The changes to the parser to connect the operator !! to double_factorial(z)
195 ;(def-mheader |$!!| (%double_factorial))
197 ;(def-led (|$!!| 160.) (op left)
198 ; (list '$expr
199 ; (mheader '$!!)
200 ; (convert left '$expr)))
202 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
203 ;;;
204 ;;; The implementation of the function Double factorial
205 ;;;
206 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
211 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
213 ;;; Set properties to give full support to the parser and display
221 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
223 ;;; Double factorial is a simplifying function
227 ;;; Double factorial distributes over bags
231 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
233 ;;; Double factorial has mirror symmetry
237 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
239 ;;; Differentiation of Double factorial
255 grad)
257 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
264 ;; Positive Integer less then $factlim or $factlim is -1. Call gfact.
270 ;; Even negative integer or real representation. Not defined.
281 ;; Odd negative integer, float or complex float.
288 ;; bigfloat or complex bigfloat.
292 ;; double_factorial(inf) -> inf
302 ;; Special case double_factorial(n-1)
303 ;; Not sure if this simplification is useful.
307 ;; Special case double_factorial(2*k+n), k integer
320 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
321 ;;; Double factorial for a complex float argument. The result is a CL complex.
322 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
326 (*
333 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
334 ;;; Double factorial for a bigfloat or complex bigfloat argument
335 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
352 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
353 ;;;
354 ;;; The implementation of the Incomplete Gamma function
355 ;;;
356 ;;; A&S 6.5.3:
357 ;;;
358 ;;; inf
359 ;;; /
360 ;;; [ a - 1 - t
361 ;;; gamma_incomplete(a, z) = I t %e dt
362 ;;; ]
363 ;;; /
364 ;;; z
365 ;;;
366 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
373 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
375 ;;; Set properties to give full support to the parser and display
383 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
385 ;;; Incomplete Gamma function is a simplifying function
389 ;;; Incomplete Gamma distributes over bags
393 ;;; Incomplete Gamma function has not mirror symmetry for z on the negative
394 ;;; real axis. We support a conjugate-function which test this case.
401 ;; Definitely not on the negative real axis for z. Mirror symmetry.
408 ;; On the negative real axis or no information. Unsimplified.
413 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
415 ;;; Derivative of the Incomplete Gamma function
421 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2
422 ;; function and the Generalized Incomplete Gamma function
423 ;; (functions.wolfram.com), only for a>0.
439 ;; No derivative. Maxima generates a noun form.
440 nil)))
441 ;; The derivative wrt z
447 ;;; Integral of the Incomplete Gamma function
451 nil
455 integral)
457 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
459 ;;; We support a simplim%function. The function is looked up in simplimit and
460 ;;; handles specific values of the function.
465 ;; Look for the limit of the arguments.
476 ;; Handle an argument 0 at this place.
483 ;; Call the simplifier of the function.
486 ;; All other cases are handled by the simplifier of the function.
489 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
491 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
492 ;;;
493 ;;; Implementation of the Lower Incomplete gamma function:
494 ;;;
495 ;;; A&S 6.5.2
496 ;;;
497 ;;; z
498 ;;; /
499 ;;; [ a - 1 - t
500 ;;; gamma_incomplete_lower(a, z) = I t %e dt
501 ;;; ]
502 ;;; /
503 ;;; 0
504 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
508 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
516 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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.
548 ;; Try to express gamma_incomplete_lower(a,z) in simpler terms for
549 ;; special values of the order a. If we can't, return NIL to say so.
552 ;; gamma_incomplete_lower(0, x) = 1-exp(x)
555 ;; gamma_incomplete_lower(a,z) can be simplified if a is a
556 ;; positive integer. Since gamma_incomplete_lower(a,z) =
557 ;; gamma(a) - gamma_incomplete(a,z), use gamma_incomplete to
558 ;; get the result.
562 ;; A&S 6.5.12:
563 ;;
564 ;; gamma_incomplete_lower(1/2,x^2) = sqrt(%pi)*erf(x)
565 ;;
566 ;; gamma_incomplete_lower(1/2,z) = sqrt(%pi)*erf(sqrt(x))
567 ;;
571 ;; gamma_incomplete_lower(a+n,z), where n is an integer
576 ;; gamma_incomplete_lower(a+n,z). See DLMF 8.8.7:
577 ;;
578 ;; gamma_incomplete_lower(a+n,z)
579 ;; = pochhammer(a,n)*gamma_incomplete_lower(a,z)
580 ;; -z^a*exp(-z)*sum(gamma(a+n)gamma(a+k+1)*z^k,k,0,n-1)
599 ;; See DLMF 8.8.8. For simplicity let g(a,z) = gamma_incomplete_lower(a,z).
600 ;;
601 ;; g(a,z) = gamma(a)/gamma(a-n)*g(a-n,z)
602 ;; - z^(a-1)*exp(z)*sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
603 ;;
604 ;; Rewrite:
605 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
606 ;; + gamma(a-n)/gamma(a)*z^(a-1)*exp(-z)
607 ;; * sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
608 ;; Or
609 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
610 ;; + z^(a-1)*exp(-z)
611 ;; * sum(gamma(a-n)/gamma(a-k)*z^(-k),k,0,n-1)
632 ;; gamma_incomplete_lower of (numeric) rational order.
633 ;; Expand it out so that the resulting order is between 0 and
634 ;; 1.
637 ;; a = n + order where 0 <= order < 1.
641 ;; Nothing to do if the order is already between 0 and 1
644 ;; Use gamma_incomplete(a+n,z) above. and then substitue
645 ;; a=order. This works for n positive or negative.
650 ;; No expansion so return nil to indicate that
651 nil)))
653 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
664 ;; Check for specific values
667 ;; gamma_incomplete(v,0) is gamma(v) only if the realpart(v) >
668 ;; 0. If realpart(v) <= 0, gamma_incomplete is undefined. For
669 ;; all other cases, return the noun form.
675 a z))
684 ;; Check for numerical evaluation in Float or Bigfloat precision
689 ;; a and z are Maxima numbers, at least one has a float value
696 ;; a is zero or a negative float representing an integer.
697 ;; For these cases the numerical routines of gamma-incomplete
698 ;; do not work. Call the numerical routine for the Exponential
699 ;; Integral E(n,z). The routine is called with a positive integer!.
709 ;; a and z are Maxima numbers, at least one is a complex value and
710 ;; we have at least one float part
718 ;; Call expintegral-e. See comment above.
733 ;; Call bfloat-expintegral-e. See comment above.
754 ;; Call bfloat-expintegral-e. See comment above.
764 ;; Check for transformations and argument simplification
767 ;; Integer or a symbol declared to be an integer. Expand in a series.
789 ;; Simplify gamma, but do not expand to avoid division
790 ;; by zero.
818 ;; We have a half integral order and $gamma_expand is not NIL.
819 ;; We expand in a series with the Erfc function
869 ;; Handle gamma_incomplete(a+n, z), where n is a numerical integer
874 ;; See DLMF 8.8.9: https://dlmf.nist.gov/8.8.E9
875 ;;
876 ;; gamma_incomplete(a+n,z) = pochhammer(a,n)*gamma_incomplete(a,z)
877 ;; + z^a*exp(-z)*sum(gamma(a+n)/gamma(a+k+1)*z^k,k,0,n-1)
895 ;; See http://functions.wolfram.com/06.06.17.0004.01
896 ;;
897 ;; gamma_incomplate(a,z) =
898 ;; (-1)^n*pochhammer(1-a,n)
899 ;; *[gamma_incomplete(a-n,z)
900 ;; + z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)]
901 ;;
902 ;; Rerarrange this in terms of gamma_incomplete(a-n,z):
903 ;;
904 ;; gamma_incomplete(a-n,z) =
905 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
906 ;; -z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)
907 ;;
908 ;; Change the summation index to go from k = 0 to n-1:
909 ;;
910 ;; z^(a-n-1)*sum(z^k/pochhammer(a-n,k),k,1,n)
911 ;; = z^(a-n-1)*sum(z^(k+1)/pochhammer(a-n,k+1),k,0,n-1)
912 ;; = z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
913 ;;
914 ;; Thuus:
915 ;; gamma_incomplete(a-n,z) =
916 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
917 ;; -z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
919 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
928 ;; sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
938 ;; gamma_incomplete of (numeric) rational order. Expand it out
939 ;; so that the resulting order is between 0 and 1.
942 ;; a = n + order where 0 <= order < 1.
946 ;; Nothing to do if the order is already between 0 and 1
949 ;; Use gamma_incomplete(a+n,z) above. and then substitue
950 ;; a=order. This works for n positive or negative.
957 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
958 ;;; Numerical evaluation of the Incomplete Gamma function
959 ;;;
960 ;;; gamma-incomplete (a,z) - real and complex double float
961 ;;; bfloat-gamma-incomplete (a z) - bigfloat
962 ;;; complex-bfloat-gamma-incomplete (a z) - complex bigfloat
963 ;;;
964 ;;; Expansion in a power series for abs(x) < R, where R is
965 ;;; *gamma-radius* + real(a) if real(a) > 0 or *gamma-radius*
966 ;;; otherwise.
967 ;;;
968 ;;; (A&S 6.5.29):
969 ;;;
970 ;;; inf
971 ;;; ===
972 ;;; \ gamma(a)
973 ;;; gamma(a,z) = exp(-x)*z^a * > ------------ * z^n
974 ;;; / gamma(a+1+n)
975 ;;; ===
976 ;;; n=0
977 ;;;
978 ;;; This expansion does not work for an integer a<=0, because the Gamma function
979 ;;; in the denominator is not defined for a=0 and negative integers. For this
980 ;;; case we use the Exponential Integral E for numerically evaluation. The
981 ;;; Incomplete Gamma function and the Exponential integral are connected by
982 ;;;
983 ;;; gamma(a,z) = z^a * expintegral_e(1-a,z)
984 ;;;
985 ;;; When the series is not used, two forms of the continued fraction
986 ;;; are used. When z is not near the negative real axis use the
987 ;;; continued fractions (A&S 6.5.31):
988 ;;;
989 ;;; 1 1-a 1 2-a 2
990 ;;; gamma(a,z) = exp(-z) z^a *( -- --- --- --- --- ... )
991 ;;; z+ 1+ z+ 1+ z+
992 ;;;
993 ;;; The accuracy is controlled by *gamma-incomplete-eps* for double float
994 ;;; precision. For bigfloat precision epsilon is 10^(-$fpprec). The expansions
995 ;;; in a power series or continued fractions stops if *gamma-incomplete-maxit*
996 ;;; is exceeded and an Maxima error is thrown.
997 ;;;
998 ;;; The above fraction does not converge on the negative real axis and
999 ;;; converges very slowly near the axis. In this case, use the
1000 ;;; relationship
1001 ;;;
1002 ;;; gamma(a,z) = gamma(a) - gamma_lower(a,z)
1003 ;;;
1004 ;;; The continued fraction for gamma_incomplete_lower(a,z) is
1005 ;;; (http://functions.wolfram.com/06.06.10.0009.01):
1006 ;;;
1007 ;;; gamma_lower(a,z) = exp(-z) * z^a / cf(a,z)
1008 ;;;
1009 ;;; where
1010 ;;;
1011 ;;; -a*z z (a+1)*z 2*z (a+2)*z
1012 ;;; cf(a,z) = a + ---- ---- ------- ---- -------
1013 ;;; a+1+ a+2- a+3+ a+4- a+5+
1014 ;;;
1015 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1025 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1027 ;;; The numerical evaluation for CL float or complex values a and x
1028 ;;; When the flag regularized is T, the result is divided by gamma(a) and
1029 ;;; Maxima returns the numercial result for gamma_incomplete_regularized
1035 ;; Compute the factor needed to scale the series or continued
1036 ;; fraction. This is x^a*exp(-x) or x^a*exp(-x)/gamma(a)
1037 ;; depending on whether we want a non-regularized or
1038 ;; regularized form. We want to compute the factor carefully
1039 ;; to avoid unnecessary overflow if possible.
1043 ;; gammafloat is more accurate for real
1044 ;; values of a.
1051 ;; Easy way failed. Use logs to compute the
1052 ;; result. This loses some precision, especially
1053 ;; for large values of a and/or x because we use
1054 ;; many bits to hold the exponent part.
1056 x
1058 a
1064 ;; Easy way failed, so use the log form.
1066 x)))))))
1070 ;; %gamma-incomplete compute the lower incomplete gamma
1071 ;; function, so we need to subtract that from gamma(a),
1072 ;; more or less.
1080 ;; Continued fraction used. Just multiply by the factor
1081 ;; to get the final result.
1084 ;; Compute the key part of the gamma incomplete function using either
1085 ;; a series expression or a continued fraction expression. Two values
1086 ;; are returned: the value itself and a boolean, indicating what the
1087 ;; computed value is. If the boolean non-NIL, then the computed value
1088 ;; is the lower incomplete gamma function.
1096 ;; The series expansion is done for x within a circle of radius
1097 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1098 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1099 ;; continued fraction is used.
1105 ;; For x near the negative real axis, use the
1106 ;; relationship gamma_incomplete(a,z) = gamma(a) -
1107 ;; gamma_incomplete_lower(a,z), where
1108 ;; gamma_incomplete_lower(a,z) is the lower poart of the
1109 ;; incomplete gamma function. We can evaluate that
1110 ;; using a continued fraction from
1111 ;; http://functions.wolfram.com/06.06.10.0009.01. (Note
1112 ;; that the alternative fraction,
1113 ;; http://functions.wolfram.com/06.06.10.0007.01,
1114 ;; appears to be less accurate.)
1115 ;;
1116 ;; Also note that this appears to be valid for all
1117 ;; values of x (real or complex), but we don't want to
1118 ;; use this everywhere for gamma_incomplete. Consider
1119 ;; what happens for large real x. gamma_incomplete(a,x)
1120 ;; is small, but gamma_incomplete(a,x) = gamma(x) - cf
1121 ;; will have large roundoff errors.
1134 t)))
1136 ;; Expansion in continued fractions for gamma_incomplete.
1147 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1156 ;; Return nil to indicate we used the continued fraction.
1159 ;; Expansion in a series
1167 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1170 ;; Return T to indicate we used the series and the series
1171 ;; is for the integral from 0 to x, not x to inf.
1174 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1176 ;;; This function is called for a and x real
1184 ;; The series expansion is done for x within a circle of radius
1185 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1186 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1187 ;; continued fraction is used.
1194 ;; Expansion in continued fractions of the Incomplete Gamma function
1203 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1224 ;; Expand to multiply everything out.
1225 ($expand
1226 ;; Expansion in continued fraction for the lower incomplete gamma.
1228 ;; NOTE: We want (power x a) instead of bigfloat:expt
1229 ;; because this preserves how maxima computes x^a when
1230 ;; x is negative and a is rational. For, example
1231 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1232 ;; principal value.
1245 x))))))))))))))
1248 ;; Series expansion of the Incomplete Gamma function
1254 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1268 ($rectform
1273 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1285 ;; The series expansion is done for x within a circle of radius
1286 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1287 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1288 ;; continued fraction is used.
1300 ;; Expansion in continued fractions of the Incomplete Gamma function
1311 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1322 gm-eps))
1331 ;; Expand to multiply everything out.
1332 ($expand
1333 ;; Expansion in continued fraction for the lower incomplete gamma.
1335 ;; NOTE: We want (power x a) instead of bigfloat:expt
1336 ;; because this preserves how maxima computes x^a when
1337 ;; x is negative and a is rational. For, example
1338 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1339 ;; principal value.
1352 x))))))))))))))
1354 ;; Series expansion of the Incomplete Gamma function
1362 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1375 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1376 ;;;
1377 ;;; Implementation of the Generalized Incomplete Gamma function
1378 ;;;
1379 ;;; z2
1380 ;;; /
1381 ;;; [ a - 1 - t
1382 ;;; gamma_incomplete_generalized(a, z1, z2) = I t %e dt
1383 ;;; ]
1384 ;;; /
1385 ;;; z1
1386 ;;;
1387 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1392 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1394 ;;; Set the properties alias, reversealias, noun and verb
1400 $gamma_incomplete_generalized reversealias)
1402 $gamma_incomplete_generalized noun)
1404 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1405 ;;; on the negative real axis.
1406 ;;; We support a conjugate-function which test this case.
1409 conjugate-gamma-incomplete-generalized conjugate-function)
1414 ;; z1 and z2 definitely not on the negative real axis.
1415 ;; Mirror symmetry.
1423 ;; On the negative real axis or no information. Unsimplified.
1428 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1430 ;;; Generalized Incomplete Gamma function is a simplifying function
1433 simp-gamma-incomplete-generalized operators)
1435 ;;; Generalized Incomplete Gamma distributes over bags
1439 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1441 ;;; Differentiation of Generalized Incomplete Gamma function
1445 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1446 ;; and the Generalized Incomplete Gamma function (functions.wolfram.com)
1468 ;; The derivative wrt z1
1472 ;; The derivative wrt z2
1476 grad)
1478 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1490 ;; Check for specific values
1517 ;; Check for numerical evaluation in Float or Bigfloat precision
1518 ;; Use the numerical routines of the Incomplete Gamma function
1545 ;; Check for transformations and argument simplification
1548 ;; Expand gamma_incomplete_generalized(a+n,z1,z2) with n an integer
1601 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1602 ;;;
1603 ;;; Implementation of the Regularized Incomplete Gamma function
1604 ;;;
1605 ;;; A&S 6.5.1
1606 ;;;
1607 ;;; gamma_incomplete(a, z)
1608 ;;; gamma_incomplete_regularized(a, z) = ----------------------
1609 ;;; gamma(a)
1610 ;;;
1611 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1616 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1622 $gamma_incomplete_regularized reversealias)
1624 $gamma_incomplete_regularized noun)
1626 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1628 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1629 ;;; on the negative real axis.
1630 ;;; We support a conjugate-function which test this case.
1633 conjugate-gamma-incomplete-regularized conjugate-function)
1638 ;; z definitely not on the negative real axis. Mirror symmetry.
1645 ;; On the negative real axis or no information. Unsimplified.
1650 ;;; Regularized Incomplete Gamma function is a simplifying function
1653 simp-gamma-incomplete-regularized operators)
1655 ;;; Regularized Incomplete Gamma distributes over bags
1659 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1661 ;;; Differentiation of Regularized Incomplete Gamma function
1665 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1666 ;; and the Regularized Generalized Incomplete Gamma function
1667 ;; (functions.wolfram.com)
1681 ;; The derivative wrt z
1683 -1
1687 grad)
1689 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1701 ;; Check for specific values
1709 a z))
1716 ;; Check for numerical evaluation in Float or Bigfloat precision
1720 ;; gamma_incomplete returns a regularized result
1726 ;; gamma_incomplete returns a regularized result
1738 ($rectform
1743 ;; Check for transformations and argument simplification
1746 ;; An integer. Expand the expression.
1763 ;; We have a half integral order and $gamma_expand is not NIL.
1764 ;; We expand in a series with the Erfc function
1817 ;; We factor the second summand.
1818 ;; Some factors vanish and the result is more readable.
1819 ($factor
1834 ;; We factor the second summand.
1835 ($factor
1849 ;; gamma_incomplete_regularized of numeric rational order.
1850 ;; Expand it out so that the resulting order is between 0 and
1851 ;; 1. Use gamma_incomplete_regularized(a+n,z) to do the dirty
1852 ;; work.
1855 ;; a = n + order where 0 <= order < 1.
1859 ;; Nothing to do if the order is already between 0 and 1
1862 ;; Use gamma_incomplete_regularized(a+n,z) above. and
1863 ;; then substitue a=order. This works for n positive or
1864 ;; negative.
1871 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1872 ;;;
1873 ;;; Implementation of the Logarithm of the Gamma function
1874 ;;;
1875 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1880 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1888 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1892 ;;; Logarithm of the Gamma function distributes over bags
1896 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1901 grad)
1903 ;; integrate(log_gamma(x),x) = psi[-2](x)
1909 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1916 ;; Check for specific values
1922 ;; We have zero, a negative integer or a float or bigfloat representation.
1928 ;; Check for numerical evaluation
1945 ;; Transform to Logarithm of Factorial for integer values
1946 ;; At this point the integer value is positive and not zero.
1953 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1954 ;;; The functions log-gamma-lanczos, bfloat-log-gamma and
1955 ;;; complex-bfloat-log-gamma are modified versions of the related functions
1956 ;;; gamma-lanczos, bffac and cbffac. The functions return the Logarithm of
1957 ;;; the Gamma function. If we have to calculate the quotient of Gamma functions,
1958 ;;; e. g. for the Beta function, it is much more appropriate to use the
1959 ;;; logarithmic versions to avoid overflow.
1960 ;;;
1961 ;;; Be careful log(gamma(z)) is only for realpart(z) positive equal to
1962 ;;; log_gamma(z). For a negative realpart(z) log_gamma differ by multiple of
1963 ;;; %pi from log(gamma(z)). But we always have exp(log_gamma(z))= gamma(z).
1964 ;;; The terms to get the transformation for log_gamma(-z) are taken from
1965 ;;; functions.wolfram.com.
1966 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1972 :initial-contents
1974 57.156235665862923517
1975 -59.597960355475491248
1976 14.136097974741747174
1977 -0.49191381609762019978
1978 .33994649984811888699e-4
1979 .46523628927048575665e-4
1980 -.98374475304879564677e-4
1981 .15808870322491248884e-3
1982 -.21026444172410488319e-3
1983 .21743961811521264320e-3
1984 -.16431810653676389022e-3
1985 .84418223983852743293e-4
1986 -.26190838401581408670e-4
1992 (+
1993 (*
2001 (*
2015 sum)
2021 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2042 bigfloat%pi
2045 bigfloat%pi '$%i
2063 y)))
2069 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2090 bigfloat%pi
2093 bigfloat%pi '$%i
2112 y))))
2113 ($rectform
2119 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2120 ;;;
2121 ;;; Implementation of the Error function Erf(z)
2122 ;;;
2123 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2128 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2136 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2138 ;;; erf has mirror symmetry
2142 ;;; erf is an odd function
2146 ;;; erf is a simplifying function
2150 ;;; erf distributes over bags
2154 ;;; Derivative of the Error function erf
2161 grad)
2163 ;;; Integral of the Error function erf
2172 integral)
2174 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2177 ; See A&S 7.1.21 or http://functions.wolfram.com/06.25.26.0001.01
2179 z
2191 ;; Check for specific values
2197 ;; Check for numerical evaluation
2210 ;; Argument simplification
2219 ;; Representation through more general functions
2221 ($hypergeometric_representation
2224 ;; Transformation to Erfc or Erfi
2229 (%erfc
2231 (%erfi
2239 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2242 ;; We use the slatec routine for float values.
2245 ;; Compute erf(z) using the relationship
2246 ;;
2247 ;; erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2,z^2)/sqrt(%pi))
2248 ;;
2249 ;; When z is real sqrt(z^2)/z is signum(z). For complex z,
2250 ;; sqrt(z^2)/z = 1 if -%pi/2 < arg(z) <= %pi/2 and -1 otherwise.
2251 ;;
2252 ;; This relationship has serious round-off issues when z is small
2253 ;; because gamma_incomplete(1/2,z^2)/sqrt(%pi) is near 1.
2254 ;;
2255 ;; complex-erf is for (lisp) complex numbers; bfloat-erf is for real
2256 ;; bfloats, and complex-bfloat-erf is for complex bfloats. Care is
2257 ;; taken to return real results for real arguments and imaginary
2258 ;; results for imaginary arguments
2261 (*
2263 -1
2269 ;; Pure real argument, the result is real
2272 ;; Pure imaginary argument, the result is pure imaginary
2275 result))))
2278 ;; Warning! This has round-off problems when abs(z) is very small.
2280 ;; The argument is real, the result is real too
2281 ($realpart
2290 ;; Warning! This has round-off problems when abs(z) is very small.
2304 ;; Pure real argument, the result is real
2307 ;; Pure imaginary argument, the result is pure imaginary
2310 ;; A general complex result
2311 result))))
2315 ;; Erf(z) for all z. Z must be a CL real or complex number or a
2316 ;; BIGFLOAT or COMPLEX-BIGFLOAT object. The result will be of the
2317 ;; same type as Z.
2320 ;; Use Slatec derf, which should be faster than the series.
2323 ;; Use the series A&S 7.1.5 for small x:
2324 ;;
2325 ;; erf(z) = 2*z/sqrt(%pi) * sum((-1)^n*z^(2*n)/n!/(2*n+1), n, 0, inf)
2326 ;;
2327 ;; The threshold is approximately erf(x) = 0.5. (Doesn't
2328 ;; have to be super accurate.) This gives max accuracy when
2329 ;; using the identity erf(x) = 1 - erfc(x).
2335 k
2339 ;; The general case.
2346 (bigfloat (maxima::$bfloat (maxima::$expand (maxima::complex-bfloat-erf (maxima::to z))))))))))
2349 ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
2350 ;; near 1. Wolfram says
2351 ;;
2352 ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
2353 ;;
2354 ;; For real(z) > 0, sqrt(z^2)/z is 1 so
2355 ;;
2356 ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2357 ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
2358 ;;
2359 ;; For real(z) < 0, sqrt(z^2)/z is -1 so
2360 ;;
2361 ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2362 ;; = 2 - 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
2384 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2385 ;;;
2386 ;;; Implementation of the Generalized Error function Erf(z1,z2)
2387 ;;;
2388 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2393 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2401 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2403 ;;; Generalized Erf has mirror symmetry
2407 ;;; Generalized Erf distributes over bags
2411 ;;; Generalized Erf is antisymmetric Erf(z1,z2) = - Erf(z2,z1)
2418 ;; Let's remove built-in symbols from list for user-defined properties.
2421 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2425 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2429 ;; derivative wrt z1
2433 ;; derviative wrt z2
2437 grad)
2439 ;;; ----------------------------------------------------------------------------
2444 ;; Look for the limit of the arguments.
2448 ;; Handle infinities at this place.
2462 ;; All other cases are handled by the simplifier of the function.
2465 ;;; ----------------------------------------------------------------------------
2474 ;; Check for specific values
2492 ;; Check for numerical evaluation. Use erf(z1,z2) = erf(z2)-erf(z1)
2515 ;; Argument simplification
2520 ;; Representation through more general functions
2522 ($hypergeometric_representation
2525 ;; Transformation to Erf
2527 ($erf_representation
2533 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2534 ;;;
2535 ;;; Implementation of the Complementary Error function Erfc(z)
2536 ;;;
2537 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2542 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2550 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2554 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2558 ;;; Complementary Error function distributes over bags
2562 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2569 grad)
2571 ;;; Integral of the Error function erfc
2580 integral)
2582 ;;; ----------------------------------------------------------------------------
2587 ;; Look for the limit of the arguments.
2590 ;; Handle infinities at this place.
2607 ;; All other cases are handled by the simplifier of the function.
2610 ;;; ----------------------------------------------------------------------------
2617 ;; Check for specific values
2623 ;; Check for numerical evaluation.
2628 ;; Argument simplification
2634 ;; Representation through more general functions
2636 ($hypergeometric_representation
2639 ;; Transformation to Erf or Erfi
2644 (%erf
2646 (%erfi
2654 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2655 ;;;
2656 ;;; Implementation of the Imaginary Error function Erfi(z)
2657 ;;;
2658 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2663 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2671 ;;; erfi has mirror symmetry
2675 ;;; erfi is an odd function
2679 ;;; erfi is an simplifying function
2683 ;;; erfi distributes over bags
2687 ;;; Derivative of the Error function erfi
2694 grad)
2696 ;;; Integral of the Error function erfi
2705 integral)
2707 ;;; ----------------------------------------------------------------------------
2712 ;; Look for the limit of the arguments.
2715 ;; Handle infinities at this place.
2719 ;; All other cases are handled by the simplifier of the function.
2722 ;;; ----------------------------------------------------------------------------
2727 ;; Use the relationship erfi(z) = -%i*erf(%i*z)
2731 ;; Take care to return real results when the argument is real.
2745 ;; Check for specific values
2751 ;; Check for numerical evaluation. Use erfi(z) = -%i*erf(%i*z).
2756 ;; Argument simplification
2764 ;; Representation through more general functions
2766 ($hypergeometric_representation
2769 ;; Transformation to Erf or Erfc
2774 (%erf
2776 (%erfc
2784 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2785 ;;;
2786 ;;; The implementation of the Inverse Error function
2787 ;;;
2788 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2793 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2795 ;;; Set properties to give full support to the parser and display
2803 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2805 ;;; The Inverse Error function is a simplifying function
2809 ;;; The Inverse Error function distributes over bags
2813 ;;; inverse_erf is the inverse of the erf function
2818 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2820 ;;; Differentiation of the Inverse Error function
2828 grad)
2830 ;;; Integral of the Inverse Error function
2837 integral)
2839 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2841 ;;; We support a simplim%function. The function is looked up in simplimit and
2842 ;;; handles specific values of the function.
2847 ;; Look for the limit of the argument.
2850 ;; Handle an argument 1 at this place
2854 ;; All other cases are handled by the simplifier of the function.
2857 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2874 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2875 ;;;
2876 ;;; The implementation of the Inverse Complementary Error function
2877 ;;;
2878 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2883 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2885 ;;; Set properties to give full support to the parser and display
2893 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2895 ;;; Inverse Complementary Error function is a simplifying function
2899 ;;; Inverse Complementary Error function distributes over bags
2903 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2905 ;;; inverse_erfc is the inverse of the erfc function
2911 ;;; Differentiation of the Inverse Complementary Error function
2919 grad)
2921 ;;; Integral of the Inverse Complementary Error function
2929 integral)
2931 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2933 ;;; We support a simplim%function. The function is looked up in simplimit and
2934 ;;; handles specific values of the function.
2939 ;; Look for the limit of the argument.
2942 ;; Handle an argument 0 at this place
2949 ;; All other cases are handled by the simplifier of the function.
2952 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2969 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2971 ;;; Implementation of the Newton algorithm to calculate numerical values
2972 ;;; of the Inverse Error functions in float or bigfloat precision.
2973 ;;; The algorithm is a modified version of the routine in newton1.mac.
3010 ;; Compute inverse_erf(z) for z a real or complex number, including
3011 ;; bigfloat objects. The value is computing using a Newton iteration
3012 ;; to solve erf(x) = z.
3015 z)
3019 z))
3021 ;; inverse_erf is odd because erf is.
3026 ;; Find an approximate solution for x = inverse_erf(z).
3029 ;; For small z, inverse_erf(z) = z*sqrt(%pi)/2
3030 ;; + O(z^3). Thus, x = z*sqrt(%pi)/2 is our
3031 ;; initial starting point.
3034 ;; For |z| > 1 and realpart(z) >= 0, we have
3035 ;; the asymptotic series z = erf(x) = 1 -
3036 ;; exp(-x^2)/x/sqrt(%pi).
3037 ;;
3038 ;; Then
3039 ;; x = sqrt(-log(x*sqrt(%pi)*(1-z))
3040 ;;
3041 ;; We can use this as a fixed-point iteration
3042 ;; to find x, and we can start the iteration at
3043 ;; x = 1. Just do one more iteration. I (RLT)
3044 ;; think that's close enough to get the Newton
3045 ;; algorithm to converge.
3052 result)))
3055 ;; Try to pick a reasonable epsilon value for the
3056 ;; Newton iteration.
3067 ;; Derivative of erf(x)
3070 #'diff
3071 z
3073 eps)))))))
3079 z))
3085 ;; Find an approximate solution for x =
3086 ;; inverse_erfc(z). We have inverse_erfc(z) =
3087 ;; inverse_erf(1-z), so that's a good starting point.
3088 ;; We don't need full precision, so a float value is
3089 ;; good enough. But if 1-z is 1, inverse_erf is
3090 ;; undefined, so we need to do something else.
3101 result)))
3107 ;; Derivative of erfc(x)
3110 #'diff
3111 z
3113 eps)))))))
3115 ;; Newton iteration for solving f(x) = z, given f and the derivative
3116 ;; of f.
3128 (intl:gettext "bf-newton: failed to converge after ~M iterations: delta = ~S, x = ~S eps = ~S")
3129 count delta x eps)
3130 x))
3138 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3139 ;;;
3140 ;;; Implementation of the Fresnel Integral S(z)
3141 ;;;
3142 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3147 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3149 ;;; Set properties to give full support to the parser and display
3157 ;;; fresnel_s is a simplifying function
3161 ;;; fresnel_s distributes over bags
3165 ;;; fresnel_s has mirror symmetry
3169 ;;; fresnel_s is an odd function
3170 ;;;
3171 ;;; Maxima has two mechanism to define a reflection property
3172 ;;; 1. Declare the feature oddfun or evenfun
3173 ;;; 2. Put a reflection rule on the property list
3174 ;;;
3175 ;;; The second way is used for the trig functions. We use it here too.
3177 ;;; This would be the first way to give the property of an odd function.
3178 ;(eval-when
3179 ; #+gcl (load eval)
3180 ; #-gcl (:load-toplevel :execute)
3181 ; (let (($context '$global) (context '$global))
3182 ; (meval '(($declare) %fresnel_s $oddfun))
3183 ; ;; Let's remove built-in symbols from list for user-defined properties.
3184 ; (setq $props (remove '%fresnel_s $props))))
3188 ;;; Differentiation of the Fresnel Integral S
3193 grad)
3195 ;;; Integration of the Fresnel Integral S
3204 integral)
3206 ;;; Limits of the Fresnel Integral S
3218 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3224 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3225 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3234 ;; For very small x, we have
3235 ;; fresnel_s(x) = %pi/6*z^3
3236 ;; fresnel_c(x) = x
3254 (maxima::merror (intl:gettext "fresnel: series expansion failed for (COMPLEX-BFLOAT-FRESNEL ~:M).") z))
3273 (* sqrt-pi
3274 z)))
3276 (* sqrt-pi
3277 z)))
3288 ;; A pure imaginary argument. Use fresnel_s(%i*x)=-%i*fresnel_s(x).
3294 fs))))
3298 ;; A pure imaginary argument. Use fresnel_c(%i*x)=%i*fresnel_c(x).
3304 fc))))
3312 ;; Check for specific values
3318 ;; Check for numerical evaluation
3322 ;; Check for argument simplification
3329 ;; Representation through equivalent functions
3331 ($erf_representation
3346 ($hypergeometric_representation
3356 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3357 ;;;
3358 ;;; Implementation of the Fresnel Integral C(z)
3359 ;;;
3360 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3365 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3367 ;;; Set properties to give full support to the parser and display
3375 ;;; fresnel_c is a simplifying function
3379 ;;; fresnel_c distributes over bags
3383 ;;; fresnel_c has mirror symmetry
3387 ;;; fresnel_c is an odd function
3388 ;;; Maxima has two mechanism to define a reflection property
3389 ;;; 1. Declare the feature oddfun or evenfun
3390 ;;; 2. Put a reflection rule on the property list
3391 ;;;
3392 ;;; The second way is used for the trig functions. We use it here too.
3396 ;;; Differentiation of the Fresnel Integral C
3401 grad)
3403 ;;; Integration of the Fresnel Integral C
3412 integral)
3414 ;;; Limits of the Fresnel Integral C
3426 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3433 ;; Check for specific values
3439 ;; Check for numerical evaluation
3444 ;; Check for argument simplification
3451 ;; Representation through equivalent functions
3453 ($erf_representation
3468 ($hypergeometric_representation
3479 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3480 ;;;
3481 ;;; Implementation of the Beta function
3482 ;;;
3483 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3485 ;;; The code for the implementation of the beta function is in the files
3486 ;;; csimp2.lisp, simp.lisp and mactex.lisp.
3487 ;;; At this place we only implement the operator property SYMMETRIC.
3489 ;;; Beta is symmetric beta(a,b) = beta(b,a)
3496 ;; Let's remove built-in symbols from list for user-defined properties.
3499 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3500 ;;;
3501 ;;; Implementation of the Incomplete Beta function
3502 ;;;
3503 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3508 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3513 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3521 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3525 ;;; beta_incomplete distributes over bags
3529 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3533 ;; Derivative wrt a
3541 z)
3543 ;; Derivative wrt b
3561 ;; The derivative wrt z
3565 grad)
3567 ;;; Integral of the Incomplete Beta function
3571 nil
3572 nil
3576 integral)
3578 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3594 ;; Check for specific values
3602 a b z))
3604 z)
3609 ;; z=1 and realpart(b)>0. Simplify to a Beta function.
3610 ;; If we have to evaluate numerically preserve the type.
3625 ;; The argument a is zero or a is negative and the argument b is
3626 ;; not in a valid range. Beta incomplete is undefined.
3627 ;; It would be more correct to return Complex infinity.
3632 ;; Check for numerical evaluation in Float or Bigfloat precision
3640 ;; Negative integer a and b is in a valid range. Expand.
3641 ($rectform
3652 ;; Negative integer a and b is in a valid range. Expand.
3653 ($rectform
3657 ;; Argument simplifications and transformations
3663 ;; Expand for b a positive integer and a not a negative integer.
3677 ;; Expand for a a positive integer.
3693 ;; Expand for a a negative integer and b an integer with b <= -a.
3722 index))
3725 index)))
3751 index)))
3770 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3777 ($rectform
3806 (merror (intl:gettext "beta_incomplete: continued fractions failed for beta_incomplete(~:M, ~:M, ~:M).") a b z))
3844 result))))))
3846 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3847 ;;;
3848 ;;; Implementation of the Generalized Incomplete Beta function
3849 ;;;
3850 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3855 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3863 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3866 simp-beta-incomplete-generalized operators)
3868 ;;; beta_incomplete_generalized distributes over bags
3872 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3874 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
3875 ;;; but not on the negative real axis and for z1 or z2 real and > 1.
3876 ;;; We support a conjugate-function which test these cases.
3879 conjugate-beta-incomplete-generalized conjugate-function)
3903 a b z1 z2)))))))
3905 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3909 ;; Derivative wrt a
3921 z1)
3927 z2)
3930 ;; Derivative wrt b
3953 ;; The derivative wrt z1
3959 ;; The derivative wrt z2
3965 grad)
3967 ;;; Integral of the Incomplete Beta function
3971 nil
3972 nil
3973 ;; Integral wrt z1
3977 ;; Integral wrt z2
3982 integral)
3984 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3996 ;; Check for specific values
4004 a b z1 z2))
4017 a b z1 z2))
4042 ;; Check for numerical evaluation in Float or Bigfloat precision
4055 ;; Check for argument simplifications and transformations
4124 index))
4127 index)))
4156 index)))
4167 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4168 ;;;
4169 ;;; Implementation of the Regularized Incomplete Beta function
4170 ;;;
4171 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4176 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4184 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4187 simp-beta-incomplete-regularized operators)
4189 ;;; beta_incomplete_regularized distributes over bags
4193 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4197 ;; Derivative wrt a
4213 ;; Derivative wrt b
4230 ;; The derivative wrt z
4235 grad)
4237 ;;; Integral of the Generalized Incomplete Beta function
4241 nil
4242 nil
4243 ;; Integral wrt z
4249 integral)
4251 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4262 ;; Check for specific values
4270 a b z))
4290 ;; Problem: for b a negative integer the Regularized Incomplete
4291 ;; Beta function is defined to be zero. BUT: When we calculate
4292 ;; e.g. beta_incomplete(1.0,-2.0,1/2)/beta(1.0,-2.0) we get the
4293 ;; result -3.0, because beta_incomplete and beta are defined for
4294 ;; for this case. How do we get a consistent behaviour?
4300 ;; TODO: The following line does not work for bigfloats.
4302 ;; Does $beta_incomplete or simpbeta underflow in this case?
4308 ;; Check for numerical evaluation in Float or Bigfloat precision
4312 beta ibeta )
4316 ;; In case of underflow (see bug #2999) or precision loss use bigfloats
4317 ;; and emporarily give some extra precision but avoid fpprec dependency.
4318 ;; Is this workaround correct for complex values?
4327 ($rectform
4331 ;; Check for argument simplifications and transformations
4355 index))
4358 index)))
4378 index)))
4386 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;