概要

# 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

説明

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.

精度

                      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

概要

# 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

説明

If

z = x + iy,

then

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

精度

                      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.

概要

# 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

説明

If

z = x + iy,

then

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

精度

                      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

概要

# 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

説明

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.

精度

                      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.

概要

# 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

説明

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.

精度

                      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.

概要

# 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

説明

Inverse complex sine:

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

精度

                      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.

概要

# 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

説明

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

精度

                      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

概要

# 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

説明

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.

精度

                      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

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. See also clog().

概要

# 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

説明

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

精度

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

概要

# 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) 

説明

casinh z = -i casin iz .

精度

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

概要

# 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

説明

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

精度

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

概要

# 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 

説明

acosh z = i acos z .

精度

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

概要

# 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  

説明

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

精度

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

概要

# 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  

説明

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

精度

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

概要

# 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  

説明

Raises complex A to the complex Zth power. Definition is per AMS55 # 4.2.8, analytically equivalent to cpow(a,z) = cexp(z clog(a)).

精度

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

概要

# 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

説明

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

精度

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
    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

概要

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

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

説明

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.

精度

                      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

概要

# 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

説明

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 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.

精度

                      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.

概要

extern double nameofconstant;

説明

This file contains a number 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.

概要

# double x, y, cosh();

$y = cosh( $x );

説明

Returns hyperbolic cosine of argument in the range MINLOG to MAXLOG. cosh(x) = ( exp(x) + exp(-x) )/2.

精度

                      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

エラーメッセージ

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

概要

# double x, y, dawsn();

$y = dawsn( $x );

説明

Approximates the 積分

                             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.

精度

                      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

概要

# double y, drand();

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

説明

Yields a random number 1.0 <= y < 2.0. 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. Versions invoked by the different arithmetic compile time options DEC, IBMPC, and MIEEE, 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. It may be used on all computers. However, the low order bits of a double precision number may not be adequately random, and may vary due to arithmetic implementation details on different computers. The other compile options 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.

概要

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

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

説明

Approximates the 積分

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

of amplitude phi and modulus m, using the arithmetic - geometric mean algorithm.

精度

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

概要

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

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

説明

Approximates the 積分

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

of amplitude phi and modulus m, using the arithmetic - geometric mean algorithm.

精度

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

概要

# double m1, y, ellpe();

$y = ellpe( $m1 );

説明

Approximates the 積分

            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.

精度

                      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

エラーメッセージ

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

概要

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

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

説明

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. These functions are periodic, with quarter-period on the real axis equal to the complete elliptic 積分 ellpk(1.0-m). Relation to 不完全 elliptic 積分: If u = ellik(phi,m), then sn(u|m) = sin(phi), and cn(u|m) = cos(phi). Phi is called the amplitude of u. 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.

精度

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 不完全 elliptic 積分. Accuracy deteriorates when u is large.

概要

# double m1, y, ellpk();

$y = ellpk( $m1 );

説明

Approximates the 積分

            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.

K(0) = pi/2.

精度

                      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

エラーメッセージ

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

概要

# 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.

概要

# double x, y, exp();

$y = exp( $x );

説明

Returns e (2.71828...) raised to the x power. 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(x**2)/( Q(x**2) - P(x**2) ) of degree 2/3 is used to approximate exp(f) in the basic interval [-0.5, 0.5].

精度

                      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

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.

エラーメッセージ

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

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

概要

# double x, y, exp10();

$y = exp10( $x );

説明

Returns 10 raised to the x power. Range reduction is accomplished by expressing the argument as 10**x = 2**n 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.

精度

                      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

エラーメッセージ

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

DEC arithmetic: MAXL10 = 38.230809449325611792. IEEE arithmetic: MAXL10 = 308.2547155599167.

概要

# double x, y, exp2();

$y = exp2( $x );

説明

Returns 2 raised to the x power. 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].

精度

                      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.

エラーメッセージ

   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.

概要

#double x, y, ei();

$y = ei( $x );

説明

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

Not defined for x <= 0. See also expn.c.

精度

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

概要

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

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

説明

Evaluates the exponential 積分

                 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.

精度

                      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

概要

# double x, y;

$y = fabs( $x );

説明

Returns the absolute value of the argument.

概要

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

$y = fac( $i );

説明

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.

精度

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

概要

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

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

説明

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 分布s with df1 and df2 degrees of freedom, respectively. The 不完全 beta 積分 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.

精度

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.

エラーメッセージ

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

概要

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

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

