| blob | 74a74db63882dfe25ed7fe8a7114856d6a591ad7 |
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 rectform"
129 ;; Always save the result from trisplit. But for backward
130 ;; compatibility, only set the flag if any item is a float.
134 ;; Return the values in the same order as the args!
137 ;;; Test for numerically evaluation in bigfloat precision
148 ;;; Test for numerically evaluation in complex bigfloat precision
151 "Determine if ARGS consists of numerical values by determining if
152 the real and imaginary parts of each arg are nuemrical (including
153 bigfloats). A non-NIL result is returned if at least one of args is
154 a floating-point value or if numer is true. If the result is
155 non-NIL, it is a list of the arguments reduced via rectform."
164 ;; Always save the result from trisplit. But for backward
165 ;; compatibility, only set the flag if any item is a bfloat.
170 ;; Return the values in the same order as the args!
173 ;;; Test for numerical evaluation in any precision, real or complex.
180 ;;; Check for an integer or a float or bigfloat representation. When we
181 ;;; have a float or bigfloat representation return the integer value.
190 val)
193 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
195 ;;; The changes to the parser to connect the operator !! to double_factorial(z)
197 ;(def-mheader |$!!| (%double_factorial))
199 ;(def-led (|$!!| 160.) (op left)
200 ; (list '$expr
201 ; (mheader '$!!)
202 ; (convert left '$expr)))
204 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
205 ;;;
206 ;;; The implementation of the function Double factorial
207 ;;;
208 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
213 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
215 ;;; Set properties to give full support to the parser and display
223 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
225 ;;; Double factorial is a simplifying function
229 ;;; Double factorial distributes over bags
233 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
235 ;;; Double factorial has mirror symmetry
239 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
241 ;;; Differentiation of Double factorial
257 grad)
259 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
266 ;; Positive Integer less then $factlim or $factlim is -1. Call gfact.
272 ;; Even negative integer or real representation. Not defined.
283 ;; Odd negative integer, float or complex float.
290 ;; bigfloat or complex bigfloat.
294 ;; double_factorial(inf) -> inf
304 ;; Special case double_factorial(n-1)
305 ;; Not sure if this simplification is useful.
309 ;; Special case double_factorial(2*k+n), k integer
322 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
323 ;;; Double factorial for a complex float argument. The result is a CL complex.
324 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
328 (*
335 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
336 ;;; Double factorial for a bigfloat or complex bigfloat argument
337 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
354 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
355 ;;;
356 ;;; The implementation of the Incomplete Gamma function
357 ;;;
358 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
365 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
367 ;;; Set properties to give full support to the parser and display
375 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
377 ;;; Incomplete Gamma function is a simplifying function
381 ;;; Incomplete Gamma distributes over bags
385 ;;; Incomplete Gamma function has not mirror symmetry for z on the negative
386 ;;; real axis. We support a conjugate-function which test this case.
393 ;; Definitely not on the negative real axis for z. Mirror symmetry.
400 ;; On the negative real axis or no information. Unsimplified.
405 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
407 ;;; Derivative of the Incomplete Gamma function
413 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2
414 ;; function and the Generalized Incomplete Gamma function
415 ;; (functions.wolfram.com), only for a>0.
431 ;; No derivative. Maxima generates a noun form.
432 nil)))
433 ;; The derivative wrt z
439 ;;; Integral of the Incomplete Gamma function
443 nil
447 integral)
449 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
451 ;;; We support a simplim%function. The function is looked up in simplimit and
452 ;;; handles specific values of the function.
457 ;; Look for the limit of the arguments.
468 ;; Handle an argument 0 at this place.
475 ;; Call the simplifier of the function.
478 ;; All other cases are handled by the simplifier of the function.
481 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
486 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
494 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
498 ;;; distribute over bags (aggregates)
502 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
504 ;; (defprop %gamma_incomplete_lower ??? grad) WHAT TO PUT HERE ??
506 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
523 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
534 ;; Check for specific values
537 ;; gamma_incomplete(v,0) is gamma(v) only if the realpart(v) >
538 ;; 0. If realpart(v) <= 0, gamma_incomplete is undefined. For
539 ;; all other cases, return the noun form.
545 a z))
554 ;; Check for numerical evaluation in Float or Bigfloat precision
559 ;; a and z are Maxima numbers, at least one has a float value
566 ;; a is zero or a negative float representing an integer.
567 ;; For these cases the numerical routines of gamma-incomplete
568 ;; do not work. Call the numerical routine for the Exponential
569 ;; Integral E(n,z). The routine is called with a positive integer!.
579 ;; a and z are Maxima numbers, at least one is a complex value and
580 ;; we have at least one float part
588 ;; Call expintegral-e. See comment above.
603 ;; Call bfloat-expintegral-e. See comment above.
624 ;; Call bfloat-expintegral-e. See comment above.
634 ;; Check for transformations and argument simplification
637 ;; Integer or a symbol declared to be an integer. Expand in a series.
659 ;; Simplify gamma, but do not expand to avoid division
660 ;; by zero.
688 ;; We have a half integral order and $gamma_expand is not NIL.
689 ;; We expand in a series with the Erfc function
778 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
779 ;;; Numerical evaluation of the Incomplete Gamma function
780 ;;;
781 ;;; gamma-incomplete (a,z) - real and complex double float
782 ;;; bfloat-gamma-incomplete (a z) - bigfloat
783 ;;; complex-bfloat-gamma-incomplete (a z) - complex bigfloat
784 ;;;
785 ;;; Expansion in a power series for abs(x) < R, where R is
786 ;;; *gamma-radius* + real(a) if real(a) > 0 or *gamma-radius*
787 ;;; otherwise.
788 ;;;
789 ;;; (A&S 6.5.29):
790 ;;;
791 ;;; inf
792 ;;; ===
793 ;;; \ gamma(a)
794 ;;; gamma(a,z) = exp(-x)*z^a * > ------------ * z^n
795 ;;; / gamma(a+1+n)
796 ;;; ===
797 ;;; n=0
798 ;;;
799 ;;; This expansion does not work for an integer a<=0, because the Gamma function
800 ;;; in the denominator is not defined for a=0 and negative integers. For this
801 ;;; case we use the Exponential Integral E for numerically evaluation. The
802 ;;; Incomplete Gamma function and the Exponential integral are connected by
803 ;;;
804 ;;; gamma(a,z) = z^a * expintegral_e(1-a,z)
805 ;;;
806 ;;; When the series is not used, two forms of the continued fraction
807 ;;; are used. When z is not near the negative real axis use the
808 ;;; continued fractions (A&S 6.5.31):
809 ;;;
810 ;;; 1 1-a 1 2-a 2
811 ;;; gamma(a,z) = exp(-z) z^a *( -- --- --- --- --- ... )
812 ;;; z+ 1+ z+ 1+ z+
813 ;;;
814 ;;; The accuracy is controlled by *gamma-incomplete-eps* for double float
815 ;;; precision. For bigfloat precision epsilon is 10^(-$fpprec). The expansions
816 ;;; in a power series or continued fractions stops if *gamma-incomplete-maxit*
817 ;;; is exceeded and an Maxima error is thrown.
818 ;;;
819 ;;; The above fraction does not converge on the negative real axis and
820 ;;; converges very slowly near the axis. In this case, use the
821 ;;; relationship
822 ;;;
823 ;;; gamma(a,z) = gamma(a) - gamma_lower(a,z)
824 ;;;
825 ;;; The continued fraction for gamma_incomplete_lower(a,z) is
826 ;;; (http://functions.wolfram.com/06.06.10.0009.01):
827 ;;;
828 ;;; gamma_lower(a,z) = exp(-z) * z^a / cf(a,z)
829 ;;;
830 ;;; where
831 ;;;
832 ;;; -a*z z (a+1)*z 2*z (a+2)*z
833 ;;; cf(a,z) = a + ---- ---- ------- ---- -------
834 ;;; a+1+ a+2- a+3+ a+4- a+5+
835 ;;;
836 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
846 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
848 ;;; The numerical evaluation for CL float or complex values a and x
849 ;;; When the flag regularized is T, the result is divided by gamma(a) and
850 ;;; Maxima returns the numercial result for gamma_incomplete_regularized
856 ;; Compute the factor needed to scale the series or continued
857 ;; fraction. This is x^a*exp(-x) or x^a*exp(-x)/gamma(a)
858 ;; depending on whether we want a non-regularized or
859 ;; regularized form. We want to compute the factor carefully
860 ;; to avoid unnecessary overflow if possible.
864 ;; gammafloat is more accurate for real
865 ;; values of a.
872 ;; Easy way failed. Use logs to compute the
873 ;; result. This loses some precision, especially
874 ;; for large values of a and/or x because we use
875 ;; many bits to hold the exponent part.
877 x
879 a
885 ;; Easy way failed, so use the log form.
887 x)))))))
891 ;; %gamma-incomplete compute the lower incomplete gamma
892 ;; function, so we need to subtract that from gamma(a),
893 ;; more or less.
901 ;; Continued fraction used. Just multiply by the factor
902 ;; to get the final result.
905 ;; Compute the key part of the gamma incomplete function using either
906 ;; a series expression or a continued fraction expression. Two values
907 ;; are returned: the value itself and a boolean, indicating what the
908 ;; computed value is. If the boolean non-NIL, then the computed value
909 ;; is the lower incomplete gamma function.
917 ;; The series expansion is done for x within a circle of radius
918 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
919 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
920 ;; continued fraction is used.
926 ;; For x near the negative real axis, use the
927 ;; relationship gamma_incomplete(a,z) = gamma(a) -
928 ;; gamma_incomplete_lower(a,z), where
929 ;; gamma_incomplete_lower(a,z) is the lower poart of the
930 ;; incomplete gamma function. We can evaluate that
931 ;; using a continued fraction from
932 ;; http://functions.wolfram.com/06.06.10.0009.01. (Note
933 ;; that the alternative fraction,
934 ;; http://functions.wolfram.com/06.06.10.0007.01,
935 ;; appears to be less accurate.)
936 ;;
937 ;; Also note that this appears to be valid for all
938 ;; values of x (real or complex), but we don't want to
939 ;; use this everywhere for gamma_incomplete. Consider
940 ;; what happens for large real x. gamma_incomplete(a,x)
941 ;; is small, but gamma_incomplete(a,x) = gamma(x) - cf
942 ;; will have large roundoff errors.
955 t)))
957 ;; Expansion in continued fractions for gamma_incomplete.
968 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
977 ;; Return nil to indicate we used the continued fraction.
980 ;; Expansion in a series
988 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
991 ;; Return T to indicate we used the series and the series
992 ;; is for the integral from 0 to x, not x to inf.
995 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
997 ;;; This function is called for a and x real
1005 ;; The series expansion is done for x within a circle of radius
1006 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1007 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1008 ;; continued fraction is used.
1015 ;; Expansion in continued fractions of the Incomplete Gamma function
1024 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1045 ;; Expand to multiply everything out.
1046 ($expand
1047 ;; Expansion in continued fraction for the lower incomplete gamma.
1049 ;; NOTE: We want (power x a) instead of bigfloat:expt
1050 ;; because this preserves how maxima computes x^a when
1051 ;; x is negative and a is rational. For, example
1052 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1053 ;; principal value.
1066 x))))))))))))))
1069 ;; Series expansion of the Incomplete Gamma function
1075 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1089 ($rectform
1094 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1106 ;; The series expansion is done for x within a circle of radius
1107 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1108 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1109 ;; continued fraction is used.
1121 ;; Expansion in continued fractions of the Incomplete Gamma function
1132 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1143 gm-eps))
1152 ;; Expand to multiply everything out.
1153 ($expand
1154 ;; Expansion in continued fraction for the lower incomplete gamma.
1156 ;; NOTE: We want (power x a) instead of bigfloat:expt
1157 ;; because this preserves how maxima computes x^a when
1158 ;; x is negative and a is rational. For, example
1159 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1160 ;; principal value.
1173 x))))))))))))))
1175 ;; Series expansion of the Incomplete Gamma function
1183 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1196 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1197 ;;;
1198 ;;; Implementation of the Generalized Incomplete Gamma function
1199 ;;;
1200 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1205 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1207 ;;; Set the properties alias, reversealias, noun and verb
1213 $gamma_incomplete_generalized reversealias)
1215 $gamma_incomplete_generalized noun)
1217 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1218 ;;; on the negative real axis.
1219 ;;; We support a conjugate-function which test this case.
1222 conjugate-gamma-incomplete-generalized conjugate-function)
1227 ;; z1 and z2 definitely not on the negative real axis.
1228 ;; Mirror symmetry.
1236 ;; On the negative real axis or no information. Unsimplified.
1241 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1243 ;;; Generalized Incomplete Gamma function is a simplifying function
1246 simp-gamma-incomplete-generalized operators)
1248 ;;; Generalized Incomplete Gamma distributes over bags
1252 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1254 ;;; Differentiation of Generalized Incomplete Gamma function
1258 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1259 ;; and the Generalized Incomplete Gamma function (functions.wolfram.com)
1281 ;; The derivative wrt z1
1285 ;; The derivative wrt z2
1289 grad)
1291 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1303 ;; Check for specific values
1330 ;; Check for numerical evaluation in Float or Bigfloat precision
1331 ;; Use the numerical routines of the Incomplete Gamma function
1358 ;; Check for transformations and argument simplification
1361 ;; Expand gamma_incomplete_generalized(a+n,z1,z2) with n an integer
1414 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1415 ;;;
1416 ;;; Implementation of the Regularized Incomplete Gamma function
1417 ;;;
1418 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1423 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1429 $gamma_incomplete_regularized reversealias)
1431 $gamma_incomplete_regularized noun)
1433 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1435 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1436 ;;; on the negative real axis.
1437 ;;; We support a conjugate-function which test this case.
1440 conjugate-gamma-incomplete-regularized conjugate-function)
1445 ;; z definitely not on the negative real axis. Mirror symmetry.
1452 ;; On the negative real axis or no information. Unsimplified.
1457 ;;; Regularized Incomplete Gamma function is a simplifying function
1460 simp-gamma-incomplete-regularized operators)
1462 ;;; Regularized Incomplete Gamma distributes over bags
1466 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1468 ;;; Differentiation of Regularized Incomplete Gamma function
1472 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1473 ;; and the Regularized Generalized Incomplete Gamma function
1474 ;; (functions.wolfram.com)
1488 ;; The derivative wrt z
1490 -1
1494 grad)
1496 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1508 ;; Check for specific values
1516 a z))
1523 ;; Check for numerical evaluation in Float or Bigfloat precision
1527 ;; gamma_incomplete returns a regularized result
1533 ;; gamma_incomplete returns a regularized result
1545 ($rectform
1550 ;; Check for transformations and argument simplification
1553 ;; An integer. Expand the expression.
1570 ;; We have a half integral order and $gamma_expand is not NIL.
1571 ;; We expand in a series with the Erfc function
1624 ;; We factor the second summand.
1625 ;; Some factors vanish and the result is more readable.
1626 ($factor
1641 ;; We factor the second summand.
1642 ($factor
1656 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1657 ;;;
1658 ;;; Implementation of the Logarithm of the Gamma function
1659 ;;;
1660 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1665 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1673 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1677 ;;; Logarithm of the Gamma function distributes over bags
1681 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1686 grad)
1688 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1695 ;; Check for specific values
1701 ;; We have zero, a negative integer or a float or bigfloat representation.
1707 ;; Check for numerical evaluation
1724 ;; Transform to Logarithm of Factorial for integer values
1725 ;; At this point the integer value is positive and not zero.
1732 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1733 ;;; The functions log-gamma-lanczos, bfloat-log-gamma and
1734 ;;; complex-bfloat-log-gamma are modified versions of the related functions
1735 ;;; gamma-lanczos, bffac and cbffac. The functions return the Logarithm of
1736 ;;; the Gamma function. If we have to calculate the quotient of Gamma functions,
1737 ;;; e. g. for the Beta function, it is much more appropriate to use the
1738 ;;; logarithmic versions to avoid overflow.
1739 ;;;
1740 ;;; Be careful log(gamma(z)) is only for realpart(z) positive equal to
1741 ;;; log_gamma(z). For a negative realpart(z) log_gamma differ by multiple of
1742 ;;; %pi from log(gamma(z)). But we always have exp(log_gamma(z))= gamma(z).
1743 ;;; The terms to get the transformation for log_gamma(-z) are taken from
1744 ;;; functions.wolfram.com.
1745 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1751 :initial-contents
1753 57.156235665862923517
1754 -59.597960355475491248
1755 14.136097974741747174
1756 -0.49191381609762019978
1757 .33994649984811888699e-4
1758 .46523628927048575665e-4
1759 -.98374475304879564677e-4
1760 .15808870322491248884e-3
1761 -.21026444172410488319e-3
1762 .21743961811521264320e-3
1763 -.16431810653676389022e-3
1764 .84418223983852743293e-4
1765 -.26190838401581408670e-4
1771 (+
1772 (*
1780 (*
1794 sum)
1800 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1821 bigfloat%pi
1824 bigfloat%pi '$%i
1842 y)))
1848 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1869 bigfloat%pi
1872 bigfloat%pi '$%i
1891 y))))
1892 ($rectform
1898 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1899 ;;;
1900 ;;; Implementation of the Error function Erf(z)
1901 ;;;
1902 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1907 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1915 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1917 ;;; erf has mirror symmetry
1921 ;;; erf is an odd function
1925 ;;; erf is a simplifying function
1929 ;;; erf distributes over bags
1933 ;;; Derivative of the Error function erf
1940 grad)
1942 ;;; Integral of the Error function erf
1951 integral)
1953 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1960 ;; Check for specific values
1966 ;; Check for numerical evaluation
1979 ;; Argument simplification
1988 ;; Representation through equivalent functions
1990 ($hypergeometric_representation
1998 ;; Transformation to Erfc or Erfi
2003 (%erfc
2005 (%erfi
2013 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2016 ;; We use the slatec routine for float values.
2019 ;; Compute erf(z) using the relationship
2020 ;;
2021 ;; erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2,z^2)/sqrt(%pi))
2022 ;;
2023 ;; When z is real sqrt(z^2)/z is signum(z). For complex z,
2024 ;; sqrt(z^2)/z = 1 if -%pi/2 < arg(z) <= %pi/2 and -1 otherwise.
2025 ;;
2026 ;; This relationship has serious round-off issues when z is small
2027 ;; because gamma_incomplete(1/2,z^2)/sqrt(%pi) is near 1.
2028 ;;
2029 ;; complex-erf is for (lisp) complex numbers; bfloat-erf is for real
2030 ;; bfloats, and complex-bfloat-erf is for complex bfloats. Care is
2031 ;; taken to return real results for real arguments and imaginary
2032 ;; results for imaginary arguments
2035 (*
2037 -1
2043 ;; Pure real argument, the result is real
2046 ;; Pure imaginary argument, the result is pure imaginary
2049 result))))
2052 ;; Warning! This has round-off problems when abs(z) is very small.
2054 ;; The argument is real, the result is real too
2055 ($realpart
2064 ;; Warning! This has round-off problems when abs(z) is very small.
2078 ;; Pure real argument, the result is real
2081 ;; Pure imaginary argument, the result is pure imaginary
2084 ;; A general complex result
2085 result))))
2089 ;; Erf(z) for all z. Z must be a CL real or complex number or a
2090 ;; BIGFLOAT or COMPLEX-BIGFLOAT object. The result will be of the
2091 ;; same type as Z.
2094 ;; Use Slatec derf, which should be faster than the series.
2097 ;; Use the series A&S 7.1.5 for small x:
2098 ;;
2099 ;; erf(z) = 2*z/sqrt(%pi) * sum((-1)^n*z^(2*n)/n!/(2*n+1), n, 0, inf)
2100 ;;
2101 ;; The threshold is approximately erf(x) = 0.5. (Doesn't
2102 ;; have to be super accurate.) This gives max accuracy when
2103 ;; using the identity erf(x) = 1 - erfc(x).
2109 k
2113 ;; The general case.
2120 (bigfloat (maxima::$bfloat (maxima::$expand (maxima::complex-bfloat-erf (maxima::to z))))))))))
2123 ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
2124 ;; near 1. Wolfram says
2125 ;;
2126 ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
2127 ;;
2128 ;; For real(z) > 0, sqrt(z^2)/z is 1 so
2129 ;;
2130 ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2131 ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
2132 ;;
2133 ;; For real(z) < 0, sqrt(z^2)/z is -1 so
2134 ;;
2135 ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2136 ;; = 2 - 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
2158 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2159 ;;;
2160 ;;; Implementation of the Generalized Error function Erf(z1,z2)
2161 ;;;
2162 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2167 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2175 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2177 ;;; Generalized Erf has mirror symmetry
2181 ;;; Generalized Erf distributes over bags
2185 ;;; Generalized Erf is antisymmetric Erf(z1,z2) = - Erf(z2,z1)
2192 ;; Let's remove built-in symbols from list for user-defined properties.
2195 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2199 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2203 ;; derivative wrt z1
2207 ;; derviative wrt z2
2211 grad)
2213 ;;; ----------------------------------------------------------------------------
2218 ;; Look for the limit of the arguments.
2222 ;; Handle infinities at this place.
2236 ;; All other cases are handled by the simplifier of the function.
2239 ;;; ----------------------------------------------------------------------------
2248 ;; Check for specific values
2266 ;; Check for numerical evaluation. Use erf(z1,z2) = erf(z2)-erf(z1)
2289 ;; Argument simplification
2294 ;; Transformation to Erf
2296 ($erf_representation
2302 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2303 ;;;
2304 ;;; Implementation of the Complementary Error function Erfc(z)
2305 ;;;
2306 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2311 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2319 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2323 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2327 ;;; Complementary Error function distributes over bags
2331 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2338 grad)
2340 ;;; Integral of the Error function erfc
2349 integral)
2351 ;;; ----------------------------------------------------------------------------
2356 ;; Look for the limit of the arguments.
2359 ;; Handle infinities at this place.
2376 ;; All other cases are handled by the simplifier of the function.
2379 ;;; ----------------------------------------------------------------------------
2386 ;; Check for specific values
2392 ;; Check for numerical evaluation.
2397 ;; Argument simplification
2403 ;; Representation through equivalent functions
2405 ($hypergeometric_representation
2414 ;; Transformation to Erf or Erfi
2419 (%erf
2421 (%erfi
2429 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2430 ;;;
2431 ;;; Implementation of the Imaginary Error function Erfi(z)
2432 ;;;
2433 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2438 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2446 ;;; erfi has mirror symmetry
2450 ;;; erfi is an odd function
2454 ;;; erfi is an simplifying function
2458 ;;; erfi distributes over bags
2462 ;;; Derivative of the Error function erfi
2469 grad)
2471 ;;; Integral of the Error function erfi
2480 integral)
2482 ;;; ----------------------------------------------------------------------------
2487 ;; Look for the limit of the arguments.
2490 ;; Handle infinities at this place.
2494 ;; All other cases are handled by the simplifier of the function.
2497 ;;; ----------------------------------------------------------------------------
2502 ;; Use the relationship erfi(z) = -%i*erf(%i*z)
2506 ;; Take care to return real results when the argument is real.
2520 ;; Check for specific values
2526 ;; Check for numerical evaluation. Use erfi(z) = -%i*erf(%i*z).
2531 ;; Argument simplification
2539 ;; Representation through equivalent functions
2541 ($hypergeometric_representation
2549 ;; Transformation to Erf or Erfc
2554 (%erf
2556 (%erfc
2564 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2565 ;;;
2566 ;;; The implementation of the Inverse Error function
2567 ;;;
2568 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2573 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2575 ;;; Set properties to give full support to the parser and display
2583 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2585 ;;; The Inverse Error function is a simplifying function
2589 ;;; The Inverse Error function distributes over bags
2593 ;;; inverse_erf is the inverse of the erf function
2598 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2600 ;;; Differentiation of the Inverse Error function
2608 grad)
2610 ;;; Integral of the Inverse Error function
2617 integral)
2619 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2621 ;;; We support a simplim%function. The function is looked up in simplimit and
2622 ;;; handles specific values of the function.
2627 ;; Look for the limit of the argument.
2630 ;; Handle an argument 1 at this place
2634 ;; All other cases are handled by the simplifier of the function.
2637 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2654 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2655 ;;;
2656 ;;; The implementation of the Inverse Complementary Error function
2657 ;;;
2658 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2663 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2665 ;;; Set properties to give full support to the parser and display
2673 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2675 ;;; Inverse Complementary Error function is a simplifying function
2679 ;;; Inverse Complementary Error function distributes over bags
2683 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2685 ;;; inverse_erfc is the inverse of the erfc function
2691 ;;; Differentiation of the Inverse Complementary Error function
2699 grad)
2701 ;;; Integral of the Inverse Complementary Error function
2709 integral)
2711 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2713 ;;; We support a simplim%function. The function is looked up in simplimit and
2714 ;;; handles specific values of the function.
2719 ;; Look for the limit of the argument.
2722 ;; Handle an argument 0 at this place
2729 ;; All other cases are handled by the simplifier of the function.
2732 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2749 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2751 ;;; Implementation of the Newton algorithm to calculate numerical values
2752 ;;; of the Inverse Error functions in float or bigfloat precision.
2753 ;;; The algorithm is a modified version of the routine in newton1.mac.
2790 ;; Compute inverse_erf(z) for z a real or complex number, including
2791 ;; bigfloat objects. The value is computing using a Newton iteration
2792 ;; to solve erf(x) = z.
2795 z)
2799 z))
2801 ;; inverse_erf is odd because erf is.
2806 ;; Find an approximate solution for x = inverse_erf(z).
2809 ;; For small z, inverse_erf(z) = z*sqrt(%pi)/2
2810 ;; + O(z^3). Thus, x = z*sqrt(%pi)/2 is our
2811 ;; initial starting point.
2814 ;; For |z| > 1 and realpart(z) >= 0, we have
2815 ;; the asymptotic series z = erf(x) = 1 -
2816 ;; exp(-x^2)/x/sqrt(%pi).
2817 ;;
2818 ;; Then
2819 ;; x = sqrt(-log(x*sqrt(%pi)*(1-z))
2820 ;;
2821 ;; We can use this as a fixed-point iteration
2822 ;; to find x, and we can start the iteration at
2823 ;; x = 1. Just do one more iteration. I (RLT)
2824 ;; think that's close enough to get the Newton
2825 ;; algorithm to converge.
2832 result)))
2835 ;; Try to pick a reasonable epsilon value for the
2836 ;; Newton iteration.
2847 ;; Derivative of erf(x)
2850 #'diff
2851 z
2853 eps)))))))
2859 z))
2865 ;; Find an approximate solution for x =
2866 ;; inverse_erfc(z). We have inverse_erfc(z) =
2867 ;; inverse_erf(1-z), so that's a good starting point.
2868 ;; We don't need full precision, so a float value is
2869 ;; good enough. But if 1-z is 1, inverse_erf is
2870 ;; undefined, so we need to do something else.
2881 result)))
2887 ;; Derivative of erfc(x)
2890 #'diff
2891 z
2893 eps)))))))
2895 ;; Newton iteration for solving f(x) = z, given f and the derivative
2896 ;; of f.
2908 (intl:gettext "bf-newton: failed to converge after ~M iterations: delta = ~S, x = ~S eps = ~S")
2909 count delta x eps)
2910 x))
2918 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2919 ;;;
2920 ;;; Implementation of the Fresnel Integral S(z)
2921 ;;;
2922 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2927 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2929 ;;; Set properties to give full support to the parser and display
2937 ;;; fresnel_s is a simplifying function
2941 ;;; fresnel_s distributes over bags
2945 ;;; fresnel_s has mirror symmetry
2949 ;;; fresnel_s is an odd function
2950 ;;;
2951 ;;; Maxima has two mechanism to define a reflection property
2952 ;;; 1. Declare the feature oddfun or evenfun
2953 ;;; 2. Put a reflection rule on the property list
2954 ;;;
2955 ;;; The second way is used for the trig functions. We use it here too.
2957 ;;; This would be the first way to give the property of an odd function.
2958 ;(eval-when
2959 ; #+gcl (load eval)
2960 ; #-gcl (:load-toplevel :execute)
2961 ; (let (($context '$global) (context '$global))
2962 ; (meval '(($declare) %fresnel_s $oddfun))
2963 ; ;; Let's remove built-in symbols from list for user-defined properties.
2964 ; (setq $props (remove '%fresnel_s $props))))
2968 ;;; Differentiation of the Fresnel Integral S
2973 grad)
2975 ;;; Integration of the Fresnel Integral S
2984 integral)
2986 ;;; Limits of the Fresnel Integral S
2998 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3004 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3005 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3014 ;; For very small x, we have
3015 ;; fresnel_s(x) = %pi/6*z^3
3016 ;; fresnel_c(x) = x
3034 (maxima::merror (intl:gettext "fresnel: series expansion failed for (COMPLEX-BFLOAT-FRESNEL ~:M).") z))
3053 (* sqrt-pi
3054 z)))
3056 (* sqrt-pi
3057 z)))
3068 ;; A pure imaginary argument. Use fresnel_s(%i*x)=-%i*fresnel_s(x).
3074 fs))))
3078 ;; A pure imaginary argument. Use fresnel_c(%i*x)=%i*fresnel_c(x).
3084 fc))))
3092 ;; Check for specific values
3098 ;; Check for numerical evaluation
3102 ;; Check for argument simplification
3109 ;; Representation through equivalent functions
3111 ($erf_representation
3126 ($hypergeometric_representation
3136 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3137 ;;;
3138 ;;; Implementation of the Fresnel Integral C(z)
3139 ;;;
3140 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3145 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3147 ;;; Set properties to give full support to the parser and display
3155 ;;; fresnel_c is a simplifying function
3159 ;;; fresnel_c distributes over bags
3163 ;;; fresnel_c has mirror symmetry
3167 ;;; fresnel_c is an odd function
3168 ;;; Maxima has two mechanism to define a reflection property
3169 ;;; 1. Declare the feature oddfun or evenfun
3170 ;;; 2. Put a reflection rule on the property list
3171 ;;;
3172 ;;; The second way is used for the trig functions. We use it here too.
3176 ;;; Differentiation of the Fresnel Integral C
3181 grad)
3183 ;;; Integration of the Fresnel Integral C
3192 integral)
3194 ;;; Limits of the Fresnel Integral C
3206 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3213 ;; Check for specific values
3219 ;; Check for numerical evaluation
3224 ;; Check for argument simplification
3231 ;; Representation through equivalent functions
3233 ($erf_representation
3248 ($hypergeometric_representation
3259 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3260 ;;;
3261 ;;; Implementation of the Beta function
3262 ;;;
3263 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3265 ;;; The code for the implementation of the beta function is in the files
3266 ;;; csimp2.lisp, simp.lisp and mactex.lisp.
3267 ;;; At this place we only implement the operator property SYMMETRIC.
3269 ;;; Beta is symmetric beta(a,b) = beta(b,a)
3276 ;; Let's remove built-in symbols from list for user-defined properties.
3279 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3280 ;;;
3281 ;;; Implementation of the Incomplete Beta function
3282 ;;;
3283 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3288 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3293 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3301 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3305 ;;; beta_incomplete distributes over bags
3309 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3313 ;; Derivative wrt a
3321 z)
3323 ;; Derivative wrt b
3341 ;; The derivative wrt z
3345 grad)
3347 ;;; Integral of the Incomplete Beta function
3351 nil
3352 nil
3356 integral)
3358 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3374 ;; Check for specific values
3382 a b z))
3384 z)
3389 ;; z=1 and realpart(b)>0. Simplify to a Beta function.
3390 ;; If we have to evaluate numerically preserve the type.
3405 ;; The argument a is zero or a is negative and the argument b is
3406 ;; not in a valid range. Beta incomplete is undefined.
3407 ;; It would be more correct to return Complex infinity.
3412 ;; Check for numerical evaluation in Float or Bigfloat precision
3420 ;; Negative integer a and b is in a valid range. Expand.
3421 ($rectform
3432 ;; Negative integer a and b is in a valid range. Expand.
3433 ($rectform
3437 ;; Argument simplifications and transformations
3443 ;; Expand for b a positive integer and a not a negative integer.
3457 ;; Expand for a a positive integer.
3473 ;; Expand for a a negative integer and b an integer with b <= -a.
3502 index))
3505 index)))
3531 index)))
3550 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3557 ($rectform
3586 (merror (intl:gettext "beta_incomplete: continued fractions failed for beta_incomplete(~:M, ~:M, ~:M).") a b z))
3624 result))))))
3626 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3627 ;;;
3628 ;;; Implementation of the Generalized Incomplete Beta function
3629 ;;;
3630 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3635 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3643 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3646 simp-beta-incomplete-generalized operators)
3648 ;;; beta_incomplete_generalized distributes over bags
3652 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3654 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
3655 ;;; but not on the negative real axis and for z1 or z2 real and > 1.
3656 ;;; We support a conjugate-function which test these cases.
3659 conjugate-beta-incomplete-generalized conjugate-function)
3683 a b z1 z2)))))))
3685 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3689 ;; Derivative wrt a
3701 z1)
3707 z2)
3710 ;; Derivative wrt b
3733 ;; The derivative wrt z1
3739 ;; The derivative wrt z2
3745 grad)
3747 ;;; Integral of the Incomplete Beta function
3751 nil
3752 nil
3753 ;; Integral wrt z1
3757 ;; Integral wrt z2
3762 integral)
3764 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3776 ;; Check for specific values
3784 a b z1 z2))
3797 a b z1 z2))
3822 ;; Check for numerical evaluation in Float or Bigfloat precision
3835 ;; Check for argument simplifications and transformations
3904 index))
3907 index)))
3936 index)))
3947 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3948 ;;;
3949 ;;; Implementation of the Regularized Incomplete Beta function
3950 ;;;
3951 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3956 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3964 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3967 simp-beta-incomplete-regularized operators)
3969 ;;; beta_incomplete_regularized distributes over bags
3973 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3977 ;; Derivative wrt a
3993 ;; Derivative wrt b
4010 ;; The derivative wrt z
4015 grad)
4017 ;;; Integral of the Generalized Incomplete Beta function
4021 nil
4022 nil
4023 ;; Integral wrt z
4029 integral)
4031 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
4042 ;; Check for specific values
4050 a b z))
4070 ;; Problem: for b a negative integer the Regularized Incomplete
4071 ;; Beta function is defined to be zero. BUT: When we calculate
4072 ;; e.g. beta_incomplete(1.0,-2.0,1/2)/beta(1.0,-2.0) we get the
4073 ;; result -3.0, because beta_incomplete and beta are defined for
4074 ;; for this case. How do we get a consistent behaviour?
4080 ;; TODO: The following line does not work for bigfloats.
4082 ;; Does $beta_incomplete or simpbeta underflow in this case?
4088 ;; Check for numerical evaluation in Float or Bigfloat precision
4092 beta ibeta )
4096 ;; In case of underflow (see bug #2999) or precision loss use bigfloats
4097 ;; and emporarily give some extra precision but avoid fpprec dependency.
4098 ;; Is this workaround correct for complex values?
4107 ($rectform
4111 ;; Check for argument simplifications and transformations
4135 index))
4138 index)))
4158 index)))
4166 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;