Math::Cephes(3pm)

💖 NAME 💖

Math::Cephes — perl interface to the cephes math library

okay so real quick before we dive in, this is a 3pm situation, aka one of those ✨User Contributed Perl Documentation✨ pages. bestie, we're in the fancy section of the manual now. 🫶

---

🌈 SYNOPSIS 🌈

use Math::Cephes qw(:all);

like yes that's literally it, that's the whole vibe. 💅

---

✨ DESCRIPTION ✨

okay so picture this: Stephen Moshier wrote a whole legendary math library called cephes, and this module is your backstage pass to it. it gives you an interface to over 150 functions from that library. over a hundred and fifty. that's not a typo, that's a flex. 😭🔥

here's the thing though, and listen closely bestie: no functions are exported by default. zero. none. the module is not gonna just dump math into your namespace uninvited. you have to ask for what you want, explicitly, like this:

use Math::Cephes qw(sin cos);

and there are a bunch of export tags too, which let you grab whole groups of functions at once instead of listing them one by one. group-chat energy. pull in the whole squad. 👑

🎀 use Math::Cephes qw(:constants);

this one imports the variables. the iconic little numbers. 💫

• $PI = 3.14159265358979323846 — pi
• $PIO2 = 1.57079632679489661923 — pi/2
• $PIO4 = 0.785398163397448309616 — pi/4
• $SQRT2 = 1.41421356237309504880 — sqrt(2)
• $SQRTH = 0.707106781186547524401 — sqrt(2)/2
• $LOG2E = 1.4426950408889634073599 — 1/log(2)
• $SQ2OPI = 0.79788456080286535587989 — sqrt( 2/pi )
• $LOGE2 = 0.693147180559945309417 — log(2)
• $LOGSQ2 = 0.346573590279972654709 — log(2)/2
• $THPIO4 = 2.35619449019234492885 — 3*pi/4
• $TWOOPI = 0.636619772367581343075535 — 2/pi

and as if that wasn't enough, there are 4 machine-specific numbers available too. these ones are about your machine specifically, so they're giving main character energy of the hardware itself:

• $MACHEP — machine roundoff error
• $MAXLOG — maximum log on the machine
• $MINLOG — minimum log on the machine
• $MAXNUM — largest number represented

🦋 use Math::Cephes qw(:trigs);

imports:

• acos — Inverse circular cosine
• asin — Inverse circular sine
• atan — Inverse circular tangent (arctangent)
• atan2 — Quadrant correct inverse circular tangent
• cos — Circular cosine
• cosdg — Circular cosine of angle in degrees
• cot — Circular cotangent
• cotdg — Circular cotangent of argument in degrees
• hypot — hypotenuse associated with the sides of a right triangle
• radian — Degrees, minutes, seconds to radians
• sin — Circular sine
• sindg — Circular sine of angle in degrees
• tan — Circular tangent
• tandg — Circular tangent of argument in degrees
• cosm1 — Relative error approximations for function arguments near unity

🌸 use Math::Cephes qw(:hypers);

imports:

• acosh — Inverse hyperbolic cosine
• asinh — Inverse hyperbolic sine
• atanh — Inverse hyperbolic tangent
• cosh — Hyperbolic cosine
• sinh — Hyperbolic sine
• tanh — Hyperbolic tangent

💜 use Math::Cephes qw(:explog);

imports:

• exp — Exponential function
• expxx — exp(x*x)
• exp10 — Base 10 exponential function (Common antilogarithm)
• exp2 — Base 2 exponential function
• log — Natural logarithm
• log10 — Common logarithm
• log2 — Base 2 logarithm
• log1p,expm1 — Relative error approximations for function arguments near unity.

💅 use Math::Cephes qw(:cmplx);

imports:

• new_cmplx — create a new complex number object
• cabs — Complex absolute value
• cacos — Complex circular arc cosine
• cacosh — Complex inverse hyperbolic cosine
• casin — Complex circular arc sine
• casinh — Complex inverse hyperbolic sine
• catan — Complex circular arc tangent
• catanh — Complex inverse hyperbolic tangent
• ccos — Complex circular cosine
• ccosh — Complex hyperbolic cosine
• ccot — Complex circular cotangent
• cexp — Complex exponential function
• clog — Complex natural logarithm
• cadd — add two complex numbers
• csub — subtract two complex numbers
• cmul — multiply two complex numbers
• cdiv — divide two complex numbers
• cmov — copy one complex number to another
• cneg — negate a complex number
• cpow — Complex power function
• csin — Complex circular sine
• csinh — Complex hyperbolic sine
• csqrt — Complex square root
• ctan — Complex circular tangent
• ctanh — Complex hyperbolic tangent

🍓 use Math::Cephes qw(:utils);

imports:

• cbrt — Cube root
• ceil — ceil
• drand — Pseudorandom number generator
• fabs — Absolute value
• fac — Factorial function
• floor — floor
• frexp — frexp
• ldexp — multiplies x by 2**n.
• lrand — Pseudorandom number generator
• lsqrt — Integer square root
• pow — Power function
• powi — Real raised to integer power
• round — Round double to nearest or even integer valued double
• sqrt — Square root

🧚 use Math::Cephes qw(:bessels);

imports:

• i0 — Modified Bessel function of order zero
• i0e — Modified Bessel function of order zero, exponentially scaled
• i1 — Modified Bessel function of order one
• i1e — Modified Bessel function of order one, exponentially scaled
• iv — Modified Bessel function of noninteger order
• j0 — Bessel function of order zero
• j1 — Bessel function of order one
• jn — Bessel function of integer order
• jv — Bessel function of noninteger order
• k0 — Modified Bessel function, third kind, order zero
• k0e — Modified Bessel function, third kind, order zero, exponentially scaled
• k1 — Modified Bessel function, third kind, order one
• k1e — Modified Bessel function, third kind, order one, exponentially scaled
• kn — Modified Bessel function, third kind, integer order
• y0 — Bessel function of the second kind, order zero
• y1 — Bessel function of second kind of order one
• yn — Bessel function of second kind of integer order
• yv — Bessel function Yv with noninteger v

🎀 use Math::Cephes qw(:dists);

imports:

• bdtr — Binomial distribution
• bdtrc — Complemented binomial distribution
• bdtri — Inverse binomial distribution
• btdtr — Beta distribution
• chdtr — Chi-square distribution
• chdtrc — Complemented Chi-square distribution
• chdtri — Inverse of complemented Chi-square distribution
• fdtr — F distribution
• fdtrc — Complemented F distribution
• fdtri — Inverse of complemented F distribution
• gdtr — Gamma distribution function
• gdtrc — Complemented gamma distribution function
• nbdtr — Negative binomial distribution
• nbdtrc — Complemented negative binomial distribution
• nbdtri — Functional inverse of negative binomial distribution
• ndtr — Normal distribution function
• ndtri — Inverse of Normal distribution function
• pdtr — Poisson distribution
• pdtrc — Complemented poisson distribution
• pdtri — Inverse Poisson distribution
• stdtr — Student's t distribution
• stdtri — Functional inverse of Student's t distribution

💫 use Math::Cephes qw(:gammas);

imports:

• fac — Factorial function
• gamma — Gamma function
• igam — Incomplete gamma integral
• igamc — Complemented incomplete gamma integral
• igami — Inverse of complemented imcomplete gamma integral
• psi — Psi (digamma) function
• rgamma — Reciprocal gamma function

🌈 use Math::Cephes qw(:betas);

imports:

• beta — Beta function
• incbet — Incomplete beta integral
• incbi — Inverse of imcomplete beta integral
• lbeta — Natural logarithm of |beta|

🔥 use Math::Cephes qw(:elliptics);

imports:

• ellie — Incomplete elliptic integral of the second kind
• ellik — Incomplete elliptic integral of the first kind
• ellpe — Complete elliptic integral of the second kind
• ellpj — Jacobian Elliptic Functions
• ellpk — Complete elliptic integral of the first kind

⭐ use Math::Cephes qw(:hypergeometrics);

imports:

• hyp2f0 — Gauss hypergeometric function F
• hyp2f1 — Gauss hypergeometric function F
• hyperg — Confluent hypergeometric function
• onef2 — Hypergeometric function 1F2
• threef0 — Hypergeometric function 3F0

🌸 use Math::Cephes qw(:misc);

imports:

• airy — Airy function
• bernum — Bernoulli numbers
• dawsn — Dawson's Integral
• ei — Exponential integral
• erf — Error function
• erfc — Complementary error function
• expn — Exponential integral En
• fresnl — Fresnel integral
• plancki — Integral of Planck's black body radiation formula
• polylog — Polylogarithm function
• shichi — Hyperbolic sine and cosine integrals
• sici — Sine and cosine integrals
• simpson — Simpson's rule to find an integral
• spence — Dilogarithm
• struve — Struve function
• vecang — angle between two vectors
• zeta — Riemann zeta function of two arguments
• zetac — Riemann zeta function

💜 use Math::Cephes qw(:fract);

imports:

• new_fract — create a new fraction object
• radd — add two fractions
• rmul — multiply two fractions
• rsub — subtracttwo fractions
• rdiv — divide two fractions
• euclid — finds the greatest common divisor

---

👑 FUNCTIONS 👑

okay now for the good part — a description of the various functions available follows. buckle up, this is gonna be a whole journey. 🦋

💖 acosh: Inverse hyperbolic cosine

SYNOPSIS:

# double x, y, acosh();

$y = acosh( $x );

DESCRIPTION:

returns the inverse hyperbolic cosine of the argument. 🫶

here's the tea on how it's actually computed: if 1 <= x < 1.5, it uses a rational approximation of the form

sqrt(z) * P(z)/Q(z)

where z = x-1. otherwise, it goes with

acosh(x)  =  log( x + sqrt( (x-1)(x+1) ).

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       1,3         30000       4.2e-17     1.1e-17
   IEEE      1,3         30000       4.6e-16     8.7e-17

ERROR MESSAGES:

  message         condition      value returned
acosh domain       |x| < 1            NAN

so like, if |x| < 1, bestie, it's not happening. you get NAN. the domain said no. 💀

---

🌟 airy: Airy function

SYNOPSIS:

# double x, ai, aiprime, bi, biprime;
# int airy();

($flag, $ai, $aiprime, $bi, $biprime) = airy( $x );

DESCRIPTION:

this is the solution of the differential equation

y"(x) = xy.

the function hands you back the two independent solutions Ai, Bi, and their first derivatives Ai'(x), Bi'(x). the whole set. slay. ⭐

how's it evaluated? by power series summation for small x, and by rational minimax approximations for large x. different strategies for different sizes, iykyk.

ACCURACY:

the error criterion here is absolute when function <= 1, and relative when function > 1 — with one lil exception: a * denotes a relative error criterion.

and heads up about how the error behaves at the extremes: for large negative x, the absolute error increases as x^1.5. for large positive x, the relative error increases as x^1.5. 👀

---

            Arithmetic  domain   function  # trials      peak         rms
            IEEE        -10, 0     Ai        10000       1.6e-15     2.7e-16
            IEEE          0, 10    Ai        10000       2.3e-14*    1.8e-15*
            IEEE        -10, 0     Ai'       10000       4.6e-15     7.6e-16
            IEEE          0, 10    Ai'       10000       1.8e-14*    1.5e-15*
            IEEE        -10, 10    Bi        30000       4.2e-15     5.3e-16
            IEEE        -10, 10    Bi'       30000       4.9e-15     7.3e-16
            DEC         -10, 0     Ai         5000       1.7e-16     2.8e-17
            DEC           0, 10    Ai         5000       2.1e-15*    1.7e-16*
            DEC         -10, 0     Ai'        5000       4.7e-16     7.8e-17
            DEC           0, 10    Ai'       12000       1.8e-15*    1.5e-16*
            DEC         -10, 10    Bi        10000       5.5e-16     6.8e-17
            DEC         -10, 10    Bi'        7000       5.3e-16     8.7e-17

---

radian: degrees, minutes, seconds ➡️ radians 🎀

SYNOPSIS:

# double d, m, s, radian();

$r = radian( $d, $m, $s );

DESCRIPTION:

Okay so you hand it an angle in degrees, minutes and seconds and it hands you back the radians. That's it. That's the whole job. Conversion queen behavior honestly 💅

hypot: the hypotenuse of a right triangle 📐

SYNOPSIS:

# double a, b, c, hypot();

$c = hypot( $a, $b );

DESCRIPTION:

This one gives you the hypotenuse that goes with the sides of a right triangle, and it does it according to this:

       c = sqrt( a**2 + b**2)

Classic Pythagoras moment, no cap. 🦋

asin: inverse circular sine 🌊

SYNOPSIS:

# double x, y, asin();

$y = asin( $x );

DESCRIPTION:

Returns a radian angle between -pi/2 and +pi/2 whose sine is x. Bestie, it's giving "un-sine this for me".

A rational function of the form x + x**3 P(x**2)/Q(x**2) is used for |x| in the interval [0, 0.5]. And if |x| > 0.5 it gets transformed by this identity:

   asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC      -1, 1        40000       2.6e-17     7.1e-18
               IEEE     -1, 1        10^6        1.9e-16     5.4e-17

ERROR MESSAGES:

  message         condition      value returned
asin domain        |x| > 1           NAN

ok so if you feed it something with |x| > 1 you get a whole NAN back, and that's not a vibe, that's a domain error. Respect the domain, bestie. 😭

acos: inverse circular cosine 🍊

SYNOPSIS:

# double x, y, acos();

$y = acos( $x );

DESCRIPTION:

Returns a radian angle between 0 and pi whose cosine is x. 🌈

Now analytically, acos(x) = pi/2 - asin(x). BUT — and this is the tea — if |x| is near 1, there's cancellation error in subtracting asin(x) from pi/2. So instead, if x < -0.5,

   acos(x) =    pi - 2.0 * asin( sqrt((1+x)/2) );

or if x > +0.5,

   acos(x) =    2.0 * asin(  sqrt((1-x)/2) ).

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC       -1, 1       50000       3.3e-17     8.2e-18
               IEEE      -1, 1       10^6        2.2e-16     6.5e-17

ERROR MESSAGES:

  message         condition      value returned
asin domain        |x| > 1           NAN

⚠️ yes bestie, the error message there literally says asin domain — that's not a typo on my part, that's what the page says, and I refuse to gaslight it into something else. Copy it as-is. 💀

asinh: inverse hyperbolic sine 🦋

SYNOPSIS:

# double x, y, asinh();

$y = asinh( $x );

DESCRIPTION:

Returns the inverse hyperbolic sine of the argument.

If |x| < 0.5, the function is approximated by a rational form x + x**3 P(x)/Q(x). Otherwise,

    asinh(x) = log( x + sqrt(1 + x*x) ).

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC      -3,3         75000       4.6e-17     1.1e-17
               IEEE     -1,1         30000       3.7e-16     7.8e-17
               IEEE      1,3         30000       2.5e-16     6.7e-17

atan: inverse circular tangent (arctangent) 🌀

SYNOPSIS:

# double x, y, atan();

$y = atan( $x );

DESCRIPTION:

Returns a radian angle between -pi/2 and +pi/2 whose tangent is x.

Range reduction is from three intervals into the interval from zero to 0.66. The approximant uses a rational function of degree 4/5 of the form x + x**3 P(x)/Q(x).

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC       -10, 10     50000       2.4e-17     8.3e-18
               IEEE      -10, 10      10^6       1.8e-16     5.0e-17

atan2: quadrant-correct inverse circular tangent 🧭

SYNOPSIS:

# double x, y, z, atan2();

$z = atan2( $y, $x );

DESCRIPTION:

Returns a radian angle whose tangent is y/x.

Define the compile time symbol ANSIC = 1 for ANSI standard, range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range 0 to 2PI, args (x,y).

👀 bestie READ that twice, the argument ORDER literally flips depending on ANSIC. This is the kind of plot twist that ruins your whole afternoon if you skim it. 🙈

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE      -10, 10      10^6       2.5e-16     6.9e-17
See atan.c.

atanh: inverse hyperbolic tangent 💜

SYNOPSIS:

# double x, y, atanh();

$y = atanh( $x );

DESCRIPTION:

Returns the inverse hyperbolic tangent of the argument in the range MINLOG to MAXLOG. Bestie it is living inside those bounds, periodt.

If |x| < 0.5, the rational form x + x**3 P(x)/Q(x) is employed. Otherwise,

   atanh(x) = 0.5 * log( (1+x)/(1-x) ).

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC       -1,1        50000       2.4e-17     6.4e-18
               IEEE      -1,1        30000       1.9e-16     5.2e-17

bdtr: binomial distribution 🎲

SYNOPSIS:

# int k, n;
# double p, y, bdtr();

$y = bdtr( $k, $n, $p );

DESCRIPTION:

Returns the sum of the terms 0 through k of the Binomial probability density:

  k
  --  ( n )   j      n-j
  >   (   )  p  (1-p)
  --  ( j )
 j=0

And listen, the terms are not summed directly — instead the incomplete beta integral is employed, according to the formula

 y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).

