| blob | 448f4d162559e6e7a662cdd0590d994249031a87 |
1 /*
2 * special_functions - special functions (e.g.: gamma, zeta, psi)
3 *
4 * Copyright (C) 2013 Christoph Zurnieden
5 *
6 * special_functions is open software; you can redistribute it and/or modify
7 * it under the terms of the version 2.1 of the GNU Lesser General Public
8 * License as published by the Free Software Foundation.
9 *
10 * special_functions is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
13 * General Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL. You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 *
20 * @(#) $Revision: 30.4 $
21 * @(#) $Id: specialfunctions.cal,v 30.4 2013/08/11 08:41:38 chongo Exp $
22 * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/specialfunctions.cal,v $
23 *
24 * Under source code control: 2013/08/11 01:31:28
25 * File existed as early as: 2013
26 */
29 /*
30 * hide internal function from resource debugging
31 */
32 static resource_debug_level;
33 resource_debug_level = config("resource_debug", 0);
36 /*
37 get dependencies
38 */
39 read -once zeta2;
42 /*
43 zeta2(x,s) is in the extra file "zeta2.cal" because of a different license:
44 GPL instead of LGPL.
45 */
46 define zeta(s)
47 {
48 /* Not the best way, I'm afraid, but a way. */
49 return hurwitzzeta(s,1);
50 }
52 define psi0(z){
53 local i k m x y eps_digits eps ret;
55 /*
56 One can use the Stirling series, too, which might be faster for some
57 values. The series used here converges very fast but needs a lot of
58 bernoulli numbers which are quite expensive to compute.
59 */
61 eps = epsilon();
62 epsilon(eps*1e-3);
63 if(isint(z) && z<=0)return newerror("psi(z); z is a negative integer or 0");
64 if(re(z) < 0){
65 return ( pi() * cot(pi() * (0-z) ) ) + psi0(1-z);
66 }
67 eps_digits = digits(1/epsilon());
68 /*
69 R.W. Gosper has r = .41 as the relation, empirical tests showed that
70 for d<100 r = .375 is sufficient and r = .396 for d<200.
71 It does not save much time but for a long series even a little
72 can grow large.
73 */
74 if(eps_digits < 100)
75 k = 2 * ceil(.375 * eps_digits);
76 else if(eps_digits < 200)
77 k = 2 * ceil(.395 * eps_digits);
78 else
79 k = 2 * ceil(11/27 * eps_digits);
80 m = 0;
81 y = (z+k)^2;
82 x = 0.0;
83 /*
84 There is a chance to speed up the first partial sum with binary splitting.
85 The second partial sum is dominated by the calculation of the Bernoulli
86 numbers but can profit from binary splitting when the Bernoulli numbers are
87 already cached.
88 */
89 for(i = 1; i <= (k/2); i++){
90 m = 1 / ( z + 2*i - 1) + 1 / ( z + 2*i -2 ) + m;
91 x = ( x + bernoulli(k-2*i+2) / (k-2*i+2) ) / y ;
92 }
93 ret = ln(z+k) - 1 / ( 2 * (z+k) ) - x - m ;
94 epsilon(eps);
95 return ret;
96 }
98 define psi(z)
99 {
100 return psi0(z);
101 }
103 define polygamma(m,z)
104 {
105 /*
106 TODO:
107 http://functions.wolfram.com/GammaBetaErf/PolyGamma2/16/01/01/0002/
108 */
109 if(!isint(m))return newerror("polygamma(m,z); m not an integer");
110 if( m<0 )return newerror("polygamma(m,z); m is < 0");
111 /*
112 Reflection formula not implemented yet, needs cot-differentiation
113 http://functions.wolfram.com/ElementaryFunctions/Cot/20/02/0003/
114 which is not implemented yet.
115 */
116 if( m == 0){
117 return psi0(z);
118 }
119 /*
120 use factorial for m a (small) integer?
121 use lngamma for m large?
122 */
123 if(isodd(m+1)){
124 return (-1)* gamma(m+1) * hurwitzzeta(m+1,z)
125 }
126 return gamma(m+1) * hurwitzzeta(m+1,z);
127 }
129 /*
130 Cache for the variable independent coefficients in the sum for the
131 Gamma-function.
132 */
133 static __CZ__Ck;
134 /*
135 Log-Gamma function for Re(z) > 0.
136 */
137 define __CZ__lngammarp(z)
138 {
139 local epsilon accuracy a factrl sum n ret holds_enough term;
141 epsilon = epsilon();
142 accuracy = digits(int(1/epsilon)) + 3;
144 epsilon(1e-18);
145 a = ceil(1.252850440912568095810522965 * accuracy);
147 epsilon(epsilon*10^(-(digits(1/epsilon)//2)));
149 holds_enough = 0;
151 if(size( __CZ__Ck) != a) {
152 __CZ__Ck = mat[a];
153 holds_enough = 1;
154 }
156 factrl = 1.0;
158 __CZ__Ck[0] = sqrt(2*pi()); /* c_0*/
159 for(n = 1; n < a; n++){
160 if(holds_enough == 1){
161 __CZ__Ck[n] = (a - n)^(n - 0.5) * exp(a - n) /factrl;
162 factrl *= -n
163 }
164 }
165 sum = __CZ__Ck[0];
166 for (n = 1; n < a; n++){
167 sum += __CZ__Ck[n]/(z+n);
168 }
170 ret = ln(sum)+(-(z+a)) + ln(z+a)*( z+0.5);
171 ret = ret-ln(z);
173 /*
174 Will take some time for large im(z) but almost all time is spend above
175 in that case.
176 */
177 if(im(ret))
178 ret = re(ret) + ln( exp( im(ret) *1i ) );
180 epsilon(epsilon);
181 return ret;
182 }
184 /* Simple lngamma with low precision*/
185 define __CZ__lngamma_lanczos(z){
186 local a k g zghalf lanczos;
187 mat lanczos[15] = {
188 9.9999999999999709182e-1,
189 5.7156235665862923516e1,
190 -5.9597960355475491248e1,
191 1.4136097974741747173e1,
192 -4.9191381609762019978e-1,
193 3.3994649984811888638e-5,
194 4.6523628927048576010e-5,
195 -9.8374475304879566105e-5,
196 1.5808870322491249322e-4,
197 -2.1026444172410489480e-4,
198 2.1743961811521265523e-4,
199 -1.6431810653676390482e-4,
200 8.4418223983852751308e-5,
201 -2.6190838401581411237e-5,
202 3.6899182659531626821e-6
203 };
204 g = 607/128;
205 z = z-1;
206 a = 0;
207 for(k = 12; k >= 1; k--){
208 a += lanczos[k]/(z+k);
209 }
210 a += lanczos[0];
211 zghalf = z + g + .5;
212 return ( ln(sqrt(2*pi())) + ln(a) -zghalf ) + (z+.5)*ln( zghalf );
213 }
215 /* Prints the Spouge coefficients actually in use. */
216 define __CZ__print__CZ__Ck(){
217 local k;
218 if(size(__CZ__Ck) <=1){
219 __CZ__lngammarp(2.2-2.2i);
220 }
221 for(k=0;k<size(__CZ__Ck);k++){
222 print __CZ__Ck[k];
223 }
224 }
226 /*Kronecker delta function */
227 define kroneckerdelta(i,j){
228 if(isnull(j)){
229 if(i==0) return 1;
230 else return 0;
231 }
232 if(i!=j) return 0;
233 else return 1;
234 }
236 /* (Pre-)Computed terms of the Stirling series */
237 static __CZ__precomp_stirling;
238 /* Number of terms in mat[0,0], precision in mat[0,1] with |z| >= mat[0,2]*/
239 static __CZ__stirling_params;
241 define __CZ__lngstirling(z,n){
242 local k head sum z2 bernterm zz;
243 head = (z-1/2)*ln(z)-z+(ln(2*pi())/2);
244 sum = 0;
245 bernterm=0;
246 zz = z;
247 z2 = z^2;
249 if(size(__CZ__precomp_stirling)<n){
250 __CZ__precomp_stirling = mat[n];
251 for(k=1;k<=n;k++){
252 bernterm = bernoulli(2*k)/( 2*k*(2*k-1 ) );
253 __CZ__precomp_stirling[k-1] = bernterm;
254 sum += bernterm /zz;
255 zz *= z2;
256 }
257 }
258 else{
259 for(k=1;k<=n;k++){
260 bernterm = __CZ__precomp_stirling[k-1];
261 sum += bernterm /zz;
262 zz *= z2;
263 }
264 }
265 return head + sum;
266 }
268 /* Compute error for Stirling series for "z" with "k" terms*/
269 define R(z,k){
270 local a b ab stirlingterm;
271 stirlingterm = bernoulli(2*k)/( 2*k*(2*k-1)*z^(2*k));
272 a = abs( stirlingterm );
273 b = (cos(.5*arg(z)^(2*k)));
274 ab = a/b;
275 /* a/b might round to zero */
276 if( abs(ab) < epsilon()){
277 return digits(1/epsilon());
278 }
279 return digits(1/(a/b));
280 ##return a/b;
281 }
283 /*Low precision lngamma(z) for testing the branch correction */
284 define lngammasmall(z){
285 local ret eps flag corr;
286 flag = 0;
287 if(isint(z)){
288 if(z <=0)
289 return newerror("gamma(z): z is a negative integer");
290 else{
291 /* may hold up accuracy a bit longer, but YMMV */
292 if(z < 20)
293 return ln((z-1)!);
294 else
295 return __CZ__lngamma_lanczos(z);
296 }
297 }
298 else{
299 if(re(z)<0){
300 if(im(z) && im(z)<0){
301 z = conj(z);
302 flag = 1;
303 }
305 ret = ln( pi() / sin(pi() *z ) ) -__CZ__lngamma_lanczos(1-z);
307 if(!im(z)){
308 if(abs(z) <= 1/2)
309 ret = re(ret) - pi()*1i;
310 else if( isodd(floor(abs(re(z)))) ){
311 ret = re(ret) + ( ceil(abs(z)) * pi() * 1i);
312 }
313 else if( iseven(floor(abs(re(z)))) ){
314 /* < n+1/2 */
315 if(iseven(floor(abs(z)))){
316 ret = re(ret) + ( ceil(abs(z)-1/2 -1 ) * pi() * 1i);
317 }
318 else{
319 ret = re(ret) + ( ceil(abs(z) -1/2 ) * pi() * 1i);
320 }
321 }
322 }
323 else{
324 corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);
325 ret = ret +2*(corr *pi() )*1i;
326 }
327 if(flag == 1)
328 ret = conj(ret);
329 return ret;
330 }
331 ret = (__CZ__lngamma_lanczos(z));
332 return ret;
333 }
334 }
337 /*
338 The logarithmic gamma function lngamma(z)
339 It has a different principal branch than log(gamma(z))
340 */
341 define lngamma(z)
342 {
343 local ret eps flag increasedby Z k tmp tmp2 decdigits;
344 local corr d37 d52 termcount;
345 /* z \in \mathbf{Z}*/
346 if(isint(z)){
347 /*The gamma-function has poles at zero and the negative integers*/
348 if(z <=0)
349 return newerror("lngamma(z): z is a negative integer");
350 else{
351 /* may hold up accuracy a bit longer, but YMMV */
352 if(z < 20)
353 return ln( (z-1)! );
354 else
355 return __CZ__lngammarp(z);
356 }
357 }
358 else{/*if(isint(z))*/
359 if(re(z) > 0){
360 flag = 0;
361 eps = epsilon();
362 epsilon(eps*1e-3);
364 /* Compute values on the real line with Spouge's algorithm*/
365 if(!im(z)){#
366 ret = __CZ__lngammarp(z);
367 epsilon(eps);
368 return ret;
369 }
370 /* Do the rest with the Stirling series.*/
371 /* This code repeats down under. Booh! */
372 /* Make it a positive im(z) */
373 if(im(z)<0){
374 z = conj(z);
375 flag = 1;
376 }
377 /* Evaluate the number of terms needed */
378 decdigits = floor( digits(1/eps) );
379 /* 20 dec. digits is the default in calc(?)*/
380 if(decdigits <= 21){
381 /* set 20 as the minimum */
382 epsilon(1e-22);
383 increasedby = 0;
384 /* inflate z */
385 Z=z;
386 while(abs(z) < 10){
387 z++;
388 increasedby++;
389 }
391 ret = __CZ__lngstirling(z,11);
392 /* deflate z */
393 if(increasedby > 1){
394 for(k=0;k<increasedby;k++){
395 ret = ret - ln(Z+(k));
396 }
397 }
398 ## ret = ret + 2*pi()*1i*ceil(re(x+1)/2-3/4 );
399 /* Undo conjugate if one has been applied */
400 if(flag == 1)
401 ret = conj(ret);
402 epsilon(eps);
403 return ret;
404 }/* if(decdigits <= 20) */
405 else{
406 /*
407 Compute the number of terms for the Stirling-sum first.
408 The error is bound by
409 |T(k)|
410 |R(k,z)| <= ----------------------------------------
411 2 k 2 k
412 2 k (2 k - 1) z cos(0.5 arg(z) )
413 with
414 1 - 2 k
415 B(2 k) z
416 T(k) = ----------------------
417 2 k (2 k - 1)
418 the terms of the Stirling-sum and B(n) the n-th Bernoulli number.
420 We do know the current precision and the result of the helper
421 function __CZ__R(z,k) returns the same number as "decdigits" if
422 fed with the correct "z" and "n".
423 The largest error is near the real line, so we try Z = re(1-z).
424 If Z is small, say < 20 we increment it until it is.
425 Computing Bernoulli numbers is expensive but they are small on
426 average, get cached and because of the method used, computing
427 any Bernoulli number > 0 produces the rest as a by-product.
428 Nevertheless we could use a faster asymptotic approximation.
430 We increment the number of terms "k" in R(z,k) until the result
431 quits incrementing (it may even get smaller!). If the result is
432 still smaller than the current precision we increment "z" with
433 fixed "k" untill the result quits incrementing.
434 The results, the current precision, abs(re(z)) and "k" are kept.
436 BTW: incrementing the number of terms might be more costly than
437 incrementing "z" -- computing large Bernoulli numbers vs.
438 computing a large number of complex logarithms is a fight with
439 a hard to know result -- and that the series isn't convergent
440 is of not much help either. E.g:
441 R(25,68) = 71 max
442 R(50,55) = 101
443 R(50,145) = 140 max
444 R(60,170) = 167 max
445 R(70,209) = 195 max
446 R(75,173) = 200 max
447 R(80,147) = 200 max
448 R(90,124) = 200 max
449 R(100,111) = 200 max
450 Bernoulli(222) has a denominator of 9388 with a 254 digit
451 numerator. Computing up to 100 complex logarithms on the
452 other side ...
454 D.E.G. Hare has found the bounds
455 |im(z)| > .37d or re(z) >= 0 and |z| > .52d
456 to be usefull to compute "z" to d digits accuracy. The numbers
457 correspond to the table above.
459 To avoid repeated expensive computation, the result is cached
460 together with the current precision. It might be a good idea
461 to keep it more permanently in a config-file?
462 */
463 d37 = decdigits * .37;
464 d52 = decdigits * .52;
465 termcount = ceil(d52);
466 if(abs(z) >= d52){
467 if(abs(im(z))>= d37 )
468 termcount = ceil(d37);
469 else
470 termcount = ceil(d52);
471 }
473 Z=z;
474 increasedby = 0;
475 /* inflate z */
476 if( abs(im(z))>= d37){
477 while(abs(z) < d52+1){
478 z++;
479 increasedby++;
480 }
481 }
482 else{
483 tmp = R(z,termcount);
484 tmp2 = tmp;
485 while(tmp2 < decdigits){
486 z++;
487 increasedby++;
488 tmp2 = R(z,termcount);
489 if(tmp2 < tmp)
490 return newerror("lngamma(1): something happend that "
491 "should not have happend");
492 }
493 }
495 corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);
497 ret = __CZ__lngstirling(z,termcount);
499 /* deflate z */
500 if(increasedby > 1){
501 for(k=0;k<increasedby;k++){
502 ret = ret - ln(Z+(k));
503 }
504 }
505 /* Undo conjugate if one has been applied */
506 if(flag == 1)
507 ret = conj(ret);
508 epsilon(eps);
509 return ret;
510 }/* if(decdigits <= 20){...} else */
511 }/* if(re(z) > 0) */
512 else{/* re(z)<0 */
513 eps = epsilon();
514 epsilon(eps*1e-3);
516 /* Use Spouge's approximation on the real line */
517 if(!im(z)){
518 /* reflection */
519 ret = ln( pi() / sin(pi() *z ) ) - __CZ__lngammarp(1-z);
520 /* it is log(gamma) and im(log(even(-x))) = k\pi, therefore: */
521 if(abs(z) <= 1/2)
522 ret = re(ret) - pi()*1i;
523 else if( isodd(floor(abs(re(z)))) ){
524 ret = re(ret) + ( ceil(abs(z)) * pi() * 1i);
525 }
526 else if( iseven(floor(abs(re(z)))) ){
527 /* < n+1/2 */
528 if(iseven(floor(abs(z)))){
529 ret = re(ret) + ( ceil(abs(z)-1/2 -1 ) * pi() * 1i);
530 }
531 else{
532 ret = re(ret) + ( ceil(abs(z) -1/2 ) * pi() * 1i);
533 }
534 }
535 epsilon(eps);
536 return ret;
537 }/*if(!im(z))*/
538 /* Use Stirlinsg approximation for the rest of the complex plane */
539 else{
540 /* Make it a positive im(z) */
541 if(im(z)<0){
542 z = conj(z);
543 flag = 1;
544 }
545 /* Evaluate the number of terms needed */
546 decdigits = floor( digits(1/eps) );
547 /*
548 Evaluate the correction term for the imaginary part needed because
549 of the reflection.
550 See
551 D. E. G. Hare, "Computing the Principal Branch of log-Gamma",
552 Journal of Algorithms 25(2):221-236 (1997)
553 http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.84.2063
554 */
556 /* 20 dec. digits is the default in calc(?)*/
557 if(decdigits <= 21){
558 /* set 20 as the minimum */
559 epsilon(1e-22);
560 termcount = 11;
561 increasedby = 0;
562 Z=z;
563 /* inflate z */
564 if( im(z)>= digits(1/epsilon()) * .37){
565 while(abs(1-z) < 10){
566 /* making z more negative makes 1-z more positive. */
567 z--;
568 increasedby++;
569 }
570 }
571 else{
572 tmp = R(1-z,termcount);
573 tmp2 = tmp;
574 while(tmp2 < 21){
575 z--;
576 increasedby++;
577 tmp2 = R(1-z,termcount);
578 if(tmp2 < tmp)
579 return newerror("lngamma(1): something happend "
580 "that should not have happend");
581 }
582 }
584 corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);
586 /* reflection */
587 ret = ln( pi() / sin(pi() * z ) ) - __CZ__lngstirling(1-z,termcount);
589 /* deflate z */
590 if(increasedby > 0){
591 for(k=0;k<increasedby;k++){
592 ret = ret + ln(z+(k));
593 }
594 }
595 /* Apply correction term */
596 ret = ret +2*(corr *pi() )*1i;
597 /* Undo conjugate if one has been applied */
598 if(flag == 1)
599 ret = conj(ret);
601 epsilon(eps);
602 return ret;
603 }/* if(decdigits <= 20) */
604 else{
605 d37 = decdigits * .37;
606 d52 = decdigits * .52;
607 termcount = ceil(d52);
608 if(abs(z) >= d52){
609 if(abs(im(z))>= d37 )
610 termcount = ceil(d37);
611 else
612 termcount = ceil(d52);
613 }
614 increasedby = 0;
615 ##print "Z 1: ",z;
616 Z=z;
617 /* inflate z */
618 if( abs(im(z))>= d37){
619 while(abs(1-z) < d52+1){
620 /* making z more negative makes 1-z more positive. */
621 z--;
622 increasedby++;
623 }
624 }
625 else{
626 tmp = R(1-z,termcount);
627 tmp2 = tmp;
628 while(tmp2 < decdigits){
629 z--;
630 increasedby++;
631 tmp2 = R(1-z,termcount);
632 if(tmp2 < tmp)
633 return newerror("lngamma(1): something happend that "
634 "should not have happend");
635 }
636 }
637 corr = ceil( re(z)/2 -3/4 - kroneckerdelta(im(z))/4);
638 /* reflection */
639 ret = ln( pi() / sin(pi() *(z) ) ) -
640 __CZ__lngstirling(1-z,termcount);
641 /* deflate z */
642 if(increasedby > 0){
643 for(k=0;k<increasedby;k++){
644 ret = ret + ln(z+(k));
645 }
646 }
647 /* Apply correction term */
648 ret = ret +2*((corr) *pi() )*1i;
649 /* Undo conjugate if one has been applied */
650 if(flag == 1)
651 ret = conj(ret);
652 epsilon(eps);
653 return ret;
654 }/* if(decdigits <= 20){...} else */
655 }/*if(!im(z)){...} else*/
656 }/* if(re(z) > 0){...} else */
657 }/*if(isint(z)){...} else*/
658 }
660 /* Warning about large values? */
661 define gamma(z)
662 {
664 /* exp(log(gamma(z))) = exp(lngamma(z)), so use Spouge here?*/
665 local ret eps;
666 if(isint(z)){
667 if(z <=0)
668 return newerror("gamma(z): z is a negative integer");
669 else{
670 /* may hold up accuracy a bit longer, but YMMV */
671 if(z < 20)
672 return (z-1)!*1.0;
673 else{
674 eps = epsilon(epsilon()*1e-2);
675 ret = exp(lngamma(z));
676 epsilon(eps);
677 return ret;
678 }
679 }
680 }
681 else{
682 eps = epsilon(epsilon()*1e-2);
683 ret = exp(lngamma(z));
684 epsilon(eps);
685 return ret;
686 }
687 }
689 define __CZ__harmonicF( a, b ,s)
690 {
691 local c;
692 if( b == a) return s ;
693 if( b-a > 1){
694 c= (b + a) >> 1;
695 return( __CZ__harmonicF(a, c,1/a) + __CZ__harmonicF(c+1, b,1/b));
696 }
697 return (1/a+1/b);
698 }
700 define harmonic(limit)
701 {
702 if( !isint(limit) )
703 return newerror("harmonic(limit): limit is not an integer");
704 if( limit <= 0 )
705 return newerror("harmonic(limit): limit is <=0");
706 /* The binary splitting algorithm returns 0 for f(1,1,0) */
707 if( limit == 1 ) return 1;
708 return __CZ__harmonicF( 1, limit ,0);
709 }
711 /* lower incomplete gamma function */
713 /* lower
714 for z <= 1.1
715 */
716 define __CZ__gammainc_series_lower(a,z){
717 local k ret tmp eps fact;
718 ret = 0;
719 k=0;
720 tmp=1;
721 fact = 1;
722 eps = epsilon();
723 while(abs(tmp-ret) > eps){
724 tmp = ret;
725 ret += (z^(k+a))/( (a+k)*fact );
726 k++;
727 fact *= -k;
728 }
729 return gamma(a)-ret;
730 }
732 /* lower
733 for z > 1.1
734 */
735 define __CZ__gammainc_cf(a,z,n){
736 local ret k num1 denom1 num2 denom2 ;
737 ret = 0;
738 for(k=n+1;k>1;k--){
739 ret = ((1-k)*(k-1-a))/(2*k-1+z-a+ret);
740 }
741 return ((z^a*exp(-z))/(1+z-a+ret));
742 }
744 /* G(n,z) lower*/
745 define __CZ__gammainc_integer_a(a,z){
746 local k sum fact zz;
747 for(k=0;k<a;k++){
748 sum += (z^k)/(k!);
749 }
750 return (a-1)!*exp(-z)*sum;
751 }
753 /*
754 P = precision in dec digits
755 n = 1,2,3...
756 a,z => G(a,z)
757 */
758 define __CZ__endcf(n,a,z,P){
759 local ret;
761 ret = P*ln(10)+ln(4*pi()*sqrt(n)) + re(z+(3/2-a)*ln(z)-lngamma(1-a));
762 ret = ret /( sqrt( 8*(abs(z)+re(z)) ) );
763 return ret^2;
765 }
767 /* lower incomplete gamma function */ define gammainc(a,z){
768 local ret nterms eps epsilon tmp_before tmp_after n A B C sum k;
769 if(z == 0)
770 return 1;
771 if(isint(a)){
772 if(a>0){
773 if(a==1) return exp(-z); return __CZ__gammainc_integer_a(a,z);
774 } else{
775 if(a==0){
776 return -expoint(-z)+1/2*( ln(-z) - ln(-1/z) ) -ln(z);
777 } else if(a==-1){
778 return expoint(-z)+1/2*(ln(-1/z)-ln(-z))+ln(z)+(exp(-z)/z);
779 } else{
780 A = (-1)^((-a)-1)/((-a)!); B =
781 expoint(-z)-1/2*(ln(-z)-ln(1/-z))+ln(z); C = exp(-z); sum =0;
782 for(k=1;k<-a;k++){
783 sum += (z^(k-a-1))/( fallingfactorial(-a,k) );
784 } return A * B - C *sum;
785 }
786 }
787 } if(re(z)<=1.1|| re(z) < a+1){##print "series";
788 eps = epsilon(epsilon()*1e-10); ret =
789 __CZ__gammainc_series_lower(a,z); epsilon(eps); return ret;
790 } else{##print "cf";
791 eps = epsilon(epsilon()*1e-10); if(abs(exp(-z)) <= eps) return 0;
792 tmp_before = 0; tmp_after = 1; n = 1; while(ceil(tmp_before) !=
793 ceil(tmp_after)){
794 tmp_before = tmp_after; tmp_after =
795 __CZ__endcf(n++,a,z,digits(1/eps)); /* a still quite arbitrary
796 limit */ if(n > 10){
797 return newerror("gammainc: evaluating limit for continued "
798 "fraction does not converge");
799 }
800 } ret = __CZ__gammainc_cf(a,z,ceil(tmp_after)); epsilon(eps);
801 return ret;
802 }
803 }
805 define heavisidestep(x){
806 return (1+sign(x))/2;
807 }
809 define NUMBER_POSITIVE_INFINITY(){ return 1/epsilon();}
811 define NUMBER_NEGATIVE_INFINITY(){ return -(1/epsilon());}
813 static TRUE = 1;
814 static FALSE = 0;
816 define g(prec){
817 local eps ret;
818 if(!isnull(prec)){
819 eps = epsilon(prec);
820 ret = -psi(1);
821 epsilon(eps);
822 return ret;
823 }
824 return -psi(1);
825 }
827 define __CZ__series_converged(new,old,max){
828 local eps;
829 if(isnull(max))
830 eps = epsilon();
831 else
832 eps = max;
833 if( abs(re(new - old)) <= eps * abs(re(new))
834 && abs(im(new - old)) <= eps * abs(im(new)))
835 return TRUE;
836 return FALSE;
837 }
839 define __CZ__ei_power(z){
840 local tmp ei k old;
841 ei = g() + ln(abs(z)) + sgn(im(z)) * 1i * abs(arg(z));
842 ##ei = g() + ln(z) -1i*pi()*floor( (arg(z)+pi()) / (2*pi()) );
843 tmp = 1;
844 k = 1;
845 while(k){
846 tmp *= z / k;
847 old = ei;
848 ei += tmp / k;
849 if (__CZ__series_converged(ei,old)) break;
850 k++;
851 }
852 return ei;
853 }
855 define __CZ__ei_asymp(z){
856 local ei old tmp k;
857 ei = sgn(im(z)) * 1i * pi();
858 tmp = exp(z) / z;
859 for(k=1; k<=floor(abs(z))+1; k++){
860 old = ei;
861 ei += tmp;
862 if (__CZ__series_converged(ei, old)) return ei;
863 tmp *= k / z;
864 }
865 return newerror("expoint: asymptotic series does not converge");
866 }
868 define __CZ__ei_cf(z){
869 local ei c d k old;
870 ei = sgn(im(z)) * 1i * pi();
871 if(ei != 0){
872 c = 1 / ei;
873 d = 0;
874 c = 1 / (1 - z - exp(z) * c);
875 d = 1 / (1 - z - exp(z) * d);
876 ei *= d / c;
877 }
878 else{
879 c = NUMBER_POSITIVE_INFINITY();
880 d = 0;
881 c = 0;
882 d = 1 / (1 - z - exp(z) * d);
883 ei = d * (- exp(z));
884 }
885 k = 1;
886 while(1){
887 c = 1 / (2 * k + 1 - z - k * k * c);
888 d = 1 / (2 * k + 1 - z - k * k * d);
889 old = ei;
890 ei *= d / c;
891 if (__CZ__series_converged(ei, old)) break;
892 k++;
893 }
894 return ei;
895 }
897 define expoint(z){
898 local ei eps ret;
899 eps=epsilon(epsilon()*1e-5);
900 if(abs(z) >= NUMBER_POSITIVE_INFINITY()){
901 if (config("user_debug") > 0) {
902 print "expoint: abs(z) > +inf";
903 }
904 ret = sgn(im(z)) * 1i * pi() + exp(z) / z;
905 epsilon(eps);
906 return ret;
907 }
908 if(abs(z) > 2 - 1.035 * log(epsilon())){
909 if (config("user_debug") > 0) {
910 print "expoint: using asymptotic series";
911 }
912 ei = __CZ__ei_asymp(z);
913 if (!iserror(ei)){
914 ret = ei;
915 epsilon(eps);
916 return ret;
917 }
918 }
919 if(abs(z) > 1 && (re(z) < 0 || abs(im(z)) > 1)){
920 if (config("user_debug") > 0) {
921 print "expoint: using continued fraction";
922 }
923 ret = __CZ__ei_cf(z);
924 epsilon(eps);
925 return ret;
926 }
927 if(abs(z) > 0){
928 if (config("user_debug") > 0) {
929 print "expoint: using power series";
930 }
931 ret = __CZ__ei_power(z);
932 epsilon(eps);
933 return ret;
934 }
935 if(abs(z) == 0){
936 if (config("user_debug") > 0) {
937 print "expoint: abs(z) = zero ";
938 }
939 epsilon(eps);
940 return NUMBER_NEGATIVE_INFINITY();
941 }
942 }
944 define erf(z){
945 return sqrt(z^2)/z * ( 1-1/sqrt(pi()) *gammainc(1/2,z^2) );
946 }
948 define erfc(z){
949 return 1-erf(z);
950 }
952 define erfi(z){
953 return -1i*erf(1i*z);
954 }
956 define faddeeva(z){
957 return exp(-z^2)*erfc(-1i*z);
958 }
960 define gammap(a,z){
961 return gammainc(a,z)/gamma(a);
962 }
964 define gammaq(a,z){
965 return 1-gammap(a,z);
966 }
968 define lnbeta(a,b){
969 local ret eps;
970 eps=epsilon(epsilon()*1e-3);
971 ret = (lngamma(a)+lngamma(b))-lngamma(a+b);
972 epsilon(eps);
973 return ret;
974 }
976 define beta(a,b){
977 return exp(lnbeta(a,b));
978 }
980 define __CZ__ibetacf_a(a,b,z,n){
981 local A B m places;
982 if(n==1) return 1;
983 m=n-1;
984 places = highbit(1 + int(1/epsilon())) + 1;
985 A = bround((a+m-1) * (a+b+m-1) * m * (b-m) * z^2,places++);
986 B = bround((a+2*(m)-1)^2,places++);
987 return A/B;
988 }
990 define __CZ__ibetacf_b(a,b,z,n){
991 local A B m places;
992 places = highbit(1 + int(1/epsilon())) + 1;
993 m=n-1;
994 A = bround((m*(b-m)*z)/(a+2*m-1),places++);
995 B = bround(( (a+m) * (a-(a+b)*z+1+m*(2-z)) )/(a+2*m+1),places++);
996 return m+A+B;
997 }
999 /* Didonato-Morris */
1000 define __CZ__ibeta_cf_var_dm(a,b,z,max){
1001 local m f c d check del h qab qam qap eps places;
1003 eps= epsilon();
1005 if(isnull(max)) max = 100;
1006 places = highbit(1 + int(1/epsilon())) + 1;
1007 f = eps;
1008 c = f;
1009 d = 0;
1010 for(m=1;m<=max;m++){
1011 d =bround( __CZ__ibetacf_b(a,b,z,m)+__CZ__ibetacf_a(a,b,z,m)*d,places++);
1012 if(abs(d)<eps ) d=eps;
1013 c =bround( __CZ__ibetacf_b(a,b,z,m)+__CZ__ibetacf_a(a,b,z,m)/c,places++);
1014 if(abs(c)<eps ) c=eps;
1015 d=1/d;
1016 check = c*d;
1017 f = f*check;
1018 if(abs(check-1)<eps)break;
1019 }
1020 if(m > max) return newerror("ibeta: continous fraction does not converge");
1021 return f;
1022 }
1024 define betainc_complex(z,a,b){
1025 local factor ret eps cf sum k N places tmp tmp2;
1027 if(z == 0){
1028 if(re(a) > 0) return 0;
1029 if(re(a) < 0) return newerror("betainc_complex: z == 0 and re(a) < 0");
1030 }
1031 if(z == 1){
1032 if(re(b)>0) return 1;
1033 else return newerror("betainc_complex: z == 1 and re(b) < 0");
1034 }
1035 if(b<=0){
1036 if(isint(b)) return 0;
1037 else return newerror("betainc_complex: b <= 0");
1038 }
1039 if(z==1/2 && (a==b)){return 1/2;
1040 }
1041 ##if(2==1){
1042 if(isint(a) && isint(b)){
1043 eps=epsilon(epsilon()*1e-10);
1044 N = a+b-1;
1045 sum = 0;
1046 for(k=a;k<=N;k++){
1047 tmp = ln(z)*k+ln(1-z)*(N-k);
1048 tmp2 = exp(ln(comb(N,k))+tmp);
1049 sum += tmp2;
1050 }
1051 epsilon(eps);
1052 return sum
1053 }
1054 else if(re(z) <= re((a+1)/(a+b+2))){
1055 eps=epsilon(epsilon()*1e-10);
1056 places = highbit(1 + int(1/epsilon())) + 1;
1057 factor = bround(( ln(z^a*(1-z)^b ) - lnbeta(a,b) ),places);
1058 cf =bround( __CZ__ibeta_cf_var_dm(a,b,z),places);
1059 ret = factor + ln(cf);
1060 if(abs(ret//ln(2)) >= places)
1061 ret = 0;
1062 else
1063 ret = bround(exp(factor + ln(cf)),places);
1064 epsilon(eps);
1065 return ret;
1066 }
1067 else if( re(z) > re( (a+1)/(a+b+2) ) || re(1-z) < re( (b+1)/(a+b+2) ) ){
1068 ret = 1 - betainc_complex(1-z,b,a);
1069 }
1070 return ret;
1071 }
1074 /******************************************************************************/
1075 /*
1076 Purpose:
1078 __CZ__ibetaas63 computes the incomplete Beta function ratio.
1080 Licensing:
1082 This code is distributed under the GNU LGPL license.
1084 Modified:
1086 2013-08-03 20:52:05 +0000
1088 Author:
1090 Original FORTRAN77 version by KL Majumder, GP Bhattacharjee.
1091 C version by John Burkardt
1092 Calc version by Christoph Zurnieden
1094 Reference:
1096 KL Majumder, GP Bhattacharjee,
1097 Algorithm AS 63:
1098 The incomplete Beta Integral,
1099 Applied Statistics,
1100 Volume 22, Number 3, 1973, pages 409-411.
1102 Parameters:
1104 Input, x, the argument, between 0 and 1.
1106 Input, a, b, the parameters, which
1107 must be positive.
1110 Output, the value of the incomplete
1111 Beta function ratio.
1112 */
1113 define __CZ__ibetaas63(x, a, b,beta){
1114 local ai betain cx indx ns aa asb bb rx temp term value xx acu places;
1115 acu = epsilon();
1117 value = x;
1118 /* inverse incbeta calculates it already */
1119 if(isnull(beta))
1120 beta = lnbeta(a,b);
1122 if ( a <= 0.0 || b <= 0.0 ){
1123 return newerror("betainc: domain error: a < 0 and/or b < 0");
1124 }
1125 if ( x < 0.0 || 1.0 < x ){
1126 return newerror("betainc: domain error: x<0 or x>1");
1127 }
1128 if ( x == 0.0 || x == 1.0 ){
1129 return value;
1130 }
1131 asb = a + b;
1132 cx = 1.0 - x;
1134 if ( a < asb * x ){
1135 xx = cx;
1136 cx = x;
1137 aa = b;
1138 bb = a;
1139 indx = 1;
1140 }
1141 else{
1142 xx = x;
1143 aa = a;
1144 bb = b;
1145 indx = 0;
1146 }
1148 term = 1.0;
1149 ai = 1.0;
1150 value = 1.0;
1151 ns = floor( bb + cx * asb );
1153 rx = xx / cx;
1154 temp = bb - ai;
1155 if ( ns == 0 ){
1156 rx = xx;
1157 }
1158 places = highbit(1 + int(1/acu)) + 1;
1159 while(1){
1160 term = bround(term * temp * rx / ( aa + ai ),places++);
1161 value = value + term;;
1162 temp = abs ( term );
1164 if ( temp <= acu && temp <= abs(acu * value) ){
1165 value = value * exp ( aa * ln ( xx )
1166 + ( bb - 1.0 ) * ln ( cx ) - beta ) / aa;
1168 if ( indx ){
1169 value = 1.0 - value;
1170 }
1171 break;
1172 }
1174 ai = ai + 1.0;
1175 ns = ns - 1;
1177 if ( 0 <= ns ) {
1178 temp = bb - ai;
1179 if ( ns == 0 ) {
1180 rx = xx;
1181 }
1182 }
1183 else {
1184 temp = asb;
1185 asb = asb + 1.0;
1186 }
1187 }
1188 epsilon(acu);
1189 return value;
1190 }
1192 /*
1193 z
1194 /
1195 [ b - 1 a - 1
1196 1/beta(a,b) * I (1 - t) t dt
1197 ]
1198 /
1199 0
1201 */
1203 define betainc(z,a,b){
1204 local factor ret eps cf sum k N places tmp tmp2;
1206 if(im(z) || im(a) || im(b))
1207 return betainc_complex(z,a,b);
1209 if(z == 0){
1210 if(re(a) > 0) return 0;
1211 if(re(a) < 0) return newerror("betainc: z == 0 and re(a) < 0");
1212 }
1213 if(z == 1){
1214 if(re(b)>0) return 1;
1215 else return newerror("betainc: z == 1 and re(b) < 0");
1216 }
1217 if(b<=0){
1218 if(isint(b)) return 0;
1219 else return newerror("betainc: b <= 0");
1220 }
1221 if(z==1/2 && a==b){
1222 return 1/2;
1223 }
1224 return __CZ__ibetaas63(z,a,b);
1226 }
1228 define __CZ__erfinvapprox(x){
1229 local a;
1230 a =0.147;
1231 return sgn(x)*sqrt(sqrt((2/(pi()*a)+(ln(1-x^2))/2)^2-(ln(1-x^2))/a)-
1232 (2/(pi()*a)+(ln(1-x^2))/2));
1233 }
1235 /* complementary inverse errror function, faster at about x < 1-.91
1236 Henry E. Fettis. "A stable algorithm for computing the inverse error function
1237 in the 'tail-end' region" Math. Comp., 28:585-587, 1974.
1238 */
1239 define __CZ__inverffettis(x,n){
1240 local y sqrtpi oldy k places;
1241 if (isnull(n))
1242 n = 205;
1243 y = erfinvapprox(1-x);
1244 places = highbit(1 + int(1/epsilon())) + 1;
1245 sqrtpi = sqrt(pi());
1246 do
1247 {
1248 oldy = y;k++;
1249 y = bround((ln( __CZ__fettiscf(y,n) / (sqrtpi * x)))^(1/2),places);
1250 } while( abs(y - oldy)/y > epsilon());
1251 return y;
1252 }
1254 /* cf for erfc() */
1255 define __CZ__fettiscf(y,n){
1256 local k t tt r a b ;
1257 t = 1/y;
1258 tt = t^2/2;
1259 for (k=n;k> 0;k--){
1260 a = 1;
1261 b = k*tt;
1262 r = b / (a + r);
1263 }
1264 return t / (1+r);
1265 }
1267 /* inverse errror function, faster at about x<=.91*/
1268 define __CZ__inverfbin(x){
1269 local places approx flow fhigh eps high low mid fmid epsilon;
1270 approx = erfinvapprox(x);
1271 epsilon = epsilon();
1272 high = approx + 1e-4;
1273 low = -1;
1274 places = highbit(1 + int(1/epsilon)) + 1;
1275 fhigh = x-erf(high);
1276 flow = x-erf(low);
1277 while(1){
1278 mid = bround(high - fhigh * (high - low) / (fhigh - flow), places);
1279 if ((mid == low) || (mid == high))
1280 places++;
1281 fmid = x-erf(mid);
1282 if (abs(fmid) < epsilon)
1283 return mid;
1284 if (sgn(fmid) == sgn(flow)) {
1285 low = mid;
1286 flow = fmid;
1287 }
1288 else {
1289 high = mid;
1290 fhigh = fmid;
1291 }
1292 }
1293 }
1295 define erfinv(x){
1296 local ret approx a eps y old places errfunc sqrtpihalf flag k;
1297 if(x<-1 || x > 1) return newerror("erfinv: input out of domain (-1<=x<=1)");
1298 if(x == 0) return 0;
1299 if(x == -1) return NUMBER_NEGATIVE_INFINITY();
1300 if(x == +1) return NUMBER_POSITIVE_INFINITY();
1302 if(x<0){
1303 x = -x;
1304 flag = 1;
1305 }
1306 /* No need for full pecision */
1307 eps=epsilon(1e-20);
1308 if(eps >= 1e-40){
1309 /* Winitzki, Sergei (6 February 2008). "A handy approximation for the error
1310 function and its inverse"*/
1311 a = 0.147;
1312 y = sgn(x)*sqrt(sqrt((2/(pi()*a)
1313 +(ln(1-x^2))/2)^2
1314 -(ln(1-x^2))/a)
1315 -(2/(pi()*a)+(ln(1-x^2))/2));
1317 }
1318 else {
1319 /* 20 digits instead of 5 */
1320 if(x <= .91)
1321 y = __CZ__inverfbin(x);
1322 else
1323 y = __CZ__inverffettis(1-x);
1325 if(eps <= 1e-20){
1326 epsilon(eps);
1327 return y;
1328 }
1329 }
1330 epsilon(eps);
1331 /* binary digits in number (here: number = epsilon()) */
1332 places = highbit(1 + int(1/eps)) + 1;
1333 sqrtpihalf = 2/sqrt(pi());
1334 k = 0;
1335 /*
1336 Do some Newton-Raphson steps to reach final accuracy.
1337 Only a couple of steps are necessary but calculating the error function at
1338 higher precision is quite costly;
1339 */
1340 do{
1341 old = y;
1342 errfunc = bround( erf(y),places);
1343 if( abs(errfunc-x) <= eps ) break;
1344 y = bround(y-( errfunc -x) / ( sqrtpihalf * exp(-y^2)),places);
1345 k++;
1346 }while(1);
1347 /*
1348 This is not really necessary but e.g:
1349 ; epsilon(1e-50)
1350 0.00000000000000000000000000000000000000000000000001
1351 ; erfinv(.9999999999999999999999999999999)
1352 8.28769266865549025938
1353 ; erfinv(.999999999999999999999999999999)
1354 8.14861622316986460738453487549552168842204512959346
1355 ; erf(8.28769266865549025938)
1356 0.99999999999999999999999999999990000000000000000000
1357 ; erf(8.14861622316986460738453487549552168842204512959346)
1358 0.99999999999999999999999999999900000000000000000000
1359 The precision "looks too short".
1360 */
1361 if(k == 0)
1362 y = bround(y-( errfunc -x) / ( sqrtpihalf * exp(-y^2)),places);
1363 if(flag == 1)
1364 y = -y;
1365 return y;
1366 }
1369 /*
1370 * restore internal function from resource debugging
1371 */
1372 config("resource_debug", resource_debug_level),;
1373 if (config("resource_debug") & 3) {
1374 print "zeta(z)";
1375 print "psi(z)";
1376 print "polygamma(m,z)";
1377 print "lngamma(z)";
1378 print "gamma(z)";
1379 print "harmonic(limit)";
1380 print "gammainc(a,z)";
1381 print "heavisidestep(x)";
1382 print "expoint(z)";
1383 print "erf(z)";
1384 print "erfinv(x)";
1385 print "erfc(z)";
1386 print "erfi(z)";
1387 print "erfinv(x)";
1388 print "faddeeva(z)";
1389 print "gammap(a,z)";
1390 print "gammaq(a,z)";
1391 print "beta(a,b)";
1392 print "lnbeta(a,b)";
1393 print "betainc(z,a,b)";
1394 }