| blob | bb86c3cedce6ed90d74c0a461759a00804d444a4 |
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
84 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
85 ;;; The following functions test if numerical evaluation has to be done.
86 ;;; The functions should help to test for numerical evaluation more consistent
87 ;;; and without complicated conditional tests including more than one or two
88 ;;; arguments.
89 ;;;
90 ;;; The functions take a list of arguments. All arguments have to be a CL or
91 ;;; Maxima number. If all arguments are numbers we have two cases:
92 ;;; 1. $numer is T we return T. The function has to be evaluated numerically.
93 ;;; 2. One of the args is a float or a bigfloat. Evaluate numerically.
94 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
96 ;;; Test for numerically evaluation in float precision
106 ;;; Test for numerically evaluation in complex float precision
109 "Determine if ARGS consists of numerical values by determining if
110 the real and imaginary parts of each arg are nuemrical (but not
111 bigfloats). A non-NIL result is returned if at least one of args is
112 a floating-point value or if numer is true. If the result is
113 non-NIL, it is a list of the arguments reduced via COMPLEX-NUMBER-P"
120 ;; Always save the result from complex-number-p. But for backward
121 ;; compatibility, only set the flag if any item is a float.
125 ;; Return the values in the same order as the args!
128 ;;; Test for numerically evaluation in bigfloat precision
139 ;;; Test for numerically evaluation in complex bigfloat precision
142 "Determine if ARGS consists of numerical values by determining if
143 the real and imaginary parts of each arg are nuemrical (including
144 bigfloats). A non-NIL result is returned if at least one of args is
145 a floating-point value or if numer is true. If the result is
146 non-NIL, it is a list of the arguments reduced via COMPLEX-NUMBER-P."
154 ;; Always save the result from complex-number-p. But for backward
155 ;; compatibility, only set the flag if any item is a bfloat.
160 ;; Return the values in the same order as the args!
163 ;;; Test for numerical evaluation in any precision, real or complex.
170 ;;; Check for an integer or a float or bigfloat representation. When we
171 ;;; have a float or bigfloat representation return the integer value.
180 val)
183 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
185 ;;; The changes to the parser to connect the operator !! to double_factorial(z)
187 ;(def-mheader |$!!| (%double_factorial))
189 ;(def-led (|$!!| 160.) (op left)
190 ; (list '$expr
191 ; (mheader '$!!)
192 ; (convert left '$expr)))
194 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
195 ;;;
196 ;;; The implementation of the function Double factorial
197 ;;;
198 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
200 ;;; Double factorial distributes over bags
204 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
206 ;;; Double factorial has mirror symmetry
210 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
212 ;;; Differentiation of Double factorial
228 grad)
230 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
232 ;;; Double factorial is a simplifying function
237 ;; Positive Integer less then $factlim or $factlim is -1. Call gfact.
243 ;; Even negative integer or real representation. Not defined.
254 ;; Odd negative integer, float or complex float.
261 ;; bigfloat or complex bigfloat.
265 ;; double_factorial(inf) -> inf
275 ;; Special case double_factorial(n-1)
276 ;; Not sure if this simplification is useful.
280 ;; Special case double_factorial(2*k+n), k integer
293 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
294 ;;; Double factorial for a complex float argument. The result is a CL complex.
295 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
299 (*
306 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
307 ;;; Double factorial for a bigfloat or complex bigfloat argument
308 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
325 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
326 ;;;
327 ;;; The implementation of the Incomplete Gamma function
328 ;;;
329 ;;; A&S 6.5.3:
330 ;;;
331 ;;; inf
332 ;;; /
333 ;;; [ a - 1 - t
334 ;;; gamma_incomplete(a, z) = I t %e dt
335 ;;; ]
336 ;;; /
337 ;;; z
338 ;;;
339 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
343 ;; TODO: This is currently called by integrate-exp-special in sin.lisp.
344 ;; Need to fix that before this can be removed.
348 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
351 ;;; Incomplete Gamma distributes over bags
355 ;;; Incomplete Gamma function has not mirror symmetry for z on the negative
356 ;;; real axis. We support a conjugate-function which test this case.
363 ;; Definitely not on the negative real axis for z. Mirror symmetry.
370 ;; On the negative real axis or no information. Unsimplified.
375 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
377 ;;; Derivative of the Incomplete Gamma function
383 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2
384 ;; function and the Generalized Incomplete Gamma function
385 ;; (functions.wolfram.com), only for a>0.
401 ;; No derivative. Maxima generates a noun form.
402 nil)))
403 ;; The derivative wrt z
409 ;;; Integral of the Incomplete Gamma function
413 nil
417 integral)
419 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
421 ;;; We support a simplim%function. The function is looked up in simplimit and
422 ;;; handles specific values of the function.
427 ;; Look for the limit of the arguments.
438 ;; Handle an argument 0 at this place.
445 ;; Call the simplifier of the function.
448 ;; All other cases are handled by the simplifier of the function.
451 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
453 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
454 ;;;
455 ;;; Implementation of the Lower Incomplete gamma function:
456 ;;;
457 ;;; A&S 6.5.2
458 ;;;
459 ;;; z
460 ;;; /
461 ;;; [ a - 1 - t
462 ;;; gamma_incomplete_lower(a, z) = I t %e dt
463 ;;; ]
464 ;;; /
465 ;;; 0
466 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
467 ;;; distribute over bags (aggregates)
471 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
473 ;; (defprop %gamma_incomplete_lower ??? grad) WHAT TO PUT HERE ??
475 ;;
476 ;; Handles some special cases for the order a and simplifies it to an
477 ;; equivalent form, possibly involving erf and gamma_incomplete_lower
478 ;; to a lower order.
491 ;; Try to express gamma_incomplete_lower(a,z) in simpler terms for
492 ;; special values of the order a. If we can't, return NIL to say so.
495 ;; gamma_incomplete_lower(0, x) = 1-exp(x)
498 ;; gamma_incomplete_lower(a,z) can be simplified if a is a
499 ;; positive integer. Since gamma_incomplete_lower(a,z) =
500 ;; gamma(a) - gamma_incomplete(a,z), use gamma_incomplete to
501 ;; get the result.
505 ;; A&S 6.5.12:
506 ;;
507 ;; gamma_incomplete_lower(1/2,x^2) = sqrt(%pi)*erf(x)
508 ;;
509 ;; gamma_incomplete_lower(1/2,z) = sqrt(%pi)*erf(sqrt(x))
510 ;;
514 ;; gamma_incomplete_lower(a+n,z), where n is an integer
519 ;; gamma_incomplete_lower(a+n,z). See DLMF 8.8.7:
520 ;;
521 ;; gamma_incomplete_lower(a+n,z)
522 ;; = pochhammer(a,n)*gamma_incomplete_lower(a,z)
523 ;; -z^a*exp(-z)*sum(gamma(a+n)gamma(a+k+1)*z^k,k,0,n-1)
542 ;; See DLMF 8.8.8. For simplicity let g(a,z) = gamma_incomplete_lower(a,z).
543 ;;
544 ;; g(a,z) = gamma(a)/gamma(a-n)*g(a-n,z)
545 ;; - z^(a-1)*exp(z)*sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
546 ;;
547 ;; Rewrite:
548 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
549 ;; + gamma(a-n)/gamma(a)*z^(a-1)*exp(-z)
550 ;; * sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
551 ;; Or
552 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
553 ;; + z^(a-1)*exp(-z)
554 ;; * sum(gamma(a-n)/gamma(a-k)*z^(-k),k,0,n-1)
575 ;; gamma_incomplete_lower of (numeric) rational order.
576 ;; Expand it out so that the resulting order is between 0 and
577 ;; 1.
580 ;; a = n + order where 0 <= order < 1.
584 ;; Nothing to do if the order is already between 0 and 1
587 ;; Use gamma_incomplete(a+n,z) above. and then substitute
588 ;; a=order. This works for n positive or negative.
593 ;; No expansion so return nil to indicate that
594 nil)))
596 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
598 ;;; Incomplete Gamma function is a simplifying function
605 ;; Check for specific values
608 ;; gamma_incomplete(v,0) is gamma(v) only if the realpart(v) >
609 ;; 0. If realpart(v) <= 0, gamma_incomplete is undefined. For
610 ;; all other cases, return the noun form.
616 a z))
625 ;; Check for numerical evaluation in Float or Bigfloat precision
630 ;; a and z are Maxima numbers, at least one has a float value
637 ;; a is zero or a negative float representing an integer.
638 ;; For these cases the numerical routines of gamma-incomplete
639 ;; do not work. Call the numerical routine for the Exponential
640 ;; Integral E(n,z). The routine is called with a positive integer!.
650 ;; a and z are Maxima numbers, at least one is a complex value and
651 ;; we have at least one float part
659 ;; Call expintegral-e. See comment above.
674 ;; Call bfloat-expintegral-e. See comment above.
695 ;; Call bfloat-expintegral-e. See comment above.
705 ;; Check for transformations and argument simplification
708 ;; Integer or a symbol declared to be an integer. Expand in a series.
730 ;; Simplify gamma, but do not expand to avoid division
731 ;; by zero.
759 ;; We have a half integral order and $gamma_expand is not NIL.
760 ;; We expand in a series with the Erfc function
810 ;; Handle gamma_incomplete(a+n, z), where n is a numerical integer
815 ;; See DLMF 8.8.9: https://dlmf.nist.gov/8.8.E9
816 ;;
817 ;; gamma_incomplete(a+n,z) = pochhammer(a,n)*gamma_incomplete(a,z)
818 ;; + z^a*exp(-z)*sum(gamma(a+n)/gamma(a+k+1)*z^k,k,0,n-1)
836 ;; See http://functions.wolfram.com/06.06.17.0004.01
837 ;;
838 ;; gamma_incomplate(a,z) =
839 ;; (-1)^n*pochhammer(1-a,n)
840 ;; *[gamma_incomplete(a-n,z)
841 ;; + z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)]
842 ;;
843 ;; Rerarrange this in terms of gamma_incomplete(a-n,z):
844 ;;
845 ;; gamma_incomplete(a-n,z) =
846 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
847 ;; -z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)
848 ;;
849 ;; Change the summation index to go from k = 0 to n-1:
850 ;;
851 ;; z^(a-n-1)*sum(z^k/pochhammer(a-n,k),k,1,n)
852 ;; = z^(a-n-1)*sum(z^(k+1)/pochhammer(a-n,k+1),k,0,n-1)
853 ;; = z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
854 ;;
855 ;; Thuus:
856 ;; gamma_incomplete(a-n,z) =
857 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
858 ;; -z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
860 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
869 ;; sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
879 ;; gamma_incomplete of (numeric) rational order. Expand it out
880 ;; so that the resulting order is between 0 and 1.
883 ;; a = n + order where 0 <= order < 1.
887 ;; Nothing to do if the order is already between 0 and 1
890 ;; Use gamma_incomplete(a+n,z) above. and then substitute
891 ;; a=order. This works for n positive or negative.
898 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
899 ;;; Numerical evaluation of the Incomplete Gamma function
900 ;;;
901 ;;; gamma-incomplete (a,z) - real and complex double float
902 ;;; bfloat-gamma-incomplete (a z) - bigfloat
903 ;;; complex-bfloat-gamma-incomplete (a z) - complex bigfloat
904 ;;;
905 ;;; Expansion in a power series for abs(x) < R, where R is
906 ;;; *gamma-radius* + real(a) if real(a) > 0 or *gamma-radius*
907 ;;; otherwise.
908 ;;;
909 ;;; (A&S 6.5.29):
910 ;;;
911 ;;; inf
912 ;;; ===
913 ;;; \ gamma(a)
914 ;;; gamma(a,z) = exp(-x)*z^a * > ------------ * z^n
915 ;;; / gamma(a+1+n)
916 ;;; ===
917 ;;; n=0
918 ;;;
919 ;;; This expansion does not work for an integer a<=0, because the Gamma function
920 ;;; in the denominator is not defined for a=0 and negative integers. For this
921 ;;; case we use the Exponential Integral E for numerically evaluation. The
922 ;;; Incomplete Gamma function and the Exponential integral are connected by
923 ;;;
924 ;;; gamma(a,z) = z^a * expintegral_e(1-a,z)
925 ;;;
926 ;;; When the series is not used, two forms of the continued fraction
927 ;;; are used. When z is not near the negative real axis use the
928 ;;; continued fractions (A&S 6.5.31):
929 ;;;
930 ;;; 1 1-a 1 2-a 2
931 ;;; gamma(a,z) = exp(-z) z^a *( -- --- --- --- --- ... )
932 ;;; z+ 1+ z+ 1+ z+
933 ;;;
934 ;;; The accuracy is controlled by *gamma-incomplete-eps* for double float
935 ;;; precision. For bigfloat precision epsilon is 10^(-$fpprec). The expansions
936 ;;; in a power series or continued fractions stops if *gamma-incomplete-maxit*
937 ;;; is exceeded and an Maxima error is thrown.
938 ;;;
939 ;;; The above fraction does not converge on the negative real axis and
940 ;;; converges very slowly near the axis. In this case, use the
941 ;;; relationship
942 ;;;
943 ;;; gamma(a,z) = gamma(a) - gamma_lower(a,z)
944 ;;;
945 ;;; The continued fraction for gamma_incomplete_lower(a,z) is
946 ;;; (http://functions.wolfram.com/06.06.10.0009.01):
947 ;;;
948 ;;; gamma_lower(a,z) = exp(-z) * z^a / cf(a,z)
949 ;;;
950 ;;; where
951 ;;;
952 ;;; -a*z z (a+1)*z 2*z (a+2)*z
953 ;;; cf(a,z) = a + ---- ---- ------- ---- -------
954 ;;; a+1+ a+2- a+3+ a+4- a+5+
955 ;;;
956 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
966 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
968 ;;; The numerical evaluation for CL float or complex values a and x
969 ;;; When the flag regularized is T, the result is divided by gamma(a) and
970 ;;; Maxima returns the numercial result for gamma_incomplete_regularized
976 ;; Compute the factor needed to scale the series or continued
977 ;; fraction. This is x^a*exp(-x) or x^a*exp(-x)/gamma(a)
978 ;; depending on whether we want a non-regularized or
979 ;; regularized form. We want to compute the factor carefully
980 ;; to avoid unnecessary overflow if possible.
984 ;; gammafloat is more accurate for real
985 ;; values of a.
992 ;; Easy way failed. Use logs to compute the
993 ;; result. This loses some precision, especially
994 ;; for large values of a and/or x because we use
995 ;; many bits to hold the exponent part.
997 x
999 a
1005 ;; Easy way failed, so use the log form.
1007 x)))))))
1011 ;; %gamma-incomplete compute the lower incomplete gamma
1012 ;; function, so we need to subtract that from gamma(a),
1013 ;; more or less.
1021 ;; Continued fraction used. Just multiply by the factor
1022 ;; to get the final result.
1025 ;; Compute the key part of the gamma incomplete function using either
1026 ;; a series expression or a continued fraction expression. Two values
1027 ;; are returned: the value itself and a boolean, indicating what the
1028 ;; computed value is. If the boolean non-NIL, then the computed value
1029 ;; is the lower incomplete gamma function.
1037 ;; The series expansion is done for x within a circle of radius
1038 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1039 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1040 ;; continued fraction is used.
1046 ;; For x near the negative real axis, use the
1047 ;; relationship gamma_incomplete(a,z) = gamma(a) -
1048 ;; gamma_incomplete_lower(a,z), where
1049 ;; gamma_incomplete_lower(a,z) is the lower poart of the
1050 ;; incomplete gamma function. We can evaluate that
1051 ;; using a continued fraction from
1052 ;; http://functions.wolfram.com/06.06.10.0009.01. (Note
1053 ;; that the alternative fraction,
1054 ;; http://functions.wolfram.com/06.06.10.0007.01,
1055 ;; appears to be less accurate.)
1056 ;;
1057 ;; Also note that this appears to be valid for all
1058 ;; values of x (real or complex), but we don't want to
1059 ;; use this everywhere for gamma_incomplete. Consider
1060 ;; what happens for large real x. gamma_incomplete(a,x)
1061 ;; is small, but gamma_incomplete(a,x) = gamma(x) - cf
1062 ;; will have large roundoff errors.
1075 t)))
1077 ;; Expansion in continued fractions for gamma_incomplete.
1088 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1097 ;; Return nil to indicate we used the continued fraction.
1100 ;; Expansion in a series
1108 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1111 ;; Return T to indicate we used the series and the series
1112 ;; is for the integral from 0 to x, not x to inf.
1115 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1117 ;;; This function is called for a and x real
1125 ;; The series expansion is done for x within a circle of radius
1126 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1127 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1128 ;; continued fraction is used.
1135 ;; Expansion in continued fractions of the Incomplete Gamma function
1144 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1165 ;; Expand to multiply everything out.
1166 ($expand
1167 ;; Expansion in continued fraction for the lower incomplete gamma.
1169 ;; NOTE: We want (power x a) instead of bigfloat:expt
1170 ;; because this preserves how maxima computes x^a when
1171 ;; x is negative and a is rational. For, example
1172 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1173 ;; principal value.
1186 x))))))))))))))
1189 ;; Series expansion of the Incomplete Gamma function
1195 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1209 ($rectform
1214 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1226 ;; The series expansion is done for x within a circle of radius
1227 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1228 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1229 ;; continued fraction is used.
1241 ;; Expansion in continued fractions of the Incomplete Gamma function
1252 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1263 gm-eps))
1272 ;; Expand to multiply everything out.
1273 ($expand
1274 ;; Expansion in continued fraction for the lower incomplete gamma.
1276 ;; NOTE: We want (power x a) instead of bigfloat:expt
1277 ;; because this preserves how maxima computes x^a when
1278 ;; x is negative and a is rational. For, example
1279 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1280 ;; principal value.
1293 x))))))))))))))
1295 ;; Series expansion of the Incomplete Gamma function
1303 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1316 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1317 ;;;
1318 ;;; Implementation of the Generalized Incomplete Gamma function
1319 ;;;
1320 ;;; z2
1321 ;;; /
1322 ;;; [ a - 1 - t
1323 ;;; gamma_incomplete_generalized(a, z1, z2) = I t %e dt
1324 ;;; ]
1325 ;;; /
1326 ;;; z1
1327 ;;;
1328 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1330 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1331 ;;; on the negative real axis.
1332 ;;; We support a conjugate-function which test this case.
1335 conjugate-gamma-incomplete-generalized conjugate-function)
1340 ;; z1 and z2 definitely not on the negative real axis.
1341 ;; Mirror symmetry.
1349 ;; On the negative real axis or no information. Unsimplified.
1354 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1356 ;;; Generalized Incomplete Gamma distributes over bags
1360 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1362 ;;; Differentiation of Generalized Incomplete Gamma function
1366 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1367 ;; and the Generalized Incomplete Gamma function (functions.wolfram.com)
1389 ;; The derivative wrt z1
1393 ;; The derivative wrt z2
1397 grad)
1399 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1401 ;;; Generalized Incomplete Gamma function is a simplifying function
1407 ;; Check for specific values
1434 ;; Check for numerical evaluation in Float or Bigfloat precision
1435 ;; Use the numerical routines of the Incomplete Gamma function
1462 ;; Check for transformations and argument simplification
1465 ;; Expand gamma_incomplete_generalized(a+n,z1,z2) with n an integer
1518 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1519 ;;;
1520 ;;; Implementation of the Regularized Incomplete Gamma function
1521 ;;;
1522 ;;; A&S 6.5.1
1523 ;;;
1524 ;;; gamma_incomplete(a, z)
1525 ;;; gamma_incomplete_regularized(a, z) = ----------------------
1526 ;;; gamma(a)
1527 ;;;
1528 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1530 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1532 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1533 ;;; on the negative real axis.
1534 ;;; We support a conjugate-function which test this case.
1537 conjugate-gamma-incomplete-regularized conjugate-function)
1542 ;; z definitely not on the negative real axis. Mirror symmetry.
1549 ;; On the negative real axis or no information. Unsimplified.
1554 ;;; Regularized Incomplete Gamma distributes over bags
1558 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1560 ;;; Differentiation of Regularized Incomplete Gamma function
1564 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1565 ;; and the Regularized Generalized Incomplete Gamma function
1566 ;; (functions.wolfram.com)
1580 ;; The derivative wrt z
1582 -1
1586 grad)
1588 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1590 ;;; Regularized Incomplete Gamma function is a simplifying function
1598 ;; Check for specific values
1606 a z))
1613 ;; Check for numerical evaluation in Float or Bigfloat precision
1617 ;; gamma_incomplete returns a regularized result
1623 ;; gamma_incomplete returns a regularized result
1635 ($rectform
1640 ;; Check for transformations and argument simplification
1643 ;; An integer. Expand the expression.
1660 ;; We have a half integral order and $gamma_expand is not NIL.
1661 ;; We expand in a series with the Erfc function
1714 ;; We factor the second summand.
1715 ;; Some factors vanish and the result is more readable.
1716 ($factor
1731 ;; We factor the second summand.
1732 ($factor
1746 ;; gamma_incomplete_regularized of numeric rational order.
1747 ;; Expand it out so that the resulting order is between 0 and
1748 ;; 1. Use gamma_incomplete_regularized(a+n,z) to do the dirty
1749 ;; work.
1752 ;; a = n + order where 0 <= order < 1.
1756 ;; Nothing to do if the order is already between 0 and 1
1759 ;; Use gamma_incomplete_regularized(a+n,z) above. and
1760 ;; then substitute a=order. This works for n positive or
1761 ;; negative.
1768 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1769 ;;;
1770 ;;; Implementation of the Logarithm of the Gamma function
1771 ;;;
1772 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1777 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1785 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1789 ;;; Logarithm of the Gamma function distributes over bags
1793 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1798 grad)
1800 ;; integrate(log_gamma(x),x) = psi[-2](x)
1806 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1813 ;; Check for specific values
1819 ;; We have zero, a negative integer or a float or bigfloat representation.
1825 ;; Check for numerical evaluation
1842 ;; Transform to Logarithm of Factorial for integer values
1843 ;; At this point the integer value is positive and not zero.
1850 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1851 ;;; The functions log-gamma-lanczos, bfloat-log-gamma and
1852 ;;; complex-bfloat-log-gamma are modified versions of the related functions
1853 ;;; gamma-lanczos, bffac and cbffac. The functions return the Logarithm of
1854 ;;; the Gamma function. If we have to calculate the quotient of Gamma functions,
1855 ;;; e. g. for the Beta function, it is much more appropriate to use the
1856 ;;; logarithmic versions to avoid overflow.
1857 ;;;
1858 ;;; Be careful log(gamma(z)) is only for realpart(z) positive equal to
1859 ;;; log_gamma(z). For a negative realpart(z) log_gamma differ by multiple of
1860 ;;; %pi from log(gamma(z)). But we always have exp(log_gamma(z))= gamma(z).
1861 ;;; The terms to get the transformation for log_gamma(-z) are taken from
1862 ;;; functions.wolfram.com.
1863 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1869 :initial-contents
1871 57.156235665862923517
1872 -59.597960355475491248
1873 14.136097974741747174
1874 -0.49191381609762019978
1875 .33994649984811888699e-4
1876 .46523628927048575665e-4
1877 -.98374475304879564677e-4
1878 .15808870322491248884e-3
1879 -.21026444172410488319e-3
1880 .21743961811521264320e-3
1881 -.16431810653676389022e-3
1882 .84418223983852743293e-4
1883 -.26190838401581408670e-4
1889 (+
1890 (*
1898 (*
1912 sum)
1918 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1939 bigfloat%pi
1942 bigfloat%pi '$%i
1960 y)))
1966 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1987 bigfloat%pi
1990 bigfloat%pi '$%i
2009 y))))
2010 ($rectform
2016 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2017 ;;;
2018 ;;; Implementation of the Error function Erf(z)
2019 ;;;
2020 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2022 ;; TODO: This is called by simp-expintegral-e in expintegral.lisp.
2023 ;; Can't remove this until that is fixed.
2027 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2029 ;;; erf has mirror symmetry
2033 ;;; erf is an odd function
2037 ;;; erf distributes over bags
2041 ;;; Derivative of the Error function erf
2048 grad)
2050 ;;; Integral of the Error function erf
2059 integral)
2061 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2064 ; See A&S 7.1.21 or http://functions.wolfram.com/06.25.26.0001.01
2066 z
2073 ;;; erf is a simplifying function
2078 ;; Check for specific values
2084 ;; Check for numerical evaluation
2097 ;; Argument simplification
2106 ;; Representation through more general functions
2108 ($hypergeometric_representation
2111 ;; Transformation to Erfc or Erfi
2116 (%erfc
2118 (%erfi
2126 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2129 ;; We use the slatec routine for float values.
2132 ;; Compute erf(z) using the relationship
2133 ;;
2134 ;; erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2,z^2)/sqrt(%pi))
2135 ;;
2136 ;; When z is real sqrt(z^2)/z is signum(z). For complex z,
2137 ;; sqrt(z^2)/z = 1 if -%pi/2 < arg(z) <= %pi/2 and -1 otherwise.
2138 ;;
2139 ;; This relationship has serious round-off issues when z is small
2140 ;; because gamma_incomplete(1/2,z^2)/sqrt(%pi) is near 1.
2141 ;;
2142 ;; complex-erf is for (lisp) complex numbers; bfloat-erf is for real
2143 ;; bfloats, and complex-bfloat-erf is for complex bfloats. Care is
2144 ;; taken to return real results for real arguments and imaginary
2145 ;; results for imaginary arguments
2148 (*
2150 -1
2156 ;; Pure real argument, the result is real
2159 ;; Pure imaginary argument, the result is pure imaginary
2162 result))))
2165 ;; Warning! This has round-off problems when abs(z) is very small.
2167 ;; The argument is real, the result is real too
2168 ($realpart
2177 ;; Warning! This has round-off problems when abs(z) is very small.
2191 ;; Pure real argument, the result is real
2194 ;; Pure imaginary argument, the result is pure imaginary
2197 ;; A general complex result
2198 result))))
2202 ;; Erf(z) for all z. Z must be a CL real or complex number or a
2203 ;; BIGFLOAT or COMPLEX-BIGFLOAT object. The result will be of the
2204 ;; same type as Z.
2207 ;; Use Slatec derf, which should be faster than the series.
2210 ;; Use the series A&S 7.1.5 for small x:
2211 ;;
2212 ;; erf(z) = 2*z/sqrt(%pi) * sum((-1)^n*z^(2*n)/n!/(2*n+1), n, 0, inf)
2213 ;;
2214 ;; The threshold is approximately erf(x) = 0.5. (Doesn't
2215 ;; have to be super accurate.) This gives max accuracy when
2216 ;; using the identity erf(x) = 1 - erfc(x).
2222 k
2226 ;; The general case.
2233 (bigfloat (maxima::$bfloat (maxima::$expand (maxima::complex-bfloat-erf (maxima::to z))))))))))
2236 ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
2237 ;; near 1. Wolfram says
2238 ;;
2239 ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
2240 ;;
2241 ;; For real(z) > 0, sqrt(z^2)/z is 1 so
2242 ;;
2243 ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2244 ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
2245 ;;
2246 ;; For real(z) < 0, sqrt(z^2)/z is -1 so
2247 ;;
2248 ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2249 ;; = 2 - 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
2271 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2272 ;;;
2273 ;;; Implementation of the Generalized Error function Erf(z1,z2)
2274 ;;;
2275 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2277 ;;; Generalized Erf has mirror symmetry
2281 ;;; Generalized Erf distributes over bags
2285 ;;; Generalized Erf is antisymmetric Erf(z1,z2) = - Erf(z2,z1)
2291 ;; Let's remove built-in symbols from list for user-defined properties.
2294 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2298 ;; derivative wrt z1
2302 ;; derviative wrt z2
2306 grad)
2308 ;;; ----------------------------------------------------------------------------
2313 ;; Look for the limit of the arguments.
2317 ;; Handle infinities at this place.
2331 ;; All other cases are handled by the simplifier of the function.
2334 ;;; ----------------------------------------------------------------------------
2339 ;; Check for specific values
2357 ;; Check for numerical evaluation. Use erf(z1,z2) = erf(z2)-erf(z1)
2380 ;; Argument simplification
2385 ;; Representation through more general functions
2387 ($hypergeometric_representation
2390 ;; Transformation to Erf
2392 ($erf_representation
2398 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2399 ;;;
2400 ;;; Implementation of the Complementary Error function Erfc(z)
2401 ;;;
2402 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2404 ;; TODO: called by simp-expintegral-e in expintegral.lisp. Need to
2405 ;; keep this around until that is fixed.
2409 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2413 ;;; Complementary Error function distributes over bags
2417 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2424 grad)
2426 ;;; Integral of the Error function erfc
2435 integral)
2437 ;;; ----------------------------------------------------------------------------
2442 ;; Look for the limit of the arguments.
2445 ;; Handle infinities at this place.
2462 ;; All other cases are handled by the simplifier of the function.
2465 ;;; ----------------------------------------------------------------------------
2470 ;; Check for specific values
2476 ;; Check for numerical evaluation.
2481 ;; Argument simplification
2487 ;; Representation through more general functions
2489 ($hypergeometric_representation
2492 ;; Transformation to Erf or Erfi
2497 (%erf
2499 (%erfi
2507 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2508 ;;;
2509 ;;; Implementation of the Imaginary Error function Erfi(z)
2510 ;;;
2511 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2513 ;; TODO: Called by integrate-exp-special in sin.lisp. That needs to
2514 ;; be fixed before this can be removed.
2518 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2520 ;;; erfi has mirror symmetry
2524 ;;; erfi is an odd function
2528 ;;; erfi distributes over bags
2532 ;;; Derivative of the Error function erfi
2539 grad)
2541 ;;; Integral of the Error function erfi
2550 integral)
2552 ;;; ----------------------------------------------------------------------------
2557 ;; Look for the limit of the arguments.
2560 ;; Handle infinities at this place.
2564 ;; All other cases are handled by the simplifier of the function.
2567 ;;; ----------------------------------------------------------------------------
2572 ;; Use the relationship erfi(z) = -%i*erf(%i*z)
2576 ;; Take care to return real results when the argument is real.
2585 ;;; erfi is an simplifying function
2590 ;; Check for specific values
2596 ;; Check for numerical evaluation. Use erfi(z) = -%i*erf(%i*z).
2601 ;; Argument simplification
2609 ;; Representation through more general functions
2611 ($hypergeometric_representation
2614 ;; Transformation to Erf or Erfc
2619 (%erf
2621 (%erfc
2629 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2630 ;;;
2631 ;;; The implementation of the Inverse Error function
2632 ;;;
2633 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2635 ;;; The Inverse Error function distributes over bags
2639 ;;; inverse_erf is the inverse of the erf function
2644 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2646 ;;; Differentiation of the Inverse Error function
2654 grad)
2656 ;;; Integral of the Inverse Error function
2663 integral)
2665 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2667 ;;; We support a simplim%function. The function is looked up in simplimit and
2668 ;;; handles specific values of the function.
2673 ;; Look for the limit of the argument.
2676 ;; Handle an argument 1 at this place
2680 ;; All other cases are handled by the simplifier of the function.
2683 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2685 ;;; The Inverse Error function is a simplifying function
2700 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2701 ;;;
2702 ;;; The implementation of the Inverse Complementary Error function
2703 ;;;
2704 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2706 ;;; Inverse Complementary Error function distributes over bags
2710 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2712 ;;; inverse_erfc is the inverse of the erfc function
2718 ;;; Differentiation of the Inverse Complementary Error function
2726 grad)
2728 ;;; Integral of the Inverse Complementary Error function
2736 integral)
2738 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2740 ;;; We support a simplim%function. The function is looked up in simplimit and
2741 ;;; handles specific values of the function.
2746 ;; Look for the limit of the argument.
2749 ;; Handle an argument 0 at this place
2756 ;; All other cases are handled by the simplifier of the function.
2759 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2761 ;;; Inverse Complementary Error function is a simplifying function
2775 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2777 ;;; Implementation of the Newton algorithm to calculate numerical values
2778 ;;; of the Inverse Error functions in float or bigfloat precision.
2779 ;;; The algorithm is a modified version of the routine in newton1.mac.
2816 ;; Compute inverse_erf(z) for z a real or complex number, including
2817 ;; bigfloat objects. The value is computing using a Newton iteration
2818 ;; to solve erf(x) = z.
2821 z)
2825 z))
2827 ;; inverse_erf is odd because erf is.
2832 ;; Find an approximate solution for x = inverse_erf(z).
2835 ;; For small z, inverse_erf(z) = z*sqrt(%pi)/2
2836 ;; + O(z^3). Thus, x = z*sqrt(%pi)/2 is our
2837 ;; initial starting point.
2840 ;; For |z| > 1 and realpart(z) >= 0, we have
2841 ;; the asymptotic series z = erf(x) = 1 -
2842 ;; exp(-x^2)/x/sqrt(%pi).
2843 ;;
2844 ;; Then
2845 ;; x = sqrt(-log(x*sqrt(%pi)*(1-z))
2846 ;;
2847 ;; We can use this as a fixed-point iteration
2848 ;; to find x, and we can start the iteration at
2849 ;; x = 1. Just do one more iteration. I (RLT)
2850 ;; think that's close enough to get the Newton
2851 ;; algorithm to converge.
2858 result)))
2861 ;; Try to pick a reasonable epsilon value for the
2862 ;; Newton iteration.
2873 ;; Derivative of erf(x)
2876 #'diff
2877 z
2879 eps)))))))
2885 z))
2891 ;; Find an approximate solution for x =
2892 ;; inverse_erfc(z). We have inverse_erfc(z) =
2893 ;; inverse_erf(1-z), so that's a good starting point.
2894 ;; We don't need full precision, so a float value is
2895 ;; good enough. But if 1-z is 1, inverse_erf is
2896 ;; undefined, so we need to do something else.
2907 result)))
2913 ;; Derivative of erfc(x)
2916 #'diff
2917 z
2919 eps)))))))
2921 ;; Newton iteration for solving f(x) = z, given f and the derivative
2922 ;; of f.
2934 (intl:gettext "bf-newton: failed to converge after ~M iterations: delta = ~S, x = ~S eps = ~S")
2935 count delta x eps)
2936 x))
2944 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2945 ;;;
2946 ;;; Implementation of the Fresnel Integral S(z)
2947 ;;;
2948 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2950 ;;; fresnel_s distributes over bags
2954 ;;; fresnel_s has mirror symmetry
2958 ;;; fresnel_s is an odd function
2959 ;;;
2960 ;;; Maxima has two mechanism to define a reflection property
2961 ;;; 1. Declare the feature oddfun or evenfun
2962 ;;; 2. Put a reflection rule on the property list
2963 ;;;
2964 ;;; The second way is used for the trig functions. We use it here too.
2966 ;;; This would be the first way to give the property of an odd function.
2967 ;(eval-when
2968 ; #+gcl (load eval)
2969 ; #-gcl (:load-toplevel :execute)
2970 ; (let (($context '$global) (context '$global))
2971 ; (meval '(($declare) %fresnel_s $oddfun))
2972 ; ;; Let's remove built-in symbols from list for user-defined properties.
2973 ; (setq $props (remove '%fresnel_s $props))))
2977 ;;; Differentiation of the Fresnel Integral S
2982 grad)
2984 ;;; Integration of the Fresnel Integral S
2993 integral)
2995 ;;; Limits of the Fresnel Integral S
3007 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3013 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3014 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3023 ;; For very small x, we have
3024 ;; fresnel_s(x) = %pi/6*z^3
3025 ;; fresnel_c(x) = x
3043 (maxima::merror (intl:gettext "fresnel: series expansion failed for (COMPLEX-BFLOAT-FRESNEL ~:M).") z))
3062 (* sqrt-pi
3063 z)))
3065 (* sqrt-pi
3066 z)))
3077 ;; A pure imaginary argument. Use fresnel_s(%i*x)=-%i*fresnel_s(x).
3083 fs))))
3087 ;; A pure imaginary argument. Use fresnel_c(%i*x)=%i*fresnel_c(x).
3093 fc))))
3097 ;;; fresnel_s is a simplifying function
3101 ;; Check for specific values
3107 ;; Check for numerical evaluation
3111 ;; Check for argument simplification
3118 ;; Representation through equivalent functions
3120 ($erf_representation
3135 ($hypergeometric_representation
3144 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3145 ;;;
3146 ;;; Implementation of the Fresnel Integral C(z)
3147 ;;;
3148 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3150 ;;; fresnel_c distributes over bags
3154 ;;; fresnel_c has mirror symmetry
3158 ;;; fresnel_c is an odd function
3159 ;;; Maxima has two mechanism to define a reflection property
3160 ;;; 1. Declare the feature oddfun or evenfun
3161 ;;; 2. Put a reflection rule on the property list
3162 ;;;
3163 ;;; The second way is used for the trig functions. We use it here too.
3167 ;;; Differentiation of the Fresnel Integral C
3172 grad)
3174 ;;; Integration of the Fresnel Integral C
3183 integral)
3185 ;;; Limits of the Fresnel Integral C
3197 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3199 ;;; fresnel_c is a simplifying function
3203 ;; Check for specific values
3209 ;; Check for numerical evaluation
3214 ;; Check for argument simplification
3221 ;; Representation through equivalent functions
3223 ($erf_representation
3238 ($hypergeometric_representation
3247 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3248 ;;;
3249 ;;; Implementation of the Beta function
3250 ;;;
3251 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3253 ;;; The code for the implementation of the beta function is in the files
3254 ;;; csimp2.lisp, simp.lisp and mactex.lisp.
3255 ;;; At this place we only implement the operator property SYMMETRIC.
3257 ;;; Beta is symmetric beta(a,b) = beta(b,a)
3263 ;; Let's remove built-in symbols from list for user-defined properties.
3266 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3267 ;;;
3268 ;;; Implementation of the Incomplete Beta function
3269 ;;;
3270 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3275 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3277 ;;; beta_incomplete distributes over bags
3281 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3285 ;; Derivative wrt a
3293 z)
3295 ;; Derivative wrt b
3312 ;; The derivative wrt z
3316 grad)
3318 ;;; Integral of the Incomplete Beta function
3322 nil
3323 nil
3327 integral)
3329 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3340 ;; Check for specific values
3348 a b z))
3350 z)
3355 ;; z=1 and realpart(b)>0. Simplify to a Beta function.
3356 ;; If we have to evaluate numerically preserve the type.
3371 ;; The argument a is zero or a is negative and the argument b is
3372 ;; not in a valid range. Beta incomplete is undefined.
3373 ;; It would be more correct to return Complex infinity.
3378 ;; Check for numerical evaluation in Float or Bigfloat precision
3386 ;; Negative integer a and b is in a valid range. Expand.
3387 ($rectform
3398 ;; Negative integer a and b is in a valid range. Expand.
3399 ($rectform
3403 ;; Argument simplifications and transformations
3409 ;; Expand for b a positive integer and a not a negative integer.
3423 ;; Expand for a a positive integer.
3439 ;; Expand for a a negative integer and b an integer with b <= -a.
3468 index))
3471 index)))
3497 index)))
3516 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3523 ($rectform
3552 (merror (intl:gettext "beta_incomplete: continued fractions failed for beta_incomplete(~:M, ~:M, ~:M).") a b z))
3590 result))))))
3592 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3593 ;;;
3594 ;;; Implementation of the Generalized Incomplete Beta function
3595 ;;;
3596 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3598 ;;; beta_incomplete_generalized distributes over bags
3602 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3604 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
3605 ;;; but not on the negative real axis and for z1 or z2 real and > 1.
3606 ;;; We support a conjugate-function which test these cases.
3609 conjugate-beta-incomplete-generalized conjugate-function)
3633 a b z1 z2)))))))
3635 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3639 ;; Derivative wrt a
3651 z1)
3657 z2)
3660 ;; Derivative wrt b
3683 ;; The derivative wrt z1
3689 ;; The derivative wrt z2
3695 grad)
3697 ;;; Integral of the Incomplete Beta function
3701 nil
3702 nil
3703 ;; Integral wrt z1
3707 ;; Integral wrt z2
3712 integral)
3714 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3720 ;; Check for specific values
3728 a b z1 z2))
3740 a b z1 z2))
3762 ;; Check for numerical evaluation in Float or Bigfloat precision
3775 ;; Check for argument simplifications and transformations
3844 index))
3847 index)))
3876 index)))
3887 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3888 ;;;
3889 ;;; Implementation of the Regularized Incomplete Beta function
3890 ;;;
3891 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3893 ;;; beta_incomplete_regularized distributes over bags
3897 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3901 ;; Derivative wrt a
3917 ;; Derivative wrt b
3934 ;; The derivative wrt z
3939 grad)
3941 ;;; Integral of the Generalized Incomplete Beta function
3945 nil
3946 nil
3947 ;; Integral wrt z
3953 integral)
3955 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3961 ;; Check for specific values
3969 a b z))
3987 ;; Problem: for b a negative integer the Regularized Incomplete
3988 ;; Beta function is defined to be zero. BUT: When we calculate
3989 ;; e.g. beta_incomplete(1.0,-2.0,1/2)/beta(1.0,-2.0) we get the
3990 ;; result -3.0, because beta_incomplete and beta are defined for
3991 ;; for this case. How do we get a consistent behaviour?
3997 ;; TODO: The following line does not work for bigfloats.
3999 ;; Does $beta_incomplete or simpbeta underflow in this case?
4005 ;; Check for numerical evaluation in Float or Bigfloat precision
4009 beta ibeta )
4013 ;; In case of underflow (see bug #2999) or precision loss use bigfloats
4014 ;; and emporarily give some extra precision but avoid fpprec dependency.
4015 ;; Is this workaround correct for complex values?
4024 ($rectform
4028 ;; Check for argument simplifications and transformations
4052 index))
4055 index)))
4075 index)))
4083 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;