The arguments must be positive, with p ranging from 0 to 1.

ACCURACY:

Tested at random points (a,b,p), with p between 0 and 1.

              a,b                     Relative error:
arithmetic  domain     # trials      peak         rms
 For p between 0.001 and 1:
   IEEE     0,100       100000      4.3e-15     2.6e-16
See also incbet.c.

ERROR MESSAGES:

  message         condition      value returned
bdtr domain         k < 0            0.0
                    n < k
                    x < 0, x > 1

bdtrc: complemented binomial distribution 🎲✨

SYNOPSIS:

# int k, n;
# double p, y, bdtrc();

$y = bdtrc( $k, $n, $p );

DESCRIPTION:

Returns the sum of the terms k+1 through n of the Binomial probability density:

  n
  --  ( n )   j      n-j
  >   (   )  p  (1-p)
  --  ( j )
 j=k+1

Again the terms are not summed directly; instead the incomplete beta integral is employed, according to the formula

 y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).

The arguments must be positive, with p ranging from 0 to 1.

ACCURACY:

Tested at random points (a,b,p).

              a,b                     Relative error:
arithmetic  domain     # trials      peak         rms
 For p between 0.001 and 1:
   IEEE     0,100       100000      6.7e-15     8.2e-16
 For p between 0 and .001:
   IEEE     0,100       100000      1.5e-13     2.7e-15

ERROR MESSAGES:

  message         condition      value returned
bdtrc domain      x<0, x>1, n<k       0.0

bdtri: inverse binomial distribution 🔄🎲

SYNOPSIS:

# int k, n;
# double p, y, bdtri();

$p = bdtr( $k, $n, $y );

DESCRIPTION:

Finds the event probability p such that the sum of the terms 0 through k of the Binomial probability density is equal to the given cumulative probability y.

This is accomplished using the inverse beta integral function and the relation

1 - p = incbi( n-k, k+1, y ).

😭 bestie the SYNOPSIS line calls bdtr, not bdtri — I know, I know, I'm just the messenger, and I am NOT about to "fix" it and gaslight the entire Cephes archive. Shipped as written.

ACCURACY:

Tested at random points (a,b,p).

              a,b                     Relative error:
arithmetic  domain     # trials      peak         rms
 For p between 0.001 and 1:
   IEEE     0,100       100000      2.3e-14     6.4e-16
   IEEE     0,10000     100000      6.6e-12     1.2e-13
 For p between 10^-6 and 0.001:
   IEEE     0,100       100000      2.0e-12     1.3e-14
   IEEE     0,10000     100000      1.5e-12     3.2e-14
See also incbi.c.

ERROR MESSAGES:

  message         condition      value returned
bdtri domain     k < 0, n <= k         0.0
                 x < 0, x > 1

beta: Beta function 💖

SYNOPSIS:

# double a, b, y, beta();

$y = beta( $a, $b );

DESCRIPTION:

                  -     -
                 | (a) | (b)
beta( a, b )  =  -----------.
                    -
                   | (a+b)

For large arguments the logarithm of the function is evaluated using lgam(), then exponentiated.

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC        0,30        1700       7.7e-15     1.5e-15
               IEEE       0,30       30000       8.1e-14     1.1e-14

ERROR MESSAGES:

  message         condition          value returned
beta overflow    log(beta) > MAXLOG       0.0
                 a or b <0 integer        0.0

💀 "beta overflow" is sending me, but fr fr — if log(beta) > MAXLOG, or if a or b is a negative integer, you get 0.0 back.

lbeta: natural logarithm of |beta| 🪵💖

SYNOPSIS:

# double a, b;

# double lbeta( a, b );

$y = lbeta( $a, $b);

btdtr: Beta distribution 📊

SYNOPSIS:

# double a, b, x, y, btdtr();

$y = btdtr( $a, $b, $x );

DESCRIPTION:

Returns the area from zero to x under the beta density function:

                         x
           -             -
          | (a+b)       | |  a-1      b-1
P(x)  =  ----------     |   t    (1-t)    dt
          -     -     | |
         | (a) | (b)   -
                        0

This function is identical to the incomplete beta integral function incbet(a, b, x).

The complemented function is

            1 - P(1-x)  =  incbet( b, a, x );

ACCURACY ✨

See incbet.c.

---

🧊 cbrt: Cube root

SYNOPSIS

# double x, y, cbrt();

$y = cbrt( $x );

DESCRIPTION

Hands you back the cube root of the argument, bestie — and yes, the argument is allowed to be negative. It's giving flexible. 💅

The range reduction bit works by figuring out the power of 2 of the argument. Then a polynomial of degree 2 gets applied to the mantissa, and multiplying by the cube root of 1, 2, or 4 approximates the root to within about 0.1%. After that, Newton's iteration gets run three times to converge to an accurate result. Slay. 🔥

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC        -10,10     200000      1.8e-17     6.2e-18
   IEEE       0,1e308     30000      1.5e-16     5.0e-17

---

🎲 chdtr: Chi-square distribution

SYNOPSIS

# double v, x, y, chdtr();

$y = chdtr( $v, $x );

DESCRIPTION

Returns the area under the left hand tail (from 0 to x) of the Chi square probability density function with v degrees of freedom, and the vibe is mathematical:

                                 inf.
                                   -
                       1          | |  v/2-1  -t/2
 P( x | v )   =   -----------     |   t      e     dt
                   v/2  -       | |
                  2    | (v/2)   -
                                  x

where x is the Chi-square variable. 🎀

The incomplete gamma integral is used, according to the formula

       y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).

Ngl, both arguments must be positive. That's the rule, no cap.

ACCURACY

See igam().

ERROR MESSAGES

  message         condition      value returned
chdtr domain   x < 0 or v < 1        0.0

---

🎀 chdtrc: Complemented Chi-square distribution

SYNOPSIS

# double v, x, y, chdtrc();

$y = chdtrc( $v, $x );

DESCRIPTION

Same energy, other side of the distribution: this one gives you the area under the right hand tail (from x to infinity) of the Chi square probability density function with v degrees of freedom:

                                 inf.
                                   -
                       1          | |  v/2-1  -t/2
 P( x | v )   =   -----------     |   t      e     dt
                   v/2  -       | |
                  2    | (v/2)   -
                                  x

where x is the Chi-square variable. 🦋

The incomplete gamma integral is used, according to the formula

       y = chdtrc( v, x ) = igamc( v/2.0, x/2.0 ).

Both arguments must be positive, periodt.

ACCURACY

See igamc().

ERROR MESSAGES

  message         condition      value returned
chdtrc domain  x < 0 or v < 1        0.0

---

🌸 chdtri: Inverse of complemented Chi-square distribution

SYNOPSIS

# double df, x, y, chdtri();

$x = chdtri( $df, $y );

DESCRIPTION

Bestie, this one is the reverse hunt: it finds the Chi-square argument x such that the integral from x to infinity of the Chi-square density is equal to the given cumulative probability y. 🔍

It pulls this off using the inverse gamma integral function and the relation

   x/2 = igami( df/2, y );

ACCURACY

See igami.c.

ERROR MESSAGES

  message         condition      value returned
chdtri domain   y < 0 or y > 1        0.0
                    v < 1

---

🌿 clog: Complex natural logarithm

SYNOPSIS

