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