説明

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 不完全 beta 積分 is used, according to the formula

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

精度

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.

エラーメッセージ

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

概要

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

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

説明

Finds the F density argument x such that the 積分 from x to infinity of the F density is equal to the given probability p. This is accomplished using the inverse beta 積分 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 分布:

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

精度

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.

エラーメッセージ

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

概要

ceil() returns the smallest integer greater than or equal to x. It truncates toward plus infinity. SYNOPSIS:

# double x, y, ceil();

$y = ceil( $x );

概要

floor() returns the largest integer less than or equal to x. It truncates toward minus infinity. SYNOPSIS:

# double x, y, floor();

$y = floor( $x );

概要

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 );

概要

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

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

概要

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

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

説明

Evaluates the Fresnel 積分s

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

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

The 積分s 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 )

精度

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

概要

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

$y = gamma( $x );

説明

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.

精度

                      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

Error for arguments outside the test range will be larger owing to error amplification by the exponential function.

概要

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

$y = lgam( $x );

説明

Returns the base e (2.718...) logarithm of the absolute value of the gamma function of the argument. The sign (+1 or -1) of the gamma function is returned in a global (extern) variable named sgngam. For arguments greater than 13, the logarithm of the gamma function is approximated by the logarithmic version of Stirling's formula using a polynomial approximation of degree 4. Arguments between -33 and +33 are reduced by recurrence to the interval [2,3] of a rational approximation. The cosecant reflection formula is employed for arguments less than -33. Arguments greater than MAXLGM return MAXNUM and an error message. MAXLGM = 2.035093e36 for DEC arithmetic or 2.556348e305 for IEEE arithmetic.

精度

 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 following test used the relative error criterion, though at certain points the relative error could be much higher than indicated.

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

概要

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

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

説明

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

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

The 不完全 gamma 積分 is used, according to the relation

y = igam( b, ax ).

精度

See igam().

エラーメッセージ

   message         condition      value returned
 gdtr domain         x < 0            0.0

概要

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

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

説明

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

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

The 不完全 gamma 積分 is used, according to the relation

y = igamc( b, ax ).

精度

See igamc().

エラーメッセージ

   message         condition      value returned
 gdtrc domain         x < 0            0.0

概要

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

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

概要

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

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

説明

  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

Cases addressed are Tests and escapes for negative integer a, b, or c Linear transformation if c - a or c - b 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 integer Conditionally, a recurrence on c to make c-a-b > 0 |x| > 1 is rejected. The parameters a, b, c are considered to be integer valued if they are within 1.0e-14 of the nearest integer (1.0e-13 for IEEE arithmetic).

精度

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

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

エラーメッセージ

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

概要

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

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

説明

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. 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 error due to roundoff, cancellation, and nonconvergence is estimated. The result with smaller estimated error is returned.

精度

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 are not in the region of near convergence of the asymptotic series. An error message is printed if the self-estimated relative error is greater than 1.0e-12.

概要

# double x, y, i0();

$y = i0( $x );

説明

Returns modified ベッセル関数 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.

精度

                      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

概要

# double x, y, i0e();

$y = i0e( $x );

説明

Returns exponentially scaled modified ベッセル関数 of order zero of the argument. The function is defined as i0e(x) = exp(-|x|) j0( ix ).

精度

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

See i0().

概要

# double x, y, i1();

$y = i1( $x );

説明

Returns modified ベッセル関数 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.

精度

                      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

概要

# double x, y, i1e();

$y = i1e( $x );

説明

Returns exponentially scaled modified ベッセル関数 of order one of the argument. The function is defined as i1(x) = -i exp(-|x|) j1( ix ).

精度

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

See i1().

概要

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

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

説明

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 積分 is evaluated by either a power series or continued fraction expansion, depending on the relative values of a and x.

精度

                      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

概要

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

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

説明

The function is defined by

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

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

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

精度

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

概要

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

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

説明

Given p, the function finds x such that

igamc( a, x ) = p.

Starting with the approximate value

         3
  x = a t

where

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

and

d = 1/9a,

the routine performs up to 10 Newton iterations to find the root of igamc(a,x) - p = 0.

精度

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

概要

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

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

説明

Returns 不完全 beta 積分 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. The 積分 from x to 1 may be obtained by the symmetry relation

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

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

精度

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.

エラーメッセージ

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

概要

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

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

説明

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.

精度

                      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

概要

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

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

説明

Returns modified ベッセル関数 of order v of the argument. If x is negative, v must be integer valued. 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.