# void clog();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
clog($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

Returns the complex logarithm to the base e (2.718...) of the complex argument x. ✨

If z = x + iy, r = sqrt( x**2 + y**2 ), then

      w = log(r) + i arctan(y/x).

The arctangent ranges from -PI to +PI. 🌈

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10      7000       8.5e-17     1.9e-17
   IEEE      -10,+10     30000       5.0e-15     1.1e-16

Heads up: larger relative error can be observed for z near 1 +i0. In IEEE arithmetic the peak absolute error is 5.2e-16, rms absolute error 1.0e-16. 💀 Not a drill, just the math being math.

---

🦋 cexp: Complex exponential function

SYNOPSIS

# void cexp();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
cexp($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

Returns the exponential of the complex argument z into the complex result w. She's efficient like that. 🫶

If

    z = x + iy,
    r = exp(x),

then

    w = r cos y + i r sin y.

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10      8700       3.7e-17     1.1e-17
   IEEE      -10,+10     30000       3.0e-16     8.7e-17

---

🧚 csin: Complex circular sine

SYNOPSIS

# void csin();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
csin($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

If

    z = x + iy,

then

    w = sin x  cosh y  +  i cos x sinh y.

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10      8400       5.3e-17     1.3e-17
   IEEE      -10,+10     30000       3.8e-16     1.0e-16

Also tested by csin(casin(z)) = z. 🔁

---

🍓 ccos: Complex circular cosine

SYNOPSIS

# void ccos();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
ccos($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

If

    z = x + iy,

then

    w = cos x  cosh y  -  i sin x sinh y.

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10      8400       4.5e-17     1.3e-17
   IEEE      -10,+10     30000       3.8e-16     1.0e-16

---

🔥 ctan: Complex circular tangent

SYNOPSIS

# void ctan();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
ctan($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

If

    z = x + iy,

then

          sin 2x  +  i sinh 2y
    w  =  --------------------.
           cos 2x  +  cosh 2y

Okay so on the real axis the denominator is zero at odd multiples of PI/2. The denominator is evaluated by its Taylor series near these points, so it doesn't just explode on us. Crisis averted. 💅

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10      5200       7.1e-17     1.6e-17
   IEEE      -10,+10     30000       7.2e-16     1.2e-16

Also tested by ctan * ccot = 1 and catan(ctan(z)) = z. 🔁

---

💜 ccot: Complex circular cotangent

SYNOPSIS

# void ccot();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
ccot($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

If

    z = x + iy,

then

          sin 2x  -  i sinh 2y
    w  =  --------------------.
           cosh 2y  -  cos 2x

On the real axis, the denominator has zeros at even multiples of PI/2. Near these points it is evaluated by a Taylor series. Same trick, different trig girly. ✨

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10      3000       6.5e-17     1.6e-17
   IEEE      -10,+10     30000       9.2e-16     1.2e-16

Also tested by ctan * ccot = 1 + i0. 🔁

---

🌈 casin: Complex circular arc sine

SYNOPSIS

# void casin();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
casin($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

Inverse complex sine, and here's the whole plot:

                                      2
w = -i clog( iz + csqrt( 1 - z ) ).

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10     10100       2.1e-15     3.4e-16
   IEEE      -10,+10     30000       2.2e-14     2.7e-15

Larger relative error can be observed for z near zero. 🙈 Also tested by csin(casin(z)) = z. 💖

---

⭐ cacos: Complex circular arc cosine

SYNOPSIS

# void cacos();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
cacos($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

Simply:

w = arccos z  =  PI/2 - arcsin z.

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10      5200      1.6e-15      2.8e-16
   IEEE      -10,+10     30000      1.8e-14      2.2e-15

---

👑 catan: Complex circular arc tangent

SYNOPSIS

# void catan();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
catan($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

DESCRIPTION

If

    z = x + iy,

then

         1       (    2x     )
Re w  =  - arctan(-----------)  +  k PI
         2       (     2    2)
                 (1 - x  - y )

              ( 2         2)
         1    (x  +  (y+1) )
Im w  =  - log(------------)
         4    ( 2         2)
              (x  +  (y-1) )

Where k is an arbitrary integer. 🎀

ACCURACY

                             Relative error:
        arithmetic   domain     # trials      peak         rms
           DEC       -10,+10      5900       1.3e-16     7.8e-18
           IEEE      -10,+10     30000       2.3e-15     8.5e-17

okay so the check catan( ctan(z) ) = z, with |x| and |y| < PI/2, had peak relative error 1.5e-16, rms relative error 2.9e-17. bestie also go peek at clog(). 🫶

---

💜 csinh: Complex hyperbolic sine

✨ SYNOPSIS:

# void csinh();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
csinh($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

📖 DESCRIPTION:

csinh z = (cexp(z) - cexp(-z))/2
        = sinh x * cos y  +  i cosh x * sin y .

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      -10,+10     30000       3.1e-16     8.2e-17

---

🦋 casinh: Complex inverse hyperbolic sine

✨ SYNOPSIS:

# void casinh();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
casinh($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w
print_new_cmplx($w);                 # prints $w as Re($w) + i Im($w)

📖 DESCRIPTION:

casinh z = -i casin iz .

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      -10,+10     30000       1.8e-14     2.6e-15

---

🎀 ccosh: Complex hyperbolic cosine

✨ SYNOPSIS:

# void ccosh();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
ccosh($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

📖 DESCRIPTION:

ccosh(z) = cosh x  cos y + i sinh x sin y .

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      -10,+10     30000       2.9e-16     8.1e-17

---

🌸 cacosh: Complex inverse hyperbolic cosine

✨ SYNOPSIS:

# void cacosh();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
cacosh($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

📖 DESCRIPTION:

acosh z = i acos z .

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      -10,+10     30000       1.6e-14     2.1e-15

---

💅 ctanh: Complex hyperbolic tangent

✨ SYNOPSIS:

# void ctanh();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
ctanh($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

📖 DESCRIPTION:

tanh z = (sinh 2x  +  i sin 2y) / (cosh 2x + cos 2y) .

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      -10,+10     30000       1.7e-14     2.4e-16

---

⭐ catanh: Complex inverse hyperbolic tangent

✨ SYNOPSIS:

# void catanh();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
catanh($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

📖 DESCRIPTION:

Inverse tanh, equal to -i catan (iz); 🌀

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      -10,+10     30000       2.3e-16     6.2e-17

---

👑 cpow: Complex power function

✨ SYNOPSIS:

# void cpow();
# cmplx a, z, w;

$a = new_cmplx(5, 6);    # $z = 5 + 6 i
$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
cpow($a, $z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

📖 DESCRIPTION:

Raises complex A to the complex Zth power. 💥 The definition is per AMS55 # 4.2.8, and it's analytically equivalent to cpow(a,z) = cexp(z clog(a)).

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      -10,+10     30000       9.4e-15     1.5e-15

---

🧚 cmplx: Complex number arithmetic

✨ SYNOPSIS:

# typedef struct {
#     double r;     real part
#     double i;     imaginary part
#    }cmplx;

# cmplx *a, *b, *c;

$a = new_cmplx(3, 5);   # $a = 3 + 5 i
$b = new_cmplx(2, 3);   # $b = 2 + 3 i
$c = new_cmplx();

cadd( $a, $b, $c );  #   c = b + a
csub( $a, $b, $c );  #   c = b - a
cmul( $a, $b, $c );  #   c = b * a
cdiv( $a, $b, $c );  #   c = b / a
cneg( $c );          #   c = -c
cmov( $b, $c );      #   c = b

print $c->{r}, '  ', $c->{i};   # prints real and imaginary parts of $c

📖 DESCRIPTION:

Addition: 🫶

c.r  =  b.r + a.r
c.i  =  b.i + a.i

Subtraction:

c.r  =  b.r - a.r
c.i  =  b.i - a.i

Multiplication:

c.r  =  b.r * a.r  -  b.i * a.i
c.i  =  b.r * a.i  +  b.i * a.r

Division:

d    =  a.r * a.r  +  a.i * a.i
c.r  = (b.r * a.r  + b.i * a.i)/d
c.i  = (b.i * a.r  -  b.r * a.i)/d

📏 ACCURACY:

So in DEC arithmetic, the test (1/z) * z = 1 had peak relative error 3.1e-17, rms 1.2e-17. And the test (y/z) * (z/y) = 1 had peak relative error 8.3e-17, rms 2.1e-17. no cap. 💅

Tests in the rectangle {-10,+10}:

                     Relative error:
arithmetic   function  # trials      peak         rms
   DEC        cadd       10000       1.4e-17     3.4e-18
   IEEE       cadd      100000       1.1e-16     2.7e-17
   DEC        csub       10000       1.4e-17     4.5e-18
   IEEE       csub      100000       1.1e-16     3.4e-17
   DEC        cmul        3000       2.3e-17     8.7e-18
   IEEE       cmul      100000       2.1e-16     6.9e-17
   DEC        cdiv       18000       4.9e-17     1.3e-17
   IEEE       cdiv      100000       3.7e-16     1.1e-16

---

🌈 cabs: Complex absolute value

✨ SYNOPSIS:

# double a, cabs();
# cmplx z;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$a = cabs( $z );

📖 DESCRIPTION:

If z = x + iy

then

      a = sqrt( x**2 + y**2 ).

Overflow and underflow are avoided by testing the magnitudes of x and y before squaring. If either is outside half of the floating point full scale range, both are rescaled. 🧠✨

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -30,+30     30000       3.2e-17     9.2e-18
   IEEE      -10,+10    100000       2.7e-16     6.9e-17

---

🍓 csqrt: Complex square root

✨ SYNOPSIS:

# void csqrt();
# cmplx z, w;

$z = new_cmplx(2, 3);    # $z = 2 + 3 i
$w = new_cmplx();
csqrt($z, $w );
print $w->{r}, '  ', $w->{i};  # prints real and imaginary parts of $w

📖 DESCRIPTION:

If z = x + iy, r = |z|, then

                      1/2
Im w  =  [ (r - x)/2 ]   ,

Re w  =  y / 2 Im w.

Heads up bestie: -w is also a square root of z. The root chosen is always in the upper half plane. 👑

Because of the potential for cancellation error in r - x, the result is sharpened by doing a Heron iteration (see sqrt.c) in complex arithmetic.

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -10,+10     25000       3.2e-17     9.6e-18
   IEEE      -10,+10    100000       3.2e-16     7.7e-17

                       2
Also tested by csqrt( z ) = z, and tested by arguments
close to the real axis.

---

💖 machconst: Globally declared constants

✨ SYNOPSIS:

extern double nameofconstant;

📖 DESCRIPTION:

This file contains a bunch of mathematical constants and also some needed size parameters of the computer arithmetic. 🧮 The values are supplied as arrays of hexadecimal integers for IEEE arithmetic; arrays of octal constants for DEC arithmetic; and in a normal decimal scientific notation for other machines. The particular notation used is determined by a symbol (DEC, IBMPC, or UNK) defined in the include file mconf.h.

The default size parameters are as follows. 🫶

For DEC and UNK modes:

MACHEP =  1.38777878078144567553E-17       2**-56
MAXLOG =  8.8029691931113054295988E1       log(2**127)
MINLOG = -8.872283911167299960540E1        log(2**-128)
MAXNUM =  1.701411834604692317316873e38    2**127

For IEEE arithmetic (IBMPC):

MACHEP =  1.11022302462515654042E-16       2**-53
MAXLOG =  7.09782712893383996843E2         log(2**1024)
MINLOG = -7.08396418532264106224E2         log(2**-1022)
MAXNUM =  1.7976931348623158E308           2**1024

These lists are subject to change. 🙈

---

🌸 cosh: Hyperbolic cosine

✨ SYNOPSIS:

# double x, y, cosh();

$y = cosh( $x );

📖 DESCRIPTION:

Returns hyperbolic cosine of argument in the range MINLOG to MAXLOG.

cosh(x)  =  ( exp(x) + exp(-x) )/2.

📏 ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       +- 88       50000       4.0e-17     7.7e-18
   IEEE     +-MAXLOG     30000       2.6e-16     5.7e-17

🚨 ERROR MESSAGES:

  message         condition      value returned
cosh overflow    |x| > MAXLOG       MAXNUM

---

💀 dawsn: Dawson's Integral

✨ SYNOPSIS:

# double x, y, dawsn();

$y = dawsn( $x );

📖 DESCRIPTION:

Approximates the integral

                            x
                            -
                     2     | |        2
 dawsn(x)  =  exp( -x  )   |    exp( t  ) dt
                         | |
                          -
                          0

Three different rational approximations are employed, for the intervals 0 to 3.25; 3.25 to 6.25; and 6.25 up. 🔥

📏 ACCURACY:

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE      0,10        10000       6.9e-16     1.0e-16
               DEC       0,10         6000       7.4e-17     1.4e-17

---

🎲 drand — pseudorandom number generator 🎲

SYNOPSIS:

# double y, drand();

($flag, $y) = drand( );

DESCRIPTION:

Bestie. It hands you a random number and it lives in the range 1.0 <= y < 2.0. That's it. That's the whole vibe. 🫶

The algorithm doing the heavy lifting is the three-generator congruential algorithm by Brian Wichmann and David Hill (BYTE magazine, March, 1987, pp 127-8), and the period they give for it is 6953607871644. Long period. Iconic. 💅

Now here's the tea on the different arithmetic compile time options, DEC, IBMPC, and MIEEE — they produce approximately the same sequences, differing only in the least significant bits of the numbers. The UNK option implements the algorithm as recommended in the BYTE article, and it may be used on all computers. But — and this is the catch, bestie — the low order bits of a double precision number may not be adequately random, and they may vary due to arithmetic implementation details on different computers. 😭

The other compile options, though? They generate an additional random integer that overwrites the low order bits of the double precision number ✨. This reduces the period by a factor of two but tends to overcome the problems mentioned. Trade-off queen behaviour tbh 👑

---

🥧 ellie — incomplete elliptic integral of the second kind 🥧

SYNOPSIS:

# double phi, m, y, ellie();

$y = ellie( $phi, $m );

DESCRIPTION:

It approximates the integral

               phi
                -
               | |
               |                   2
E(phi_\m)  =    |    sqrt( 1 - m sin t ) dt
               |
             | |
              -
               0

of amplitude phi and modulus m, using the arithmetic - geometric mean algorithm. Math girlies stay winning 🔥

ACCURACY:

Tested at random arguments with phi in [-10, 10] and m in [0, 1].

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC        0,2         2000       1.9e-16     3.4e-17
               IEEE     -10,10      150000       3.3e-15     1.4e-16

---

🍰 ellik — incomplete elliptic integral of the first kind 🍰

SYNOPSIS:

# double phi, m, y, ellik();

$y = ellik( $phi, $m );

DESCRIPTION:

It approximates the integral

               phi
                -
               | |
               |           dt
F(phi_\m)  =    |    ------------------
               |                   2
             | |    sqrt( 1 - m sin t )
              -
               0

of amplitude phi and modulus m, using the arithmetic - geometric mean algorithm. Same algorithm, different integral, she's consistent 💅

ACCURACY:

Tested at random points with m in [0, 1] and phi as indicated.

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE     -10,10       200000      7.4e-16     1.0e-16

---

🌸 ellpe — complete elliptic integral of the second kind 🌸

SYNOPSIS:

# double m1, y, ellpe();

$y = ellpe( $m1 );

DESCRIPTION:

It approximates the integral

           pi/2
            -
           | |                 2
E(m)  =    |    sqrt( 1 - m sin t ) dt
         | |
          -
           0

Where m = 1 - m1, using the approximation

     P(x)  -  x log x Q(x).

Now, listen — though there are no singularities, the argument m1 is used rather than m for compatibility with ellpk(). She's keeping the friend group consistent 🫶

Also, for the record: E(1) = 1; E(0) = pi/2.

ACCURACY:

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC        0, 1       13000       3.1e-17     9.4e-18
               IEEE       0, 1       10000       2.1e-16     7.3e-17

ERROR MESSAGES:

  message         condition      value returned
ellpe domain      x<0, x>1            0.0

---

🦋 ellpj — Jacobian Elliptic Functions 🦋

SYNOPSIS:

# double u, m, sn, cn, dn, phi;
# int ellpj();

($flag, $sn, $cn, $dn, $phi) = ellpj( $u, $m );

DESCRIPTION:

Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m), and dn(u|m) of parameter m between 0 and 1, and real argument u. Look at all those return values. She's an overachiever and we love her for it 💖

These functions are periodic, with quarter-period on the real axis equal to the complete elliptic integral ellpk(1.0-m).

Relation to incomplete elliptic integral: if u = ellik(phi,m), then sn(u|m) = sin(phi), and cn(u|m) = cos(phi). Phi is called the amplitude of u. iykyk ✨

Computation is by means of the arithmetic-geometric mean algorithm, except when m is within 1e-9 of 0 or 1. In the latter case with m close to 1, the approximation applies only for phi < pi/2. Heads up on that one 🙈

ACCURACY:

Tested at random points with u between 0 and 10, m between 0 and 1.

           Absolute error (* = relative error):
arithmetic   function   # trials      peak         rms
   DEC       sn           1800       4.5e-16     8.7e-17
   IEEE      phi         10000       9.2e-16*    1.4e-16*
   IEEE      sn          50000       4.1e-15     4.6e-16
   IEEE      cn          40000       3.6e-15     4.4e-16
   IEEE      dn          10000       1.3e-12     1.8e-14

Peak error observed in consistency check using addition theorem for sn(u+v) was 4e-16 (absolute). Also tested by the above relation to the incomplete elliptic integral. Accuracy deteriorates when u is large. So keep your u modest, bestie 😭

---

👑 ellpk — Complete elliptic integral of the first kind 👑

SYNOPSIS:

# double m1, y, ellpk();

$y = ellpk( $m1 );

DESCRIPTION:

It approximates the integral

           pi/2
            -
           | |
           |           dt
K(m)  =    |    ------------------
           |                   2
         | |    sqrt( 1 - m sin t )
          -
           0

where m = 1 - m1, using the approximation

    P(x)  -  log x Q(x).

The argument m1 is used rather than m so that the logarithmic singularity at m = 1 will be shifted to the origin — this preserves maximum accuracy. It's giving strategic relocation 💅

Also: K(0) = pi/2.

ACCURACY:

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC        0,1        16000       3.5e-17     1.1e-17
               IEEE       0,1        30000       2.5e-16     6.8e-17

ERROR MESSAGES:

  message         condition      value returned
ellpk domain       x<0, x>1           0.0

---

🍓 euclid — Rational arithmetic routines 🍓

SYNOPSIS:

# typedef struct
#     {
#     double n;  numerator
#     double d;  denominator
#     }fract;

$a = new_fract(3, 4);  # a = 3 / 4
$b = new_fract(2, 3);  # b = 2 / 3
$c = new_fract();
radd( $a, $b, $c ); #     c = b + a
rsub( $a, $b, $c ); #     c = b - a
rmul( $a, $b, $c ); #     c = b * a
rdiv( $a, $b, $c ); #     c = b / a
print $c->{n}, ' ', $c->{d};  # prints numerator and denominator of $c

($gcd, $m_reduced, $n_reduced) = euclid($m, $n);
# returns the greatest common divisor of $m and $n, as well as
# the result of reducing $m and $n by $gcd

Arguments of the routines are pointers to the structures. The double precision numbers are assumed, without checking, to be integer valued ⚠️. Overflow conditions are reported. So like, she trusts you, but she's also watching. 🫶

---

🧚 exp — Exponential function 🧚

SYNOPSIS:

# double x, y, exp();

$y = exp( $x );

DESCRIPTION:

Returns e (2.71828...) raised to the x power. The classic. The moment. 🌈

Range reduction is accomplished by separating the argument into an integer k and fraction f such that

    x    k  f
   e  = 2  e.

A Pade' form 1 + 2x P(x2)/( Q(x2) - P(x**2) ) of degree 2/3 is used to approximate exp(f) in the basic interval [-0.5, 0.5]. She said Pade' form and meant it 💅

ACCURACY:

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC       +- 88       50000       2.8e-17     7.0e-18
               IEEE      +- 708      40000       2.0e-16     5.6e-17

⚠️ Okay so this next bit is important and I need you to actually read it: error amplification in the exponential function can be a serious matter 😭. The error propagation involves

exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),

which shows that a 1 lsb error in representing X produces a relative error of X times 1 lsb in the function. While the routine gives an accurate result for arguments that are exactly represented by a double precision computer number, the result contains amplified roundoff error for large arguments not exactly represented. Not me warning you twice 💀

ERROR MESSAGES:

  message         condition      value returned
exp underflow    x < MINLOG         0.0
exp overflow     x > MAXLOG         INFINITY

---

⭐ expxx — exp(x*x) ⭐

#  double x, y, expxx();
# int sign;

  $y = expxx( $x, $sign );

DESCRIPTION:

Computes y = exp(x*x) while suppressing the error amplification that would ordinarily arise from the inexactness of the exponential argument x*x. She's doing damage control before the damage even happens. Queen behaviour 👑

If sign < 0, exp(-x*x) is returned. If sign > 0, or omitted, exp(x*x) is returned.

ACCURACY:

                              Relative error:
arithmetic    domain     # trials      peak         rms
   IEEE      -26.6, 26.6    10^7       3.9e-16     8.9e-17

---

🔥 exp10 — Base 10 exponential function (Common antilogarithm) 🔥

SYNOPSIS:

# double x, y, exp10();

$y = exp10( $x );

DESCRIPTION:

Returns 10 raised to the x power. Ten out of ten, literally 💅

Range reduction is accomplished by expressing the argument as 10x = 2n 10**f, with |f| < 0.5 log10(2). The Pade' form

   1 + 2x P(x**2)/( Q(x**2) - P(x**2) )

is used to approximate 10**f. Same Pade' energy as exp, different base. Sisters fr fr 🫶

ACCURACY:

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE     -307,+307    30000       2.2e-16     5.5e-17
            Test result from an earlier version (2.1):
               DEC       -38,+38     70000       3.1e-17     7.0e-18

ERROR MESSAGES:

              message         condition      value returned
            exp10 underflow    x < -MAXL10        0.0
            exp10 overflow     x > MAXL10       MAXNUM

For DEC arithmetic: MAXL10 = 38.230809449325611792.

For IEEE arithmetic: MAXL10 = 308.2547155599167.

---

exp2: Base 2 exponential function 🌸

SYNOPSIS:

# double x, y, exp2();

$y = exp2( $x );

DESCRIPTION:

Returns 2 raised to the x power. Yes bestie, it's the glow-up of exp but for base 2 specifically 💅

Range reduction is accomplished by separating the argument into an integer k and fraction f such that

    x    k  f
   2  = 2  2.

A Pade' form

1 + 2x P(x**2) / (Q(x**2) - x P(x**2) )

approximates 2**x in the basic range [-0.5, 0.5].

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE    -1022,+1024   30000       1.8e-16     5.4e-17

See exp.c for comments on error amplification.

ERROR MESSAGES:

  message         condition      value returned
exp underflow    x < -MAXL2        0.0
exp overflow     x > MAXL2         MAXNUM

• For DEC arithmetic, MAXL2 = 127.
• For IEEE arithmetic, MAXL2 = 1024.

---

ei: Exponential integral ✨

SYNOPSIS:

#double x, y, ei();

$y = ei( $x );

DESCRIPTION:

              x
               -     t
              | |   e
   Ei(x) =   -|-   ---  dt .
            | |     t
             -
            -inf

Not defined for x <= 0. Like, at all. Don't even try it.

See also expn.c.

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE       0,100       50000      8.6e-16     1.3e-16

---

expn: Exponential integral En 💖

SYNOPSIS:

# int n;
# double x, y, expn();

$y = expn( $n, $x );

DESCRIPTION:

Evaluates the exponential integral

                inf.
                  -
                 | |   -xt
                 |    e
     E (x)  =    |    ----  dt.
      n          |      n
               | |     t
                -
                 1

Both n and x must be nonnegative 👑

The routine employs either a power series, a continued fraction, or an asymptotic formula depending on the relative values of n and x. It picks its own vibe based on the situation, slay.

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       0, 30        5000       2.0e-16     4.6e-17
   IEEE      0, 30       10000       1.7e-15     3.6e-16

---

fabs: Absolute value 🦋

SYNOPSIS:

# double x, y;

$y = fabs( $x );

DESCRIPTION:

Returns the absolute value of the argument. Negative energy cancelled, we only do positive vibes here 💅

---

fac: Factorial function 🍓

SYNOPSIS:

# double y, fac();
# int i;

$y = fac( $i );

DESCRIPTION:

Returns factorial of i = 1 * 2 * 3 * ... * i.

fac(0) = 1.0.

Due to machine arithmetic bounds the largest value of i accepted is 33 in DEC arithmetic or 170 in IEEE arithmetic. Greater values, or negative ones, produce an error message and return MAXNUM.

ACCURACY:

For i < 34 the values are simply tabulated, and have full machine accuracy. If i > 55, fac(i) = gamma(i+1); see gamma.c.

                     Relative error:
arithmetic   domain      peak
   IEEE      0, 170    1.4e-15
   DEC       0, 33      1.4e-17

---

fdtr: F distribution ✨

SYNOPSIS:

# int df1, df2;
# double x, y, fdtr();

$y = fdtr( $df1, $df2, $x );

DESCRIPTION:

Returns the area from zero to x under the F density function (also known as Snedcor's density or the variance ratio density). This is the density of x = (u1/df1)/(u2/df2), where u1 and u2 are random variables having Chi square distributions with df1 and df2 degrees of freedom, respectively.

The incomplete beta integral is used, according to the formula

       P(x) = incbet( df1/2, df2/2, df1*x/(df2 + df1*x) ).

The arguments a and b are greater than zero, and x is nonnegative.

ACCURACY:

Tested at random points (a,b,x).

               x     a,b                     Relative error:
arithmetic  domain  domain     # trials      peak         rms
   IEEE      0,1    0,100       100000      9.8e-15     1.7e-15
   IEEE      1,5    0,100       100000      6.5e-15     3.5e-16
   IEEE      0,1    1,10000     100000      2.2e-11     3.3e-12
   IEEE      1,5    1,10000     100000      1.1e-11     1.7e-13

See also incbet.c.

ERROR MESSAGES:

  message         condition      value returned
fdtr domain     a<0, b<0, x<0         0.0

---

fdtrc: Complemented F distribution 💜

SYNOPSIS:

# int df1, df2;
# double x, y, fdtrc();

$y = fdtrc( $df1, $df2, $x );

DESCRIPTION:

Returns the area from x to infinity under the F density function (also known as Snedcor's density or the variance ratio density).

                     inf.
                      -
             1       | |  a-1      b-1
1-P(x)  =  ------    |   t    (1-t)    dt
           B(a,b)  | |
                    -
                     x

The incomplete beta integral is used, according to the formula

       P(x) = incbet( df2/2, df1/2, df2/(df2 + df1*x) ).

ACCURACY:

Tested at random points (a,b,x) in the indicated intervals.

               x     a,b                     Relative error:
arithmetic  domain  domain     # trials      peak         rms
   IEEE      0,1    1,100       100000      3.7e-14     5.9e-16
   IEEE      1,5    1,100       100000      8.0e-15     1.6e-15
   IEEE      0,1    1,10000     100000      1.8e-11     3.5e-13
   IEEE      1,5    1,10000     100000      2.0e-11     3.0e-12

See also incbet.c.

ERROR MESSAGES:

  message         condition      value returned
fdtrc domain    a<0, b<0, x<0         0.0

---

fdtri: Inverse of complemented F distribution 🎀

SYNOPSIS:

# int df1, df2;
# double x, p, fdtri();

$x = fdtri( $df1, $df2, $p );

DESCRIPTION:

Finds the F density argument x such that the integral from x to infinity of the F density is equal to the given probability p.

This is accomplished using the inverse beta integral function and the relations

     z = incbi( df2/2, df1/2, p )
     x = df2 (1-z) / (df1 z).

Note: the following relations hold for the inverse of the uncomplemented F distribution:

     z = incbi( df1/2, df2/2, p )
     x = df2 z / (df1 (1-z)).

ACCURACY:

Tested at random points (a,b,p).

             a,b                     Relative error:
arithmetic  domain     # trials      peak         rms
 For p between .001 and 1:
   IEEE     1,100       100000      8.3e-15     4.7e-16
   IEEE     1,10000     100000      2.1e-11     1.4e-13
 For p between 10^-6 and 10^-3:
   IEEE     1,100        50000      1.3e-12     8.4e-15
   IEEE     1,10000      50000      3.0e-12     4.8e-14

See also fdtrc.c.

ERROR MESSAGES:

  message         condition      value returned
fdtri domain   p <= 0 or p > 1       0.0
                    v < 1

---

ceil: ceil 🌈

ceil() returns the smallest integer greater than or equal to x. It truncates toward plus infinity. ⬆️ bestie only goes UP

SYNOPSIS:

# double x, y, ceil();

$y = ceil( $x );

---

floor: floor 🌸

floor() returns the largest integer less than or equal to x. It truncates toward minus infinity. ⬇️ and this one only goes DOWN

SYNOPSIS:

# double x, y, floor();

$y = floor( $x );

---

frexp: frexp 💫

frexp() extracts the exponent from x. It returns an integer power of two to expnt and the significand between 0.5 and 1 to y. Thus x = y * 2**expn.

SYNOPSIS:

# double x, y, frexp();
# int expnt;

($y, $expnt)  = frexp( $x );

---

ldexp: multiplies x by 2**n. 🔥

SYNOPSIS:

# double x, y, ldexp();
# int n;

$y = ldexp( $x, $n );

---

fresnl: Fresnel integral ✨

SYNOPSIS:

# double x, S, C;
# void fresnl();

($flag, $S, $C) = fresnl( $x );

DESCRIPTION:

Evaluates the Fresnel integrals

          x
          -
         | |
C(x) =   |   cos(pi/2 t**2) dt,
       | |
        -
         0

          x
          -
         | |
S(x) =   |   sin(pi/2 t**2) dt.
       | |
        -
         0

The integrals are evaluated by a power series for x < 1. For x >= 1 auxiliary functions f(x) and g(x) are employed such that

C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )

ACCURACY:

Relative error.

Arithmetic  function   domain     # trials      peak         rms
  IEEE       S(x)      0, 10       10000       2.0e-15     3.2e-16
  IEEE       C(x)      0, 10       10000       1.8e-15     3.3e-16
  DEC        S(x)      0, 10        6000       2.2e-16     3.9e-17
  DEC        C(x)      0, 10        5000       2.3e-16     3.9e-17

---

gamma: Gamma function 👑

SYNOPSIS:

# double x, y, gamma();
# extern int sgngam;

$y = gamma( $x );

DESCRIPTION:

Returns gamma function of the argument. The result is correctly signed, and the sign (+1 or -1) is also returned in a global (extern) variable named sgngam. This variable is also filled in by the logarithmic gamma function lgam().

Arguments |x| <= 34 are reduced by recurrence and the function approximated by a rational function of degree 6/7 in the interval (2,3). Large arguments are handled by Stirling's formula. Large negative arguments are made positive using a reflection formula.

ACCURACY:

Relative error:

            arithmetic   domain     # trials      peak         rms
               DEC      -34, 34      10000       1.3e-16     2.5e-17
               IEEE    -170,-33      20000       2.3e-15     3.3e-16
               IEEE     -33,  33     20000       9.4e-16     2.2e-16
               IEEE      33, 171.6   20000       2.3e-15     3.2e-16

and bestie, listen — if you feed it arguments outside the test range, the error gets bigger. Why? Because the exponential function amplifies it. That's just how it is, no cap. 😭

---

✨ lgam — Natural logarithm of gamma function ✨

SYNOPSIS:

# double x, y, lgam();
# extern int sgngam;

$y = lgam( $x );

DESCRIPTION:

Okay so this one gives you the base e (that's 2.718...) logarithm of the absolute value of the gamma function of your argument. 💜

And here's the sneaky part: the sign of the gamma function — +1 or -1 — doesn't come back as a return value, it gets dropped into a global (extern) variable called sgngam. So check that too, fr fr. 👑

For arguments greater than 13, the logarithm of the gamma function gets approximated by the logarithmic version of Stirling's formula, using a polynomial approximation of degree 4. Arguments between -33 and +33 get reduced by recurrence to the interval [2,3] of a rational approximation. And the cosecant reflection formula handles arguments less than -33.

Arguments greater than MAXLGM return MAXNUM and print an error message. 🌈 And MAXLGM = 2.035093e36 for DEC arithmetic or 2.556348e305 for IEEE arithmetic.

ACCURACY:

arithmetic      domain        # trials     peak         rms
   DEC     0, 3                  7000     5.2e-17     1.3e-17
   DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
   IEEE    0, 3                 28000     5.4e-16     1.1e-16
   IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17

The error criterion was relative when the function magnitude was greater than one, but absolute when it was less than one. 💅

The next test used the relative error criterion — though at certain points the relative error could be way higher than what's shown:

   IEEE    -200, -4             10000     4.8e-16     1.3e-16

---

🦋 gdtr — Gamma distribution function 🦋

SYNOPSIS:

# double a, b, x, y, gdtr();

$y = gdtr( $a, $b, $x );

DESCRIPTION:

Returns the integral from zero to x of the gamma probability density function:

               x
       b       -
      a       | |   b-1  -at
y =  -----    |    t    e    dt
      -     | |
     | (b)   -
              0

The incomplete gamma integral is used, according to the relation

y = igam( b, ax ).

ACCURACY:

See igam(). 👀

ERROR MESSAGES:

  message         condition      value returned
gdtr domain         x < 0            0.0

---

🌸 gdtrc — Complemented gamma distribution function 🌸

SYNOPSIS:

# double a, b, x, y, gdtrc();

$y = gdtrc( $a, $b, $x );

DESCRIPTION:

Returns the integral from x to infinity of the gamma probability density function:

              inf.
       b       -
      a       | |   b-1  -at
y =  -----    |    t    e    dt
      -     | |
     | (b)   -
              x

The incomplete gamma integral is used, according to the relation

y = igamc( b, ax ).

ACCURACY:

See igamc(). 💖

ERROR MESSAGES:

  message         condition      value returned
gdtrc domain         x < 0            0.0

---

🎀 hyp2f0 — Gauss hypergeometric function 2F0 🎀

SYNOPSIS:

# double a, b, x, value, *err;
# int type;    /* determines what converging factor to use */

($value, $err) =  hyp2f0( $a, $b, $x, $type )

---

💅 hyp2f1 — Gauss hypergeometric function 2F1 💅

SYNOPSIS:

# double a, b, c, x, y, hyp2f1();

$y = hyp2f1( $a, $b, $c, $x );

DESCRIPTION:

 hyp2f1( a, b, c, x )  =   F ( a, b; c; x )
                          2 1

          inf.
           -   a(a+1)...(a+k) b(b+1)...(b+k)   k+1
  =  1 +   >   -----------------------------  x   .
           -         c(c+1)...(c+k) (k+1)!
         k = 0

The cases it handles are:

• Tests and escapes for negative integer a, b, or c
• Linear transformation if c - a or c - b is a negative integer
• Special case c = a or c = b
• Linear transformation for x near +1
• Transformation for x < -0.5
• Psi function expansion if x > 0.5 and c - a - b is integer
• Conditionally, a recurrence on c to make c-a-b > 0

|x| > 1 is rejected. 🚫 Not playing with that.

The parameters a, b, c are considered to be integer valued if they're within 1.0e-14 of the nearest integer (1.0e-13 for IEEE arithmetic). 👑

ACCURACY:

              Relative error (-1 < x < 1):
arithmetic   domain     # trials      peak         rms
   IEEE      -1,7        230000      1.2e-11     5.2e-14

Several special cases were also tested with a, b, c in the range -7 to 7. ✨

ERROR MESSAGES:

A "partial loss of precision" message gets printed if the internally estimated relative error exceeds 1^-12. A "singularity" message gets printed on overflow or in cases not addressed (like x < -1). 😭

---

🔥 hyperg — Confluent hypergeometric function 🔥

SYNOPSIS:

# double a, b, x, y, hyperg();

$y = hyperg( $a, $b, $x );

DESCRIPTION:

Computes the confluent hypergeometric function

                         1           2
                      a x    a(a+1) x
  F ( a,b;x )  =  1 + ---- + --------- + ...
 1 1                  b 1!   b(b+1) 2!

Many higher transcendental functions are special cases of this power series. it's giving iconic source material. 🌟

As is evident from the formula, b must not be a negative integer or zero — unless a is an integer with 0 >= a > b.

The routine attempts both a direct summation of the series and an asymptotic expansion. In each case it estimates the error from roundoff, cancellation, and nonconvergence. Whichever result has the smaller estimated error gets returned. 🫶

ACCURACY:

Tested at random points (a, b, x), all three variables ranging from 0 to 30.

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       0,30         2000       1.2e-15     1.3e-16
   IEEE      0,30        30000       1.8e-14     1.1e-15

Larger errors can be observed when b is near a negative integer or zero. Certain combinations of arguments yield serious cancellation error in the power series summation and also aren't in the region of near convergence of the asymptotic series. 💀 An error message gets printed if the self-estimated relative error is greater than 1.0e-12.

---

⭐ i0 — Modified Bessel function of order zero ⭐

SYNOPSIS:

# double x, y, i0();

$y = i0( $x );

DESCRIPTION:

Returns the modified Bessel function of order zero of the argument.

The function is defined as i0(x) = j0( ix ).

The range is partitioned into the two intervals [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval. 🎀

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       0,30         6000       8.2e-17     1.9e-17
   IEEE      0,30        30000       5.8e-16     1.4e-16

---

🌈 i0e — Modified Bessel function of order zero, exponentially scaled 🌈

SYNOPSIS:

# double x, y, i0e();

$y = i0e( $x );

DESCRIPTION:

Returns the exponentially scaled modified Bessel function of order zero of the argument.

The function is defined as i0e(x) = exp(-|x|) j0( ix ).

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      0,30        30000       5.4e-16     1.2e-16

See i0(). 💖

---

🦋 i1 — Modified Bessel function of order one 🦋

SYNOPSIS:

# double x, y, i1();

$y = i1( $x );

DESCRIPTION:

Returns the modified Bessel function of order one of the argument.

The function is defined as i1(x) = -i j1( ix ).

The range is partitioned into the two intervals [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval. ✨

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       0, 30        3400       1.2e-16     2.3e-17
   IEEE      0, 30       30000       1.9e-15     2.1e-16

---

💜 i1e — Modified Bessel function of order one, exponentially scaled 💜

SYNOPSIS:

# double x, y, i1e();

$y = i1e( $x );

DESCRIPTION:

Returns the exponentially scaled modified Bessel function of order one of the argument.

The function is defined as i1(x) = -i exp(-|x|) j1( ix ).

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      0, 30       30000       2.0e-15     2.0e-16

See i1(). 🌸

---

👑 igam — Incomplete gamma integral 👑

SYNOPSIS:

# double a, x, y, igam();

$y = igam( $a, $x );

DESCRIPTION:

The function is defined by

                          x
                           -
                  1       | |  -t  a-1
 igam(a,x)  =   -----     |   e   t   dt.
                 -      | |
                | (a)    -
                          0

In this implementation both arguments must be positive. 🫶 The integral is evaluated by either a power series or continued fraction expansion, depending on the relative values of a and x.

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      0,30       200000       3.6e-14     2.9e-15
   IEEE      0,100      300000       9.9e-14     1.5e-14

---

🍓 igamc — Complemented incomplete gamma integral 🍓

SYNOPSIS:

# double a, x, y, igamc();

$y = igamc( $a, $x );

DESCRIPTION:

The function is defined by

 igamc(a,x)   =   1 - igam(a,x)

And then, the other side of the coin — the upper incomplete gamma, computed from x all the way out to infinity, a.k.a. the complemented version:

                                       inf.
                                         -
                                1       | |  -t  a-1
                          =   -----     |   e   t   dt.
                               -      | |
                              | (a)    -
                                        x

Bestie, listen: in this implementation both arguments must be positive. The integral gets evaluated either by a power series or by a continued fraction expansion, depending on the relative values of a and x. She picks whichever one slays harder for the numbers you handed her. 💅

ACCURACY:

Tested at random a, x.

               a         x                      Relative error:
arithmetic   domain   domain     # trials      peak         rms
   IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
   IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15

---

igami: Inverse of complemented imcomplete gamma integral

SYNOPSIS:

# double a, x, p, igami();

$x = igami( $a, $p );

DESCRIPTION:

You give her p, and she finds you the x such that

 igamc( a, x ) = p.

It's valid in the right-hand tail of the distribution, so p < 0.5. (Stay on the right side of the tail, bb, that's her zone. ✨)

She starts with the approximate value

        3
 x = a t

where

 t = 1 - d - ndtri(p) sqrt(d)

and

 d = 1/9a,

and then the routine performs up to 10 Newton iterations to find the root of igamc(a,x) - p = 0. Yes, she's iterating. Ten times. Iconic problem-solving behavior. 👑

ACCURACY:

Tested at random a, p in the intervals indicated.

               a        p                      Relative error:
arithmetic   domain   domain     # trials      peak         rms
   IEEE     0.5,100   0,0.5       100000       1.0e-14     1.7e-15
   IEEE     0.01,0.5  0,0.5       100000       9.0e-14     3.4e-15
   IEEE    0.5,10000  0,0.5        20000       2.3e-13     3.8e-14

---

incbet: Incomplete beta integral

SYNOPSIS:

# double a, b, x, y, incbet();

$y = incbet( $a, $b, $x );

DESCRIPTION:

Returns the incomplete beta integral of the arguments, evaluated from zero to x. The function is defined as

                 x
    -            -
   | (a+b)      | |  a-1     b-1
 -----------    |   t   (1-t)   dt.
  -     -     | |
 | (a) | (b)   -
                0

The domain of definition is 0 <= x <= 1. In this implementation a and b are restricted to positive values. ✨

Want the integral from x to 1 instead? That's literally free, just use the symmetry relation:

   1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).

The integral is evaluated by a continued fraction expansion or, when b*x is small, by a power series.

ACCURACY:

Tested at uniformly distributed random points (a,b,x) with a and b in "domain" and x between 0 and 1.

                                           Relative error
arithmetic   domain     # trials      peak         rms
   IEEE      0,5         10000       6.9e-15     4.5e-16
   IEEE      0,85       250000       2.2e-13     1.7e-14
   IEEE      0,1000      30000       5.3e-12     6.3e-13
   IEEE      0,10000    250000       9.3e-11     7.1e-12
   IEEE      0,100000    10000       8.7e-10     4.8e-11

Outputs smaller than the IEEE gradual underflow threshold were excluded from these statistics. (Not me being picky about data hygiene, it's just science 💅)

ERROR MESSAGES:

  message         condition      value returned
incbet domain      x<0, x>1          0.0
incbet underflow                     0.0

---

incbi: Inverse of imcomplete beta integral

SYNOPSIS:

# double a, b, x, y, incbi();

$x = incbi( $a, $b, $y );

DESCRIPTION:

Given y, the function finds x such that

 incbet( a, b, x ) = y .

The routine performs interval halving or Newton iterations to find the root of incbet(a,b,x) - y = 0. Two strategies because she contains multitudes. 🦋

ACCURACY:

                     Relative error:
               x     a,b
arithmetic   domain  domain  # trials    peak       rms
   IEEE      0,1    .5,10000   50000    5.8e-12   1.3e-13
   IEEE      0,1   .25,100    100000    1.8e-13   3.9e-15
   IEEE      0,1     0,5       50000    1.1e-12   5.5e-15
   VAX       0,1    .5,100     25000    3.5e-14   1.1e-15

With a and b constrained to half-integer or integer values:

   IEEE      0,1    .5,10000   50000    5.8e-12   1.1e-13
   IEEE      0,1    .5,100    100000    1.7e-14   7.9e-16

With a = .5, b constrained to half-integer or integer values:

   IEEE      0,1    .5,10000   10000    8.3e-11   1.0e-11

---

iv: Modified Bessel function of noninteger order

SYNOPSIS:

# double v, x, y, iv();

$y = iv( $v, $x );

DESCRIPTION:

Returns the modified Bessel function of order v of the argument. Heads up bestie: if x is negative, v must be integer valued. No exceptions, that's the vibe check. 🚨

The function is defined as Iv(x) = Jv( ix ). It is here computed in terms of the confluent hypergeometric function, according to the formula

             v  -x
Iv(x) = (x/2)  e   hyperg( v+0.5, 2v+1, 2x ) / gamma(v+1)

If v is a negative integer, then v is replaced by -v.

ACCURACY:

Tested at random points (v, x), with v between 0 and 30, x between 0 and 28.

                             Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       0,30          2000      3.1e-15     5.4e-16
   IEEE      0,30         10000      1.7e-14     2.7e-15

Accuracy is diminished if v is near a negative integer. 😭

See also hyperg.c.

---

j0: Bessel function of order zero

SYNOPSIS:

# double x, y, j0();

$y = j0( $x );

DESCRIPTION:

Returns the Bessel function of order zero of the argument.

The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval the following rational approximation is used:

       2         2
(w - r  ) (w - r  ) P (w) / Q (w)
      1         2    3       8

           2
where w = x  and the two r's are zeros of the function.

In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7. She's got range. 💅

ACCURACY:

                     Absolute error:
arithmetic   domain     # trials      peak         rms
   DEC       0, 30       10000       4.4e-17     6.3e-18
   IEEE      0, 30       60000       4.2e-16     1.1e-16

---

y0: Bessel function of the second kind, order zero

SYNOPSIS:

# double x, y, y0();

$y = y0( $x );

DESCRIPTION:

Returns the Bessel function of the second kind, of order zero, of the argument.

The domain is divided into the intervals [0, 5] and (5, infinity). In the first interval a rational approximation R(x) is employed to compute

  y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.

Thus a call to j0() is required. (Yes, she's calling her friend. Collaboration era. 🫶)

In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7.

ACCURACY:

Absolute error, when y0(x) < 1; else relative error:

arithmetic   domain     # trials      peak         rms
   DEC       0, 30        9400       7.0e-17     7.9e-18
   IEEE      0, 30       30000       1.3e-15     1.6e-16

---

j1: Bessel function of order one

SYNOPSIS:

# double x, y, j1();

$y = j1( $x );

DESCRIPTION:

Returns the Bessel function of order one of the argument.

The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 24 term Chebyshev expansion is used. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5.

ACCURACY:

                     Absolute error:
arithmetic   domain      # trials      peak         rms
   DEC       0, 30       10000       4.0e-17     1.1e-17
   IEEE      0, 30       30000       2.6e-16     1.1e-16

---

y1: Bessel function of second kind of order one

SYNOPSIS:

# double x, y, y1();

$y = y1( $x );

DESCRIPTION:

Returns the Bessel function of the second kind of order one of the argument.

The domain is divided into the intervals [0, 8] and (8, infinity). In the first interval a 25 term Chebyshev expansion is used, and a call to j1() is required. In the second, the asymptotic trigonometric representation is employed using two rational functions of degree 5/5.

ACCURACY:

                     Absolute error:
arithmetic   domain      # trials      peak         rms
   DEC       0, 30       10000       8.6e-17     1.3e-17
   IEEE      0, 30       30000       1.0e-15     1.3e-16

(error criterion relative when |y1| > 1).

---

jn: Bessel function of integer order

SYNOPSIS:

# int n;
# double x, y, jn();

$y = jn( $n, $x );

DESCRIPTION:

Returns the Bessel function of order n, where n is a (possibly negative) integer.

The ratio of jn(x) to j0(x) is computed by backward recurrence. First the ratio jn/jn-1 is found by a continued fraction expansion. Then the recurrence relating successive orders is applied until j0 or j1 is reached. She's walking backwards through the orders like it's nothing. Main character energy fr fr. 🚶‍♀️✨

If n = 0 or 1 the routine for j0 or j1 is called directly.

ACCURACY:

                     Absolute error:
arithmetic   range      # trials      peak         rms
   DEC       0, 30        5500       6.9e-17     9.3e-18
   IEEE      0, 30        5000       4.4e-16     7.9e-17

Not suitable for large n or x. Use jv() instead.

---

jv: Bessel function of noninteger order

SYNOPSIS:

# double v, x, y, jv();

$y = jv( $v, $x );

DESCRIPTION:

Returns the Bessel function of order v of the argument, where v is real. Negative x is allowed if v is an integer.

Several expansions are included: the ascending power series, the Hankel expansion, and two transitional expansions for large v. If v is not too large, it is reduced by recurrence to a region of best accuracy. The transitional expansions give 12D accuracy for v > 500. Twelve digits, bestie. She ate. 🔥

ACCURACY:

Results for integer v are indicated by *, where x and v both vary from -125 to +125. Otherwise, x ranges from 0 to 125, v ranges as indicated by "domain." Error criterion is absolute, except relative when |jv()| > 1.

…picking up RIGHT where we left off, here's the rest of that accuracy table bestie, no cap 💅

            arithmetic  v domain  x domain    # trials      peak       rms
               IEEE      0,125     0,125      100000      4.6e-15    2.2e-16
               IEEE   -125,0       0,125       40000      5.4e-11    3.7e-13
               IEEE      0,500     0,500       20000      4.4e-15    4.0e-16
            Integer v:
               IEEE   -125,125   -125,125      50000      3.5e-15*   1.9e-16*

---

🎀 k0: Modified Bessel function, third kind, order zero 🎀

SYNOPSIS

# double x, y, k0();

$y = k0( $x );

DESCRIPTION

Okay so this one hands you back the modified Bessel function of the third kind, order zero, of whatever argument you feed it. ✨

The range gets sliced into two lil intervals — [0,8] and (8, infinity) — and Chebyshev polynomial expansions do the work in each one. Iconic behavior honestly.

ACCURACY

It was tested at 2000 random points between 0 and 8. Peak absolute error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15. Slay.

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC       0, 30        3100       1.3e-16     2.1e-17
               IEEE      0, 30       30000       1.2e-15     1.6e-16

ERROR MESSAGES

  message         condition      value returned
 K0 domain          x <= 0          MAXNUM

So if x is <= 0? Yeah that's a domain issue and you get MAXNUM back. Not the vibe you wanted, but them's the rules. 😭

---

🦋 k0e: Modified Bessel function, third kind, order zero, exponentially scaled 🦋

SYNOPSIS

# double x, y, k0e();

$y = k0e( $x );

DESCRIPTION

Returns the exponentially scaled modified Bessel function of the third kind of order zero of the argument. And like, here's the formula so you know exactly what "exponentially scaled" means:

     k0e(x) = exp(x) * k0(x).

ACCURACY

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE      0, 30       30000       1.4e-15     1.4e-16

See k0(). 👑

---

💜 k1: Modified Bessel function, third kind, order one 💜

SYNOPSIS

# double x, y, k1();

$y = k1( $x );

DESCRIPTION

Computes the modified Bessel function of the third kind of order one of the argument.

The range is partitioned into the two intervals [0,2] and (2, infinity), and yeah, Chebyshev polynomial expansions are employed in each interval. Same energy as k0, different order. 💅

ACCURACY

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC       0, 30        3300       8.9e-17     2.2e-17
               IEEE      0, 30       30000       1.2e-15     1.6e-16

ERROR MESSAGES

  message         condition      value returned
k1 domain          x <= 0          MAXNUM

x <= 0 again? Yeah bestie, domain error, MAXNUM returned. Note the lowercase k1 in the message this time, that's literally how it's written, don't @ me. 🙈

---

🌸 k1e: Modified Bessel function, third kind, order one, exponentially scaled 🌸

SYNOPSIS

# double x, y, k1e();

$y = k1e( $x );

DESCRIPTION

Returns the exponentially scaled modified Bessel function of the third kind of order one of the argument:

     k1e(x) = exp(x) * k1(x).

ACCURACY

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE      0, 30       30000       7.8e-16     1.2e-16

See k1(). 🫶

---

⭐ kn: Modified Bessel function, third kind, integer order ⭐

SYNOPSIS

# double x, y, kn();
# int n;

$y = kn( $n, $x );

DESCRIPTION

Returns the modified Bessel function of the third kind of order n of the argument.

The range is partitioned into the two intervals [0,9.55] and (9.55, infinity). An ascending power series is used in the low range, and an asymptotic expansion in the high range. Two strategies, one function, main character energy. ✨

ACCURACY

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC       0,30         3000       1.3e-9      5.8e-11
               IEEE      0,30        90000       1.8e-8      3.0e-10

Heads up: the error is high only near the crossover point x = 9.55 between the two expansions used. So like, right at the handoff it gets a little messy. It's giving seam. 💀

---

🌈 log: Natural logarithm 🌈

SYNOPSIS

# double x, y, log();

$y = log( $x );

DESCRIPTION

Returns the base e (2.718...) logarithm of x.

The argument is separated into its exponent and fractional parts. If the exponent is between -1 and +1, the logarithm of the fraction is approximated by

    log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).

Otherwise, setting z = 2(x-1)/x+1),

    log(x) = z + z**3 P(z)/Q(z).

ACCURACY

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE      0.5, 2.0    150000      1.44e-16    5.06e-17
               IEEE      +-MAXNUM    30000       1.20e-16    4.78e-17
               DEC       0, 10       170000      1.8e-17     6.3e-18

In the tests over the interval [+-MAXNUM], the logarithms of the random arguments were uniformly distributed over [0, MAXLOG]. Fancy. 💖

ERROR MESSAGES

• log singularity: x = 0; returns -INFINITY
• log domain: x < 0; returns NAN

Okay so zero is a singularity and negatives are a domain error. Negative logs aren't real (literally 💀). You get NAN and you deserve to know why.

---

🍓 log10: Common logarithm 🍓

SYNOPSIS

# double x, y, log10();

$y = log10( $x );

DESCRIPTION

Returns the logarithm to the base 10 of x.

The argument is separated into its exponent and fractional parts. The logarithm of the fraction is approximated by

    log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).

ACCURACY

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE      0.5, 2.0     30000      1.5e-16     5.0e-17
               IEEE      0, MAXNUM    30000      1.4e-16     4.8e-17
               DEC       1, MAXNUM    50000      2.5e-17     6.0e-18

In the tests over the interval [1, MAXNUM], the logarithms of the random arguments were uniformly distributed over [0, MAXLOG].

ERROR MESSAGES

• log10 singularity: x = 0; returns -INFINITY
• log10 domain: x < 0; returns NAN

---

🧚 log2: Base 2 logarithm 🧚

SYNOPSIS

# double x, y, log2();

$y = log2( $x );

DESCRIPTION

Returns the base 2 logarithm of x. Yes queen, binary logs. 👑

The argument is separated into its exponent and fractional parts. If the exponent is between -1 and +1, the base e logarithm of the fraction is approximated by

    log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).

Otherwise, setting z = 2(x-1)/x+1),

    log(x) = z + z**3 P(z)/Q(z).

ACCURACY

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE      0.5, 2.0    30000       2.0e-16     5.5e-17
               IEEE      exp(+-700)  40000       1.3e-16     4.6e-17

In the tests over the interval [exp(+-700)], the logarithms of the random arguments were uniformly distributed.

ERROR MESSAGES

• log2 singularity: x = 0; returns -INFINITY
• log2 domain: x < 0; returns NAN

---

💀 lrand: Pseudorandom number generator 💀

SYNOPSIS

long y, lrand();

$y = lrand( );

DESCRIPTION

Yields a long integer random number. RNG girlies rise up. 🎲

The three-generator congruential algorithm by Brian Wichmann and David Hill (BYTE magazine, March, 1987, pp 127-8) is used. The period, given by them, is 6953607871644. That's a big number and I respect it.

---

✨ lsqrt: Integer square root ✨

SYNOPSIS

long x, y;
long lsqrt();

$y = lsqrt( $x );

DESCRIPTION

Returns a long integer square root of the long integer argument. The computation is by binary long division.

The largest possible result is lsqrt(2,147,483,647) = 46341.

If x < 0, the square root of |x| is returned, and an error message is printed.

Not me taking the absolute value of your negative input and then complaining about it out loud. That's literally what it does. 🙈

ACCURACY

An extra, roundoff, bit is computed; hence the result is the nearest integer to the actual square root.

NOTE: only DEC arithmetic is currently supported. Only DEC. That's it. That's the list. 😭

---

🫶 mtherr: Library common error handling routine 🫶

SYNOPSIS

char *fctnam;
# int code;
# int mtherr();

mtherr( $fctnam, $code );

DESCRIPTION

This routine may be called to report one of the following error conditions (in the include file mconf.h):

  Mnemonic        Value          Significance

   DOMAIN            1       argument domain error
   SING              2       function singularity
   OVERFLOW          3       overflow range error
   UNDERFLOW         4       underflow range error
   TLOSS             5       total loss of precision
   PLOSS             6       partial loss of precision
   EDOM             33       Unix domain error code
   ERANGE           34       Unix range error code

The default version of the file prints the function name, passed to it by the pointer fctnam, followed by the error condition. The display is directed to the standard output device. The routine then returns to the calling program. Users may wish to modify the program to abort by calling exit() under severe error conditions such as domain errors.

Since all error conditions pass control to this function, the display may be easily changed, eliminated, or directed to an error logging device. Customize your error UI, bestie, it's your world. 💅

SEE ALSO: mconf.h

---

🎀 nbdtr: Negative binomial distribution 🎀

SYNOPSIS

# int k, n;
# double p, y, nbdtr();

$y = nbdtr( $k, $n, $p );

DESCRIPTION

Returns the sum of the terms 0 through k of the negative binomial distribution:

  k
  --  ( n+j-1 )   n      j
  >   (       )  p  (1-p)
  --  (   j   )
 j=0

In a sequence of Bernoulli trials, this is the probability that k or fewer failures precede the *n*th success. Cute lil probability moment. ✨

The terms are not computed individually; instead the incomplete beta integral is employed, according to the formula

y = nbdtr( k, n, p ) = incbet( n, k+1, p ).

The arguments must be positive, with p ranging from 0 to 1.

ACCURACY

Tested at random points (a,b,p), with p between 0 and 1.

              a,b                     Relative error:
arithmetic  domain     # trials      peak         rms
   IEEE     0,100       100000      1.7e-13     8.8e-15

See also incbet.c.

---

🦋 nbdtrc: Complemented negative binomial distribution 🦋

SYNOPSIS

# int k, n;
# double p, y, nbdtrc();

$y = nbdtrc( $k, $n, $p );

DESCRIPTION

Returns the sum of the terms k+1 to infinity of the negative binomial distribution:

              inf
              --  ( n+j-1 )   n      j
              >   (       )  p  (1-p)
              --  (   j   )
             j=k+1

And bestie, here's the thing — we do NOT sit there adding up every single term like some kind of feral for-loop. Nope! We use the incomplete beta integral instead, per this formula right here:

y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).

Both arguments have to be positive, and p lives between 0 and 1. That's non-negotiable 💅

🎀 ACCURACY

Tested at random points (a,b,p), with p between 0 and 1.

              a,b                     Relative error:
arithmetic  domain     # trials      peak         rms
   IEEE     0,100       100000      1.7e-13     8.8e-15

See also incbet.c.

---

🦋 nbdtrc: Complemented negative binomial distribution

SYNOPSIS

# int k, n;
# double p, y, nbdtrc();

$y = nbdtrc( $k, $n, $p );

DESCRIPTION

Takes you straight to the sum of terms k+1 all the way out to infinity of the negative binomial distribution. Infinity!! Like, the whole tail end of it:

  inf
  --  ( n+j-1 )   n      j
  >   (       )  p  (1-p)
  --  (   j   )
 j=k+1

And again — no individual term crunching. We let the incomplete beta integral do the heavy lifting, according to the formula

y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).

Arguments: positive, always. And p ranges from 0 to 1. Periodt.

ACCURACY

See incbet.c.

---

✨ nbdtri: Functional inverse of negative binomial distribution

SYNOPSIS

# int k, n;
# double p, y, nbdtri();

$p = nbdtri( $k, $n, $y );

DESCRIPTION

It finds the argument p such that nbdtr(k,n,p) is equal to y. That's it, that's the whole vibe — reverse engineering the p you'd need.

ACCURACY

Tested at random points (a,b,y), with y between 0 and 1.

              a,b                     Relative error:
arithmetic  domain     # trials      peak         rms
   IEEE     0,100       100000      1.5e-14     8.5e-16

See also incbi.c.

---

🌸 ndtr: Normal distribution function

SYNOPSIS

# double x, y, ndtr();

$y = ndtr( $x );

DESCRIPTION

Returns the area under the Gaussian probability density function, integrated from minus infinity to x:

                           x
                            -
                  1        | |          2
   ndtr(x)  = ---------    |    exp( - t /2 ) dt
              sqrt(2pi)  | |
                          -
                         -inf.

            =  ( 1 + erf(z) ) / 2

where z = x/sqrt(2). Computed via the functions erf and erfc — teamwork fr fr 🫶

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC      -13,0         8000       2.1e-15     4.8e-16
   IEEE     -13,0        30000       3.4e-14     6.7e-15

ERROR MESSAGES

  message         condition         value returned
erfc underflow    x > 37.519379347       0.0

---

💜 erf: Error function

SYNOPSIS

# double x, y, erf();

$y = erf( $x );

DESCRIPTION

The integral is

                          x
                           -
                2         | |          2
  erf(x)  =  --------     |    exp( - t  ) dt.
             sqrt(pi)   | |
                         -
                          0

The magnitude of x is limited to 9.231948545 for DEC arithmetic; 1 or -1 is returned outside this range. (So DEC girlies, know your limits 😭)

For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise erf(x) = 1 - erfc(x).

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       0,1         14000       4.7e-17     1.5e-17
   IEEE      0,1         30000       3.7e-16     1.0e-16

---

🌈 erfc: Complementary error function

SYNOPSIS

# double x, y, erfc();

$y = erfc( $x );

DESCRIPTION

 1 - erf(x) =

                                  inf.
                                    -
                         2         | |          2
          erfc(x)  =  --------     |    exp( - t  ) dt
                      sqrt(pi)   | |
                                  -
                                   x

For small x, erfc(x) = 1 - erf(x); otherwise rational approximations are computed. She adapts 💅

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       0, 9.2319   12000       5.1e-16     1.2e-16
   IEEE      0,26.6417   30000       5.7e-14     1.5e-14

ERROR MESSAGES

  message         condition              value returned
erfc underflow    x > 9.231948545 (DEC)       0.0

---

⭐ ndtri: Inverse of Normal distribution function

SYNOPSIS

# double x, y, ndtri();

$x = ndtri( $y );

DESCRIPTION

Returns the argument, x, for which the area under the Gaussian probability density function (integrated from minus infinity to x) is equal to y. It's the glow-up reverse of the last one 🔥

For small arguments 0 < y < exp(-2), the program computes

z = sqrt( -2.0 * log(y) );

then the approximation is

x = z - log(z)/z  - (1/z) P(1/z) / Q(1/z).

There are two rational functions P/Q, one for 0 < y < exp(-32) and the other for y up to exp(-2). For larger arguments,

w = y - 0.5, and  x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).

ACCURACY

                     Relative error:
arithmetic   domain        # trials      peak         rms
   DEC      0.125, 1         5500       9.5e-17     2.1e-17
   DEC      6e-39, 0.135     3500       5.7e-17     1.3e-17
   IEEE     0.125, 1        20000       7.2e-16     1.3e-16
   IEEE     3e-308, 0.135   50000       4.6e-16     9.8e-17

ERROR MESSAGES

  message         condition    value returned
ndtri domain       x <= 0        -MAXNUM
ndtri domain       x >= 1         MAXNUM

---

🍓 pdtr: Poisson distribution

SYNOPSIS

# int k;
# double m, y, pdtr();

$y = pdtr( $k, $m );

DESCRIPTION

Returns the sum of the first k terms of the Poisson distribution:

  k         j
  --   -m  m
  >   e    --
  --       j!
 j=0

The terms are not summed directly; instead the incomplete gamma integral is employed, according to the relation

y = pdtr( k, m ) = igamc( k+1, m ).

The arguments must both be positive. Both. Positive. Every time 🙏

ACCURACY

See igamc().

---

🦄 pdtrc: Complemented poisson distribution

SYNOPSIS

# int k;
# double m, y, pdtrc();

$y = pdtrc( $k, $m );

DESCRIPTION

Returns the sum of the terms k+1 to infinity of the Poisson distribution:

 inf.       j
  --   -m  m
  >   e    --
  --       j!
 j=k+1

The terms are not summed directly; instead the incomplete gamma integral is employed, according to the formula

y = pdtrc( k, m ) = igam( k+1, m ).

The arguments must both be positive, bestie, I cannot stress this enough.

ACCURACY

See igam.c.

---

💖 pdtri: Inverse Poisson distribution

SYNOPSIS

# int k;
# double m, y, pdtr();

$m = pdtri( $k, $y );

DESCRIPTION

Finds the Poisson variable x such that the integral from 0 to x of the Poisson density is equal to the given probability y.

This is accomplished using the inverse gamma integral function and the relation

   m = igami( k+1, y ).

ACCURACY

See igami.c.

ERROR MESSAGES

  message         condition      value returned
pdtri domain    y < 0 or y >= 1       0.0
                    k < 0

---

🧚 pow: Power function

SYNOPSIS

# double x, y, z, pow();

$z = pow( $x, $y );

DESCRIPTION

Computes x raised to the yth power. Analytically,

     x**y  =  exp( y log(x) ).

Following Cody and Waite, this program uses a lookup table of 2**-i/16 and pseudo extended precision arithmetic to obtain an extra three bits of accuracy in both the logarithm and the exponential. Extra bits? She ate 👑

ACCURACY

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE     -26,26       30000      4.2e-16      7.7e-17
   DEC      -26,26       60000      4.8e-17      9.1e-18
1/26 < x < 26, with log(x) uniformly distributed.
-26 < y < 26, y uniformly distributed.
   IEEE     0,8700       30000      1.5e-14      2.1e-15
0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.

ERROR MESSAGES

  message         condition      value returned
pow overflow     x**y > MAXNUM      INFINITY
pow underflow   x**y < 1/MAXNUM       0.0
pow domain      x<0 and y noninteger  0.0

---

🔥 powi: Real raised to integer power

SYNOPSIS

# double x, y, powi();
# int n;

$y = powi( $x, $n );

DESCRIPTION

Returns argument x raised to the nth power. The routine efficiently decomposes n as a sum of powers of two. The desired power is a product of two-to-the-kth powers of x. Thus to compute the 32767 power of x requires 28 multiplications instead of 32767 multiplications.

Twenty. Eight. Multiplications. Instead of thirty-two thousand seven hundred sixty-seven. The efficiency is giving main character energy, no cap 💅✨

ACCURACY

                     Relative error:
arithmetic   x domain   n domain  # trials      peak         rms
   DEC       .04,26     -26,26    100000       2.7e-16     4.3e-17
   IEEE      .04,26     -26,26     50000       2.0e-15     3.8e-16
   IEEE        1,2    -1022,1023   50000       8.6e-14     1.6e-14

Returns MAXNUM on overflow, zero on underflow.

---

👑 psi: Psi (digamma) function

SYNOPSIS

# double x, y, psi();

$y = psi( $x );

DESCRIPTION

             d      -
  psi(x)  =  -- ln | (x)
             dx

is the logarithmic derivative of the gamma function. For integer x,

                          n-1
                           -
        psi(n) = -EUL  +   >  1/k.
                           -
                          k=1

This formula is used for 0 < n <= 10. If x is negative, it is transformed to a positive argument by the reflection formula

psi(1-x) = psi(x) + pi cot(pi x).

For general positive x, the argument is made greater than 10 using the recurrence

psi(x+1) = psi(x) + 1/x.

Then the following asymptotic expansion is applied:

...

                                      inf.   B
                                       -      2k
            psi(x) = log(x) - 1/2x -   >   -------
                                       -        2k
                                      k=1   2k x

and those B2k cuties are the Bernoulli numbers, as always 💖

ACCURACY 📊

okay so heads up: this one is relative error, except it's absolute when |psi| < 1. literally read the fine print bestie.

• DEC, domain 0,30, 2500 trials — peak 1.7e-16, rms 2.0e-17
• IEEE, domain 0,30, 30000 trials — peak 1.3e-15, rms 1.4e-16
• IEEE, domain -30,0, 40000 trials — peak 1.5e-15, rms 2.2e-16

ERROR MESSAGES 💀

• message psi singularity — condition: x integer <=0 — value returned: MAXNUM

---

🦋 rgamma — reciprocal gamma function 🦋

SYNOPSIS ✨

# double x, y, rgamma();

$y = rgamma( $x );

DESCRIPTION 💖

It hands you one divided by the gamma function of the argument. that's the whole vibe.

So how does it do it? Chebyshev expansion living in the interval [0,1], that's how. Then range reduction by recurrence for arguments between -34.034 and +34.84425627277176174 (yes that precision, they were not playing 😭).

Go outside that range with a positive argument? You get 1/MAXNUM back. Arguments less than -34.034 get the cosecant reflection formula instead, and logarithms are employed so we don't get unnecessary overflow. iconic, honestly.

The reciprocal gamma function has no singularities — queen never breaks — but overflow and underflow can still happen for large arguments, and those conditions return either MAXNUM or 1/MAXNUM, with the appropriate sign. no cap.

ACCURACY 🔥

Relative error, always:

• DEC, domain -30,+30, 4000 trials — peak 1.2e-16, rms 1.8e-17
• IEEE, domain -30,+30, 30000 trials — peak 1.1e-15, rms 2.0e-16

For arguments less than -34.034, the peak error is on the order of 5e-15 (DEC), unless there's overflow or underflow. 💅

---

🎀 round — round double to nearest or even integer valued double 🎀

SYNOPSIS ✨

# double x, y, round();

$y = round( $x );

DESCRIPTION 💖

Returns the nearest integer to x as a double precision floating point result. and if x ends in 0.5 exactly? the nearest even integer is chosen. banker's rounding ate, periodt.

ACCURACY 📊

Real talk: if x is greater than 1/(2*MACHEP), its closest machine representation is already an integer, so the rounding just... doesn't change anything. it's giving nothing to do here 💀

---

🌸 shichi — hyperbolic sine and cosine integrals 🌸

SYNOPSIS ✨

# double x, Chi, Shi, shichi();

($flag, $Shi, $Chi) = shichi( $x );

DESCRIPTION 💖

Approximates these integrals:

                           x
                           -
                          | |   cosh t - 1
  Chi(x) = eul + ln x +   |    -----------  dt,
                        | |          t
                         -
                         0

              x
              -
             | |  sinh t
  Shi(x) =   |    ------  dt
           | |       t
            -
            0

where eul = 0.57721566490153286061 — that's Euler's constant, bestie, write it down.

The integrals are evaluated by power series for x < 8, and by Chebyshev expansions for x between 8 and 88. For large x, both functions approach exp(x)/2x. Arguments greater than 88 in magnitude? You get MAXNUM. bye ✋

ACCURACY 📊

Test interval 0 to 88.

Relative error:

• DEC, function Shi, 3000 trials — peak 9.1e-17
• IEEE, function Shi, 30000 trials — peak 6.9e-16, rms 1.6e-16

Absolute error, except relative when |Chi| > 1:

• DEC, function Chi, 2500 trials — peak 9.3e-17
• IEEE, function Chi, 30000 trials — peak 8.4e-16, rms 1.4e-16

---

💜 sici — sine and cosine integrals 💜

SYNOPSIS ✨

# double x, Ci, Si, sici();

($flag, $Si, $Ci) = sici( $x );

DESCRIPTION 💖

Evaluates:

                             x
                             -
                            |  cos t - 1
  Ci(x) = eul + ln x +      |  --------- dt,
                            |      t
                           -
                            0
                x
                -
               |  sin t
  Si(x) =      |  ----- dt
               |    t
              -
               0

where eul = 0.57721566490153286061 — Euler's constant, same icon as always 👑

The integrals are approximated by rational functions. And for x > 8, two auxiliary functions f(x) and g(x) get pulled in such that:

Ci(x) = f(x) sin(x) - g(x) cos(x)
Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)

ACCURACY 🔥

Test interval = [0,50].

Absolute error, except relative when > 1:

• IEEE, function Si, 30000 trials — peak 4.4e-16, rms 7.3e-17
• IEEE, function Ci, 30000 trials — peak 6.9e-16, rms 5.1e-17
• DEC, function Si, 5000 trials — peak 4.4e-17, rms 9.0e-18
• DEC, function Ci, 5300 trials — peak 7.9e-17, rms 5.2e-18

---

⭐ sin — circular sine ⭐

SYNOPSIS ✨

# double x, y, sin();

$y = sin( $x );

DESCRIPTION 💖

Range reduction goes into intervals of pi/4. The reduction error is nearly eliminated by contriving an extended precision modular arithmetic — yeah, they really said that, and honestly? slay.

Two polynomial approximating functions are employed:

• between 0 and pi/4, the sine is approximated by x + x**3 P(x**2).
• between pi/4 and pi/2, the cosine is represented as 1 - x**2 Q(x**2).

ACCURACY 📊

Relative error:

• DEC, domain 0, 10, 150000 trials — peak 3.0e-17, rms 7.8e-18
• IEEE, domain -1.07e9,+1.07e9, 130000 trials — peak 2.1e-16, rms 5.4e-17

ERROR MESSAGES 💀

• message sin total loss — condition x > 1.073741824e9 — value returned 0.0

⚠️ bestie listen. Partial loss of accuracy starts creeping in at x = 2**30 = 1.074e9. And it's not gradual — it jumps suddenly to about 1 part in 10e7. Results may be meaningless for x > 2**49 = 5.6e14. The routine as implemented flags a TLOSS error for x > 2**30 and returns 0.0.

---

🌈 cos — circular cosine 🌈

SYNOPSIS ✨

# double x, y, cos();

$y = cos( $x );

DESCRIPTION 💖

Range reduction goes into intervals of pi/4, and again, the reduction error is nearly eliminated by contriving an extended precision modular arithmetic. obsessed with the technique fr fr.

Two polynomial approximating functions are employed:

• between 0 and pi/4, the cosine is approximated by 1 - x**2 Q(x**2).
• between pi/4 and pi/2, the sine is represented as x + x**3 P(x**2).

ACCURACY 📊

Relative error:

• IEEE, domain -1.07e9,+1.07e9, 130000 trials — peak 2.1e-16, rms 5.4e-17
• DEC, domain 0,+1.07e9, 17000 trials — peak 3.0e-17, rms 7.2e-18

---

🍓 sindg — circular sine of angle in degrees 🍓

SYNOPSIS ✨

# double x, y, sindg();

$y = sindg( $x );

DESCRIPTION 💖

Range reduction is into intervals of 45 degrees. (degrees!! not radians!! she's a different girl ok 😌)

Two polynomial approximating functions are employed:

• between 0 and pi/4, the sine is approximated by x + x**3 P(x**2).
• between pi/4 and pi/2, the cosine is represented as 1 - x**2 P(x**2).

ACCURACY 📊

Relative error:

• DEC, domain +-1000, 3100 trials — peak 3.3e-17, rms 9.0e-18
• IEEE, domain +-1000, 30000 trials — peak 2.3e-16, rms 5.6e-17

ERROR MESSAGES 💀

• message sindg total loss — condition x > 8.0e14 (DEC) or x > 1.0e14 (IEEE) — value returned 0.0

---

🧚 cosdg — circular cosine of angle in degrees 🧚

SYNOPSIS ✨

# double x, y, cosdg();

$y = cosdg( $x );

DESCRIPTION 💖

Range reduction is into intervals of 45 degrees, same energy as her sister up there.

Two polynomial approximating functions are employed:

• between 0 and pi/4, the cosine is approximated by 1 - x**2 P(x**2).
• between pi/4 and pi/2, the sine is represented as x + x**3 P(x**2).

ACCURACY 📊

Relative error:

• DEC, domain +-1000, 3400 trials — peak 3.5e-17, rms 9.1e-18
• IEEE, domain +-1000, 30000 trials — peak 2.1e-16, rms 5.7e-17

See also sin(). 👀

---

💖 sinh — hyperbolic sine 💖

SYNOPSIS ✨

# double x, y, sinh();

$y = sinh( $x );

DESCRIPTION 💖

Returns the hyperbolic sine of an argument in the range MINLOG to MAXLOG.

The range is partitioned into two segments, because one size never fits all:

• if |x| <= 1, a rational function of the form x + x**3 P(x)/Q(x) is employed.
• otherwise, the calculation is sinh(x) = ( exp(x) - exp(-x) )/2.

ACCURACY 📊

Relative error:

• DEC, domain +- 88, 50000 trials — peak 4.0e-17, rms 7.7e-18
• IEEE, domain +-MAXLOG, 30000 trials — peak 2.6e-16, rms 5.7e-17

---

👑 spence — dilogarithm 👑

SYNOPSIS ✨

# double x, y, spence();

$y = spence( $x );

DESCRIPTION 💖

Computes this integral:

                   x
                   -
                  | | log t
spence(x)  =  -   |   ----- dt
                | |   t - 1
                 -
                 1

for x >= 0. A rational approximation gives the integral in the interval (0.5, 1.5), and transformation formulas for 1/x and 1-x are employed outside the basic expansion range. shapeshifting queen 💅

ACCURACY 📊

Relative error:

• IEEE, domain 0,4, 30000 trials — peak 3.9e-15, rms 5.4e-16
• DEC, domain 0,4, 3000 trials — peak 2.5e-16, rms 4.5e-17

---

🌟 sqrt — square root 🌟

SYNOPSIS ✨

# double x, y, sqrt();

$y = sqrt( $x );

DESCRIPTION 💖

Returns the square root of x. classic, beloved, no notes.

Range reduction involves isolating the power of two of the argument, then using a polynomial approximation to grab a rough value for the square root. After that, Heron's iteration is used three times to converge to an accurate value. three times and she's done, that's the confidence we stan 🫶

ACCURACY ✨

okay bestie, the sqrt accuracy receipts are coming in hot right below 👇

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC       0, 10       60000       2.1e-17     7.9e-18
               IEEE      0,1.7e308   30000       1.7e-16     6.3e-17

ERROR MESSAGES:

  message         condition      value returned
sqrt domain        x < 0            0.0

okay so that's the tail end of that one — moving on!! 💅

---

stdtr</p> — Student's t distribution 🎀

SYNOPSIS:

# double t, stdtr();
short k;

$y = stdtr( $k, $t );

DESCRIPTION:

So it computes the integral from minus infinity to t of the Student t distribution with integer k > 0 degrees of freedom, and bestie, the formula is giving:

                                                 t
                                                 -
                                                | |
                         -                      |         2   -(k+1)/2
                        | ( (k+1)/2 )           |  (     x   )
                  ----------------------        |  ( 1 + --- )        dx
                                -               |  (      k  )
                  sqrt( k pi ) | ( k/2 )        |
                                              | |
                                               -
                                              -inf.

Relation to incomplete beta integral:

       1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
where
       z = k/(k + t**2).

For t < -2, this is the method of computation — that's how we do it. For higher t, a direct method is derived from integration by parts. Since the function is symmetric about t=0, the area under the right tail of the density is found by calling the function with -t instead of t. lowkey elegant ✨

ACCURACY:

Tested at random 1 <= k <= 25. The "domain" refers to t.

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE     -100,-2      50000       5.9e-15     1.4e-15
               IEEE     -2,100      500000       2.7e-15     4.9e-17

---

stdtri — Functional inverse of Student's t distribution 🦋

SYNOPSIS:

# double p, t, stdtri();
# int k;

$t = stdtri( $k, $p );

DESCRIPTION:

Given probability p, it finds the argument t such that stdtr(k,t) is equal to p. That's it, that's the whole vibe. 🫶

ACCURACY:

Tested at random 1 <= k <= 100. The "domain" refers to p:

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE    .001,.999     25000       5.7e-15     8.0e-16
               IEEE    10^-6,.001    25000       2.0e-12     2.9e-14

---

struve — Struve function 🌸

SYNOPSIS:

# double v, x, y, struve();

$y = struve( $v, $x );

DESCRIPTION:

Computes the Struve function Hv(x) of order v, argument x. Heads up though: negative x is rejected unless v is an integer. no exceptions bestie.

ACCURACY:

Not accurately characterized, but spot checked against tables. 💀

---

plancki — Integral of Planck's black body radiation formula 🔥

SYNOPSIS:

# double lambda, T, y, plancki()

$y = plancki( $lambda, $T );

DESCRIPTION:

Evaluates the definite integral, from wavelength 0 to lambda, of Planck's radiation formula:

                                  -5
                        c1  lambda
                 E =  ------------------
                        c2/(lambda T)
                       e             - 1

Physical constants c1 = 3.7417749e-16 and c2 = 0.01438769 are built in to the function program. They are scaled to provide a result in watts per square meter. Argument T represents temperature in degrees Kelvin; lambda is wavelength in meters.

The integral is expressed in closed form, in terms of polylogarithms (see polylog.c). ✨

The total area under the curve is

     (-1/8) (42 zeta(4) - 12 pi^2 zeta(2) + pi^4 ) c1 (T/c2)^4
      = (pi^4 / 15)  c1 (T/c2)^4
      =  5.6705032e-8 T^4

where sigma = 5.6705032e-8 W m^2 K^-4 is the Stefan-Boltzmann constant. periodt. 👑

ACCURACY:

The left tail of the function experiences some relative error amplification in computing the dominant term exp(-c2/(lambda T)). For the right-hand tail see planckc, below.

                             Relative error.
          The domain refers to lambda T / c2.
            arithmetic   domain     # trials      peak         rms
               IEEE      0.1, 10      50000      7.1e-15     5.4e-16

---

polylog — polylogarithm function 🌈

SYNOPSIS:

# double x, y, polylog();
# int n;

    $y = polylog( $n, $x );

The polylogarithm of order n is defined by the series

              inf   k
               -   x
  Li (x)  =    >   ---  .
    n          -     n
              k=1   k

For x = 1,

               inf
                -    1
   Li (1)  =    >   ---   =  Riemann zeta function (n)  .
     n          -     n
               k=1   k

When n = 2, the function is the dilogarithm, related to Spence's integral:

                 x                      1-x
                 -                        -
                | |  -ln(1-t)            | |  ln t
   Li (x)  =    |    -------- dt    =    |    ------ dt    =   spence(1-x) .
     2        | |       t              | |    1 - t
               -                        -
                0                        1

ACCURACY:

                      Relative error:
 arithmetic   domain   n   # trials      peak         rms
    IEEE      0, 1     2     50000      6.2e-16     8.0e-17
    IEEE      0, 1     3    100000      2.5e-16     6.6e-17
    IEEE      0, 1     4     30000      1.7e-16     4.9e-17
    IEEE      0, 1     5     30000      5.1e-16     7.8e-17

---

bernum</bernum> — Bernoulli numbers 🍓

SYNOPSIS:

($num, $den) = bernum( $n);
($num_array, $den_array) = bernum();

DESCRIPTION:

This calculates the Bernoulli numbers, up to 30th order. If called with an integer argument, the numerator and denominator of that Bernoulli number is returned; if called with no argument, two array references representing the numerator and denominators of the first 30 Bernoulli numbers are returned. either way you ate 💅

---

simpson — Simpson's rule to find an integral 🧚

SYNOPSIS:

$result = simpson(\&fun, $a, $b, $abs_err, $rel_err, $nmax);

sub fun {
   my $x = shift;
   return cos($x)*exp($x);
}

DESCRIPTION:

This evaluates the area under the graph of a function, represented in a subroutine, from $a to $b, using an 8-point Newton-Cotes formula. The routine divides up the interval into equal segments, evaluates the integral, then compares that to the result with double the number of segments. If the two results agree, to within an absolute error $abs_err or a relative error $rel_err, the result is returned; otherwise, the number of segments is doubled again, and the results compared. This continues until the desired accuracy is attained, or until the maximum number of iterations $nmax is reached. it literally just keeps going until it's happy or it gives up 😭

---

vecang — angle between two vectors 💖

SYNOPSIS:

# double p[3], q[3], vecang();

   $y = vecang( $p, $q );

DESCRIPTION:

For two vectors p, q, the angle A between them is given by

     p.q / (|p| |q|)  = cos A  .

where "." represents inner product, "|x|" the length of vector x. If the angle is small, an expression in sin A is preferred. Set r = q - p. Then

    p.q = p.p + p.r ,

    |p|^2 = p.p ,

    |q|^2 = p.p + 2 p.r + r.r ,

                 p.p^2 + 2 p.p p.r + p.r^2
    cos^2 A  =  ----------------------------
                   p.p (p.p + 2 p.r + r.r)

                 p.p + 2 p.r + p.r^2 / p.p
             =  --------------------------- ,
                    p.p + 2 p.r + r.r

    sin^2 A  =  1 - cos^2 A

                  r.r - p.r^2 / p.p
             =  --------------------
                 p.p + 2 p.r + r.r

             =   (r.r - p.r^2 / p.p) / q.q  .

no cap, that algebra is iconic. 🔥

ACCURACY:

                                 Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE      -1, 1        10^6       1.7e-16     4.2e-17

---

onef2 — Hypergeometric function 1F2 ⭐

SYNOPSIS:

# double a, b, c, x, value;

# double *err;

($value, $err) = onef2( $a, $b, $c, $x)

ACCURACY:

Not accurately characterized, but spot checked against tables. 🙈

---

threef0 — Hypergeometric function 3F0 🌸

SYNOPSIS:

# double a, b, c, x, value;

# double *err;

($value, $err) = threef0( $a, $b, $c, $x )

ACCURACY:

Not accurately characterized, but spot checked against tables. (yes, same energy as the last one, they're twins fr fr)

---

yv — Bessel function Yv with noninteger v 💜

SYNOPSIS:

# double v, x;

# double yv( v, x );

$y = yv( $v, $x );

ACCURACY:

Not accurately characterized, but spot checked against tables. 🦋

---

tan — Circular tangent 🎀

SYNOPSIS:

# double x, y, tan();

$y = tan( $x );

DESCRIPTION:

Returns the circular tangent of the radian argument x.

Range reduction is modulo pi/4. A rational function

      x + x**3 P(x**2)/Q(x**2)

is employed in the basic interval [0, pi/4].

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               DEC      +-1.07e9      44000      4.1e-17     1.0e-17
               IEEE     +-1.07e9      30000      2.9e-16     8.1e-17

ERROR MESSAGES:

  message         condition          value returned
tan total loss   x > 1.073741824e9     0.0

bestie, if you're past x > 1.073741824e9 you get 0.0 back and that is the code going "total loss, i cannot". rip 😭

---

cot — Circular cotangent 🌊

SYNOPSIS:

# double x, y, cot();

$y = cot( $x );

DESCRIPTION:

Returns the circular cotangent of the radian argument x.

Range reduction is modulo pi/4. A rational function

      x + x**3 P(x**2)/Q(x**2)

is employed in the basic interval [0, pi/4].

ACCURACY:

                             Relative error:
            arithmetic   domain     # trials      peak         rms
               IEEE     +-1.07e9      30000      2.9e-16     8.2e-17

ERROR MESSAGES:

  message         condition          value returned
cot total loss   x > 1.073741824e9       0.0
cot singularity  x = 0                  INFINITY

two ways to go wrong here, babes: too far out and you get 0.0 (total loss), or land exactly on x = 0 and it's a singularity handing you INFINITY. 💀

---

tandg — Circular tangent of argument in degrees 👑

SYNOPSIS:

# double x, y, tandg();

🌸 tandg — circular tangent but make it degrees

# double x, y, tandg();

$y = tandg( $x );

Bestie, this one hands you the circular tangent of the argument x, except x is in degrees 💅. Yes, degrees. We're not doing radians today, we're doing degrees, and we're doing them gorgeously.

The range reduction is modulo pi/4, and then in the basic interval [0, pi/4] it flexes a rational function that looks like this:

x + x**3 P(x**2)/Q(x**2)

ACCURACY — the receipts, the way the numbers peak and rms and just serve:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC      0,10          8000      3.4e-17      1.2e-17
   IEEE     0,10         30000      3.2e-16      8.4e-17

ERROR MESSAGES — she will yell at you, and here's exactly how:

  message         condition          value returned
tandg total loss   x > 8.0e14 (DEC)      0.0
                   x > 1.0e14 (IEEE)
tandg singularity  x = 180 k  +  90     MAXNUM

So if x is a singularity at 180 k + 90, you get MAXNUM. And if you push past the "total loss" threshold you just get a whole 0.0 — like, thanks babe, very helpful 🙈.

---

🎀 cotdg — circular cotangent, still in degrees

# double x, y, cotdg();

$y = cotdg( $x );

Same energy, opposite vibe: this returns the circular cotangent of the argument x in degrees.

Range reduction is modulo pi/4 (yes, again, consistent queen), and the rational function used in the basic interval [0, pi/4] is:

x + x**3 P(x**2)/Q(x**2)

ERROR MESSAGES — no accuracy table for this one, she's mysterious like that ✨:

  message         condition          value returned
cotdg total loss   x > 8.0e14 (DEC)      0.0
                   x > 1.0e14 (IEEE)
cotdg singularity  x = 180 k            MAXNUM

Note the singularity is x = 180 k here — no + 90, because cotangent is a different girl 💀.

---

🦋 tanh — hyperbolic tangent, the glow-up

# double x, y, tanh();

$y = tanh( $x );

Returns the hyperbolic tangent of the argument, over the range MINLOG to MAXLOG.

How it's actually computed — because ikyk y'all wanna know:

• A rational function is used for |x| < 0.625.
• The form x + x**3 P(x)/Q(x) of Cody _& Waite is employed. (Yes, an underscore in the name. Not my typo, that's just how it's written. We respect the lore 📜.)
• Otherwise:

tanh(x) = sinh(x)/cosh(x) = 1  -  2/(exp(2x) + 1).

ACCURACY:

                     Relative error:
arithmetic   domain     # trials      peak         rms
   DEC       -2,2        50000       3.3e-17     6.4e-18
   IEEE      -2,2        30000       2.5e-16     5.8e-17

---

💖 unity — relative error approximations for arguments that are close to one

This is the "we're near 1 and floating point is about to embarrass us" support group 🫶. Three functions, all up in here:

#    log1p(x) = log(1+x)

$y = log1p( $x );

#    expm1(x) = exp(x) - 1

$y = expm1( $x );

#    cosm1(x) = cos(x) - 1

$y = cosm1( $x );

That's it, that's the whole entry. No description paragraph, no accuracy table, she just showed up, dropped three exact function signatures, and left 👑.

---

🌈 yn — Bessel function of the second kind, integer order

# double x, y, yn();
# int n;

$y = yn( $n, $x );

Returns the Bessel function of order n, where n is an integer that's possibly negative (because we're inclusive, we don't gatekeep 💅).

How she gets there:

• The function is evaluated by forward recurrence on n, starting with values computed by the routines y0() and y1().
• If n = 0 or 1, the routine for y0 or y1 is called directly — no recurrence detour, straight to the source ✨.

ACCURACY (absolute error, except relative when y > 1):

                     Absolute error, except relative
                     when y > 1:
arithmetic   domain     # trials      peak         rms
   DEC       0, 30        2200       2.9e-16     5.3e-17
   IEEE      0, 30       30000       3.4e-15     4.3e-16

ERROR MESSAGES:

  message         condition      value returned
yn singularity   x = 0              MAXNUM
yn overflow                         MAXNUM

Also: spot checked against tables for x, n between 0 and 100. That's the peer review, bestie 🔥.

---

⭐ zeta — Riemann zeta function of two arguments

# double x, q, y, zeta();

$y = zeta( $x, $q );

The definition, in all her glory:

                inf.
                 -        -x
  zeta(x,q)  =   >   (k+q)
                 -
                k=0

...where x > 1 and q is not a negative integer or zero. (Read that constraint again, it matters 👀.)

The Euler-Maclaurin summation formula is used to obtain the expansion:

               n
               -       -x
zeta(x,q)  =   >  (k+q)
               -
              k=1

          1-x                 inf.  B   x(x+1)...(x+2j)
     (n+q)           1         -     2j
 +  ---------  -  -------  +   >    --------------------
       x-1              x      -                   x+2j+1
                  2(n+q)      j=1       (2j)! (n+q)

...where the B2j are Bernoulli numbers. And note that (see zetac.c):

zeta(x,1) = zetac(x) + 1.

ACCURACY: — crickets — the table is empty, there is literally nothing under it. Not every girl gives you numbers and that's okay 🌸.

REFERENCE:

Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, Series, and Products, p. 1073; Academic Press, 1980.

---

💜 zetac — Riemann zeta function (the single-argument one)

# double x, y, zetac();

$y = zetac( $x );

The definition:

               inf.
                -    -x
  zetac(x)  =   >   k   ,   x > 1,
                -
               k=2

...and she's related to the Riemann zeta function by:

       Riemann zeta(x) = zetac(x) + 1.

Important deets, do not skip:

• Extension of the function definition for x < 1 is implemented.
• Zero is returned for x > log2(MAXNUM).
• An overflow error may occur for large negative x, due to the gamma function in the reflection formula. So like, be careful down there 😭.

ACCURACY:

Tabulated values have full machine accuracy.

                     Relative error:
arithmetic   domain     # trials      peak         rms
   IEEE      1,50        10000       9.8e-16       1.3e-16
   DEC       1,50         2000       1.1e-16     1.9e-17

---

🍓 TODO

• Include more operating systems when generating mconf.h.

Just one todo. Growth mindset tho 💅✨.

---

🧚 MAINTAINER

Shlomi Fish, <http://www.shlomifish.org/>, <https://metacpan.org/author/SHLOMIF> .

(Yes, two links, he's got range 👑.)

---

💀 BUGS

Please report any on the rt.cpan.org interface:

<https://rt.cpan.org/Dist/Display.html?Queue=Math-Cephes>

Tell them everything, bestie. They want to know 🫶.

---

🌸 VERSION CONTROL

This distribution is maintained in this GitHub repository:

<https://github.com/shlomif/Math-Cephes>.

---

🎀 SEE ALSO

Okay so the whole friend group, gather round:

• For interfaces to programs which can do symbolic manipulation, see PDL, Math::Pari, and Math::ematica.
• For a command line interface to the routines of Math::Cephes, see the included pmath script.
• For a different interface to the fraction and complex number routines, see Math::Cephes::Fraction and Math::Cephes::Complex.
• For an interface to some polynomial routines, see Math::Cephes::Polynomial.
• And for some matrix routines, see Math::Cephes::Matrix.

Squad goals fr fr 🌈.

---

💖 COPYRIGHT

The C code for the Cephes Math Library is Copyright 1984, 1987, 1989, 2002 by Stephen L. Moshier, and is available at <http://www.netlib.org/cephes/>. Direct inquiries to 30 Frost Street, Cambridge, MA 02140.

The file arrays.c included here to handle passing arrays into and out of C routines comes from the PGPLOT module of Karl Glazebrook <kgb@zzoepp.aao.gov.au>.

The perl interface is copyright 2000, 2002 by Randy Kobes. This library is free software; you can redistribute it and/or modify it under the same terms as Perl itself.

Perl interface maintained by Shlomi Fish starting from 2012. All explicit or implicit copyrights on the changes are disclaimed by him.

---

And that's the whole page, babes 💅✨. This was Math::Cephes(3pm), from perl v5.42.3, dated 2026-08-19 — so yes, go touch some grass, hydrate, and come back and call zeta whenever, because she's not going anywhere 🫶🌸