mirror of
https://github.com/lcn2/calc.git
synced 2025-08-16 01:03:29 +03:00
3c18e6e25b
Fix "Under source code control" date for new version.h file. Sorted the order of symbols printed by "make env". Test if <stdbool.h> exists and set HAVE_STDBOOL_H accordingly in have_stdbool.h. Added HAVE_STDBOOL_H to allow one to force this value. Added "bool.h" include file to support use of boolean symbols, true and false for pre-c99 C compilers. The "bool.h" include file defines TRUE as true, FALSE as false, and BOOL as bool: for backward compatibility. The sign in a ZVALUE is now of type SIGN, whcih is either SB32 when CALC2_COMPAT is defined, or a bool. Replaced in C source, TRUE with true, FALSE with false, and BOOL with bool.
1466 lines
31 KiB
Plaintext
1466 lines
31 KiB
Plaintext
/*
|
|
* specialfunctions - special functions (e.g.: gamma, zeta, psi)
|
|
*
|
|
* Copyright (C) 2013,2021,2023 Christoph Zurnieden
|
|
*
|
|
* Calc is open software; you can redistribute it and/or modify
|
|
* it under the terms of the version 2.1 of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation.
|
|
*
|
|
* Calc is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
|
|
* General Public License for more details.
|
|
*
|
|
* A copy of version 2.1 of the GNU Lesser General Public License is
|
|
* distributed with calc under the filename COPYING-LGPL. You should have
|
|
* received a copy with calc; if not, write to Free Software Foundation, Inc.
|
|
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
*
|
|
* Under source code control: 2013/08/11 01:31:28
|
|
* File existed as early as: 2013
|
|
*/
|
|
|
|
|
|
/*
|
|
* hide internal function from resource debugging
|
|
*/
|
|
static resource_debug_level;
|
|
resource_debug_level = config("resource_debug", 0);
|
|
|
|
|
|
/*
|
|
get dependencies
|
|
*/
|
|
read -once zeta2;
|
|
|
|
|
|
/*
|
|
zeta2(x,s) is in the extra file "zeta2.cal" because of a different license:
|
|
GPL instead of LGPL.
|
|
*/
|
|
define zeta(s)
|
|
{
|
|
/* Not the best way, I'm afraid, but a way. */
|
|
return hurwitzzeta(s, 1);
|
|
}
|
|
|
|
define psi0(z)
|
|
{
|
|
local i k m x y eps_digits eps ret;
|
|
|
|
/*
|
|
* One can use the Stirling series, too, which might be faster for some
|
|
* values. The series used here converges very fast but needs a lot of
|
|
* bernoulli numbers which are quite expensive to compute.
|
|
*/
|
|
|
|
eps = epsilon();
|
|
epsilon(eps * 1e-3);
|
|
if (isint(z) && z <= 0) {
|
|
return newerror("psi(z); z is a negative integer or 0");
|
|
}
|
|
if (re(z) < 0) {
|
|
return (pi() * cot(pi() * (0 - z))) + psi0(1 - z);
|
|
}
|
|
eps_digits = digits(1 / epsilon());
|
|
/*
|
|
* R.W. Gosper has r = .41 as the relation, empirical tests showed that
|
|
* for d<100 r = .375 is sufficient and r = .396 for d<200.
|
|
* It does not save much time but for a long series even a little
|
|
* can grow large.
|
|
*/
|
|
if (eps_digits < 100) {
|
|
k = 2 * ceil(.375 * eps_digits);
|
|
} else if (eps_digits < 200) {
|
|
k = 2 * ceil(.395 * eps_digits);
|
|
} else {
|
|
k = 2 * ceil(11 / 27 * eps_digits);
|
|
}
|
|
m = 0;
|
|
y = (z + k) ^ 2;
|
|
x = 0.0;
|
|
/*
|
|
* There is a chance to speed up the first partial sum with binary
|
|
* splitting. The second partial sum is dominated by the calculation
|
|
* of the Bernoulli numbers but can profit from binary splitting when
|
|
* the Bernoulli numbers are already cached.
|
|
*/
|
|
for (i = 1; i <= (k / 2); i++) {
|
|
m = 1 / (z + 2 * i - 1) + 1 / (z + 2 * i - 2) + m;
|
|
x = (x + bernoulli(k - 2 * i + 2) / (k - 2 * i + 2)) / y;
|
|
}
|
|
ret = ln(z + k) - 1 / (2 * (z + k)) - x - m;
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
|
|
define psi(z)
|
|
{
|
|
return psi0(z);
|
|
}
|
|
|
|
define polygamma(m, z)
|
|
{
|
|
/*
|
|
* TODO: http://functions.wolfram.com/GammaBetaErf/PolyGamma2/16/01/01/0002/
|
|
*/
|
|
if (!isint(m)) {
|
|
return newerror("polygamma(m,z); m not an integer");
|
|
}
|
|
if (m < 0) {
|
|
return newerror("polygamma(m,z); m is < 0");
|
|
}
|
|
/*
|
|
* Reflection formula not implemented yet, needs cot-differentiation
|
|
* http://functions.wolfram.com/ElementaryFunctions/Cot/20/02/0003/
|
|
* which is not implemented yet.
|
|
*/
|
|
if (m == 0) {
|
|
return psi0(z);
|
|
}
|
|
/*
|
|
* use factorial for m a (small) integer?
|
|
* use lngamma for m large?
|
|
*/
|
|
if (isodd(m + 1)) {
|
|
return (-1) * gamma(m + 1) * hurwitzzeta(m + 1, z)
|
|
}
|
|
return gamma(m + 1) * hurwitzzeta(m + 1, z);
|
|
}
|
|
|
|
/*
|
|
Cache for the variable independent coefficients in the sum for the
|
|
Gamma-function.
|
|
*/
|
|
static __CZ__Ck;
|
|
/*
|
|
Log-Gamma function for Re(z) > 0.
|
|
*/
|
|
define __CZ__lngammarp(z)
|
|
{
|
|
local epsilon accuracy a factrl sum n ret holds_enough term;
|
|
|
|
epsilon = epsilon();
|
|
accuracy = digits(int (1 / epsilon)) + 3;
|
|
|
|
epsilon(1e-18);
|
|
a = ceil(1.252850440912568095810522965 * accuracy);
|
|
|
|
epsilon(epsilon * 10 ^ (-(digits(1 / epsilon) //2)));
|
|
|
|
holds_enough = 0;
|
|
|
|
if (size(__CZ__Ck) != a) {
|
|
__CZ__Ck = mat[a];
|
|
holds_enough = 1;
|
|
}
|
|
|
|
factrl = 1.0;
|
|
|
|
__CZ__Ck[0] = sqrt(2 * pi()); /* c_0 */
|
|
for (n = 1; n < a; n++) {
|
|
if (holds_enough == 1) {
|
|
__CZ__Ck[n] = (a - n) ^ (n - 0.5) * exp(a - n) / factrl;
|
|
factrl *= -n}
|
|
}
|
|
sum = __CZ__Ck[0];
|
|
for (n = 1; n < a; n++) {
|
|
sum += __CZ__Ck[n] / (z + n);
|
|
}
|
|
|
|
ret = ln(sum) + (-(z + a)) + ln(z + a) * (z + 0.5);
|
|
ret = ret - ln(z);
|
|
|
|
/*
|
|
* Will take some time for large im(z) but almost all time is spend above
|
|
* in that case.
|
|
*/
|
|
if (im(ret)) {
|
|
ret = re(ret) + ln(exp(im(ret) * 1i));
|
|
}
|
|
|
|
epsilon(epsilon);
|
|
return ret;
|
|
}
|
|
|
|
/* Simple lngamma with low precision*/
|
|
define __CZ__lngamma_lanczos(z)
|
|
{
|
|
local a k g zghalf lanczos;
|
|
mat lanczos[15] = {
|
|
9.9999999999999709182e-1,
|
|
5.7156235665862923516e1,
|
|
-5.9597960355475491248e1,
|
|
1.4136097974741747173e1,
|
|
-4.9191381609762019978e-1,
|
|
3.3994649984811888638e-5,
|
|
4.6523628927048576010e-5,
|
|
-9.8374475304879566105e-5,
|
|
1.5808870322491249322e-4,
|
|
-2.1026444172410489480e-4,
|
|
2.1743961811521265523e-4,
|
|
-1.6431810653676390482e-4,
|
|
8.4418223983852751308e-5,
|
|
-2.6190838401581411237e-5,
|
|
3.6899182659531626821e-6
|
|
};
|
|
g = 607 / 128;
|
|
z = z - 1;
|
|
a = 0;
|
|
for (k = 12; k >= 1; k--) {
|
|
a += lanczos[k] / (z + k);
|
|
}
|
|
a += lanczos[0];
|
|
zghalf = z + g + .5;
|
|
return (ln(sqrt(2 * pi())) + ln(a) - zghalf) + (z + .5) * ln(zghalf);
|
|
}
|
|
|
|
/* Prints the Spouge coefficients actually in use. */
|
|
define __CZ__print__CZ__Ck()
|
|
{
|
|
local k;
|
|
if (size(__CZ__Ck) <= 1) {
|
|
__CZ__lngammarp(2.2 - 2.2i);
|
|
}
|
|
for (k = 0; k < size(__CZ__Ck); k++) {
|
|
print __CZ__Ck[k];
|
|
}
|
|
}
|
|
|
|
/*Kronecker delta function */
|
|
define kroneckerdelta(i, j)
|
|
{
|
|
if (isnull(j)) {
|
|
if (i == 0) {
|
|
return 1;
|
|
} else {
|
|
return 0;
|
|
}
|
|
}
|
|
if (i != j) {
|
|
return 0;
|
|
} else {
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
/* (Pre-)Computed terms of the Stirling series */
|
|
static __CZ__precomp_stirling;
|
|
/* Number of terms in mat[0,0], precision in mat[0,1] with |z| >= mat[0,2]*/
|
|
static __CZ__stirling_params;
|
|
|
|
define __CZ__lngstirling(z, n)
|
|
{
|
|
local k head sum z2 bernterm zz;
|
|
head = (z - 1 / 2) * ln(z) - z + (ln(2 * pi()) / 2);
|
|
sum = 0;
|
|
bernterm = 0;
|
|
zz = z;
|
|
z2 = z ^ 2;
|
|
|
|
if (size(__CZ__precomp_stirling) < n) {
|
|
__CZ__precomp_stirling = mat[n];
|
|
for (k = 1; k <= n; k++) {
|
|
bernterm = bernoulli(2 * k) / (2 * k * (2 * k - 1));
|
|
__CZ__precomp_stirling[k - 1] = bernterm;
|
|
sum += bernterm / zz;
|
|
zz *= z2;
|
|
}
|
|
} else {
|
|
for (k = 1; k <= n; k++) {
|
|
bernterm = __CZ__precomp_stirling[k - 1];
|
|
sum += bernterm / zz;
|
|
zz *= z2;
|
|
}
|
|
}
|
|
return head + sum;
|
|
}
|
|
|
|
/* Compute error for Stirling series for "z" with "k" terms*/
|
|
define R(z, k)
|
|
{
|
|
local a b ab stirlingterm;
|
|
stirlingterm = bernoulli(2 * k) / (2 * k * (2 * k - 1) * z ^ (2 * k));
|
|
a = abs(stirlingterm);
|
|
b = (cos(.5 * arg(z) ^ (2 * k)));
|
|
ab = a / b;
|
|
/* a/b might round to zero */
|
|
if (abs(ab) < epsilon()) {
|
|
return digits(1 / epsilon());
|
|
}
|
|
return digits(1 / (a / b));
|
|
}
|
|
|
|
/*Low precision lngamma(z) for testing the branch correction */
|
|
define lngammasmall(z)
|
|
{
|
|
local ret eps flag corr;
|
|
flag = 0;
|
|
if (isint(z)) {
|
|
if (z <= 0) {
|
|
return newerror("gamma(z): z is a negative integer");
|
|
} else {
|
|
/* may hold up accuracy a bit longer, but YMMV */
|
|
if (z < 20) {
|
|
return ln((z - 1) !);
|
|
} else {
|
|
return __CZ__lngamma_lanczos(z);
|
|
}
|
|
}
|
|
} else {
|
|
if (re(z) < 0) {
|
|
if (im(z) && im(z) < 0) {
|
|
z = conj(z);
|
|
flag = 1;
|
|
}
|
|
|
|
ret = ln(pi() / sin(pi() * z)) - __CZ__lngamma_lanczos(1 - z);
|
|
|
|
if (!im(z)) {
|
|
if (abs(z) <= 1 / 2) {
|
|
ret = re(ret) - pi() * 1i;
|
|
} else if (isodd(floor(abs(re(z))))) {
|
|
ret = re(ret) + (ceil(abs(z)) * pi() * 1i);
|
|
} else if (iseven(floor(abs(re(z))))) {
|
|
/* < n+1/2 */
|
|
if (iseven(floor(abs(z)))) {
|
|
ret = re(ret) + (ceil(abs(z) - 1 / 2 - 1) * pi() * 1i);
|
|
} else {
|
|
ret = re(ret) + (ceil(abs(z) - 1 / 2) * pi() * 1i);
|
|
}
|
|
}
|
|
} else {
|
|
corr = ceil(re(z) / 2 - 3 / 4 - kroneckerdelta(im(z)) / 4);
|
|
ret = ret + 2 * (corr * pi()) * 1i;
|
|
}
|
|
if (flag == 1) {
|
|
ret = conj(ret);
|
|
}
|
|
return ret;
|
|
}
|
|
ret = (__CZ__lngamma_lanczos(z));
|
|
return ret;
|
|
}
|
|
}
|
|
|
|
|
|
/*
|
|
The logarithmic gamma function lngamma(z)
|
|
It has a different principal branch than log(gamma(z))
|
|
*/
|
|
define lngamma(z)
|
|
{
|
|
local ret eps flag increasedby Z k tmp tmp2 decdigits;
|
|
local corr d37 d52 termcount;
|
|
/* z \in \mathbf{Z} */
|
|
if (isint(z)) {
|
|
/*The gamma-function has poles at zero and the negative integers */
|
|
if (z <= 0) {
|
|
return newerror("lngamma(z): z is a negative integer");
|
|
} else {
|
|
/* may hold up accuracy a bit longer, but YMMV */
|
|
if (z < 20) {
|
|
return ln((z - 1) !);
|
|
} else {
|
|
return __CZ__lngammarp(z);
|
|
}
|
|
}
|
|
} else { /*if(isint(z)) */
|
|
if (re(z) > 0) {
|
|
flag = 0;
|
|
eps = epsilon();
|
|
epsilon(eps * 1e-3);
|
|
|
|
/* Compute values on the real line with Spouge's algorithm */
|
|
if (!im(z)) {
|
|
ret = __CZ__lngammarp(z);
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
/* Do the rest with the Stirling series. */
|
|
/* This code repeats down under. Booh! */
|
|
/* Make it a positive im(z) */
|
|
if (im(z) < 0) {
|
|
z = conj(z);
|
|
flag = 1;
|
|
}
|
|
/* Evaluate the number of terms needed */
|
|
decdigits = floor(digits(1 / eps));
|
|
/* 20 dec. digits is the default in calc(?) */
|
|
if (decdigits <= 21) {
|
|
/* set 20 as the minimum */
|
|
epsilon(1e-22);
|
|
increasedby = 0;
|
|
/* inflate z */
|
|
Z = z;
|
|
while (abs(z) < 10) {
|
|
z++;
|
|
increasedby++;
|
|
}
|
|
|
|
ret = __CZ__lngstirling(z, 11);
|
|
/* deflate z */
|
|
if (increasedby > 1) {
|
|
for (k = 0; k < increasedby; k++) {
|
|
ret = ret - ln(Z + (k));
|
|
}
|
|
}
|
|
/* Undo conjugate if one has been applied */
|
|
if (flag == 1) {
|
|
ret = conj(ret);
|
|
}
|
|
epsilon(eps);
|
|
return ret;
|
|
} else {
|
|
d37 = decdigits * .37;
|
|
d52 = decdigits * .52;
|
|
termcount = ceil(d52);
|
|
if (abs(z) >= d52) {
|
|
if (abs(im(z)) >= d37) {
|
|
termcount = ceil(d37);
|
|
} else {
|
|
termcount = ceil(d52);
|
|
}
|
|
}
|
|
|
|
Z = z;
|
|
increasedby = 0;
|
|
/* inflate z */
|
|
if (abs(im(z)) >= d37) {
|
|
while (abs(z) < d52 + 1) {
|
|
z++;
|
|
increasedby++;
|
|
}
|
|
} else {
|
|
tmp = R(z, termcount);
|
|
tmp2 = tmp;
|
|
while (tmp2 < decdigits) {
|
|
z++;
|
|
increasedby++;
|
|
tmp2 = R(z, termcount);
|
|
if (tmp2 < tmp) {
|
|
return
|
|
newerror(strcat
|
|
("lngamma(1): something happened ",
|
|
"that shouldn't have happened"));
|
|
}
|
|
}
|
|
}
|
|
|
|
corr = ceil(re(z) / 2 - 3 / 4 - kroneckerdelta(im(z)) / 4);
|
|
|
|
ret = __CZ__lngstirling(z, termcount);
|
|
|
|
/* deflate z */
|
|
if (increasedby > 1) {
|
|
for (k = 0; k < increasedby; k++) {
|
|
ret = ret - ln(Z + (k));
|
|
}
|
|
}
|
|
/* Undo conjugate if one has been applied */
|
|
if (flag == 1) {
|
|
ret = conj(ret);
|
|
}
|
|
epsilon(eps);
|
|
return ret;
|
|
} /* if(decdigits <= 20){...} else */
|
|
} else { /* re(z)<0 */
|
|
eps = epsilon();
|
|
epsilon(eps * 1e-3);
|
|
|
|
/* Use Spouge's approximation on the real line */
|
|
if (!im(z)) {
|
|
/* reflection */
|
|
ret = ln(pi() / sin(pi() * z)) - __CZ__lngammarp(1 - z);
|
|
/* it is log(gamma) and im(log(even(-x))) = k\pi, therefore: */
|
|
if (abs(z) <= 1 / 2) {
|
|
ret = re(ret) - pi() * 1i;
|
|
} else if (isodd(floor(abs(re(z))))) {
|
|
ret = re(ret) + (ceil(abs(z)) * pi() * 1i);
|
|
} else if (iseven(floor(abs(re(z))))) {
|
|
/* < n+1/2 */
|
|
if (iseven(floor(abs(z)))) {
|
|
ret = re(ret) + (ceil(abs(z) - 1 / 2 - 1) * pi() * 1i);
|
|
} else {
|
|
ret = re(ret) + (ceil(abs(z) - 1 / 2) * pi() * 1i);
|
|
}
|
|
}
|
|
epsilon(eps);
|
|
return ret;
|
|
} else {
|
|
/* Make it a positive im(z) */
|
|
if (im(z) < 0) {
|
|
z = conj(z);
|
|
flag = 1;
|
|
}
|
|
/* Evaluate the number of terms needed */
|
|
decdigits = floor(digits(1 / eps));
|
|
/* 20 dec. digits is the default in calc(?) */
|
|
if (decdigits <= 21) {
|
|
/* set 20 as the minimum */
|
|
epsilon(1e-22);
|
|
termcount = 11;
|
|
increasedby = 0;
|
|
Z = z;
|
|
/* inflate z */
|
|
if (im(z) >= digits(1 / epsilon()) * .37) {
|
|
while (abs(1 - z) < 10) {
|
|
z--;
|
|
increasedby++;
|
|
}
|
|
} else {
|
|
tmp = R(1 - z, termcount);
|
|
tmp2 = tmp;
|
|
while (tmp2 < 21) {
|
|
z--;
|
|
increasedby++;
|
|
tmp2 = R(1 - z, termcount);
|
|
if (tmp2 < tmp) {
|
|
return
|
|
newerror(strcat
|
|
("lngamma(1): something happened ",
|
|
"that should not have happened"));
|
|
}
|
|
}
|
|
}
|
|
|
|
corr = ceil(re(z) / 2 - 3 / 4 - kroneckerdelta(im(z)) / 4);
|
|
|
|
/* reflection */
|
|
ret =
|
|
ln(pi() / sin(pi() * z)) - __CZ__lngstirling(1 - z,
|
|
termcount);
|
|
|
|
/* deflate z */
|
|
if (increasedby > 0) {
|
|
for (k = 0; k < increasedby; k++) {
|
|
ret = ret + ln(z + (k));
|
|
}
|
|
}
|
|
/* Apply correction term */
|
|
ret = ret + 2 * (corr * pi()) * 1i;
|
|
/* Undo conjugate if one has been applied */
|
|
if (flag == 1) {
|
|
ret = conj(ret);
|
|
}
|
|
|
|
epsilon(eps);
|
|
return ret;
|
|
} else {
|
|
d37 = decdigits * .37;
|
|
d52 = decdigits * .52;
|
|
termcount = ceil(d52);
|
|
if (abs(z) >= d52) {
|
|
if (abs(im(z)) >= d37) {
|
|
termcount = ceil(d37);
|
|
} else {
|
|
termcount = ceil(d52);
|
|
}
|
|
}
|
|
increasedby = 0;
|
|
Z = z;
|
|
/* inflate z */
|
|
if (abs(im(z)) >= d37) {
|
|
while (abs(1 - z) < d52 + 1) {
|
|
z--;
|
|
increasedby++;
|
|
}
|
|
} else {
|
|
tmp = R(1 - z, termcount);
|
|
tmp2 = tmp;
|
|
while (tmp2 < decdigits) {
|
|
z--;
|
|
increasedby++;
|
|
tmp2 = R(1 - z, termcount);
|
|
if (tmp2 < tmp) {
|
|
return
|
|
newerror(strcat
|
|
("lngamma(1): something happened ",
|
|
"that should not have happened"));
|
|
}
|
|
}
|
|
}
|
|
corr = ceil(re(z) / 2 - 3 / 4 - kroneckerdelta(im(z)) / 4);
|
|
/* reflection */
|
|
ret =
|
|
ln(pi() / sin(pi() * (z))) - __CZ__lngstirling(1 - z,
|
|
termcount);
|
|
/* deflate z */
|
|
if (increasedby > 0) {
|
|
for (k = 0; k < increasedby; k++) {
|
|
ret = ret + ln(z + (k));
|
|
}
|
|
}
|
|
/* Apply correction term */
|
|
ret = ret + 2 * ((corr) * pi()) * 1i;
|
|
/* Undo conjugate if one has been applied */
|
|
if (flag == 1) {
|
|
ret = conj(ret);
|
|
}
|
|
epsilon(eps);
|
|
return ret;
|
|
} /* if(decdigits <= 20){...} else */
|
|
} /*if(!im(z)){...} else */
|
|
} /* if(re(z) > 0){...} else */
|
|
} /*if(isint(z)){...} else */
|
|
}
|
|
|
|
/* Warning about large values? */
|
|
define gamma(z)
|
|
{
|
|
|
|
/* exp(log(gamma(z))) = exp(lngamma(z)), so use Spouge here? */
|
|
local ret eps;
|
|
if (isint(z)) {
|
|
if (z <= 0) {
|
|
return newerror("gamma(z): z is a negative integer");
|
|
} else {
|
|
/* may hold up accuracy a bit longer, but YMMV */
|
|
if (z < 20) {
|
|
return (z - 1) ! *1.0;
|
|
} else {
|
|
eps = epsilon(epsilon() * 1e-2);
|
|
ret = exp(lngamma(z));
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
}
|
|
} else {
|
|
eps = epsilon(epsilon() * 1e-2);
|
|
ret = exp(lngamma(z));
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
}
|
|
|
|
define __CZ__harmonicF(a, b, s)
|
|
{
|
|
local c;
|
|
if (b == a) {
|
|
return s;
|
|
}
|
|
if (b - a > 1) {
|
|
c = (b + a) >> 1;
|
|
return (__CZ__harmonicF(a, c, 1 / a) +
|
|
__CZ__harmonicF(c + 1, b, 1 / b));
|
|
}
|
|
return (1 / a + 1 / b);
|
|
}
|
|
|
|
define harmonic(limit)
|
|
{
|
|
if (!isint(limit)) {
|
|
return newerror("harmonic(limit): limit is not an integer");
|
|
}
|
|
if (limit <= 0) {
|
|
return newerror("harmonic(limit): limit is <=0");
|
|
}
|
|
/* The binary splitting algorithm returns 0 for f(1,1,0) */
|
|
if (limit == 1) {
|
|
return 1;
|
|
}
|
|
return __CZ__harmonicF(1, limit, 0);
|
|
}
|
|
|
|
/* lower incomplete gamma function */
|
|
|
|
/* lower
|
|
for z <= 1.1
|
|
*/
|
|
define __CZ__gammainc_series_lower(a, z)
|
|
{
|
|
local k ret tmp eps fact;
|
|
ret = 0;
|
|
k = 0;
|
|
tmp = 1;
|
|
fact = 1;
|
|
eps = epsilon();
|
|
while (abs(tmp - ret) > eps) {
|
|
tmp = ret;
|
|
ret += (z ^ (k + a)) / ((a + k) * fact);
|
|
k++;
|
|
fact *= -k;
|
|
}
|
|
return gamma(a) - ret;
|
|
}
|
|
|
|
/* lower
|
|
for z > 1.1
|
|
*/
|
|
define __CZ__gammainc_cf(a, z, n)
|
|
{
|
|
local ret k num1 denom1 num2 denom2;
|
|
ret = 0;
|
|
for (k = n + 1; k > 1; k--) {
|
|
ret = ((1 - k) * (k - 1 - a)) / (2 * k - 1 + z - a + ret);
|
|
}
|
|
return ((z ^ a * exp(-z)) / (1 + z - a + ret));
|
|
}
|
|
|
|
/* G(n,z) lower*/
|
|
define __CZ__gammainc_integer_a(a, z)
|
|
{
|
|
local k sum fact zz;
|
|
for (k = 0; k < a; k++) {
|
|
sum += (z ^ k) / (k !);
|
|
}
|
|
return (a - 1) ! *exp(-z) * sum;
|
|
}
|
|
|
|
/*
|
|
P = precision in dec digits
|
|
n = 1,2,3...
|
|
a,z => G(a,z)
|
|
*/
|
|
define __CZ__endcf(n, a, z, P)
|
|
{
|
|
local ret;
|
|
|
|
ret =
|
|
P * ln(10) + ln(4 * pi() * sqrt(n)) + re(z + (3 / 2 - a) * ln(z) -
|
|
lngamma(1 - a));
|
|
ret = ret / (sqrt(8 * (abs(z) + re(z))));
|
|
return ret ^ 2;
|
|
|
|
}
|
|
|
|
/* lower incomplete gamma function */
|
|
define gammainc(a, z)
|
|
{
|
|
local ret nterms eps epsilon tmp_before tmp_after n A B C sum k;
|
|
if (z == 0) {
|
|
return 1;
|
|
}
|
|
if (isint(a)) {
|
|
if (a > 0) {
|
|
if (a == 1) {
|
|
return exp(-z);
|
|
}
|
|
return __CZ__gammainc_integer_a(a, z);
|
|
} else {
|
|
if (a == 0) {
|
|
return -expoint(-z) + 1 / 2 * (ln(-z) - ln(-1 / z)) - ln(z);
|
|
} else if (a == -1) {
|
|
return expoint(-z) + 1 / 2 * (ln(-1 / z) - ln(-z)) + ln(z) +
|
|
(exp(-z) / z);
|
|
} else {
|
|
A = (-1) ^ ((-a) - 1) / ((-a) !);
|
|
B = expoint(-z) - 1 / 2 * (ln(-z) - ln(1 / -z)) + ln(z);
|
|
C = exp(-z);
|
|
sum = 0;
|
|
for (k = 1; k < -a; k++) {
|
|
sum += (z ^ (k - a - 1)) / (fallingfactorial(-a, k));
|
|
}
|
|
return A * B - C * sum;
|
|
}
|
|
}
|
|
}
|
|
if (re(z) <= 1.1 || re(z) < a + 1) {
|
|
eps = epsilon(epsilon() * 1e-10);
|
|
ret = __CZ__gammainc_series_lower(a, z);
|
|
epsilon(eps);
|
|
return ret;
|
|
} else {
|
|
eps = epsilon(epsilon() * 1e-10);
|
|
if (abs(exp(-z)) <= eps) {
|
|
return 0;
|
|
}
|
|
tmp_before = 0;
|
|
tmp_after = 1;
|
|
n = 1;
|
|
while (ceil(tmp_before) != ceil(tmp_after)) {
|
|
tmp_before = tmp_after;
|
|
tmp_after = __CZ__endcf(n++, a, z, digits(1 / eps));
|
|
/* a still quite arbitrary limit */
|
|
if (n > 10) {
|
|
return
|
|
newerror(strcat
|
|
("gammainc: evaluating limit for continued ",
|
|
"fraction does not converge"));
|
|
}
|
|
}
|
|
ret = __CZ__gammainc_cf(a, z, ceil(tmp_after));
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
}
|
|
|
|
|
|
define heavisidestep(x)
|
|
{
|
|
return (1 + sign(x)) / 2;
|
|
}
|
|
|
|
define NUMBER_POSITIVE_INFINITY()
|
|
{
|
|
return 1 / epsilon();
|
|
}
|
|
|
|
define NUMBER_NEGATIVE_INFINITY()
|
|
{
|
|
return -(1 / epsilon());
|
|
}
|
|
|
|
static true = 1;
|
|
static false = 0;
|
|
|
|
define g(prec)
|
|
{
|
|
local eps ret;
|
|
if (!isnull(prec)) {
|
|
eps = epsilon(prec);
|
|
ret = -psi(1);
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
return -psi(1);
|
|
}
|
|
|
|
define __CZ__series_converged(new, old, max)
|
|
{
|
|
local eps;
|
|
if (isnull(max)) {
|
|
eps = epsilon();
|
|
} else {
|
|
eps = max;
|
|
}
|
|
if (abs(re(new - old)) <= eps * abs(re(new))
|
|
&& abs(im(new - old)) <= eps * abs(im(new))) {
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
define __CZ__ei_power(z)
|
|
{
|
|
local tmp ei k old;
|
|
ei = g() + ln(abs(z)) + sgn(im(z)) * 1i * abs(arg(z));;
|
|
tmp = 1;
|
|
k = 1;
|
|
while (k) {
|
|
tmp *= z / k;
|
|
old = ei;
|
|
ei += tmp / k;
|
|
if (__CZ__series_converged(ei, old)) {
|
|
break;
|
|
}
|
|
k++;
|
|
}
|
|
return ei;
|
|
}
|
|
|
|
define __CZ__ei_asymp(z)
|
|
{
|
|
local ei old tmp k;
|
|
ei = sgn(im(z)) * 1i * pi();
|
|
tmp = exp(z) / z;
|
|
for (k = 1; k <= floor(abs(z)) + 1; k++) {
|
|
old = ei;
|
|
ei += tmp;
|
|
if (__CZ__series_converged(ei, old)) {
|
|
return ei;
|
|
}
|
|
tmp *= k / z;
|
|
}
|
|
return newerror("expoint: asymptotic series does not converge");
|
|
}
|
|
|
|
define __CZ__ei_cf(z)
|
|
{
|
|
local ei c d k old;
|
|
ei = sgn(im(z)) * 1i * pi();
|
|
if (ei != 0) {
|
|
c = 1 / ei;
|
|
d = 0;
|
|
c = 1 / (1 - z - exp(z) * c);
|
|
d = 1 / (1 - z - exp(z) * d);
|
|
ei *= d / c;
|
|
} else {
|
|
c = NUMBER_POSITIVE_INFINITY();
|
|
d = 0;
|
|
c = 0;
|
|
d = 1 / (1 - z - exp(z) * d);
|
|
ei = d * (-exp(z));
|
|
}
|
|
k = 1;
|
|
while (1) {
|
|
c = 1 / (2 * k + 1 - z - k * k * c);
|
|
d = 1 / (2 * k + 1 - z - k * k * d);
|
|
old = ei;
|
|
ei *= d / c;
|
|
if (__CZ__series_converged(ei, old)) {
|
|
break;
|
|
}
|
|
k++;
|
|
}
|
|
return ei;
|
|
}
|
|
|
|
define expoint(z)
|
|
{
|
|
local ei eps ret;
|
|
eps = epsilon(epsilon() * 1e-5);
|
|
if (abs(z) >= NUMBER_POSITIVE_INFINITY()) {
|
|
if (config("user_debug") > 0) {
|
|
print "expoint: abs(z) > +inf";
|
|
}
|
|
ret = sgn(im(z)) * 1i * pi() + exp(z) / z;
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
if (abs(z) > 2 - 1.035 * log(epsilon())) {
|
|
if (config("user_debug") > 0) {
|
|
print "expoint: using asymptotic series";
|
|
}
|
|
ei = __CZ__ei_asymp(z);
|
|
if (!iserror(ei)) {
|
|
ret = ei;
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
}
|
|
if (abs(z) > 1 && (re(z) < 0 || abs(im(z)) > 1)) {
|
|
if (config("user_debug") > 0) {
|
|
print "expoint: using continued fraction";
|
|
}
|
|
ret = __CZ__ei_cf(z);
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
if (abs(z) > 0) {
|
|
if (config("user_debug") > 0) {
|
|
print "expoint: using power series";
|
|
}
|
|
ret = __CZ__ei_power(z);
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
if (abs(z) == 0) {
|
|
if (config("user_debug") > 0) {
|
|
print "expoint: abs(z) = zero ";
|
|
}
|
|
epsilon(eps);
|
|
return NUMBER_NEGATIVE_INFINITY();
|
|
}
|
|
}
|
|
|
|
define erf(z)
|
|
{
|
|
return sqrt(z ^ 2) / z * (1 - 1 / sqrt(pi()) * gammainc(1 / 2, z ^ 2));
|
|
}
|
|
|
|
define erfc(z)
|
|
{
|
|
return 1 - erf(z);
|
|
}
|
|
|
|
define erfi(z)
|
|
{
|
|
return -1i * erf(1i * z);
|
|
}
|
|
|
|
define faddeeva(z)
|
|
{
|
|
return exp(-z ^ 2) * erfc(-1i * z);
|
|
}
|
|
|
|
define gammap(a, z)
|
|
{
|
|
return gammainc(a, z) / gamma(a);
|
|
}
|
|
|
|
define gammaq(a, z)
|
|
{
|
|
return 1 - gammap(a, z);
|
|
}
|
|
|
|
define lnbeta(a, b)
|
|
{
|
|
local ret eps;
|
|
eps = epsilon(epsilon() * 1e-3);
|
|
ret = (lngamma(a) + lngamma(b)) - lngamma(a + b);
|
|
epsilon(eps);
|
|
return ret;
|
|
}
|
|
|
|
define beta(a, b)
|
|
{
|
|
return exp(lnbeta(a, b));
|
|
}
|
|
|
|
|
|
define __CZ__ibetacf_a(a, b, z, n)
|
|
{
|
|
local A B m places;
|
|
if (n == 1) {
|
|
return 1;
|
|
}
|
|
m = n - 1;
|
|
places = highbit(1 + int (1 / epsilon())) +1;
|
|
A = bround((a + m - 1) * (a + b + m - 1) * m * (b - m) * z ^ 2, places++);
|
|
B = bround((a + 2 * (m) - 1) ^ 2, places++);
|
|
return A / B;
|
|
}
|
|
|
|
define __CZ__ibetacf_b(a, b, z, n)
|
|
{
|
|
local A B m places;
|
|
places = highbit(1 + int (1 / epsilon())) +1;
|
|
m = n - 1;
|
|
A = bround((m * (b - m) * z) / (a + 2 * m - 1), places++);
|
|
B = bround(((a + m) * (a - (a + b) * z + 1 + m * (2 - z))) / (a + 2 * m +
|
|
1), places++);
|
|
return m + A + B;
|
|
}
|
|
|
|
/* Didonato-Morris */
|
|
define __CZ__ibeta_cf_var_dm(a, b, z, max)
|
|
{
|
|
local m f c d check del h qab qam qap eps places;
|
|
|
|
eps = epsilon();
|
|
|
|
if (isnull(max)) {
|
|
max = 100;
|
|
}
|
|
places = highbit(1 + int (1 / epsilon())) + 1;
|
|
f = eps;
|
|
c = f;
|
|
d = 0;
|
|
for (m = 1; m <= max; m++) {
|
|
d = bround(__CZ__ibetacf_b(a, b, z, m) +
|
|
__CZ__ibetacf_a(a, b, z, m) * d, places++);
|
|
if (abs(d) < eps) {
|
|
d = eps;
|
|
}
|
|
c = bround(__CZ__ibetacf_b(a, b, z, m) +
|
|
__CZ__ibetacf_a(a, b, z, m) / c, places++);
|
|
if (abs(c) < eps) {
|
|
c = eps;
|
|
}
|
|
d = 1 / d;
|
|
check = c * d;
|
|
f = f * check;
|
|
if (abs(check - 1) < eps) {
|
|
break;
|
|
}
|
|
}
|
|
if (m > max) {
|
|
return newerror("ibeta: continuous fraction does not converge");
|
|
}
|
|
return f;
|
|
}
|
|
|
|
define betainc_complex(z, a, b)
|
|
{
|
|
local factor ret eps cf sum k N places tmp tmp2;
|
|
|
|
if (z == 0) {
|
|
if (re(a) > 0) {
|
|
return 0;
|
|
}
|
|
if (re(a) < 0) {
|
|
return newerror("betainc_complex: z == 0 and re(a) < 0");
|
|
}
|
|
}
|
|
if (z == 1) {
|
|
if (re(b) > 0) {
|
|
return 1;
|
|
} else {
|
|
return newerror("betainc_complex: z == 1 and re(b) < 0");
|
|
}
|
|
}
|
|
if (b <= 0) {
|
|
if (isint(b)) {
|
|
return 0;
|
|
} else {
|
|
return newerror("betainc_complex: b <= 0");
|
|
}
|
|
}
|
|
if (z == 1 / 2 && (a == b)) {
|
|
return 1 / 2;
|
|
}
|
|
if (isint(a) && isint(b)) {
|
|
eps = epsilon(epsilon() * 1e-10);
|
|
N = a + b - 1;
|
|
sum = 0;
|
|
for (k = a; k <= N; k++) {
|
|
tmp = ln(z) * k + ln(1 - z) * (N - k);
|
|
tmp2 = exp(ln(comb(N, k)) + tmp);
|
|
sum += tmp2;
|
|
}
|
|
epsilon(eps);
|
|
return sum;
|
|
} else if (re(z) <= re((a + 1) / (a + b + 2))) {
|
|
eps = epsilon(epsilon() * 1e-10);
|
|
places = highbit(1 + int (1 / epsilon())) + 1;
|
|
factor = bround((ln(z ^ a * (1 - z) ^ b) - lnbeta(a, b)), places);
|
|
cf = bround(__CZ__ibeta_cf_var_dm(a, b, z), places);
|
|
ret = factor + ln(cf);
|
|
if (abs(ret//ln(2)) >= places) {
|
|
ret = 0;
|
|
} else {
|
|
ret = bround(exp(factor + ln(cf)), places);
|
|
}
|
|
epsilon(eps);
|
|
return ret;
|
|
} else if (re(z) > re((a + 1) / (a + b + 2))
|
|
|| re(1 - z) < re((b + 1) / (a + b + 2))) {
|
|
ret = 1 - betainc_complex(1 - z, b, a);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
/******************************************************************************/
|
|
/*
|
|
Purpose:
|
|
|
|
__CZ__ibetaas63 computes the incomplete Beta function ratio.
|
|
|
|
Licensing:
|
|
|
|
This code is distributed under the GNU LGPL license.
|
|
|
|
Modified:
|
|
|
|
2013-08-03 20:52:05 +0000
|
|
|
|
Author:
|
|
|
|
Original FORTRAN77 version by KL Majumder, GP Bhattacharjee.
|
|
C version by John Burkardt
|
|
Calc version by Christoph Zurnieden
|
|
|
|
Reference:
|
|
|
|
KL Majumder, GP Bhattacharjee,
|
|
Algorithm AS 63:
|
|
The incomplete Beta Integral,
|
|
Applied Statistics,
|
|
Volume 22, Number 3, 1973, pages 409-411.
|
|
|
|
Parameters:
|
|
|
|
Input, x, the argument, between 0 and 1.
|
|
|
|
Input, a, b, the parameters, which
|
|
must be positive.
|
|
|
|
|
|
Output, the value of the incomplete
|
|
Beta function ratio.
|
|
*/
|
|
define __CZ__ibetaas63(x, a, b, beta)
|
|
{
|
|
local ai betain cx indx ns aa asb bb rx temp term value xx acu places;
|
|
acu = epsilon();
|
|
|
|
value = x;
|
|
/* inverse incbeta calculates it already */
|
|
if (isnull(beta))
|
|
beta = lnbeta(a, b);
|
|
|
|
if (a <= 0.0 || b <= 0.0) {
|
|
return newerror("betainc: domain error: a < 0 and/or b < 0");
|
|
}
|
|
if (x < 0.0 || 1.0 < x) {
|
|
return newerror("betainc: domain error: x<0 or x>1");
|
|
}
|
|
if (x == 0.0 || x == 1.0) {
|
|
return value;
|
|
}
|
|
asb = a + b;
|
|
cx = 1.0 - x;
|
|
|
|
if (a < asb * x) {
|
|
xx = cx;
|
|
cx = x;
|
|
aa = b;
|
|
bb = a;
|
|
indx = 1;
|
|
} else {
|
|
xx = x;
|
|
aa = a;
|
|
bb = b;
|
|
indx = 0;
|
|
}
|
|
|
|
term = 1.0;
|
|
ai = 1.0;
|
|
value = 1.0;
|
|
ns = floor(bb + cx * asb);
|
|
|
|
rx = xx / cx;
|
|
temp = bb - ai;
|
|
if (ns == 0) {
|
|
rx = xx;
|
|
}
|
|
places = highbit(1 + int (1 / acu)) + 1;
|
|
while (1) {
|
|
term = bround(term * temp * rx / (aa + ai), places++);
|
|
value = value + term;;
|
|
temp = abs(term);
|
|
|
|
if (temp <= acu && temp <= abs(acu * value)) {
|
|
value = value * exp(aa * ln(xx)
|
|
+ (bb - 1.0) * ln(cx) - beta) / aa;
|
|
|
|
if (indx) {
|
|
value = 1.0 - value;
|
|
}
|
|
break;
|
|
}
|
|
|
|
ai = ai + 1.0;
|
|
ns = ns - 1;
|
|
|
|
if (0 <= ns) {
|
|
temp = bb - ai;
|
|
if (ns == 0) {
|
|
rx = xx;
|
|
}
|
|
} else {
|
|
temp = asb;
|
|
asb = asb + 1.0;
|
|
}
|
|
}
|
|
epsilon(acu);
|
|
return value;
|
|
}
|
|
|
|
/*
|
|
z
|
|
/
|
|
[ b - 1 a - 1
|
|
1/beta(a,b) * I (1 - t) t dt
|
|
]
|
|
/
|
|
0
|
|
|
|
*/
|
|
|
|
define betainc(z, a, b)
|
|
{
|
|
local factor ret eps cf sum k N places tmp tmp2;
|
|
|
|
if (im(z) || im(a) || im(b))
|
|
return betainc_complex(z, a, b);
|
|
|
|
if (z == 0) {
|
|
if (re(a) > 0) {
|
|
return 0;
|
|
}
|
|
if (re(a) < 0) {
|
|
return newerror("betainc: z == 0 and re(a) < 0");
|
|
}
|
|
}
|
|
if (z == 1) {
|
|
if (re(b) > 0) {
|
|
return 1;
|
|
} else {
|
|
return newerror("betainc: z == 1 and re(b) < 0");
|
|
}
|
|
}
|
|
if (b <= 0) {
|
|
if (isint(b)) {
|
|
return 0;
|
|
} else {
|
|
return newerror("betainc: b <= 0");
|
|
}
|
|
}
|
|
if (z == 1 / 2 && a == b) {
|
|
return 1 / 2;
|
|
}
|
|
return __CZ__ibetaas63(z, a, b);
|
|
|
|
}
|
|
|
|
define __CZ__erfinvapprox(x)
|
|
{
|
|
local a;
|
|
a = 0.147;
|
|
return sgn(x) *
|
|
sqrt(sqrt
|
|
((2 / (pi() * a) + (ln(1 - x ^ 2)) / 2) ^ 2 - (ln(1 - x ^ 2)) / a)
|
|
- (2 / (pi() * a) + (ln(1 - x ^ 2)) / 2));
|
|
}
|
|
|
|
/* complementary inverse error function, faster at about x < 1-.91
|
|
Henry E. Fettis. "A stable algorithm for computing the inverse error function
|
|
in the 'tail-end' region" Math. Comp., 28:585-587, 1974.
|
|
*/
|
|
define __CZ__inverffettis(x, n)
|
|
{
|
|
local y sqrtpi oldy k places;
|
|
if (isnull(n))
|
|
n = 205;
|
|
y = erfinvapprox(1 - x);
|
|
places = highbit(1 + int (1 / epsilon())) + 1;
|
|
sqrtpi = sqrt(pi());
|
|
do {
|
|
oldy = y;
|
|
k++;
|
|
y = bround((ln(__CZ__fettiscf(y, n) / (sqrtpi * x))) ^ (1 / 2), places);
|
|
} while (abs(y - oldy) / y > epsilon());
|
|
return y;
|
|
}
|
|
|
|
/* cf for erfc() */
|
|
define __CZ__fettiscf(y, n)
|
|
{
|
|
local k t tt r a b;
|
|
t = 1 / y;
|
|
tt = t ^ 2 / 2;
|
|
for (k = n; k > 0; k--) {
|
|
a = 1;
|
|
b = k * tt;
|
|
r = b / (a + r);
|
|
}
|
|
return t / (1 + r);
|
|
}
|
|
|
|
/* inverse error function, faster at about x<=.91*/
|
|
define __CZ__inverfbin(x)
|
|
{
|
|
local places approx flow fhigh eps high low mid fmid epsilon;
|
|
approx = erfinvapprox(x);
|
|
epsilon = epsilon();
|
|
high = approx + 1e-4;
|
|
low = -1;
|
|
places = highbit(1 + int (1 / epsilon)) +1;
|
|
fhigh = x - erf(high);
|
|
flow = x - erf(low);
|
|
while (1) {
|
|
mid = bround(high - fhigh * (high - low) / (fhigh - flow), places);
|
|
if ((mid == low) || (mid == high)) {
|
|
places++;
|
|
}
|
|
fmid = x - erf(mid);
|
|
if (abs(fmid) < epsilon) {
|
|
return mid;
|
|
}
|
|
if (sgn(fmid) == sgn(flow)) {
|
|
low = mid;
|
|
flow = fmid;
|
|
} else {
|
|
high = mid;
|
|
fhigh = fmid;
|
|
}
|
|
}
|
|
}
|
|
|
|
define erfinv(x)
|
|
{
|
|
local ret approx a eps y old places errfunc sqrtpihalf flag k;
|
|
if (x < -1 || x > 1)
|
|
return newerror("erfinv: input out of domain (-1<=x<=1)");
|
|
if (x == 0)
|
|
return 0;
|
|
if (x == -1)
|
|
return NUMBER_NEGATIVE_INFINITY();
|
|
if (x == +1)
|
|
return NUMBER_POSITIVE_INFINITY();
|
|
|
|
if (x < 0) {
|
|
x = -x;
|
|
flag = 1;
|
|
}
|
|
/* No need for full precision */
|
|
eps = epsilon(1e-20);
|
|
if (eps >= 1e-40) {
|
|
/* Winitzki, Sergei (6 February 2008). "A handy approximation for the
|
|
* error function and its inverse" */
|
|
a = 0.147;
|
|
y = sgn(x) * sqrt(sqrt((2 / (pi() * a)
|
|
+ (ln(1 - x ^ 2)) / 2) ^ 2
|
|
- (ln(1 - x ^ 2)) / a)
|
|
- (2 / (pi() * a) + (ln(1 - x ^ 2)) / 2));
|
|
|
|
} else {
|
|
/* 20 digits instead of 5 */
|
|
if (x <= .91) {
|
|
y = __CZ__inverfbin(x);
|
|
} else {
|
|
y = __CZ__inverffettis(1 - x);
|
|
}
|
|
|
|
if (eps <= 1e-20) {
|
|
epsilon(eps);
|
|
return y;
|
|
}
|
|
}
|
|
epsilon(eps);
|
|
/* binary digits in number (here: number = epsilon()) */
|
|
places = highbit(1 + int (1 / eps)) +1;
|
|
sqrtpihalf = 2 / sqrt(pi());
|
|
k = 0;
|
|
/*
|
|
* Do some Newton-Raphson steps to reach final accuracy.
|
|
* Only a couple of steps are necessary but calculating the error function
|
|
* at higher precision is quite costly;
|
|
*/
|
|
do {
|
|
old = y;
|
|
errfunc = bround(erf(y), places);
|
|
if (abs(errfunc - x) <= eps) {
|
|
break;
|
|
}
|
|
y = bround(y - (errfunc - x) / (sqrtpihalf * exp(-y ^ 2)), places);
|
|
k++;
|
|
} while (1);
|
|
/*
|
|
* This is not really necessary but e.g:
|
|
* ; epsilon(1e-50)
|
|
* 0.00000000000000000000000000000000000000000000000001
|
|
* ; erfinv(.9999999999999999999999999999999)
|
|
* 8.28769266865549025938
|
|
* ; erfinv(.999999999999999999999999999999)
|
|
* 8.14861622316986460738453487549552168842204512959346
|
|
* ; erf(8.28769266865549025938)
|
|
* 0.99999999999999999999999999999990000000000000000000
|
|
* ; erf(8.14861622316986460738453487549552168842204512959346)
|
|
* 0.99999999999999999999999999999900000000000000000000
|
|
* The precision "looks too short".
|
|
*/
|
|
if (k == 0) {
|
|
y = bround(y - (errfunc - x) / (sqrtpihalf * exp(-y ^ 2)), places);
|
|
}
|
|
if (flag == 1) {
|
|
y = -y;
|
|
}
|
|
return y;
|
|
}
|
|
|
|
|
|
/*
|
|
* restore internal function from resource debugging
|
|
*/
|
|
config("resource_debug", resource_debug_level),;
|
|
if (config("resource_debug") & 3) {
|
|
print "zeta(z)";
|
|
print "psi(z)";
|
|
print "polygamma(m,z)";
|
|
print "lngamma(z)";
|
|
print "gamma(z)";
|
|
print "harmonic(limit)";
|
|
print "gammainc(a,z)";
|
|
print "heavisidestep(x)";
|
|
print "expoint(z)";
|
|
print "erf(z)";
|
|
print "erfinv(x)";
|
|
print "erfc(z)";
|
|
print "erfi(z)";
|
|
print "erfinv(x)";
|
|
print "faddeeva(z)";
|
|
print "gammap(a,z)";
|
|
print "gammaq(a,z)";
|
|
print "beta(a,b)";
|
|
print "lnbeta(a,b)";
|
|
print "betainc(z,a,b)";
|
|
}
|