精度

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.

概要

# double x, y, j0();

$y = j0( $x );

説明

Returns ベッセル関数 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.

精度

                      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

概要

# double x, y, y0();

$y = y0( $x );

説明

Returns ベッセル関数 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. In the second interval, the Hankel asymptotic expansion is employed with two rational functions of degree 6/6 and 7/7.

精度

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

概要

# double x, y, j1();

$y = j1( $x );

説明

Returns ベッセル関数 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.

精度

                      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

概要

# double x, y, y1();

$y = y1( $x );

説明

Returns ベッセル関数 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.

精度

                      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).

概要

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

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

説明

Returns ベッセル関数 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. If n = 0 or 1 the routine for j0 or j1 is called directly.

精度

                      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.

概要

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

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

説明

Returns ベッセル関数 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.

精度

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.

 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*

概要

# double x, y, k0();

$y = k0( $x );

説明

Returns modified ベッセル関数 of the third kind of order zero of the argument. The range is partitioned into the two intervals [0,8] and (8, infinity). Chebyshev polynomial expansions are employed in each interval.

精度

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

                      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

エラーメッセージ

   message         condition      value returned
  K0 domain          x <= 0          MAXNUM

概要

# double x, y, k0e();

$y = k0e( $x );

説明

Returns exponentially scaled modified ベッセル関数 of the third kind of order zero of the argument.

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

精度

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

See k0().

概要

# double x, y, k1();

$y = k1( $x );

説明

Computes the modified ベッセル関数 of the third kind of order one of the argument. The range is partitioned into the two intervals [0,2] and (2, infinity). Chebyshev polynomial expansions are employed in each interval.

精度

                      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

エラーメッセージ

   message         condition      value returned
 k1 domain          x <= 0          MAXNUM

概要

# double x, y, k1e();

$y = k1e( $x );

説明

Returns exponentially scaled modified ベッセル関数 of the third kind of order one of the argument:

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

精度

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

See k1().

概要

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

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

説明

Returns modified ベッセル関数 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.

精度

                      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

Error is high only near the crossover point x = 9.55 between the two expansions used.

概要

# double x, y, log();

$y = log( $x );

説明

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).

精度

                      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].

エラーメッセージ

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

概要

# double x, y, log10();

$y = log10( $x );

説明

Returns 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).

精度

                      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].

エラーメッセージ

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

概要

# double x, y, log2();

$y = log2( $x );

説明

Returns the base 2 logarithm of x. 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).

精度

                      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.

エラーメッセージ

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

概要

long y, lrand();

$y = lrand( );

説明

Yields a long integer random number. 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.

概要

long x, y;
long lsqrt();

$y = lsqrt( $x );

説明

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.

精度

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.

概要

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

mtherr( $fctnam, $code );

説明

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. SEE ALSO: mconf.h

概要

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

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

説明

Returns the sum of the terms 0 through k of the negative binomial 分布:

   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 nth success. The terms are not computed individually; instead the 不完全 beta 積分 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.

精度

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.

概要

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

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

説明

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

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

The terms are not computed individually; instead the 不完全 beta 積分 is employed, according to the formula

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

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

精度

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.

概要

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

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

説明

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

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

The terms are not computed individually; instead the 不完全 beta 積分 is employed, according to the formula

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

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

精度

See incbet.c.

概要

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

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

説明

Finds the argument p such that nbdtr(k,n,p) is equal to y.

精度

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.

概要

# double x, y, ndtr();

$y = ndtr( $x );

説明

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). Computation is via the functions erf and erfc.

精度

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

エラーメッセージ

   message         condition         value returned
 erfc underflow    x > 37.519379347       0.0

概要

# double x, y, erf();

$y = erf( $x );

説明

The 積分 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. For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise erf(x) = 1 - erfc(x).

精度

                      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

概要

# double x, y, erfc(); $y = erfc( $x );

説明

  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.

精度

                      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

エラーメッセージ

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

概要

# double x, y, ndtri();

$x = ndtri( $y );

説明

Returns the argument, x, for which the area under the Gaussian probability density function (integrated from minus infinity to x) is equal to y. 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)).

精度

                      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

エラーメッセージ

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

概要

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

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

説明

Returns the sum of the first k terms of the Poisson 分布:

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

The terms are not summed directly; instead the 不完全 gamma 積分 is employed, according to the relation

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

The arguments must both be positive.

精度

See igamc().

概要

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

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

説明

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

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

