File indexing completed on 2024-05-12 17:22:34
0001 /* floatincgamma.h: incomplete gamma function */ 0002 /* 0003 Copyright (C) 2009 Wolf Lammen. 0004 0005 This program is free software; you can redistribute it and/or modify 0006 it under the terms of the GNU General Public License as published by 0007 the Free Software Foundation; either version 2 of the License , or 0008 (at your option) any later version. 0009 0010 This program is distributed in the hope that it will be useful, 0011 but WITHOUT ANY WARRANTY; without even the implied warranty of 0012 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0013 GNU General Public License for more details. 0014 0015 You should have received a copy of the GNU General Public License 0016 along with this program; see the file COPYING. If not, write to: 0017 0018 The Free Software Foundation, Inc. 0019 59 Temple Place, Suite 330 0020 Boston, MA 02111-1307 USA. 0021 0022 0023 You may contact the author by: 0024 e-mail: ookami1 <at> gmx <dot> de 0025 mail: Wolf Lammen 0026 Oertzweg 45 0027 22307 Hamburg 0028 Germany 0029 0030 *************************************************************************/ 0031 0032 #include "floatincgamma.h" 0033 #include "floatconst.h" 0034 #include "floatcommon.h" 0035 #include "floathmath.h" 0036 0037 /* many ideas are from a paper of Serge Winitzki, 0038 "Computing the incomplete Gamma function to arbitrary precision" 0039 */ 0040 0041 /* The critical line of the lower incomplete gamma function lowgamma(a, z) 0042 is composed of all non-positive real 'a' along with z == 0. It cannot be 0043 continued on this line. Near this line lowgamma behaves approximately as 0044 1. z^a/a, a not close to a negative integer 0045 2. (z^a/a) + (-1)^k*z^h/h/k!, where -k is a negative integer very close to a 0046 and h = a + k very small 0047 If h << z the second summand in 2. becomes dominant, otherwise z^a/a is 0048 a rough estimate of lowgamma(a, z), diverging for fixed non-positive 'a* 0049 and z -> 0. 0050 Since z^a is true complex for real z < 0, lowgamma cannot be continued 0051 in the real number domain accross the critical line, except for 'a' a 0052 negative integer. 0053 0054 reglowgammanearpole computes 0055 f(x, k, h) = x^(k-h) * [lowgamma(-k+h, x) - (-1)^k * x^h / h / k!] 0056 where k is an integer >= 0. It should only be used for |h| <= 0.5 and |x| 0057 not much bigger than 2. 0058 f is useful for computing both the lower and upper incomplete gamma 0059 function for values a, x close to the critical line. The term 0060 (-1)^k * x^h/(k!*h) subtracts the pole Gamma(-k), and, in addition, 0061 lowgamma(a, x) is regularized by the factor x^(k-h), in order to eliminate 0062 the singularity introduced by the factor z^a/a. 0063 f is real-valued even for negative x, where lowgamma is true complex. 0064 f is not continuous at a = (-2k-1)/2, k integer > 0, where the subtracted 0065 pole changes */ 0066 static void 0067 reglowgammanearpole(floatnum x, int k, cfloatnum h, int digits){ 0068 floatstruct tmp; 0069 float_create(&tmp); 0070 if (float_iszero(x) || float_getexponent(x) < -digits-3){ 0071 // tiny x, return the first element of the series: k > 0? 1/(h-k) : -x/(1+h) 0072 if (k > 0){ 0073 float_addi(x, h, -k, digits+2); 0074 float_reciprocal(x, digits+1); 0075 } 0076 else{ 0077 float_sub(&tmp, &cMinus1, h, digits+2); 0078 float_div(x, x, &tmp, digits+1); 0079 } 0080 } 0081 else{ 0082 // evaluate the series 0083 int workprec = digits+5; 0084 int i = 0; 0085 int expx = float_getexponent(x); 0086 floatstruct summand; 0087 floatstruct sum; 0088 float_create(&summand); 0089 float_create(&sum); 0090 0091 float_copy(&tmp, &c1, EXACT); 0092 float_setzero(&sum); 0093 while (workprec > 0){ 0094 if (i != k){ 0095 float_addi(&summand, h, i-k, workprec); 0096 float_div(&summand, &tmp, &summand, workprec); 0097 float_add(&sum, &sum, &summand, digits+5); 0098 } 0099 workprec = digits - float_getexponent(&sum) 0100 + float_getexponent(&summand) + expx+5; 0101 if (workprec > 0){ 0102 float_mul(&tmp, &tmp, x, workprec); 0103 float_divi(&tmp, &tmp, -++i, workprec); 0104 } 0105 } 0106 float_move(x, &sum); // frees sum as a side-effect 0107 // float_free(&sum); 0108 float_free(&summand); 0109 } 0110 float_free(&tmp); 0111 } 0112 0113 /* computes x^a (withExp = 0) or exp(-x)*x^a (withExp = 1) 0114 x^a[*exp(-x)] is the regulizing factor of the lower incomplete gamma function */ 0115 static char 0116 regulizingfactor(floatnum x, cfloatnum a, int digits, char withExp){ 0117 char result; 0118 floatstruct tmp; 0119 float_create(&tmp); 0120 result = float_raise(&tmp, x, a, digits); 0121 if (result && withExp){ 0122 float_neg(x); 0123 result = float_exp(x, digits) 0124 && float_mul(&tmp, &tmp, x, digits); 0125 } 0126 float_move(x, &tmp); // frees tmp as a side-effect 0127 // float_free(&tmp); 0128 return result; 0129 } 0130 0131 /* adds (-1)^k * x^h / (h * k!) to lowgamma, accounting for the 0132 pole of the gamma function, h!= 0 */ 0133 static char 0134 addgammapole(floatnum lowgamma, cfloatnum x, int k, cfloatnum h, int digits){ 0135 // estimate the size of this summand, often it will not contribute 0136 // to the result 0137 char result = 1; 0138 float fh = float_getexponent(h) >= -37? float_asfloat(h) : 0; 0139 float xx; 0140 float fexp = aprxlog10fn(x); 0141 fexp *= fh; 0142 fexp -= aprxlog10fn(h); 0143 xx = aprxlngamma(k+1); 0144 xx *= 0.434294481903f; 0145 fexp += xx; 0146 fexp +=1; 0147 /* float fexp = (aprxlog2fn(x)*fh - aprxlog2fn(h)) * 0.301029995663981f 0148 - aprxlngamma(k+1) * 0.434294481903f + 1;*/ 0149 if (fexp > EXPMIN){ 0150 int exp = fexp; 0151 int explowgamma = float_getexponent(lowgamma); 0152 int workprec = digits + 3; 0153 if (exp < explowgamma-3) 0154 workprec = digits - explowgamma - exp -3; 0155 if (workprec > 0){ 0156 floatstruct pwr; 0157 floatstruct fct; 0158 float_create(&pwr); 0159 float_create(&fct); 0160 float_setinteger(&fct, k); 0161 result = float_raise(&pwr, x, h, workprec) 0162 && float_factorial(&fct, workprec) 0163 && float_div(&pwr, &pwr, &fct, workprec) 0164 && float_div(&pwr, &pwr, h, workprec); 0165 if ((k & 1) != 0) 0166 float_neg(&pwr); 0167 result = result && float_add(lowgamma, lowgamma, &pwr, digits+2); 0168 float_free(&fct); 0169 float_free(&pwr); 0170 } 0171 } 0172 return result; 0173 } 0174 0175 /* for 0.5 > a and ln Gamma(-a) / ln 10 < 2 * digits */ 0176 static char 0177 lowgammanearpole(floatnum x, cfloatnum a, int digits){ 0178 int k = float_getexponent(a); 0179 char result; 0180 floatstruct h; 0181 floatstruct lowgamma; 0182 floatstruct factor; 0183 float_create(&h); 0184 float_create(&lowgamma); 0185 float_create(&factor); 0186 float_sub(&h, a, &c1Div2, k < -1? 2 : k+3); 0187 k = -float_asinteger(&h); 0188 float_addi(&h, a, k, digits+2); 0189 float_copy(&lowgamma, x, digits+2); 0190 float_copy(&factor, x, digits+3); 0191 reglowgammanearpole(&lowgamma, k, &h, digits); 0192 result = regulizingfactor(&factor, a, digits, 0) 0193 && float_mul(&lowgamma, &lowgamma, &factor, digits+2) 0194 && addgammapole(&lowgamma, x, k, &h, digits); 0195 float_move(x, &lowgamma); // frees lowgamma 0196 float_free(&lowgamma); 0197 float_free(&h); 0198 return result; 0199 } 0200 0201 void 0202 testincgamma(floatnum x, cfloatnum a, int digits){ 0203 lowgammanearpole(x, a, digits); 0204 }