Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

1 files changed, 0 insertions, 4721 deletions

diff --git a/libm/float/README.txt b/libm/float/README.txt
deleted file mode 100644
index 30a10b083..000000000
--- a/libm/float/README.txt
+++ /dev/null
@@ -1,4721 +0,0 @@
-/*							acoshf.c
- *
- *	Inverse hyperbolic cosine
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, acoshf();
- *
- * y = acoshf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns inverse hyperbolic cosine of argument.
- *
- * If 1 <= x < 1.5, a polynomial approximation
- *
- *	sqrt(z) * P(z)
- *
- * where z = x-1, is used.  Otherwise,
- *
- * acosh(x)  =  log( x + sqrt( (x-1)(x+1) ).
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      1,3         100000      1.8e-7       3.9e-8
- *    IEEE      1,2000      100000                   3.0e-8
- *
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * acoshf domain      |x| < 1            0.0
- *
- */
-
-/*							airy.c
- *
- *	Airy function
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, ai, aip, bi, bip;
- * int airyf();
- *
- * airyf( x, _&ai, _&aip, _&bi, _&bip );
- *
- *
- *
- * DESCRIPTION:
- *
- * Solution of the differential equation
- *
- *	y"(x) = xy.
- *
- * The function returns the two independent solutions Ai, Bi
- * and their first derivatives Ai'(x), Bi'(x).
- *
- * Evaluation is by power series summation for small x,
- * by rational minimax approximations for large x.
- *
- *
- *
- * ACCURACY:
- * Error criterion is absolute when function <= 1, relative
- * when function > 1, except * denotes relative error criterion.
- * 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        50000       7.0e-7      1.2e-7
- * IEEE          0, 10    Ai        50000       9.9e-6*     6.8e-7*
- * IEEE        -10, 0     Ai'       50000       2.4e-6      3.5e-7
- * IEEE          0, 10    Ai'       50000       8.7e-6*     6.2e-7*
- * IEEE        -10, 10    Bi       100000       2.2e-6      2.6e-7
- * IEEE        -10, 10    Bi'       50000       2.2e-6      3.5e-7
- *
- */
-
-/*							asinf.c
- *
- *	Inverse circular sine
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, asinf();
- *
- * y = asinf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns radian angle between -pi/2 and +pi/2 whose sine is x.
- *
- * A polynomial of the form x + x**3 P(x**2)
- * is used for |x| in the interval [0, 0.5].  If |x| > 0.5 it is
- * transformed by the identity
- *
- *    asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE     -1, 1       100000       2.5e-7       5.0e-8
- *
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * asinf domain        |x| > 1           0.0
- *
- */
-/*							acosf()
- *
- *	Inverse circular cosine
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, acosf();
- *
- * y = acosf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns radian angle between -pi/2 and +pi/2 whose cosine
- * is x.
- *
- * Analytically, acos(x) = pi/2 - asin(x).  However if |x| is
- * near 1, there is cancellation error in subtracting asin(x)
- * from pi/2.  Hence 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
- *    IEEE      -1, 1      100000       1.4e-7      4.2e-8
- *
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * acosf domain        |x| > 1           0.0
- */
-
-/*							asinhf.c
- *
- *	Inverse hyperbolic sine
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, asinhf();
- *
- * y = asinhf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns inverse hyperbolic sine of 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
- *    IEEE     -3,3        100000       2.4e-7      4.1e-8
- *
- */
-
-/*							atanf.c
- *
- *	Inverse circular tangent
- *      (arctangent)
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, atanf();
- *
- * y = atanf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns radian angle between -pi/2 and +pi/2 whose tangent
- * is x.
- *
- * Range reduction is from four intervals into the interval
- * from zero to  tan( pi/8 ).  A polynomial approximates
- * the function in this basic interval.
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10, 10     100000      1.9e-7      4.1e-8
- *
- */
-/*							atan2f()
- *
- *	Quadrant correct inverse circular tangent
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, z, atan2f();
- *
- * z = atan2f( y, x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns radian angle whose tangent is y/x.
- * Define 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).
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10, 10     100000      1.9e-7      4.1e-8
- * See atan.c.
- *
- */
-
-/*							atanhf.c
- *
- *	Inverse hyperbolic tangent
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, atanhf();
- *
- * y = atanhf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns inverse hyperbolic tangent of argument in the range
- * MINLOGF to MAXLOGF.
- *
- * If |x| < 0.5, a polynomial approximation is used.
- * Otherwise,
- *        atanh(x) = 0.5 * log( (1+x)/(1-x) ).
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -1,1        100000      1.4e-7      3.1e-8
- *
- */
-
-/*							bdtrf.c
- *
- *	Binomial distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * int k, n;
- * float p, y, bdtrf();
- *
- * y = bdtrf( 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
- *
- * 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:
- *
- *        Relative error (p varies from 0 to 1):
- * arithmetic   domain     # trials      peak         rms
- *    IEEE       0,100       2000       6.9e-5      1.1e-5
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * bdtrf domain        k < 0            0.0
- *                     n < k
- *                     x < 0, x > 1
- *
- */
-/*							bdtrcf()
- *
- *	Complemented binomial distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * int k, n;
- * float p, y, bdtrcf();
- *
- * y = bdtrcf( 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
- *
- * 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:
- *
- *        Relative error (p varies from 0 to 1):
- * arithmetic   domain     # trials      peak         rms
- *    IEEE       0,100       2000       6.0e-5      1.2e-5
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * bdtrcf domain     x<0, x>1, n<k       0.0
- */
-/*							bdtrif()
- *
- *	Inverse binomial distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * int k, n;
- * float p, y, bdtrif();
- *
- * p = bdtrf( 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 ).
- *
- *
- *
- *
- * ACCURACY:
- *
- *        Relative error (p varies from 0 to 1):
- * arithmetic   domain     # trials      peak         rms
- *    IEEE       0,100       2000       3.5e-5      3.3e-6
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * bdtrif domain    k < 0, n <= k         0.0
- *                  x < 0, x > 1
- *
- */
-
-/*							betaf.c
- *
- *	Beta function
- *
- *
- *
- * SYNOPSIS:
- *
- * float a, b, y, betaf();
- *
- * y = betaf( 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
- *    IEEE       0,30       10000       4.0e-5      6.0e-6
- *    IEEE       -20,0      10000       4.9e-3      5.4e-5
- *
- * ERROR MESSAGES:
- *
- *   message         condition          value returned
- * betaf overflow   log(beta) > MAXLOG       0.0
- *                  a or b <0 integer        0.0
- *
- */
-
-/*							cbrtf.c
- *
- *	Cube root
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, cbrtf();
- *
- * y = cbrtf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns the cube root of the argument, which may be negative.
- *
- * Range reduction involves determining the power of 2 of
- * the argument.  A polynomial of degree 2 applied to the
- * mantissa, and multiplication by the cube root of 1, 2, or 4
- * approximates the root to within about 0.1%.  Then Newton's
- * iteration is used to converge to an accurate result.
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      0,1e38      100000      7.6e-8      2.7e-8
- *
- */
-
-/*							chbevlf.c
- *
- *	Evaluate Chebyshev series
- *
- *
- *
- * SYNOPSIS:
- *
- * int N;
- * float x, y, coef[N], chebevlf();
- *
- * y = chbevlf( x, coef, N );
- *
- *
- *
- * DESCRIPTION:
- *
- * Evaluates the series
- *
- *        N-1
- *         - '
- *  y  =   >   coef[i] T (x/2)
- *         -            i
- *        i=0
- *
- * of Chebyshev polynomials Ti at argument x/2.
- *
- * Coefficients are stored in reverse order, i.e. the zero
- * order term is last in the array.  Note N is the number of
- * coefficients, not the order.
- *
- * If coefficients are for the interval a to b, x must
- * have been transformed to x -> 2(2x - b - a)/(b-a) before
- * entering the routine.  This maps x from (a, b) to (-1, 1),
- * over which the Chebyshev polynomials are defined.
- *
- * If the coefficients are for the inverted interval, in
- * which (a, b) is mapped to (1/b, 1/a), the transformation
- * required is x -> 2(2ab/x - b - a)/(b-a).  If b is infinity,
- * this becomes x -> 4a/x - 1.
- *
- *
- *
- * SPEED:
- *
- * Taking advantage of the recurrence properties of the
- * Chebyshev polynomials, the routine requires one more
- * addition per loop than evaluating a nested polynomial of
- * the same degree.
- *
- */
-
-/*							chdtrf.c
- *
- *	Chi-square distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * float df, x, y, chdtrf();
- *
- * y = chdtrf( df, 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.
- *
- *
- *                                  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 ).
- *
- *
- * The arguments must both be positive.
- *
- *
- *
- * ACCURACY:
- *
- *        Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE       0,100       5000       3.2e-5      5.0e-6
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * chdtrf domain  x < 0 or v < 1        0.0
- */
-/*							chdtrcf()
- *
- *	Complemented Chi-square distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * float v, x, y, chdtrcf();
- *
- * y = chdtrcf( v, x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns 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 = chdtr( v, x ) = igamc( v/2.0, x/2.0 ).
- *
- *
- * The arguments must both be positive.
- *
- *
- *
- * ACCURACY:
- *
- *        Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE       0,100       5000       2.7e-5      3.2e-6
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * chdtrc domain  x < 0 or v < 1        0.0
- */
-/*							chdtrif()
- *
- *	Inverse of complemented Chi-square distribution
- *
- *
- *
- * SYNOPSIS:
- *
- * float df, x, y, chdtrif();
- *
- * x = chdtrif( df, y );
- *
- *
- *
- *
- * DESCRIPTION:
- *
- * 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.
- *
- * This is accomplished using the inverse gamma integral
- * function and the relation
- *
- *    x/2 = igami( df/2, y );
- *
- *
- *
- *
- * ACCURACY:
- *
- *        Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE       0,100       10000      2.2e-5      8.5e-7
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * chdtri domain   y < 0 or y > 1        0.0
- *                     v < 1
- *
- */
-
-/*							clogf.c
- *
- *	Complex natural logarithm
- *
- *
- *
- * SYNOPSIS:
- *
- * void clogf();
- * cmplxf z, w;
- *
- * clogf( &z, &w );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns 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
- *    IEEE      -10,+10     30000       1.9e-6       6.2e-8
- *
- * Larger relative error can be observed for z near 1 +i0.
- * In IEEE arithmetic the peak absolute error is 3.1e-7.
- *
- */
-/*							cexpf()
- *
- *	Complex exponential function
- *
- *
- *
- * SYNOPSIS:
- *
- * void cexpf();
- * cmplxf z, w;
- *
- * cexpf( &z, &w );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns the exponential of the complex argument z
- * into the complex result w.
- *
- * If
- *     z = x + iy,
- *     r = exp(x),
- *
- * then
- *
- *     w = r cos y + i r sin y.
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10,+10     30000       1.4e-7      4.5e-8
- *
- */
-/*							csinf()
- *
- *	Complex circular sine
- *
- *
- *
- * SYNOPSIS:
- *
- * void csinf();
- * cmplxf z, w;
- *
- * csinf( &z, &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
- *    IEEE      -10,+10     30000       1.9e-7      5.5e-8
- *
- */
-/*							ccosf()
- *
- *	Complex circular cosine
- *
- *
- *
- * SYNOPSIS:
- *
- * void ccosf();
- * cmplxf z, w;
- *
- * ccosf( &z, &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
- *    IEEE      -10,+10     30000       1.8e-7       5.5e-8
- */
-/*							ctanf()
- *
- *	Complex circular tangent
- *
- *
- *
- * SYNOPSIS:
- *
- * void ctanf();
- * cmplxf z, w;
- *
- * ctanf( &z, &w );
- *
- *
- *
- * DESCRIPTION:
- *
- * If
- *     z = x + iy,
- *
- * then
- *
- *           sin 2x  +  i sinh 2y
- *     w  =  --------------------.
- *            cos 2x  +  cosh 2y
- *
- * 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.
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10,+10     30000       3.3e-7       5.1e-8
- */
-/*							ccotf()
- *
- *	Complex circular cotangent
- *
- *
- *
- * SYNOPSIS:
- *
- * void ccotf();
- * cmplxf z, w;
- *
- * ccotf( &z, &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.
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10,+10     30000       3.6e-7       5.7e-8
- * Also tested by ctan * ccot = 1 + i0.
- */
-/*							casinf()
- *
- *	Complex circular arc sine
- *
- *
- *
- * SYNOPSIS:
- *
- * void casinf();
- * cmplxf z, w;
- *
- * casinf( &z, &w );
- *
- *
- *
- * DESCRIPTION:
- *
- * Inverse complex sine:
- *
- *                               2
- * w = -i clog( iz + csqrt( 1 - z ) ).
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10,+10     30000       1.1e-5      1.5e-6
- * Larger relative error can be observed for z near zero.
- *
- */
-/*							cacosf()
- *
- *	Complex circular arc cosine
- *
- *
- *
- * SYNOPSIS:
- *
- * void cacosf();
- * cmplxf z, w;
- *
- * cacosf( &z, &w );
- *
- *
- *
- * DESCRIPTION:
- *
- *
- * w = arccos z  =  PI/2 - arcsin z.
- *
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      -10,+10     30000       9.2e-6       1.2e-6
- *
- */
-/*							catan()
- *
- *	Complex circular arc tangent
- *
- *
- *
- * SYNOPSIS:
- *
- * void catan();
- * cmplxf z, w;
- *
- * catan( &z, &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
- *    IEEE      -10,+10     30000        2.3e-6      5.2e-8
- *
- */
-
-/*							cmplxf.c
- *
- *	Complex number arithmetic
- *
- *
- *
- * SYNOPSIS:
- *
- * typedef struct {
- *      float r;     real part
- *      float i;     imaginary part
- *     }cmplxf;
- *
- * cmplxf *a, *b, *c;
- *
- * caddf( a, b, c );     c = b + a
- * csubf( a, b, c );     c = b - a
- * cmulf( a, b, c );     c = b * a
- * cdivf( a, b, c );     c = b / a
- * cnegf( c );           c = -c
- * cmovf( b, c );        c = b
- *
- *
- *
- * 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:
- *
- * In DEC arithmetic, the test (1/z) * z = 1 had peak relative
- * error 3.1e-17, rms 1.2e-17.  The test (y/z) * (z/y) = 1 had
- * peak relative error 8.3e-17, rms 2.1e-17.
- *
- * Tests in the rectangle {-10,+10}:
- *                      Relative error:
- * arithmetic   function  # trials      peak         rms
- *    IEEE       cadd       30000       5.9e-8      2.6e-8
- *    IEEE       csub       30000       6.0e-8      2.6e-8
- *    IEEE       cmul       30000       1.1e-7      3.7e-8
- *    IEEE       cdiv       30000       2.1e-7      5.7e-8
- */
-
-/*							cabsf()
- *
- *	Complex absolute value
- *
- *
- *
- * SYNOPSIS:
- *
- * float cabsf();
- * cmplxf z;
- * float a;
- *
- * a = cabsf( &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
- *    IEEE      -10,+10     30000       1.2e-7      3.4e-8
- */
-/*							csqrtf()
- *
- *	Complex square root
- *
- *
- *
- * SYNOPSIS:
- *
- * void csqrtf();
- * cmplxf z, w;
- *
- * csqrtf( &z, &w );
- *
- *
- *
- * DESCRIPTION:
- *
- *
- * If z = x + iy,  r = |z|, then
- *
- *                       1/2
- * Im w  =  [ (r - x)/2 ]   ,
- *
- * Re w  =  y / 2 Im w.
- *
- *
- * Note that -w is also a square root of z.  The solution
- * reported 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
- *    IEEE      -10,+10    100000       1.8e-7       4.2e-8
- *
- */
-
-/*							coshf.c
- *
- *	Hyperbolic cosine
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, coshf();
- *
- * y = coshf( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns hyperbolic cosine of argument in the range MINLOGF to
- * MAXLOGF.
- *
- * cosh(x)  =  ( exp(x) + exp(-x) )/2.
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE     +-MAXLOGF    100000      1.2e-7      2.8e-8
- *
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * coshf overflow  |x| > MAXLOGF       MAXNUMF
- *
- *
- */
-
-/*							dawsnf.c
- *
- *	Dawson's Integral
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, dawsnf();
- *
- * y = dawsnf( 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        50000       4.4e-7      6.3e-8
- *
- *
- */
-
-/*							ellief.c
- *
- *	Incomplete elliptic integral of the second kind
- *
- *
- *
- * SYNOPSIS:
- *
- * float phi, m, y, ellief();
- *
- * y = ellief( phi, m );
- *
- *
- *
- * DESCRIPTION:
- *
- * 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.
- *
- *
- *
- * ACCURACY:
- *
- * Tested at random arguments with phi in [0, 2] and m in
- * [0, 1].
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE       0,2        10000       4.5e-7      7.4e-8
- *
- *
- */
-
-/*							ellikf.c
- *
- *	Incomplete elliptic integral of the first kind
- *
- *
- *
- * SYNOPSIS:
- *
- * float phi, m, y, ellikf();
- *
- * y = ellikf( phi, m );
- *
- *
- *
- * DESCRIPTION:
- *
- * 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.
- *
- *
- *
- *
- * ACCURACY:
- *
- * Tested at random points with phi in [0, 2] and m in
- * [0, 1].
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE      0,2         10000       2.9e-7      5.8e-8
- *
- *
- */
-
-/*							ellpef.c
- *
- *	Complete elliptic integral of the second kind
- *
- *
- *
- * SYNOPSIS:
- *
- * float m1, y, ellpef();
- *
- * y = ellpef( m1 );
- *
- *
- *
- * DESCRIPTION:
- *
- * 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).
- *
- * Though there are no singularities, the argument m1 is used
- * rather than m for compatibility with ellpk().
- *
- * E(1) = 1; E(0) = pi/2.
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    IEEE       0, 1       30000       1.1e-7      3.9e-8
- *
- *
- * ERROR MESSAGES:
- *
- *   message         condition      value returned
- * ellpef domain     x<0, x>1            0.0
- *
- */
-
-/*							ellpjf.c
- *
- *	Jacobian Elliptic Functions
- *
- *
- *
- * SYNOPSIS:
- *
- * float u, m, sn, cn, dn, phi;
- * int ellpj();
- *
- * ellpj( u, m, _&sn, _&cn, _&dn, _&phi );
- *
- *
- *
- * DESCRIPTION:
- *
- *
- * Evaluates the Jacobian ellipti