The terms are not summed directly; instead the 不完全 gamma 積分 is employed, according to the formula

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

The arguments must both be positive.

精度

See igam.c.

概要

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

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

説明

Finds the Poisson variable x such that the 積分 from 0 to x of the Poisson density is equal to the given probability y. This is accomplished using the inverse gamma 積分 function and the relation

m = igami( k+1, y ).

精度

See igami.c.

エラーメッセージ

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

概要

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

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

説明

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.

精度

                      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.

エラーメッセージ

   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

概要

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

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

説明

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.

精度

                      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.

概要

# double x, y, psi();

$y = psi( $x );

説明

              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

where the B2k are Bernoulli numbers.

精度

Relative error (except absolute when |psi| < 1):

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

エラーメッセージ

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

概要

# double x, y, rgamma();

$y = rgamma( $x );

説明

Returns one divided by the gamma function of the argument. The function is approximated by a Chebyshev expansion in the interval [0,1]. Range reduction is by recurrence for arguments between -34.034 and +34.84425627277176174. 1/MAXNUM is returned for positive arguments outside this range. For arguments less than -34.034 the cosecant reflection formula is applied; lograrithms are employed to avoid unnecessary overflow. The reciprocal gamma function has no singularities, but overflow and underflow may occur for large arguments. These conditions return either MAXNUM or 1/MAXNUM with appropriate sign.

精度

                      Relative error:
 arithmetic   domain     # trials      peak         rms
    DEC      -30,+30       4000       1.2e-16     1.8e-17
    IEEE     -30,+30      30000       1.1e-15     2.0e-16

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

概要

# double x, y, round();

$y = round( $x );

説明

Returns the nearest integer to x as a double precision floating point result. If x ends in 0.5 exactly, the nearest even integer is chosen.

精度

If x is greater than 1/(2*MACHEP), its closest machine representation is already an integer, so rounding does not change it.

概要

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

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

説明

Approximates the 積分s

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

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

where eul = 0.57721566490153286061 is Euler's constant. The 積分s 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 return MAXNUM.

精度

Test interval 0 to 88.

                      Relative error:
 arithmetic   function  # trials      peak         rms
    DEC          Shi       3000       9.1e-17
    IEEE         Shi      30000       6.9e-16     1.6e-16

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

    DEC          Chi       2500       9.3e-17
    IEEE         Chi      30000       8.4e-16     1.4e-16

概要

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

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

説明

Evaluates the 積分s

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

where eul = 0.57721566490153286061 is Euler's constant. The 積分s are approximated by rational functions. For x > 8 auxiliary functions f(x) and g(x) are employed such that

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

精度

Test interval = [0,50].

Absolute error, except relative when > 1:

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

概要

# double x, y, sin();

$y = sin( $x );

説明

Range reduction is into intervals of pi/4. The reduction error is nearly eliminated by contriving an extended precision modular arithmetic. 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).

精度

                      Relative error:
 arithmetic   domain      # trials      peak         rms
    DEC       0, 10       150000       3.0e-17     7.8e-18
    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17

エラーメッセージ

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

Partial loss of accuracy begins to occur at x = 2**30 = 1.074e9. The loss is not gradual, but 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.

概要

# double x, y, cos();

$y = cos( $x );

説明

Range reduction is into intervals of pi/4. The reduction error is nearly eliminated by contriving an extended precision modular arithmetic. 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).

精度

                      Relative error:
 arithmetic   domain      # trials      peak         rms
    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
    DEC        0,+1.07e9   17000       3.0e-17     7.2e-18

概要

# double x, y, sindg();

$y = sindg( $x );

説明

Range reduction is into intervals of 45 degrees. 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).

精度

                      Relative error:
 arithmetic   domain      # trials      peak         rms
    DEC       +-1000        3100      3.3e-17      9.0e-18
    IEEE      +-1000       30000      2.3e-16      5.6e-17

エラーメッセージ

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

概要

# double x, y, cosdg();

$y = cosdg( $x );

説明

Range reduction is into intervals of 45 degrees. 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).

精度

                      Relative error:
 arithmetic   domain      # trials      peak         rms
    DEC      +-1000         3400       3.5e-17     9.1e-18
    IEEE     +-1000        30000       2.1e-16     5.7e-17

See also sin().

概要

# double x, y, sinh();

$y = sinh( $x );

説明

Returns hyperbolic sine of argument in the range MINLOG to MAXLOG. The range is partitioned into two segments. 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.

精度

                      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

概要

# double x, y, spence();

$y = spence( $x );

