| blob | 745dcdb20b496bdb8fdf61d1cd83b81a3b87e662 |
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
78 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
92 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
93 ;;; The following functions test if numerical evaluation has to be done.
94 ;;; The functions should help to test for numerical evaluation more consitent
95 ;;; and without complicated conditional tests including more than one or two
96 ;;; arguments.
97 ;;;
98 ;;; The functions take a list of arguments. All arguments have to be a CL or
99 ;;; Maxima number. If all arguments are numbers we have two cases:
100 ;;; 1. $numer is T we return T. The function has to be evaluated numerically.
101 ;;; 2. One of the args is a float or a bigfloat. Evaluate numerically.
102 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
104 ;;; Test for numerically evaluation in float precision
114 ;;; Test for numerically evaluation in complex float precision
117 "Determine if ARGS consists of numerical values by determining if
118 the real and imaginary parts of each arg are nuemrical (but not
119 bigfloats). A non-NIL result is returned if at least one of args is
120 a floating-point value or if numer is true. If the result is
121 non-NIL, it is a list of the arguments reduced via COMPLEX-NUMBER-P"
128 ;; Always save the result from complex-number-p. But for backward
129 ;; compatibility, only set the flag if any item is a float.
133 ;; Return the values in the same order as the args!
136 ;;; Test for numerically evaluation in bigfloat precision
147 ;;; Test for numerically evaluation in complex bigfloat precision
150 "Determine if ARGS consists of numerical values by determining if
151 the real and imaginary parts of each arg are nuemrical (including
152 bigfloats). A non-NIL result is returned if at least one of args is
153 a floating-point value or if numer is true. If the result is
154 non-NIL, it is a list of the arguments reduced via COMPLEX-NUMBER-P."
162 ;; Always save the result from complex-number-p. But for backward
163 ;; compatibility, only set the flag if any item is a bfloat.
168 ;; Return the values in the same order as the args!
171 ;;; Test for numerical evaluation in any precision, real or complex.
178 ;;; Check for an integer or a float or bigfloat representation. When we
179 ;;; have a float or bigfloat representation return the integer value.
188 val)
191 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
193 ;;; The changes to the parser to connect the operator !! to double_factorial(z)
195 ;(def-mheader |$!!| (%double_factorial))
197 ;(def-led (|$!!| 160.) (op left)
198 ; (list '$expr
199 ; (mheader '$!!)
200 ; (convert left '$expr)))
202 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
203 ;;;
204 ;;; The implementation of the function Double factorial
205 ;;;
206 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
208 ;;; Double factorial distributes over bags
212 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
214 ;;; Double factorial has mirror symmetry
218 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
220 ;;; Differentiation of Double factorial
236 grad)
238 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
240 ;;; Double factorial is a simplifying function
245 ;; Positive Integer less then $factlim or $factlim is -1. Call gfact.
251 ;; Even negative integer or real representation. Not defined.
262 ;; Odd negative integer, float or complex float.
269 ;; bigfloat or complex bigfloat.
273 ;; double_factorial(inf) -> inf
283 ;; Special case double_factorial(n-1)
284 ;; Not sure if this simplification is useful.
288 ;; Special case double_factorial(2*k+n), k integer
301 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
302 ;;; Double factorial for a complex float argument. The result is a CL complex.
303 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
307 (*
314 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
315 ;;; Double factorial for a bigfloat or complex bigfloat argument
316 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
333 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
334 ;;;
335 ;;; The implementation of the Incomplete Gamma function
336 ;;;
337 ;;; A&S 6.5.3:
338 ;;;
339 ;;; inf
340 ;;; /
341 ;;; [ a - 1 - t
342 ;;; gamma_incomplete(a, z) = I t %e dt
343 ;;; ]
344 ;;; /
345 ;;; z
346 ;;;
347 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
351 ;; TODO: This is currently called by integrate-exp-special in sin.lisp.
352 ;; Need to fix that before this can be removed.
356 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
359 ;;; Incomplete Gamma distributes over bags
363 ;;; Incomplete Gamma function has not mirror symmetry for z on the negative
364 ;;; real axis. We support a conjugate-function which test this case.
371 ;; Definitely not on the negative real axis for z. Mirror symmetry.
378 ;; On the negative real axis or no information. Unsimplified.
383 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
385 ;;; Derivative of the Incomplete Gamma function
391 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2
392 ;; function and the Generalized Incomplete Gamma function
393 ;; (functions.wolfram.com), only for a>0.
409 ;; No derivative. Maxima generates a noun form.
410 nil)))
411 ;; The derivative wrt z
417 ;;; Integral of the Incomplete Gamma function
421 nil
425 integral)
427 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
429 ;;; We support a simplim%function. The function is looked up in simplimit and
430 ;;; handles specific values of the function.
435 ;; Look for the limit of the arguments.
446 ;; Handle an argument 0 at this place.
453 ;; Call the simplifier of the function.
456 ;; All other cases are handled by the simplifier of the function.
459 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
461 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
462 ;;;
463 ;;; Implementation of the Lower Incomplete gamma function:
464 ;;;
465 ;;; A&S 6.5.2
466 ;;;
467 ;;; z
468 ;;; /
469 ;;; [ a - 1 - t
470 ;;; gamma_incomplete_lower(a, z) = I t %e dt
471 ;;; ]
472 ;;; /
473 ;;; 0
474 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
475 ;;; distribute over bags (aggregates)
479 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
481 ;; (defprop %gamma_incomplete_lower ??? grad) WHAT TO PUT HERE ??
483 ;;
484 ;; Handles some special cases for the order a and simplifies it to an
485 ;; equivalent form, possibly involving erf and gamma_incomplete_lower
486 ;; to a lower order.
499 ;; Try to express gamma_incomplete_lower(a,z) in simpler terms for
500 ;; special values of the order a. If we can't, return NIL to say so.
503 ;; gamma_incomplete_lower(0, x) = 1-exp(x)
506 ;; gamma_incomplete_lower(a,z) can be simplified if a is a
507 ;; positive integer. Since gamma_incomplete_lower(a,z) =
508 ;; gamma(a) - gamma_incomplete(a,z), use gamma_incomplete to
509 ;; get the result.
513 ;; A&S 6.5.12:
514 ;;
515 ;; gamma_incomplete_lower(1/2,x^2) = sqrt(%pi)*erf(x)
516 ;;
517 ;; gamma_incomplete_lower(1/2,z) = sqrt(%pi)*erf(sqrt(x))
518 ;;
522 ;; gamma_incomplete_lower(a+n,z), where n is an integer
527 ;; gamma_incomplete_lower(a+n,z). See DLMF 8.8.7:
528 ;;
529 ;; gamma_incomplete_lower(a+n,z)
530 ;; = pochhammer(a,n)*gamma_incomplete_lower(a,z)
531 ;; -z^a*exp(-z)*sum(gamma(a+n)gamma(a+k+1)*z^k,k,0,n-1)
550 ;; See DLMF 8.8.8. For simplicity let g(a,z) = gamma_incomplete_lower(a,z).
551 ;;
552 ;; g(a,z) = gamma(a)/gamma(a-n)*g(a-n,z)
553 ;; - z^(a-1)*exp(z)*sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
554 ;;
555 ;; Rewrite:
556 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
557 ;; + gamma(a-n)/gamma(a)*z^(a-1)*exp(-z)
558 ;; * sum(gamma(a)/gamma(a-k)*z^(-k),k,0,n-1)
559 ;; Or
560 ;; g(a-n,z) = gamma(a-n)/gamma(a)*g(a,z)
561 ;; + z^(a-1)*exp(-z)
562 ;; * sum(gamma(a-n)/gamma(a-k)*z^(-k),k,0,n-1)
583 ;; gamma_incomplete_lower of (numeric) rational order.
584 ;; Expand it out so that the resulting order is between 0 and
585 ;; 1.
588 ;; a = n + order where 0 <= order < 1.
592 ;; Nothing to do if the order is already between 0 and 1
595 ;; Use gamma_incomplete(a+n,z) above. and then substitue
596 ;; a=order. This works for n positive or negative.
601 ;; No expansion so return nil to indicate that
602 nil)))
604 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
606 ;;; Incomplete Gamma function is a simplifying function
613 ;; Check for specific values
616 ;; gamma_incomplete(v,0) is gamma(v) only if the realpart(v) >
617 ;; 0. If realpart(v) <= 0, gamma_incomplete is undefined. For
618 ;; all other cases, return the noun form.
624 a z))
633 ;; Check for numerical evaluation in Float or Bigfloat precision
638 ;; a and z are Maxima numbers, at least one has a float value
645 ;; a is zero or a negative float representing an integer.
646 ;; For these cases the numerical routines of gamma-incomplete
647 ;; do not work. Call the numerical routine for the Exponential
648 ;; Integral E(n,z). The routine is called with a positive integer!.
658 ;; a and z are Maxima numbers, at least one is a complex value and
659 ;; we have at least one float part
667 ;; Call expintegral-e. See comment above.
682 ;; Call bfloat-expintegral-e. See comment above.
703 ;; Call bfloat-expintegral-e. See comment above.
713 ;; Check for transformations and argument simplification
716 ;; Integer or a symbol declared to be an integer. Expand in a series.
738 ;; Simplify gamma, but do not expand to avoid division
739 ;; by zero.
767 ;; We have a half integral order and $gamma_expand is not NIL.
768 ;; We expand in a series with the Erfc function
818 ;; Handle gamma_incomplete(a+n, z), where n is a numerical integer
823 ;; See DLMF 8.8.9: https://dlmf.nist.gov/8.8.E9
824 ;;
825 ;; gamma_incomplete(a+n,z) = pochhammer(a,n)*gamma_incomplete(a,z)
826 ;; + z^a*exp(-z)*sum(gamma(a+n)/gamma(a+k+1)*z^k,k,0,n-1)
844 ;; See http://functions.wolfram.com/06.06.17.0004.01
845 ;;
846 ;; gamma_incomplate(a,z) =
847 ;; (-1)^n*pochhammer(1-a,n)
848 ;; *[gamma_incomplete(a-n,z)
849 ;; + z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)]
850 ;;
851 ;; Rerarrange this in terms of gamma_incomplete(a-n,z):
852 ;;
853 ;; gamma_incomplete(a-n,z) =
854 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
855 ;; -z^(a-n-1)*exp(-z)*sum(z^k/pochhammer(a-n,k),k,1,n)
856 ;;
857 ;; Change the summation index to go from k = 0 to n-1:
858 ;;
859 ;; z^(a-n-1)*sum(z^k/pochhammer(a-n,k),k,1,n)
860 ;; = z^(a-n-1)*sum(z^(k+1)/pochhammer(a-n,k+1),k,0,n-1)
861 ;; = z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
862 ;;
863 ;; Thuus:
864 ;; gamma_incomplete(a-n,z) =
865 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
866 ;; -z^(a-n)*sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
868 ;; (-1)^n*gamma_incomplete(a,z)/pochhammer(1-a,n)
877 ;; sum(z^k/pochhammer(a-n,k+1),k,0,n-1)
887 ;; gamma_incomplete of (numeric) rational order. Expand it out
888 ;; so that the resulting order is between 0 and 1.
891 ;; a = n + order where 0 <= order < 1.
895 ;; Nothing to do if the order is already between 0 and 1
898 ;; Use gamma_incomplete(a+n,z) above. and then substitue
899 ;; a=order. This works for n positive or negative.
906 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
907 ;;; Numerical evaluation of the Incomplete Gamma function
908 ;;;
909 ;;; gamma-incomplete (a,z) - real and complex double float
910 ;;; bfloat-gamma-incomplete (a z) - bigfloat
911 ;;; complex-bfloat-gamma-incomplete (a z) - complex bigfloat
912 ;;;
913 ;;; Expansion in a power series for abs(x) < R, where R is
914 ;;; *gamma-radius* + real(a) if real(a) > 0 or *gamma-radius*
915 ;;; otherwise.
916 ;;;
917 ;;; (A&S 6.5.29):
918 ;;;
919 ;;; inf
920 ;;; ===
921 ;;; \ gamma(a)
922 ;;; gamma(a,z) = exp(-x)*z^a * > ------------ * z^n
923 ;;; / gamma(a+1+n)
924 ;;; ===
925 ;;; n=0
926 ;;;
927 ;;; This expansion does not work for an integer a<=0, because the Gamma function
928 ;;; in the denominator is not defined for a=0 and negative integers. For this
929 ;;; case we use the Exponential Integral E for numerically evaluation. The
930 ;;; Incomplete Gamma function and the Exponential integral are connected by
931 ;;;
932 ;;; gamma(a,z) = z^a * expintegral_e(1-a,z)
933 ;;;
934 ;;; When the series is not used, two forms of the continued fraction
935 ;;; are used. When z is not near the negative real axis use the
936 ;;; continued fractions (A&S 6.5.31):
937 ;;;
938 ;;; 1 1-a 1 2-a 2
939 ;;; gamma(a,z) = exp(-z) z^a *( -- --- --- --- --- ... )
940 ;;; z+ 1+ z+ 1+ z+
941 ;;;
942 ;;; The accuracy is controlled by *gamma-incomplete-eps* for double float
943 ;;; precision. For bigfloat precision epsilon is 10^(-$fpprec). The expansions
944 ;;; in a power series or continued fractions stops if *gamma-incomplete-maxit*
945 ;;; is exceeded and an Maxima error is thrown.
946 ;;;
947 ;;; The above fraction does not converge on the negative real axis and
948 ;;; converges very slowly near the axis. In this case, use the
949 ;;; relationship
950 ;;;
951 ;;; gamma(a,z) = gamma(a) - gamma_lower(a,z)
952 ;;;
953 ;;; The continued fraction for gamma_incomplete_lower(a,z) is
954 ;;; (http://functions.wolfram.com/06.06.10.0009.01):
955 ;;;
956 ;;; gamma_lower(a,z) = exp(-z) * z^a / cf(a,z)
957 ;;;
958 ;;; where
959 ;;;
960 ;;; -a*z z (a+1)*z 2*z (a+2)*z
961 ;;; cf(a,z) = a + ---- ---- ------- ---- -------
962 ;;; a+1+ a+2- a+3+ a+4- a+5+
963 ;;;
964 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
974 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
976 ;;; The numerical evaluation for CL float or complex values a and x
977 ;;; When the flag regularized is T, the result is divided by gamma(a) and
978 ;;; Maxima returns the numercial result for gamma_incomplete_regularized
984 ;; Compute the factor needed to scale the series or continued
985 ;; fraction. This is x^a*exp(-x) or x^a*exp(-x)/gamma(a)
986 ;; depending on whether we want a non-regularized or
987 ;; regularized form. We want to compute the factor carefully
988 ;; to avoid unnecessary overflow if possible.
992 ;; gammafloat is more accurate for real
993 ;; values of a.
1000 ;; Easy way failed. Use logs to compute the
1001 ;; result. This loses some precision, especially
1002 ;; for large values of a and/or x because we use
1003 ;; many bits to hold the exponent part.
1005 x
1007 a
1013 ;; Easy way failed, so use the log form.
1015 x)))))))
1019 ;; %gamma-incomplete compute the lower incomplete gamma
1020 ;; function, so we need to subtract that from gamma(a),
1021 ;; more or less.
1029 ;; Continued fraction used. Just multiply by the factor
1030 ;; to get the final result.
1033 ;; Compute the key part of the gamma incomplete function using either
1034 ;; a series expression or a continued fraction expression. Two values
1035 ;; are returned: the value itself and a boolean, indicating what the
1036 ;; computed value is. If the boolean non-NIL, then the computed value
1037 ;; is the lower incomplete gamma function.
1045 ;; The series expansion is done for x within a circle of radius
1046 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1047 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1048 ;; continued fraction is used.
1054 ;; For x near the negative real axis, use the
1055 ;; relationship gamma_incomplete(a,z) = gamma(a) -
1056 ;; gamma_incomplete_lower(a,z), where
1057 ;; gamma_incomplete_lower(a,z) is the lower poart of the
1058 ;; incomplete gamma function. We can evaluate that
1059 ;; using a continued fraction from
1060 ;; http://functions.wolfram.com/06.06.10.0009.01. (Note
1061 ;; that the alternative fraction,
1062 ;; http://functions.wolfram.com/06.06.10.0007.01,
1063 ;; appears to be less accurate.)
1064 ;;
1065 ;; Also note that this appears to be valid for all
1066 ;; values of x (real or complex), but we don't want to
1067 ;; use this everywhere for gamma_incomplete. Consider
1068 ;; what happens for large real x. gamma_incomplete(a,x)
1069 ;; is small, but gamma_incomplete(a,x) = gamma(x) - cf
1070 ;; will have large roundoff errors.
1083 t)))
1085 ;; Expansion in continued fractions for gamma_incomplete.
1096 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1105 ;; Return nil to indicate we used the continued fraction.
1108 ;; Expansion in a series
1116 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1119 ;; Return T to indicate we used the series and the series
1120 ;; is for the integral from 0 to x, not x to inf.
1123 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1125 ;;; This function is called for a and x real
1133 ;; The series expansion is done for x within a circle of radius
1134 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1135 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1136 ;; continued fraction is used.
1143 ;; Expansion in continued fractions of the Incomplete Gamma function
1152 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1173 ;; Expand to multiply everything out.
1174 ($expand
1175 ;; Expansion in continued fraction for the lower incomplete gamma.
1177 ;; NOTE: We want (power x a) instead of bigfloat:expt
1178 ;; because this preserves how maxima computes x^a when
1179 ;; x is negative and a is rational. For, example
1180 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1181 ;; principal value.
1194 x))))))))))))))
1197 ;; Series expansion of the Incomplete Gamma function
1203 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1217 ($rectform
1222 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1234 ;; The series expansion is done for x within a circle of radius
1235 ;; R, where R = *gamma-radius*+(realpart(a)) for realpart(a) > 0
1236 ;; and R = *gamma-radisu* for realpart(a) < 0. Otherwise a
1237 ;; continued fraction is used.
1249 ;; Expansion in continued fractions of the Incomplete Gamma function
1260 (merror (intl:gettext "gamma_incomplete: continued fractions failed for gamma_incomplete(~:M, ~:M).") a x))
1271 gm-eps))
1280 ;; Expand to multiply everything out.
1281 ($expand
1282 ;; Expansion in continued fraction for the lower incomplete gamma.
1284 ;; NOTE: We want (power x a) instead of bigfloat:expt
1285 ;; because this preserves how maxima computes x^a when
1286 ;; x is negative and a is rational. For, example
1287 ;; (-8)^(1/2) is -2. bigfloat:expt returns the
1288 ;; principal value.
1301 x))))))))))))))
1303 ;; Series expansion of the Incomplete Gamma function
1311 (merror (intl:gettext "gamma_incomplete: series expansion failed for gamma_incomplete(~:M, ~:M).") a x))
1324 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1325 ;;;
1326 ;;; Implementation of the Generalized Incomplete Gamma function
1327 ;;;
1328 ;;; z2
1329 ;;; /
1330 ;;; [ a - 1 - t
1331 ;;; gamma_incomplete_generalized(a, z1, z2) = I t %e dt
1332 ;;; ]
1333 ;;; /
1334 ;;; z1
1335 ;;;
1336 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1338 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1339 ;;; on the negative real axis.
1340 ;;; We support a conjugate-function which test this case.
1343 conjugate-gamma-incomplete-generalized conjugate-function)
1348 ;; z1 and z2 definitely not on the negative real axis.
1349 ;; Mirror symmetry.
1357 ;; On the negative real axis or no information. Unsimplified.
1362 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1364 ;;; Generalized Incomplete Gamma distributes over bags
1368 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1370 ;;; Differentiation of Generalized Incomplete Gamma function
1374 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1375 ;; and the Generalized Incomplete Gamma function (functions.wolfram.com)
1397 ;; The derivative wrt z1
1401 ;; The derivative wrt z2
1405 grad)
1407 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1409 ;;; Generalized Incomplete Gamma function is a simplifying function
1415 ;; Check for specific values
1442 ;; Check for numerical evaluation in Float or Bigfloat precision
1443 ;; Use the numerical routines of the Incomplete Gamma function
1470 ;; Check for transformations and argument simplification
1473 ;; Expand gamma_incomplete_generalized(a+n,z1,z2) with n an integer
1526 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1527 ;;;
1528 ;;; Implementation of the Regularized Incomplete Gamma function
1529 ;;;
1530 ;;; A&S 6.5.1
1531 ;;;
1532 ;;; gamma_incomplete(a, z)
1533 ;;; gamma_incomplete_regularized(a, z) = ----------------------
1534 ;;; gamma(a)
1535 ;;;
1536 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1538 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1540 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
1541 ;;; on the negative real axis.
1542 ;;; We support a conjugate-function which test this case.
1545 conjugate-gamma-incomplete-regularized conjugate-function)
1550 ;; z definitely not on the negative real axis. Mirror symmetry.
1557 ;; On the negative real axis or no information. Unsimplified.
1562 ;;; Regularized Incomplete Gamma distributes over bags
1566 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1568 ;;; Differentiation of Regularized Incomplete Gamma function
1572 ;; The derivative wrt a in terms of hypergeometric_regularized 2F2 function
1573 ;; and the Regularized Generalized Incomplete Gamma function
1574 ;; (functions.wolfram.com)
1588 ;; The derivative wrt z
1590 -1
1594 grad)
1596 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1598 ;;; Regularized Incomplete Gamma function is a simplifying function
1606 ;; Check for specific values
1614 a z))
1621 ;; Check for numerical evaluation in Float or Bigfloat precision
1625 ;; gamma_incomplete returns a regularized result
1631 ;; gamma_incomplete returns a regularized result
1643 ($rectform
1648 ;; Check for transformations and argument simplification
1651 ;; An integer. Expand the expression.
1668 ;; We have a half integral order and $gamma_expand is not NIL.
1669 ;; We expand in a series with the Erfc function
1722 ;; We factor the second summand.
1723 ;; Some factors vanish and the result is more readable.
1724 ($factor
1739 ;; We factor the second summand.
1740 ($factor
1754 ;; gamma_incomplete_regularized of numeric rational order.
1755 ;; Expand it out so that the resulting order is between 0 and
1756 ;; 1. Use gamma_incomplete_regularized(a+n,z) to do the dirty
1757 ;; work.
1760 ;; a = n + order where 0 <= order < 1.
1764 ;; Nothing to do if the order is already between 0 and 1
1767 ;; Use gamma_incomplete_regularized(a+n,z) above. and
1768 ;; then substitue a=order. This works for n positive or
1769 ;; negative.
1776 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1777 ;;;
1778 ;;; Implementation of the Logarithm of the Gamma function
1779 ;;;
1780 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1785 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1793 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1797 ;;; Logarithm of the Gamma function distributes over bags
1801 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1806 grad)
1808 ;; integrate(log_gamma(x),x) = psi[-2](x)
1814 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1821 ;; Check for specific values
1827 ;; We have zero, a negative integer or a float or bigfloat representation.
1833 ;; Check for numerical evaluation
1850 ;; Transform to Logarithm of Factorial for integer values
1851 ;; At this point the integer value is positive and not zero.
1858 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1859 ;;; The functions log-gamma-lanczos, bfloat-log-gamma and
1860 ;;; complex-bfloat-log-gamma are modified versions of the related functions
1861 ;;; gamma-lanczos, bffac and cbffac. The functions return the Logarithm of
1862 ;;; the Gamma function. If we have to calculate the quotient of Gamma functions,
1863 ;;; e. g. for the Beta function, it is much more appropriate to use the
1864 ;;; logarithmic versions to avoid overflow.
1865 ;;;
1866 ;;; Be careful log(gamma(z)) is only for realpart(z) positive equal to
1867 ;;; log_gamma(z). For a negative realpart(z) log_gamma differ by multiple of
1868 ;;; %pi from log(gamma(z)). But we always have exp(log_gamma(z))= gamma(z).
1869 ;;; The terms to get the transformation for log_gamma(-z) are taken from
1870 ;;; functions.wolfram.com.
1871 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1877 :initial-contents
1879 57.156235665862923517
1880 -59.597960355475491248
1881 14.136097974741747174
1882 -0.49191381609762019978
1883 .33994649984811888699e-4
1884 .46523628927048575665e-4
1885 -.98374475304879564677e-4
1886 .15808870322491248884e-3
1887 -.21026444172410488319e-3
1888 .21743961811521264320e-3
1889 -.16431810653676389022e-3
1890 .84418223983852743293e-4
1891 -.26190838401581408670e-4
1897 (+
1898 (*
1906 (*
1920 sum)
1926 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1947 bigfloat%pi
1950 bigfloat%pi '$%i
1968 y)))
1974 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1995 bigfloat%pi
1998 bigfloat%pi '$%i
2017 y))))
2018 ($rectform
2024 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2025 ;;;
2026 ;;; Implementation of the Error function Erf(z)
2027 ;;;
2028 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2030 ;; TODO: This is called by simp-expintegral-e in expintegral.lisp.
2031 ;; Can't remove this until that is fixed.
2035 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2037 ;;; erf has mirror symmetry
2041 ;;; erf is an odd function
2045 ;;; erf distributes over bags
2049 ;;; Derivative of the Error function erf
2056 grad)
2058 ;;; Integral of the Error function erf
2067 integral)
2069 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2072 ; See A&S 7.1.21 or http://functions.wolfram.com/06.25.26.0001.01
2074 z
2081 ;;; erf is a simplifying function
2086 ;; Check for specific values
2092 ;; Check for numerical evaluation
2105 ;; Argument simplification
2114 ;; Representation through more general functions
2116 ($hypergeometric_representation
2119 ;; Transformation to Erfc or Erfi
2124 (%erfc
2126 (%erfi
2134 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2137 ;; We use the slatec routine for float values.
2140 ;; Compute erf(z) using the relationship
2141 ;;
2142 ;; erf(z) = sqrt(z^2)/z*(1 - gamma_incomplete(1/2,z^2)/sqrt(%pi))
2143 ;;
2144 ;; When z is real sqrt(z^2)/z is signum(z). For complex z,
2145 ;; sqrt(z^2)/z = 1 if -%pi/2 < arg(z) <= %pi/2 and -1 otherwise.
2146 ;;
2147 ;; This relationship has serious round-off issues when z is small
2148 ;; because gamma_incomplete(1/2,z^2)/sqrt(%pi) is near 1.
2149 ;;
2150 ;; complex-erf is for (lisp) complex numbers; bfloat-erf is for real
2151 ;; bfloats, and complex-bfloat-erf is for complex bfloats. Care is
2152 ;; taken to return real results for real arguments and imaginary
2153 ;; results for imaginary arguments
2156 (*
2158 -1
2164 ;; Pure real argument, the result is real
2167 ;; Pure imaginary argument, the result is pure imaginary
2170 result))))
2173 ;; Warning! This has round-off problems when abs(z) is very small.
2175 ;; The argument is real, the result is real too
2176 ($realpart
2185 ;; Warning! This has round-off problems when abs(z) is very small.
2199 ;; Pure real argument, the result is real
2202 ;; Pure imaginary argument, the result is pure imaginary
2205 ;; A general complex result
2206 result))))
2210 ;; Erf(z) for all z. Z must be a CL real or complex number or a
2211 ;; BIGFLOAT or COMPLEX-BIGFLOAT object. The result will be of the
2212 ;; same type as Z.
2215 ;; Use Slatec derf, which should be faster than the series.
2218 ;; Use the series A&S 7.1.5 for small x:
2219 ;;
2220 ;; erf(z) = 2*z/sqrt(%pi) * sum((-1)^n*z^(2*n)/n!/(2*n+1), n, 0, inf)
2221 ;;
2222 ;; The threshold is approximately erf(x) = 0.5. (Doesn't
2223 ;; have to be super accurate.) This gives max accuracy when
2224 ;; using the identity erf(x) = 1 - erfc(x).
2230 k
2234 ;; The general case.
2241 (bigfloat (maxima::$bfloat (maxima::$expand (maxima::complex-bfloat-erf (maxima::to z))))))))))
2244 ;; Compute erfc(z) via 1 - erf(z) is not very accurate if erf(z) is
2245 ;; near 1. Wolfram says
2246 ;;
2247 ;; erfc(z) = 1 - sqrt(z^2)/z * (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2, z^2))
2248 ;;
2249 ;; For real(z) > 0, sqrt(z^2)/z is 1 so
2250 ;;
2251 ;; erfc(z) = 1 - (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2252 ;; = 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2)
2253 ;;
2254 ;; For real(z) < 0, sqrt(z^2)/z is -1 so
2255 ;;
2256 ;; erfc(z) = 1 + (1 - 1/sqrt(pi)*gamma_incomplete_tail(1/2,z^2))
2257 ;; = 2 - 1/sqrt(pi)*gamma_incomplete(1/2,z^2)
2279 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2280 ;;;
2281 ;;; Implementation of the Generalized Error function Erf(z1,z2)
2282 ;;;
2283 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2285 ;;; Generalized Erf has mirror symmetry
2289 ;;; Generalized Erf distributes over bags
2293 ;;; Generalized Erf is antisymmetric Erf(z1,z2) = - Erf(z2,z1)
2300 ;; Let's remove built-in symbols from list for user-defined properties.
2303 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2307 ;; derivative wrt z1
2311 ;; derviative wrt z2
2315 grad)
2317 ;;; ----------------------------------------------------------------------------
2322 ;; Look for the limit of the arguments.
2326 ;; Handle infinities at this place.
2340 ;; All other cases are handled by the simplifier of the function.
2343 ;;; ----------------------------------------------------------------------------
2348 ;; Check for specific values
2366 ;; Check for numerical evaluation. Use erf(z1,z2) = erf(z2)-erf(z1)
2389 ;; Argument simplification
2394 ;; Representation through more general functions
2396 ($hypergeometric_representation
2399 ;; Transformation to Erf
2401 ($erf_representation
2407 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2408 ;;;
2409 ;;; Implementation of the Complementary Error function Erfc(z)
2410 ;;;
2411 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2413 ;; TODO: called by simp-expintegral-e in expintegral.lisp. Need to
2414 ;; keep this around until that is fixed.
2418 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2422 ;;; Complementary Error function distributes over bags
2426 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2433 grad)
2435 ;;; Integral of the Error function erfc
2444 integral)
2446 ;;; ----------------------------------------------------------------------------
2451 ;; Look for the limit of the arguments.
2454 ;; Handle infinities at this place.
2471 ;; All other cases are handled by the simplifier of the function.
2474 ;;; ----------------------------------------------------------------------------
2479 ;; Check for specific values
2485 ;; Check for numerical evaluation.
2490 ;; Argument simplification
2496 ;; Representation through more general functions
2498 ($hypergeometric_representation
2501 ;; Transformation to Erf or Erfi
2506 (%erf
2508 (%erfi
2516 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2517 ;;;
2518 ;;; Implementation of the Imaginary Error function Erfi(z)
2519 ;;;
2520 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2522 ;; TODO: Called by integrate-exp-special in sin.lisp. That needs to
2523 ;; be fixed before this can be removed.
2527 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2529 ;;; erfi has mirror symmetry
2533 ;;; erfi is an odd function
2537 ;;; erfi distributes over bags
2541 ;;; Derivative of the Error function erfi
2548 grad)
2550 ;;; Integral of the Error function erfi
2559 integral)
2561 ;;; ----------------------------------------------------------------------------
2566 ;; Look for the limit of the arguments.
2569 ;; Handle infinities at this place.
2573 ;; All other cases are handled by the simplifier of the function.
2576 ;;; ----------------------------------------------------------------------------
2581 ;; Use the relationship erfi(z) = -%i*erf(%i*z)
2585 ;; Take care to return real results when the argument is real.
2594 ;;; erfi is an simplifying function
2599 ;; Check for specific values
2605 ;; Check for numerical evaluation. Use erfi(z) = -%i*erf(%i*z).
2610 ;; Argument simplification
2618 ;; Representation through more general functions
2620 ($hypergeometric_representation
2623 ;; Transformation to Erf or Erfc
2628 (%erf
2630 (%erfc
2638 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2639 ;;;
2640 ;;; The implementation of the Inverse Error function
2641 ;;;
2642 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2644 ;;; The Inverse Error function distributes over bags
2648 ;;; inverse_erf is the inverse of the erf function
2653 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2655 ;;; Differentiation of the Inverse Error function
2663 grad)
2665 ;;; Integral of the Inverse Error function
2672 integral)
2674 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2676 ;;; We support a simplim%function. The function is looked up in simplimit and
2677 ;;; handles specific values of the function.
2682 ;; Look for the limit of the argument.
2685 ;; Handle an argument 1 at this place
2689 ;; All other cases are handled by the simplifier of the function.
2692 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2694 ;;; The Inverse Error function is a simplifying function
2709 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2710 ;;;
2711 ;;; The implementation of the Inverse Complementary Error function
2712 ;;;
2713 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2715 ;;; Inverse Complementary Error function distributes over bags
2719 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2721 ;;; inverse_erfc is the inverse of the erfc function
2727 ;;; Differentiation of the Inverse Complementary Error function
2735 grad)
2737 ;;; Integral of the Inverse Complementary Error function
2745 integral)
2747 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2749 ;;; We support a simplim%function. The function is looked up in simplimit and
2750 ;;; handles specific values of the function.
2755 ;; Look for the limit of the argument.
2758 ;; Handle an argument 0 at this place
2765 ;; All other cases are handled by the simplifier of the function.
2768 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2770 ;;; Inverse Complementary Error function is a simplifying function
2784 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2786 ;;; Implementation of the Newton algorithm to calculate numerical values
2787 ;;; of the Inverse Error functions in float or bigfloat precision.
2788 ;;; The algorithm is a modified version of the routine in newton1.mac.
2825 ;; Compute inverse_erf(z) for z a real or complex number, including
2826 ;; bigfloat objects. The value is computing using a Newton iteration
2827 ;; to solve erf(x) = z.
2830 z)
2834 z))
2836 ;; inverse_erf is odd because erf is.
2841 ;; Find an approximate solution for x = inverse_erf(z).
2844 ;; For small z, inverse_erf(z) = z*sqrt(%pi)/2
2845 ;; + O(z^3). Thus, x = z*sqrt(%pi)/2 is our
2846 ;; initial starting point.
2849 ;; For |z| > 1 and realpart(z) >= 0, we have
2850 ;; the asymptotic series z = erf(x) = 1 -
2851 ;; exp(-x^2)/x/sqrt(%pi).
2852 ;;
2853 ;; Then
2854 ;; x = sqrt(-log(x*sqrt(%pi)*(1-z))
2855 ;;
2856 ;; We can use this as a fixed-point iteration
2857 ;; to find x, and we can start the iteration at
2858 ;; x = 1. Just do one more iteration. I (RLT)
2859 ;; think that's close enough to get the Newton
2860 ;; algorithm to converge.
2867 result)))
2870 ;; Try to pick a reasonable epsilon value for the
2871 ;; Newton iteration.
2882 ;; Derivative of erf(x)
2885 #'diff
2886 z
2888 eps)))))))
2894 z))
2900 ;; Find an approximate solution for x =
2901 ;; inverse_erfc(z). We have inverse_erfc(z) =
2902 ;; inverse_erf(1-z), so that's a good starting point.
2903 ;; We don't need full precision, so a float value is
2904 ;; good enough. But if 1-z is 1, inverse_erf is
2905 ;; undefined, so we need to do something else.
2916 result)))
2922 ;; Derivative of erfc(x)
2925 #'diff
2926 z
2928 eps)))))))
2930 ;; Newton iteration for solving f(x) = z, given f and the derivative
2931 ;; of f.
2943 (intl:gettext "bf-newton: failed to converge after ~M iterations: delta = ~S, x = ~S eps = ~S")
2944 count delta x eps)
2945 x))
2953 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2954 ;;;
2955 ;;; Implementation of the Fresnel Integral S(z)
2956 ;;;
2957 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2959 ;;; fresnel_s distributes over bags
2963 ;;; fresnel_s has mirror symmetry
2967 ;;; fresnel_s is an odd function
2968 ;;;
2969 ;;; Maxima has two mechanism to define a reflection property
2970 ;;; 1. Declare the feature oddfun or evenfun
2971 ;;; 2. Put a reflection rule on the property list
2972 ;;;
2973 ;;; The second way is used for the trig functions. We use it here too.
2975 ;;; This would be the first way to give the property of an odd function.
2976 ;(eval-when
2977 ; #+gcl (load eval)
2978 ; #-gcl (:load-toplevel :execute)
2979 ; (let (($context '$global) (context '$global))
2980 ; (meval '(($declare) %fresnel_s $oddfun))
2981 ; ;; Let's remove built-in symbols from list for user-defined properties.
2982 ; (setq $props (remove '%fresnel_s $props))))
2986 ;;; Differentiation of the Fresnel Integral S
2991 grad)
2993 ;;; Integration of the Fresnel Integral S
3002 integral)
3004 ;;; Limits of the Fresnel Integral S
3016 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3022 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3023 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3032 ;; For very small x, we have
3033 ;; fresnel_s(x) = %pi/6*z^3
3034 ;; fresnel_c(x) = x
3052 (maxima::merror (intl:gettext "fresnel: series expansion failed for (COMPLEX-BFLOAT-FRESNEL ~:M).") z))
3071 (* sqrt-pi
3072 z)))
3074 (* sqrt-pi
3075 z)))
3086 ;; A pure imaginary argument. Use fresnel_s(%i*x)=-%i*fresnel_s(x).
3092 fs))))
3096 ;; A pure imaginary argument. Use fresnel_c(%i*x)=%i*fresnel_c(x).
3102 fc))))
3106 ;;; fresnel_s is a simplifying function
3110 ;; Check for specific values
3116 ;; Check for numerical evaluation
3120 ;; Check for argument simplification
3127 ;; Representation through equivalent functions
3129 ($erf_representation
3144 ($hypergeometric_representation
3153 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3154 ;;;
3155 ;;; Implementation of the Fresnel Integral C(z)
3156 ;;;
3157 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3208 ;;; fresnel_c is a simplifying function
3212 ;; Check for specific values
3218 ;; Check for numerical evaluation
3223 ;; Check for argument simplification
3230 ;; Representation through equivalent functions
3232 ($erf_representation
3247 ($hypergeometric_representation
3256 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3257 ;;;
3258 ;;; Implementation of the Beta function
3259 ;;;
3260 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3262 ;;; The code for the implementation of the beta function is in the files
3263 ;;; csimp2.lisp, simp.lisp and mactex.lisp.
3264 ;;; At this place we only implement the operator property SYMMETRIC.
3266 ;;; Beta is symmetric beta(a,b) = beta(b,a)
3273 ;; Let's remove built-in symbols from list for user-defined properties.
3276 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3277 ;;;
3278 ;;; Implementation of the Incomplete Beta function
3279 ;;;
3280 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3285 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3287 ;;; beta_incomplete distributes over bags
3291 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3295 ;; Derivative wrt a
3303 z)
3305 ;; Derivative wrt b
3322 ;; The derivative wrt z
3326 grad)
3328 ;;; Integral of the Incomplete Beta function
3332 nil
3333 nil
3337 integral)
3339 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3350 ;; Check for specific values
3358 a b z))
3360 z)
3365 ;; z=1 and realpart(b)>0. Simplify to a Beta function.
3366 ;; If we have to evaluate numerically preserve the type.
3381 ;; The argument a is zero or a is negative and the argument b is
3382 ;; not in a valid range. Beta incomplete is undefined.
3383 ;; It would be more correct to return Complex infinity.
3388 ;; Check for numerical evaluation in Float or Bigfloat precision
3396 ;; Negative integer a and b is in a valid range. Expand.
3397 ($rectform
3408 ;; Negative integer a and b is in a valid range. Expand.
3409 ($rectform
3413 ;; Argument simplifications and transformations
3419 ;; Expand for b a positive integer and a not a negative integer.
3433 ;; Expand for a a positive integer.
3449 ;; Expand for a a negative integer and b an integer with b <= -a.
3478 index))
3481 index)))
3507 index)))
3526 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3533 ($rectform
3562 (merror (intl:gettext "beta_incomplete: continued fractions failed for beta_incomplete(~:M, ~:M, ~:M).") a b z))
3600 result))))))
3602 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3603 ;;;
3604 ;;; Implementation of the Generalized Incomplete Beta function
3605 ;;;
3606 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3608 ;;; beta_incomplete_generalized distributes over bags
3612 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3614 ;;; Generalized Incomplete Gamma function has not mirror symmetry for z1 or z2
3615 ;;; but not on the negative real axis and for z1 or z2 real and > 1.
3616 ;;; We support a conjugate-function which test these cases.
3619 conjugate-beta-incomplete-generalized conjugate-function)
3643 a b z1 z2)))))))
3645 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3649 ;; Derivative wrt a
3661 z1)
3667 z2)
3670 ;; Derivative wrt b
3693 ;; The derivative wrt z1
3699 ;; The derivative wrt z2
3705 grad)
3707 ;;; Integral of the Incomplete Beta function
3711 nil
3712 nil
3713 ;; Integral wrt z1
3717 ;; Integral wrt z2
3722 integral)
3724 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3730 ;; Check for specific values
3738 a b z1 z2))
3750 a b z1 z2))
3772 ;; Check for numerical evaluation in Float or Bigfloat precision
3785 ;; Check for argument simplifications and transformations
3854 index))
3857 index)))
3886 index)))
3897 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3898 ;;;
3899 ;;; Implementation of the Regularized Incomplete Beta function
3900 ;;;
3901 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3903 ;;; beta_incomplete_regularized distributes over bags
3907 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3911 ;; Derivative wrt a
3927 ;; Derivative wrt b
3944 ;; The derivative wrt z
3949 grad)
3951 ;;; Integral of the Generalized Incomplete Beta function
3955 nil
3956 nil
3957 ;; Integral wrt z
3963 integral)
3965 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3971 ;; Check for specific values
3979 a b z))
3997 ;; Problem: for b a negative integer the Regularized Incomplete
3998 ;; Beta function is defined to be zero. BUT: When we calculate
3999 ;; e.g. beta_incomplete(1.0,-2.0,1/2)/beta(1.0,-2.0) we get the
4000 ;; result -3.0, because beta_incomplete and beta are defined for
4001 ;; for this case. How do we get a consistent behaviour?
4007 ;; TODO: The following line does not work for bigfloats.
4009 ;; Does $beta_incomplete or simpbeta underflow in this case?
4015 ;; Check for numerical evaluation in Float or Bigfloat precision
4019 beta ibeta )
4023 ;; In case of underflow (see bug #2999) or precision loss use bigfloats
4024 ;; and emporarily give some extra precision but avoid fpprec dependency.
4025 ;; Is this workaround correct for complex values?
4034 ($rectform
4038 ;; Check for argument simplifications and transformations
4062 index))
4065 index)))
4085 index)))
4093 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;