説明

Computes the 積分

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

for x >= 0. A rational approximation gives the 積分 in the interval (0.5, 1.5). Transformation formulas for 1/x and 1-x are employed outside the basic expansion range.

精度

                      Relative error:
 arithmetic   domain     # trials      peak         rms
    IEEE      0,4         30000       3.9e-15     5.4e-16
    DEC       0,4          3000       2.5e-16     4.5e-17

概要

# double x, y, sqrt();

$y = sqrt( $x );

説明

Returns the square root of x. Range reduction involves isolating the power of two of the argument and using a polynomial approximation to obtain a rough value for the square root. Then Heron's iteration is used three times to converge to an accurate value.

精度

                      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

エラーメッセージ

   message         condition      value returned
 sqrt domain        x < 0            0.0

概要

# double t, stdtr();
short k;

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

説明

Computes the 積分 from minus infinity to t of the Student t 分布 with integer k > 0 degrees of freedom:

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

Relation to 不完全 beta 積分:

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. 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.

精度

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

概要

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

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

説明

Given probability p, finds the argument t such that stdtr(k,t) is equal to p.

精度

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

概要

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

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

説明

Computes the Struve function Hv(x) of order v, argument x. Negative x is rejected unless v is an integer.

精度

Not accurately characterized, but spot checked against tables.

概要

# double lambda, T, y, plancki()

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

説明

Evaluates the definite 積分, 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 積分 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.

精度

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

概要

# 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 積分:

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

精度

                       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

概要

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

説明

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.

概要

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

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

説明

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 積分, 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.

概要

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

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

説明

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  .

精度

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

概要

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

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

精度

Not accurately characterized, but spot checked against tables.

概要

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

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

精度

Not accurately characterized, but spot checked against tables.

概要

# double v, x;

# double yv( v, x );

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

精度

Not accurately characterized, but spot checked against tables.

概要

# double x, y, tan();

$y = tan( $x );

説明

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].

精度

                      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

エラーメッセージ

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

概要

 # double x, y, cot();

 $y = cot( $x );

説明

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].

精度

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

エラーメッセージ

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

概要

# double x, y, tandg();

$y = tandg( $x );

説明

Returns the circular tangent of the argument x in degrees. 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].

精度

                      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

エラーメッセージ

   message         condition          value returned
 tandg total loss   x > 8.0e14 (DEC)      0.0
                    x > 1.0e14 (IEEE)
 tandg singularity  x = 180 k  +  90     MAXNUM

概要

# double x, y, cotdg();

$y = cotdg( $x );

説明

Returns the circular cotangent of the argument x in degrees. 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].

エラーメッセージ

   message         condition          value returned
 cotdg total loss   x > 8.0e14 (DEC)      0.0
                    x > 1.0e14 (IEEE)
 cotdg singularity  x = 180 k            MAXNUM

概要

# double x, y, tanh();

$y = tanh( $x );

説明

Returns hyperbolic tangent of argument in the range MINLOG to MAXLOG. A rational function is used for |x| < 0.625. The form x + x**3 P(x)/Q(x) of Cody _& Waite is employed. Otherwise,

tanh(x) = sinh(x)/cosh(x) = 1  -  2/(exp(2x) + 1).

精度

                      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

概要

#    log1p(x) = log(1+x)
$y = log1p( $x );

#    expm1(x) = exp(x) - 1
$y = expm1( $x );

#    cosm1(x) = cos(x) - 1
$y = cosm1( $x );

概要

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

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

説明

Returns ベッセル関数 of order n, where n is a (possibly negative) integer. 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.

精度

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

エラーメッセージ

   message         condition      value returned
 yn singularity       x = 0          MAXNUM
 yn overflow                         MAXNUM

Spot checked against tables for x, n between 0 and 100.

概要

# double x, q, y, zeta();

$y = zeta( $x, $q );

説明

                 inf.
                  -        -x
   zeta(x,q)  =   >   (k+q)  
                  -
                 k=0

where x > 1 and q is not a negative integer or zero. 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. Note that (see zetac.c) zeta(x,1) = zetac(x) + 1.

精度

REFERENCE: Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, Series, and Products, p. 1073; Academic Press, 1980.

概要

# double x, y, zetac();

$y = zetac( $x );

説明

                inf.
                 -    -x
   zetac(x)  =   >   k   ,   x > 1,
                 -
                k=2

is related to the Riemann zeta function by

Riemann zeta(x) = zetac(x) + 1.

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.

精度

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