/* 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 */
/* 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 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 integral * ellpk(1.0-m). * * Relation to incomplete elliptic integral: * If u = ellik(phi,m), then sn(u|m) = sin(phi), * and cn(u|m) = cos(phi). Phi is called the amplitude of u. * * 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. * * ACCURACY: * * Tested at random points with u between 0 and 10, m between * 0 and 1. * * Absolute error (* = relative error): * arithmetic function # trials peak rms * IEEE sn 10000 1.7e-6 2.2e-7 * IEEE cn 10000 1.6e-6 2.2e-7 * IEEE dn 100000 3.2e-6 2.6e-7 * IEEE phi 10000 3.9e-7* 6.7e-8* * * Larger errors occur for m near 1. * Peak error observed in consistency check using addition * theorem for sn(u+v) was 4e-16 (absolute). Also tested by * the above relation to the incomplete elliptic integral. * Accuracy deteriorates when u is large. * */
/* ellpkf.c * * Complete elliptic integral of the first kind * * * * SYNOPSIS: * * float m1, y, ellpkf(); * * y = ellpkf( m1 ); * * * * DESCRIPTION: * * Approximates the integral * * * * pi/2 * - * | | * | dt * K(m) = | ------------------ * | 2 * | | sqrt( 1 - m sin t ) * - * 0 * * where m = 1 - m1, using the approximation * * P(x) - log x Q(x). * * The argument m1 is used rather than m so that the logarithmic * singularity at m = 1 will be shifted to the origin; this * preserves maximum accuracy. * * K(0) = pi/2. * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,1 30000 1.3e-7 3.4e-8 * * ERROR MESSAGES: * * message condition value returned * ellpkf domain x<0, x>1 0.0 * */
/* exp10f.c * * Base 10 exponential function * (Common antilogarithm) * * * * SYNOPSIS: * * float x, y, exp10f(); * * y = exp10f( x ); * * * * DESCRIPTION: * * 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). * A polynomial approximates 10**f. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -38,+38 100000 9.8e-8 2.8e-8 * * ERROR MESSAGES: * * message condition value returned * exp10 underflow x < -MAXL10 0.0 * exp10 overflow x > MAXL10 MAXNUM * * IEEE single arithmetic: MAXL10 = 38.230809449325611792. * */
/* exp2f.c * * Base 2 exponential function * * * * SYNOPSIS: * * float x, y, exp2f(); * * y = exp2f( x ); * * * * DESCRIPTION: * * 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 polynomial approximates 2**x in the basic range [-0.5, 0.5]. * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -127,+127 100000 1.7e-7 2.8e-8 * * * See exp.c for comments on error amplification. * * * ERROR MESSAGES: * * message condition value returned * exp underflow x < -MAXL2 0.0 * exp overflow x > MAXL2 MAXNUMF * * For IEEE arithmetic, MAXL2 = 127. */
/* expf.c * * Exponential function * * * * SYNOPSIS: * * float x, y, expf(); * * y = expf( x ); * * * * DESCRIPTION: * * 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 polynomial is used to approximate exp(f) * in the basic range [-0.5, 0.5]. * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE +- MAXLOG 100000 1.7e-7 2.8e-8 * * * 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. * * * ERROR MESSAGES: * * message condition value returned * expf underflow x < MINLOGF 0.0 * expf overflow x > MAXLOGF MAXNUMF * */
/* expnf.c * * Exponential integral En * * * * SYNOPSIS: * * int n; * float x, y, expnf(); * * y = expnf( n, x ); * * * * DESCRIPTION: * * Evaluates the exponential integral * * inf. * - * | | -xt * | e * E (x) = | ---- dt. * n | n * | | t * - * 1 * * * Both n and x must be nonnegative. * * The routine employs either a power series, a continued * fraction, or an asymptotic formula depending on the * relative values of n and x. * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0, 30 10000 5.6e-7 1.2e-7 * */
/* expx2f.c * * Exponential of squared argument * * * * SYNOPSIS: * * double x, y, expx2f(); * * y = expx2f( x ); * * * * DESCRIPTION: * * Computes y = exp(x*x) while suppressing error amplification * that would ordinarily arise from the inexactness of the argument x*x. * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -9.4, 9.4 10^7 1.7e-7 4.7e-8 * */
/* facf.c * * Factorial function * * * * SYNOPSIS: * * float y, facf(); * int i; * * y = facf( i ); * * * * DESCRIPTION: * * Returns factorial of i = 1 * 2 * 3 * ... * i. * fac(0) = 1.0. * * Due to machine arithmetic bounds the largest value of * i accepted is 33 in single precision arithmetic. * Greater values, or negative ones, * produce an error message and return MAXNUM. * * * * ACCURACY: * * For i < 34 the values are simply tabulated, and have * full machine accuracy. * */
/* fdtrf.c * * F distribution * * * * SYNOPSIS: * * int df1, df2; * float x, y, fdtrf(); * * y = fdtrf( df1, df2, x ); * * * * DESCRIPTION: * * Returns the area from zero to x under the F density * function (also known as Snedcor's density or the * variance ratio density). This is the density * of x = (u1/df1)/(u2/df2), where u1 and u2 are random * variables having Chi square distributions with df1 * and df2 degrees of freedom, respectively. * * The incomplete beta integral is used, according to the * formula * * P(x) = incbet( df1/2, df2/2, (df1*x/(df2 + df1*x) ). * * * The arguments a and b are greater than zero, and x * x is nonnegative. * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 2.2e-5 1.1e-6 * * ERROR MESSAGES: * * message condition value returned * fdtrf domain a<0, b<0, x<0 0.0 * */
/* fdtrcf() * * Complemented F distribution * * * * SYNOPSIS: * * int df1, df2; * float x, y, fdtrcf(); * * y = fdtrcf( df1, df2, x ); * * * * DESCRIPTION: * * Returns the area from x to infinity under the F density * function (also known as Snedcor's density or the * variance ratio density). * * * inf. * - * 1 | | a-1 b-1 * 1-P(x) = ------ | t (1-t) dt * B(a,b) | | * - * x * * (See fdtr.c.) * * The incomplete beta integral is used, according to the * formula * * P(x) = incbet( df2/2, df1/2, (df2/(df2 + df1*x) ). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 7.3e-5 1.2e-5 * * ERROR MESSAGES: * * message condition value returned * fdtrcf domain a<0, b<0, x<0 0.0 * */
/* fdtrif() * * Inverse of complemented F distribution * * * * SYNOPSIS: * * float df1, df2, x, y, fdtrif(); * * x = fdtrif( df1, df2, y ); * * * * * DESCRIPTION: * * Finds the F density argument x such that the integral * from x to infinity of the F density is equal to the * given probability y. * * This is accomplished using the inverse beta integral * function and the relations * * z = incbi( df2/2, df1/2, y ) * x = df2 (1-z) / (df1 z). * * Note: the following relations hold for the inverse of * the uncomplemented F distribution: * * z = incbi( df1/2, df2/2, y ) * x = df2 z / (df1 (1-z)). * * * * ACCURACY: * * arithmetic domain # trials peak rms * Absolute error: * IEEE 0,100 5000 4.0e-5 3.2e-6 * Relative error: * IEEE 0,100 5000 1.2e-3 1.8e-5 * * ERROR MESSAGES: * * message condition value returned * fdtrif domain y <= 0 or y > 1 0.0 * v < 1 * */
/* ceilf() * floorf() * frexpf() * ldexpf() * signbitf() * isnanf() * isfinitef() * * Single precision floating point numeric utilities * * * * SYNOPSIS: * * float x, y; * float ceilf(), floorf(), frexpf(), ldexpf(); * int signbit(), isnan(), isfinite(); * int expnt, n; * * y = floorf(x); * y = ceilf(x); * y = frexpf( x, &expnt ); * y = ldexpf( x, n ); * n = signbit(x); * n = isnan(x); * n = isfinite(x); * * * * DESCRIPTION: * * All four routines return a single precision floating point * result. * * sfloor() returns the largest integer less than or equal to x. * It truncates toward minus infinity. * * sceil() returns the smallest integer greater than or equal * to x. It truncates toward plus infinity. * * sfrexp() 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. * * ldexpf() multiplies x by 2**n. * * signbit(x) returns 1 if the sign bit of x is 1, else 0. * * These functions are part of the standard C run time library * for many but not all C compilers. The ones supplied are * written in C for either DEC or IEEE arithmetic. They should * be used only if your compiler library does not already have * them. * * The IEEE versions assume that denormal numbers are implemented * in the arithmetic. Some modifications will be required if * the arithmetic has abrupt rather than gradual underflow. */
/* fresnlf.c * * Fresnel integral * * * * SYNOPSIS: * * float x, S, C; * void fresnlf(); * * fresnlf( x, &S, &C ); * * * DESCRIPTION: * * Evaluates the Fresnel integrals * * x * - * | | * C(x) = | cos(pi/2 t**2) dt, * | | * - * 0 * * x * - * | | * S(x) = | sin(pi/2 t**2) dt. * | | * - * 0 * * * The integrals are evaluated by power series for small x. * For x >= 1 auxiliary functions f(x) and g(x) are employed * such that * * C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 ) * S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 ) * * * * ACCURACY: * * Relative error. * * Arithmetic function domain # trials peak rms * IEEE S(x) 0, 10 30000 1.1e-6 1.9e-7 * IEEE C(x) 0, 10 30000 1.1e-6 2.0e-7 */
/* gammaf.c * * Gamma function * * * * SYNOPSIS: * * float x, y, gammaf(); * extern int sgngamf; * * y = gammaf( x ); * * * * DESCRIPTION: * * Returns gamma function of the argument. The result is * correctly signed, and the sign (+1 or -1) is also * returned in a global (extern) variable named sgngamf. * This same variable is also filled in by the logarithmic * gamma function lgam(). * * Arguments between 0 and 10 are reduced by recurrence and the * function is approximated by a polynomial function covering * the interval (2,3). Large arguments are handled by Stirling's * formula. Negative arguments are made positive using * a reflection formula. * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,-33 100,000 5.7e-7 1.0e-7 * IEEE -33,0 100,000 6.1e-7 1.2e-7 * * */
/* lgamf() * * Natural logarithm of gamma function * * * * SYNOPSIS: * * float x, y, lgamf(); * extern int sgngamf; * * y = lgamf( x ); * * * * DESCRIPTION: * * 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 sgngamf. * * For arguments greater than 6.5, the logarithm of the gamma * function is approximated by the logarithmic version of * Stirling's formula. Arguments between 0 and +6.5 are reduced by * by recurrence to the interval [.75,1.25] or [1.5,2.5] of a rational * approximation. The cosecant reflection formula is employed for * arguments less than zero. * * Arguments greater than MAXLGM = 2.035093e36 return MAXNUM and an * error message. * * * * ACCURACY: * * * * arithmetic domain # trials peak rms * IEEE -100,+100 500,000 7.4e-7 6.8e-8 * The error criterion was relative when the function magnitude * was greater than one but absolute when it was less than one. * The routine has low relative error for positive arguments. * * The following test used the relative error criterion. * IEEE -2, +3 100000 4.0e-7 5.6e-8 * */
/* gdtrf.c * * Gamma distribution function * * * * SYNOPSIS: * * float a, b, x, y, gdtrf(); * * y = gdtrf( a, b, x ); * * * * DESCRIPTION: * * Returns the integral from zero to x of the gamma probability * density function: * * * x * b - * a | | b-1 -at * y = ----- | t e dt * - | | * | (b) - * 0 * * The incomplete gamma integral is used, according to the * relation * * y = igam( b, ax ). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 5.8e-5 3.0e-6 * * ERROR MESSAGES: * * message condition value returned * gdtrf domain x < 0 0.0 * */
/* gdtrcf.c * * Complemented gamma distribution function * * * * SYNOPSIS: * * float a, b, x, y, gdtrcf(); * * y = gdtrcf( a, b, x ); * * * * DESCRIPTION: * * Returns the integral from x to infinity of the gamma * probability density function: * * * inf. * b - * a | | b-1 -at * y = ----- | t e dt * - | | * | (b) - * x * * The incomplete gamma integral is used, according to the * relation * * y = igamc( b, ax ). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 9.1e-5 1.5e-5 * * ERROR MESSAGES: * * message condition value returned * gdtrcf domain x < 0 0.0 * */
/* hyp2f1f.c * * Gauss hypergeometric function F * 2 1 * * * SYNOPSIS: * * float a, b, c, x, y, hyp2f1f(); * * y = hyp2f1f( a, b, c, x ); * * * DESCRIPTION: * * * hyp2f1( a, b, c, x ) = F ( a, b; c; x ) * 2 1 * * inf. * - a(a+1)...(a+k) b(b+1)...(b+k) k+1 * = 1 + > ----------------------------- x . * - c(c+1)...(c+k) (k+1)! * k = 0 * * 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-6 of the nearest integer. * * ACCURACY: * * Relative error (-1 < x < 1): * arithmetic domain # trials peak rms * IEEE 0,3 30000 5.8e-4 4.3e-6 */
/* hypergf.c * * Confluent hypergeometric function * * * * SYNOPSIS: * * float a, b, x, y, hypergf(); * * y = hypergf( a, b, x ); * * * * DESCRIPTION: * * Computes the confluent hypergeometric function * * 1 2 * a x a(a+1) x * F ( a,b;x ) = 1 + ---- + --------- + ... * 1 1 b 1! b(b+1) 2! * * Many higher transcendental functions are special cases of * this power series. * * 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. * * * * ACCURACY: * * Tested at random points (a, b, x), all three variables * ranging from 0 to 30. * Relative error: * arithmetic domain # trials peak rms * IEEE 0,5 10000 6.6e-7 1.3e-7 * IEEE 0,30 30000 1.1e-5 6.5e-7 * * 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-3. * */
/* i0f.c * * Modified Bessel function of order zero * * * * SYNOPSIS: * * float x, y, i0(); * * y = i0f( x ); * * * * DESCRIPTION: * * Returns modified Bessel function of order zero of the * argument. * * The function is defined as i0(x) = j0( ix ). * * The range is partitioned into the two intervals [0,8] and * (8, infinity). Chebyshev polynomial expansions are employed * in each interval. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,30 100000 4.0e-7 7.9e-8 * */
/* i0ef.c * * Modified Bessel function of order zero, * exponentially scaled * * * * SYNOPSIS: * * float x, y, i0ef(); * * y = i0ef( x ); * * * * DESCRIPTION: * * Returns exponentially scaled modified Bessel function * of order zero of the argument. * * The function is defined as i0e(x) = exp(-|x|) j0( ix ). * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,30 100000 3.7e-7 7.0e-8 * See i0f(). * */
/* i1f.c * * Modified Bessel function of order one * * * * SYNOPSIS: * * float x, y, i1f(); * * y = i1f( x ); * * * * DESCRIPTION: * * Returns modified Bessel function of order one of the * argument. * * The function is defined as i1(x) = -i j1( ix ). * * The range is partitioned into the two intervals [0,8] and * (8, infinity). Chebyshev polynomial expansions are employed * in each interval. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0, 30 100000 1.5e-6 1.6e-7 * * */
/* i1ef.c * * Modified Bessel function of order one, * exponentially scaled * * * * SYNOPSIS: * * float x, y, i1ef(); * * y = i1ef( x ); * * * * DESCRIPTION: * * Returns exponentially scaled modified Bessel function * of order one of the argument. * * The function is defined as i1(x) = -i exp(-|x|) j1( ix ). * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0, 30 30000 1.5e-6 1.5e-7 * See i1(). * */
/* igamf.c * * Incomplete gamma integral * * * * SYNOPSIS: * * float a, x, y, igamf(); * * y = igamf( a, x ); * * * * DESCRIPTION: * * The function is defined by * * x * - * 1 | | -t a-1 * igam(a,x) = ----- | e t dt. * - | | * | (a) - * 0 * * * In this implementation both arguments must be positive. * The integral is evaluated by either a power series or * continued fraction expansion, depending on the relative * values of a and x. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,30 20000 7.8e-6 5.9e-7 * */
/* igamcf() * * Complemented incomplete gamma integral * * * * SYNOPSIS: * * float a, x, y, igamcf(); * * y = igamcf( a, x ); * * * * DESCRIPTION: * * 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 integral is evaluated by either a power series or * continued fraction expansion, depending on the relative * values of a and x. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,30 30000 7.8e-6 5.9e-7 * */
/* igamif() * * Inverse of complemented imcomplete gamma integral * * * * SYNOPSIS: * * float a, x, y, igamif(); * * x = igamif( a, y ); * * * * DESCRIPTION: * * Given y, the function finds x such that * * igamc( a, x ) = y. * * It is valid in the right-hand tail of the distribution, y < 0.5. * Starting with the approximate value * * 3 * x = a t * * where * * t = 1 - d - ndtri(y) sqrt(d) * * and * * d = 1/9a, * * the routine performs up to 10 Newton iterations to find the * root of igamc(a,x) - y = 0. * * * ACCURACY: * * Tested for a ranging from 0 to 100 and x from 0 to 1. * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 1.0e-5 1.5e-6 * */
/* incbetf.c * * Incomplete beta integral * * * SYNOPSIS: * * float a, b, x, y, incbetf(); * * y = incbetf( a, b, x ); * * * DESCRIPTION: * * Returns incomplete beta integral of the arguments, evaluated * from zero to x. The function is defined as * * x * - - * | (a+b) | | a-1 b-1 * ----------- | t (1-t) dt. * - - | | * | (a) | (b) - * 0 * * The domain of definition is 0 <= x <= 1. In this * implementation a and b are restricted to positive values. * The integral from x to 1 may be obtained by the symmetry * relation * * 1 - incbet( a, b, x ) = incbet( b, a, 1-x ). * * The integral is evaluated by a continued fraction expansion. * If a < 1, the function calls itself recursively after a * transformation to increase a to a+1. * * ACCURACY: * * Tested at random points (a,b,x) with a and b in the indicated * interval and x between 0 and 1. * * arithmetic domain # trials peak rms * Relative error: * IEEE 0,30 10000 3.7e-5 5.1e-6 * IEEE 0,100 10000 1.7e-4 2.5e-5 * The useful domain for relative error is limited by underflow * of the single precision exponential function. * Absolute error: * IEEE 0,30 100000 2.2e-5 9.6e-7 * IEEE 0,100 10000 6.5e-5 3.7e-6 * * Larger errors may occur for extreme ratios of a and b. * * ERROR MESSAGES: * message condition value returned * incbetf domain x<0, x>1 0.0 */
/* incbif() * * Inverse of imcomplete beta integral * * * * SYNOPSIS: * * float a, b, x, y, incbif(); * * x = incbif( a, b, y ); * * * * DESCRIPTION: * * Given y, the function finds x such that * * incbet( a, b, x ) = y. * * the routine performs up to 10 Newton iterations to find the * root of incbet(a,b,x) - y = 0. * * * ACCURACY: * * Relative error: * x a,b * arithmetic domain domain # trials peak rms * IEEE 0,1 0,100 5000 2.8e-4 8.3e-6 * * Overflow and larger errors may occur for one of a or b near zero * and the other large. */
/* ivf.c * * Modified Bessel function of noninteger order * * * * SYNOPSIS: * * float v, x, y, ivf(); * * y = ivf( v, x ); * * * * DESCRIPTION: * * Returns modified Bessel function 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. * * * ACCURACY: * * Tested at random points (v, x), with v between 0 and * 30, x between 0 and 28. * arithmetic domain # trials peak rms * Relative error: * IEEE 0,15 3000 4.7e-6 5.4e-7 * Absolute error (relative when function > 1) * IEEE 0,30 5000 8.5e-6 1.3e-6 * * Accuracy is diminished if v is near a negative integer. * The useful domain for relative error is limited by overflow * of the single precision exponential function. * * See also hyperg.c. * */
/* j0f.c * * Bessel function of order zero * * * * SYNOPSIS: * * float x, y, j0f(); * * y = j0f( x ); * * * * DESCRIPTION: * * Returns Bessel function of order zero of the argument. * * The domain is divided into the intervals [0, 2] and * (2, infinity). In the first interval the following polynomial * approximation is used: * * * 2 2 2 * (w - r ) (w - r ) (w - r ) P(w) * 1 2 3 * * 2 * where w = x and the three r's are zeros of the function. * * In the second interval, the modulus and phase are approximated * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x) * and Phase(x) = x + 1/x R(1/x^2) - pi/4. The function is * * j0(x) = Modulus(x) cos( Phase(x) ). * * * * ACCURACY: * * Absolute error: * arithmetic domain # trials peak rms * IEEE 0, 2 100000 1.3e-7 3.6e-8 * IEEE 2, 32 100000 1.9e-7 5.4e-8 * */
/* y0f.c * * Bessel function of the second kind, order zero * * * * SYNOPSIS: * * float x, y, y0f(); * * y = y0f( x ); * * * * DESCRIPTION: * * Returns Bessel function of the second kind, of order * zero, of the argument. * * The domain is divided into the intervals [0, 2] and * (2, infinity). In the first interval a rational approximation * R(x) is employed to compute * * 2 2 2 * y0(x) = (w - r ) (w - r ) (w - r ) R(x) + 2/pi ln(x) j0(x). * 1 2 3 * * Thus a call to j0() is required. The three zeros are removed * from R(x) to improve its numerical stability. * * In the second interval, the modulus and phase are approximated * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x) * and Phase(x) = x + 1/x S(1/x^2) - pi/4. Then the function is * * y0(x) = Modulus(x) sin( Phase(x) ). * * * * * ACCURACY: * * Absolute error, when y0(x) < 1; else relative error: * * arithmetic domain # trials peak rms * IEEE 0, 2 100000 2.4e-7 3.4e-8 * IEEE 2, 32 100000 1.8e-7 5.3e-8 * */
/* j1f.c * * Bessel function of order one * * * * SYNOPSIS: * * float x, y, j1f(); * * y = j1f( x ); * * * * DESCRIPTION: * * Returns Bessel function of order one of the argument. * * The domain is divided into the intervals [0, 2] and * (2, infinity). In the first interval a polynomial approximation * 2 * (w - r ) x P(w) * 1 * 2 * is used, where w = x and r is the first zero of the function. * * In the second interval, the modulus and phase are approximated * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x) * and Phase(x) = x + 1/x R(1/x^2) - 3pi/4. The function is * * j0(x) = Modulus(x) cos( Phase(x) ). * * * * ACCURACY: * * Absolute error: * arithmetic domain # trials peak rms * IEEE 0, 2 100000 1.2e-7 2.5e-8 * IEEE 2, 32 100000 2.0e-7 5.3e-8 * * */
/* y1 * * Bessel function of second kind of order one * * * * SYNOPSIS: * * double x, y, y1(); * * y = y1( x ); * * * * DESCRIPTION: * * Returns Bessel function of the second kind of order one * of the argument. * * The domain is divided into the intervals [0, 2] and * (2, infinity). In the first interval a rational approximation * R(x) is employed to compute * * 2 * y0(x) = (w - r ) x R(x^2) + 2/pi (ln(x) j1(x) - 1/x) . * 1 * * Thus a call to j1() is required. * * In the second interval, the modulus and phase are approximated * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x) * and Phase(x) = x + 1/x S(1/x^2) - 3pi/4. Then the function is * * y0(x) = Modulus(x) sin( Phase(x) ). * * * * * ACCURACY: * * Absolute error: * arithmetic domain # trials peak rms * IEEE 0, 2 100000 2.2e-7 4.6e-8 * IEEE 2, 32 100000 1.9e-7 5.3e-8 * * (error criterion relative when |y1| > 1). * */
/* jnf.c * * Bessel function of integer order * * * * SYNOPSIS: * * int n; * float x, y, jnf(); * * y = jnf( n, x ); * * * * DESCRIPTION: * * Returns Bessel function of order n, where n is a * (possibly negative) integer. * * The ratio of jn(x) to j0(x) is computed by backward * recurrence. First the ratio jn/jn-1 is found by a * continued fraction expansion. Then the recurrence * relating successive orders is applied until j0 or j1 is * reached. * * If n = 0 or 1 the routine for j0 or j1 is called * directly. * * * * ACCURACY: * * Absolute error: * arithmetic range # trials peak rms * IEEE 0, 15 30000 3.6e-7 3.6e-8 * * * Not suitable for large n or x. Use jvf() instead. * */
/* jvf.c * * Bessel function of noninteger order * * * * SYNOPSIS: * * float v, x, y, jvf(); * * y = jvf( v, x ); * * * * DESCRIPTION: * * Returns Bessel function of order v of the argument, * where v is real. Negative x is allowed if v is an integer. * * Several expansions are included: the ascending power * series, the Hankel expansion, and two transitional * expansions for large v. If v is not too large, it * is reduced by recurrence to a region of best accuracy. * * The single precision routine accepts negative v, but with * reduced accuracy. * * * * ACCURACY: * Results for integer v are indicated by *. * Error criterion is absolute, except relative when |jv()| > 1. * * arithmetic domain # trials peak rms * v x * IEEE 0,125 0,125 30000 2.0e-6 2.0e-7 * IEEE -17,0 0,125 30000 1.1e-5 4.0e-7 * IEEE -100,0 0,125 3000 1.5e-4 7.8e-6 */
/* k0f.c * * Modified Bessel function, third kind, order zero * * * * SYNOPSIS: * * float x, y, k0f(); * * y = k0f( x ); * * * * DESCRIPTION: * * Returns modified Bessel function 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. * * * * ACCURACY: * * 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 * IEEE 0, 30 30000 7.8e-7 8.5e-8 * * ERROR MESSAGES: * * message condition value returned * K0 domain x <= 0 MAXNUM * */
/* k0ef() * * Modified Bessel function, third kind, order zero, * exponentially scaled * * * * SYNOPSIS: * * float x, y, k0ef(); * * y = k0ef( x ); * * * * DESCRIPTION: * * Returns exponentially scaled modified Bessel function * of the third kind of order zero of the argument. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0, 30 30000 8.1e-7 7.8e-8 * See k0(). * */
/* k1f.c * * Modified Bessel function, third kind, order one * * * * SYNOPSIS: * * float x, y, k1f(); * * y = k1f( x ); * * * * DESCRIPTION: * * Computes the modified Bessel function of the third kind * of order one of the argument. * * The range is partitioned into the two intervals [0,2] and * (2, infinity). Chebyshev polynomial expansions are employed * in each interval. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0, 30 30000 4.6e-7 7.6e-8 * * ERROR MESSAGES: * * message condition value returned * k1 domain x <= 0 MAXNUM * */
/* k1ef.c * * Modified Bessel function, third kind, order one, * exponentially scaled * * * * SYNOPSIS: * * float x, y, k1ef(); * * y = k1ef( x ); * * * * DESCRIPTION: * * Returns exponentially scaled modified Bessel function * of the third kind of order one of the argument: * * k1e(x) = exp(x) * k1(x). * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0, 30 30000 4.9e-7 6.7e-8 * See k1(). * */
/* knf.c * * Modified Bessel function, third kind, integer order * * * * SYNOPSIS: * * float x, y, knf(); * int n; * * y = knf( n, x ); * * * * DESCRIPTION: * * Returns modified Bessel function of the third kind * of order n of the argument. * * The range is partitioned into the two intervals [0,9.55] and * (9.55, infinity). An ascending power series is used in the * low range, and an asymptotic expansion in the high range. * * * * ACCURACY: * * Absolute error, relative when function > 1: * arithmetic domain # trials peak rms * IEEE 0,30 10000 2.0e-4 3.8e-6 * * Error is high only near the crossover point x = 9.55 * between the two expansions used. */
/* log10f.c * * Common logarithm * * * * SYNOPSIS: * * float x, y, log10f(); * * y = log10f( x ); * * * * DESCRIPTION: * * 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). * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0.5, 2.0 100000 1.3e-7 3.4e-8 * IEEE 0, MAXNUMF 100000 1.3e-7 2.6e-8 * * In the tests over the interval [0, MAXNUM], the logarithms * of the random arguments were uniformly distributed over * [-MAXL10, MAXL10]. * * ERROR MESSAGES: * * log10f singularity: x = 0; returns -MAXL10 * log10f domain: x < 0; returns -MAXL10 * MAXL10 = 38.230809449325611792 */
/* log2f.c * * Base 2 logarithm * * * * SYNOPSIS: * * float x, y, log2f(); * * y = log2f( x ); * * * * DESCRIPTION: * * 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). * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE exp(+-88) 100000 1.1e-7 2.4e-8 * IEEE 0.5, 2.0 100000 1.1e-7 3.0e-8 * * In the tests over the interval [exp(+-88)], the logarithms * of the random arguments were uniformly distributed. * * ERROR MESSAGES: * * log singularity: x = 0; returns MINLOGF/log(2) * log domain: x < 0; returns MINLOGF/log(2) */
/* logf.c * * Natural logarithm * * * * SYNOPSIS: * * float x, y, logf(); * * y = logf( x ); * * * * DESCRIPTION: * * Returns the base e (2.718...) logarithm of x. * * The argument is separated into its exponent and fractional * parts. If the exponent is between -1 and +1, the logarithm * of the fraction is approximated by * * log(1+x) = x - 0.5 x**2 + x**3 P(x) * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0.5, 2.0 100000 7.6e-8 2.7e-8 * IEEE 1, MAXNUMF 100000 2.6e-8 * * In the tests over the interval [1, MAXNUM], the logarithms * of the random arguments were uniformly distributed over * [0, MAXLOGF]. * * ERROR MESSAGES: * * logf singularity: x = 0; returns MINLOG * logf domain: x < 0; returns MINLOG */
/* mtherrf.c * * Library common error handling routine * * * * SYNOPSIS: * * char *fctnam; * int code; * void mtherrf(); * * mtherrf( fctnam, code ); * * * * DESCRIPTION: * * This routine may be called to report one of the following * error conditions (in the include file mconf.h). * * Mnemonic Value Significance * * DOMAIN 1 argument domain error * SING 2 function singularity * OVERFLOW 3 overflow range error * UNDERFLOW 4 underflow range error * TLOSS 5 total loss of precision * PLOSS 6 partial loss of precision * EDOM 33 Unix domain error code * ERANGE 34 Unix range error code * * The default version of the file prints the function name, * passed to it by the pointer fctnam, followed by the * error condition. The display is directed to the standard * output device. The routine then returns to the calling * program. Users may wish to modify the program to abort by * calling exit() under severe error conditions such as domain * errors. * * Since all error conditions pass control to this function, * the display may be easily changed, eliminated, or directed * to an error logging device. * * SEE ALSO: * * mconf.h * */
/* nbdtrf.c * * Negative binomial distribution * * * * SYNOPSIS: * * int k, n; * float p, y, nbdtrf(); * * y = nbdtrf( k, n, p ); * * * * DESCRIPTION: * * Returns the sum of the terms 0 through k of the negative * binomial distribution: * * k * -- ( n+j-1 ) n j * > ( ) p (1-p) * -- ( j ) * j=0 * * In a sequence of Bernoulli trials, this is the probability * that k or fewer failures precede the nth success. * * The terms are not computed individually; instead the incomplete * beta integral is employed, according to the formula * * y = nbdtr( k, n, p ) = incbet( n, k+1, p ). * * The arguments must be positive, with p ranging from 0 to 1. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 1.5e-4 1.9e-5 * */
/* nbdtrcf.c * * Complemented negative binomial distribution * * * * SYNOPSIS: * * int k, n; * float p, y, nbdtrcf(); * * y = nbdtrcf( k, n, p ); * * * * DESCRIPTION: * * Returns the sum of the terms k+1 to infinity of the negative * binomial distribution: * * inf * -- ( n+j-1 ) n j * > ( ) p (1-p) * -- ( j ) * j=k+1 * * The terms are not computed individually; instead the incomplete * beta integral 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. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 1.4e-4 2.0e-5 * */
/* ndtrf.c * * Normal distribution function * * * * SYNOPSIS: * * float x, y, ndtrf(); * * y = ndtrf( x ); * * * * DESCRIPTION: * * Returns the area under the Gaussian probability density * function, integrated from minus infinity to x: * * x * - * 1 | | 2 * ndtr(x) = --------- | exp( - t /2 ) dt * sqrt(2pi) | | * - * -inf. * * = ( 1 + erf(z) ) / 2 * = erfc(z) / 2 * * where z = x/sqrt(2). Computation is via the functions * erf and erfc. * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -13,0 50000 1.5e-5 2.6e-6 * * * ERROR MESSAGES: * * See erfcf(). * */
/* erff.c * * Error function * * * * SYNOPSIS: * * float x, y, erff(); * * y = erff( x ); * * * * DESCRIPTION: * * The integral is * * x * - * 2 | | 2 * erf(x) = -------- | exp( - t ) dt. * sqrt(pi) | | * - * 0 * * The magnitude of x is limited to 9.231948545 for DEC * arithmetic; 1 or -1 is returned outside this range. * * For 0 <= |x| < 1, erf(x) = x * P(x**2); otherwise * erf(x) = 1 - erfc(x). * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -9.3,9.3 50000 1.7e-7 2.8e-8 * */
/* erfcf.c * * Complementary error function * * * * SYNOPSIS: * * float x, y, erfcf(); * * y = erfcf( x ); * * * * DESCRIPTION: * * * 1 - erf(x) = * * inf. * - * 2 | | 2 * erfc(x) = -------- | exp( - t ) dt * sqrt(pi) | | * - * x * * * For small x, erfc(x) = 1 - erf(x); otherwise polynomial * approximations 1/x P(1/x**2) are computed. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -9.3,9.3 50000 3.9e-6 7.2e-7 * * * ERROR MESSAGES: * * message condition value returned * erfcf underflow x**2 > MAXLOGF 0.0 * * */
/* ndtrif.c * * Inverse of Normal distribution function * * * * SYNOPSIS: * * float x, y, ndtrif(); * * x = ndtrif( y ); * * * * DESCRIPTION: * * Returns the argument, x, for which the area under the * Gaussian probability density function (integrated from * minus infinity to x) is equal to y. * * * For small arguments 0 < y < exp(-2), the program computes * z = sqrt( -2.0 * log(y) ); then the approximation is * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). * There are two rational functions P/Q, one for 0 < y < exp(-32) * and the other for y up to exp(-2). For larger arguments, * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 1e-38, 1 30000 3.6e-7 5.0e-8 * * * ERROR MESSAGES: * * message condition value returned * ndtrif domain x <= 0 -MAXNUM * ndtrif domain x >= 1 MAXNUM * */
/* pdtrf.c * * Poisson distribution * * * * SYNOPSIS: * * int k; * float m, y, pdtrf(); * * y = pdtrf( k, m ); * * * * DESCRIPTION: * * Returns the sum of the first k terms of the Poisson * distribution: * * k j * -- -m m * > e -- * -- j! * j=0 * * The terms are not summed directly; instead the incomplete * gamma integral is employed, according to the relation * * y = pdtr( k, m ) = igamc( k+1, m ). * * The arguments must both be positive. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 6.9e-5 8.0e-6 * */
/* pdtrcf() * * Complemented poisson distribution * * * * SYNOPSIS: * * int k; * float m, y, pdtrcf(); * * y = pdtrcf( k, m ); * * * * DESCRIPTION: * * Returns the sum of the terms k+1 to infinity of the Poisson * distribution: * * inf. j * -- -m m * > e -- * -- j! * j=k+1 * * The terms are not summed directly; instead the incomplete * gamma integral is employed, according to the formula * * y = pdtrc( k, m ) = igam( k+1, m ). * * The arguments must both be positive. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 8.4e-5 1.2e-5 * */
/* pdtrif() * * Inverse Poisson distribution * * * * SYNOPSIS: * * int k; * float m, y, pdtrf(); * * m = pdtrif( k, y ); * * * * * DESCRIPTION: * * Finds the Poisson variable x such that the integral * from 0 to x of the Poisson density is equal to the * given probability y. * * This is accomplished using the inverse gamma integral * function and the relation * * m = igami( k+1, y ). * * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,100 5000 8.7e-6 1.4e-6 * * ERROR MESSAGES: * * message condition value returned * pdtri domain y < 0 or y >= 1 0.0 * k < 0 * */
/* polevlf.c * p1evlf.c * * Evaluate polynomial * * * * SYNOPSIS: * * int N; * float x, y, coef[N+1], polevlf[]; * * y = polevlf( x, coef, N ); * * * * DESCRIPTION: * * Evaluates polynomial of degree N: * * 2 N * y = C + C x + C x +...+ C x * 0 1 2 N * * Coefficients are stored in reverse order: * * coef[0] = C , ..., coef[N] = C . * N 0 * * The function p1evl() assumes that coef[N] = 1.0 and is * omitted from the array. Its calling arguments are * otherwise the same as polevl(). * * * SPEED: * * In the interest of speed, there are no checks for out * of bounds arithmetic. This routine is used by most of * the functions in the library. Depending on available * equipment features, the user may wish to rewrite the * program in microcode or assembly language. * */
/* polynf.c * polyrf.c * Arithmetic operations on polynomials * * In the following descriptions a, b, c are polynomials of degree * na, nb, nc respectively. The degree of a polynomial cannot * exceed a run-time value MAXPOLF. An operation that attempts * to use or generate a polynomial of higher degree may produce a * result that suffers truncation at degree MAXPOL. The value of * MAXPOL is set by calling the function * * polinif( maxpol ); * * where maxpol is the desired maximum degree. This must be * done prior to calling any of the other functions in this module. * Memory for internal temporary polynomial storage is allocated * by polinif(). * * Each polynomial is represented by an array containing its * coefficients, together with a separately declared integer equal * to the degree of the polynomial. The coefficients appear in * ascending order; that is, * * 2 na * a(x) = a[0] + a[1] * x + a[2] * x + ... + a[na] * x . * * * * sum = poleva( a, na, x ); Evaluate polynomial a(t) at t = x. * polprtf( a, na, D ); Print the coefficients of a to D digits. * polclrf( a, na ); Set a identically equal to zero, up to a[na]. * polmovf( a, na, b ); Set b = a. * poladdf( a, na, b, nb, c ); c = b + a, nc = max(na,nb) * polsubf( a, na, b, nb, c ); c = b - a, nc = max(na,nb) * polmulf( a, na, b, nb, c ); c = b * a, nc = na+nb * * * Division: * * i = poldivf( a, na, b, nb, c ); c = b / a, nc = MAXPOL * * returns i = the degree of the first nonzero coefficient of a. * The computed quotient c must be divided by x^i. An error message * is printed if a is identically zero. * * * Change of variables: * If a and b are polynomials, and t = a(x), then * c(t) = b(a(x)) * is a polynomial found by substituting a(x) for t. The * subroutine call for this is * * polsbtf( a, na, b, nb, c ); * * * Notes: * poldivf() is an integer routine; polevaf() is float. * Any of the arguments a, b, c may refer to the same array. * */
/* powf.c * * Power function * * * * SYNOPSIS: * * float x, y, z, powf(); * * z = powf( x, y ); * * * * DESCRIPTION: * * Computes x raised to the yth power. Analytically, * * x**y = exp( y log(x) ). * * Following Cody and Waite, this program uses a lookup table * of 2**-i/16 and pseudo extended precision arithmetic to * obtain an extra three bits of accuracy in both the logarithm * and the exponential. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -10,10 100,000 1.4e-7 3.6e-8 * 1/10 < x < 10, x uniformly distributed. * -10 < y < 10, y uniformly distributed. * * * ERROR MESSAGES: * * message condition value returned * powf overflow x**y > MAXNUMF MAXNUMF * powf underflow x**y < 1/MAXNUMF 0.0 * powf domain x<0 and y noninteger 0.0 * */
/* powif.c * * Real raised to integer power * * * * SYNOPSIS: * * float x, y, powif(); * int n; * * y = powif( x, n ); * * * * DESCRIPTION: * * Returns argument x raised to the nth power. * The routine efficiently decomposes n as a sum of powers of * two. The desired power is a product of two-to-the-kth * powers of x. Thus to compute the 32767 power of x requires * 28 multiplications instead of 32767 multiplications. * * * * ACCURACY: * * * Relative error: * arithmetic x domain n domain # trials peak rms * IEEE .04,26 -26,26 100000 1.1e-6 2.0e-7 * IEEE 1,2 -128,128 100000 1.1e-5 1.0e-6 * * Returns MAXNUMF on overflow, zero on underflow. * */
/* psif.c * * Psi (digamma) function * * * SYNOPSIS: * * float x, y, psif(); * * y = psif( x ); * * * DESCRIPTION: * * d - * psi(x) = -- ln | (x) * dx * * is the logarithmic derivative of the gamma function. * For integer x, * n-1 * - * psi(n) = -EUL + > 1/k. * - * k=1 * * This formula is used for 0 < n <= 10. If x is negative, it * is transformed to a positive argument by the reflection * formula psi(1-x) = psi(x) + pi cot(pi x). * For general positive x, the argument is made greater than 10 * using the recurrence psi(x+1) = psi(x) + 1/x. * Then the following asymptotic expansion is applied: * * inf. B * - 2k * psi(x) = log(x) - 1/2x - > ------- * - 2k * k=1 2k x * * where the B2k are Bernoulli numbers. * * ACCURACY: * Absolute error, relative when |psi| > 1 : * arithmetic domain # trials peak rms * IEEE -33,0 30000 8.2e-7 1.2e-7 * IEEE 0,33 100000 7.3e-7 7.7e-8 * * ERROR MESSAGES: * message condition value returned * psi singularity x integer <=0 MAXNUMF */
/* rgammaf.c * * Reciprocal gamma function * * * * SYNOPSIS: * * float x, y, rgammaf(); * * y = rgammaf( x ); * * * * DESCRIPTION: * * 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/MAXNUMF is returned for positive arguments outside this * range. * * The reciprocal gamma function has no singularities, * but overflow and underflow may occur for large arguments. * These conditions return either MAXNUMF or 1/MAXNUMF with * appropriate sign. * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -34,+34 100000 8.9e-7 1.1e-7 */
/* shichif.c * * Hyperbolic sine and cosine integrals * * * * SYNOPSIS: * * float x, Chi, Shi; * * shichi( x, &Chi, &Shi ); * * * DESCRIPTION: * * Approximates the integrals * * x * - * | | cosh t - 1 * Chi(x) = eul + ln x + | ----------- dt, * | | t * - * 0 * * x * - * | | sinh t * Shi(x) = | ------ dt * | | t * - * 0 * * where eul = 0.57721566490153286061 is Euler's constant. * The integrals are evaluated by power series for x < 8 * and by Chebyshev expansions for x between 8 and 88. * For large x, both functions approach exp(x)/2x. * Arguments greater than 88 in magnitude return MAXNUM. * * * ACCURACY: * * Test interval 0 to 88. * Relative error: * arithmetic function # trials peak rms * IEEE Shi 20000 3.5e-7 7.0e-8 * Absolute error, except relative when |Chi| > 1: * IEEE Chi 20000 3.8e-7 7.6e-8 */
/* sicif.c * * Sine and cosine integrals * * * * SYNOPSIS: * * float x, Ci, Si; * * sicif( x, &Si, &Ci ); * * * DESCRIPTION: * * Evaluates the integrals * * 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 integrals 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) * * * ACCURACY: * Test interval = [0,50]. * Absolute error, except relative when > 1: * arithmetic function # trials peak rms * IEEE Si 30000 2.1e-7 4.3e-8 * IEEE Ci 30000 3.9e-7 2.2e-8 */
/* sindgf.c * * Circular sine of angle in degrees * * * * SYNOPSIS: * * float x, y, sindgf(); * * y = sindgf( x ); * * * * DESCRIPTION: * * 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 Q(x**2). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE +-3600 100,000 1.2e-7 3.0e-8 * * ERROR MESSAGES: * * message condition value returned * sin total loss x > 2^24 0.0 * */
/* cosdgf.c * * Circular cosine of angle in degrees * * * * SYNOPSIS: * * float x, y, cosdgf(); * * y = cosdgf( x ); * * * * DESCRIPTION: * * 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 Q(x**2). * Between pi/4 and pi/2 the sine is represented as * x + x**3 P(x**2). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -8192,+8192 100,000 3.0e-7 3.0e-8 */
/* sinf.c * * Circular sine * * * * SYNOPSIS: * * float x, y, sinf(); * * y = sinf( x ); * * * * DESCRIPTION: * * 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). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -4096,+4096 100,000 1.2e-7 3.0e-8 * IEEE -8192,+8192 100,000 3.0e-7 3.0e-8 * * ERROR MESSAGES: * * message condition value returned * sin total loss x > 2^24 0.0 * * Partial loss of accuracy begins to occur at x = 2^13 * = 8192. Results may be meaningless for x >= 2^24 * The routine as implemented flags a TLOSS error * for x >= 2^24 and returns 0.0. */
/* cosf.c * * Circular cosine * * * * SYNOPSIS: * * float x, y, cosf(); * * y = cosf( x ); * * * * DESCRIPTION: * * 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). * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -8192,+8192 100,000 3.0e-7 3.0e-8 */
/* sinhf.c * * Hyperbolic sine * * * * SYNOPSIS: * * float x, y, sinhf(); * * y = sinhf( x ); * * * * DESCRIPTION: * * Returns hyperbolic sine of argument in the range MINLOGF to * MAXLOGF. * * The range is partitioned into two segments. If |x| <= 1, a * polynomial approximation is used. * Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE +-MAXLOG 100000 1.1e-7 2.9e-8 * */
/* spencef.c * * Dilogarithm * * * * SYNOPSIS: * * float x, y, spencef(); * * y = spencef( x ); * * * * DESCRIPTION: * * Computes the integral * * x * - * | | log t * spence(x) = - | ----- dt * | | t - 1 * - * 1 * * for x >= 0. A rational approximation gives the integral in * the interval (0.5, 1.5). Transformation formulas for 1/x * and 1-x are employed outside the basic expansion range. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,4 30000 4.4e-7 6.3e-8 * * */
/* sqrtf.c * * Square root * * * * SYNOPSIS: * * float x, y, sqrtf(); * * y = sqrtf( x ); * * * * DESCRIPTION: * * 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. * * * * ACCURACY: * * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,1.e38 100000 8.7e-8 2.9e-8 * * * ERROR MESSAGES: * * message condition value returned * sqrtf domain x < 0 0.0 * */
/* stdtrf.c * * Student's t distribution * * * * SYNOPSIS: * * float t, stdtrf(); * short k; * * y = stdtrf( k, t ); * * * DESCRIPTION: * * Computes the integral from minus infinity to t of the Student * t distribution 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 incomplete beta integral: * * 1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z ) * where * z = k/(k + t**2). * * For t < -1, 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. * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE +/- 100 5000 2.3e-5 2.9e-6 */
/* struvef.c * * Struve function * * * * SYNOPSIS: * * float v, x, y, struvef(); * * y = struvef( v, x ); * * * * DESCRIPTION: * * Computes the Struve function Hv(x) of order v, argument x. * Negative x is rejected unless v is an integer. * * This module also contains the hypergeometric functions 1F2 * and 3F0 and a routine for the Bessel function Yv(x) with * noninteger v. * * * * ACCURACY: * * v varies from 0 to 10. * Absolute error (relative error when |Hv(x)| > 1): * arithmetic domain # trials peak rms * IEEE -10,10 100000 9.0e-5 4.0e-6 * */
/* tandgf.c * * Circular tangent of angle in degrees * * * * SYNOPSIS: * * float x, y, tandgf(); * * y = tandgf( x ); * * * * DESCRIPTION: * * Returns the circular tangent of the radian argument x. * * Range reduction is into intervals of 45 degrees. * * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE +-2^24 50000 2.4e-7 4.8e-8 * * ERROR MESSAGES: * * message condition value returned * tanf total loss x > 2^24 0.0 * */
/* cotdgf.c * * Circular cotangent of angle in degrees * * * * SYNOPSIS: * * float x, y, cotdgf(); * * y = cotdgf( x ); * * * * DESCRIPTION: * * Range reduction is into intervals of 45 degrees. * A common routine computes either the tangent or cotangent. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE +-2^24 50000 2.4e-7 4.8e-8 * * * ERROR MESSAGES: * * message condition value returned * cot total loss x > 2^24 0.0 * cot singularity x = 0 MAXNUMF * */
/* tanf.c * * Circular tangent * * * * SYNOPSIS: * * float x, y, tanf(); * * y = tanf( x ); * * * * DESCRIPTION: * * Returns the circular tangent of the radian argument x. * * Range reduction is modulo pi/4. A polynomial approximation * is employed in the basic interval [0, pi/4]. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE +-4096 100000 3.3e-7 4.5e-8 * * ERROR MESSAGES: * * message condition value returned * tanf total loss x > 2^24 0.0 * */
/* cotf.c * * Circular cotangent * * * * SYNOPSIS: * * float x, y, cotf(); * * y = cotf( x ); * * * * DESCRIPTION: * * Returns the circular cotangent of the radian argument x. * A common routine computes either the tangent or cotangent. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE +-4096 100000 3.0e-7 4.5e-8 * * * ERROR MESSAGES: * * message condition value returned * cot total loss x > 2^24 0.0 * cot singularity x = 0 MAXNUMF * */
/* tanhf.c * * Hyperbolic tangent * * * * SYNOPSIS: * * float x, y, tanhf(); * * y = tanhf( x ); * * * * DESCRIPTION: * * Returns hyperbolic tangent of argument in the range MINLOG to * MAXLOG. * * A polynomial approximation is used for |x| < 0.625. * Otherwise, * * tanh(x) = sinh(x)/cosh(x) = 1 - 2/(exp(2x) + 1). * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE -2,2 100000 1.3e-7 2.6e-8 * */
/* ynf.c * * Bessel function of second kind of integer order * * * * SYNOPSIS: * * float x, y, ynf(); * int n; * * y = ynf( n, x ); * * * * DESCRIPTION: * * Returns Bessel function 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. * * * * ACCURACY: * * * Absolute error, except relative when y > 1: * * arithmetic domain # trials peak rms * IEEE 0, 30 10000 2.3e-6 3.4e-7 * * * ERROR MESSAGES: * * message condition value returned * yn singularity x = 0 MAXNUMF * yn overflow MAXNUMF * * Spot checked against tables for x, n between 0 and 100. * */
/* zetacf.c * * Riemann zeta function * * * * SYNOPSIS: * * float x, y, zetacf(); * * y = zetacf( x ); * * * * DESCRIPTION: * * * * 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. * * ACCURACY: * * Tabulated values have full machine accuracy. * * Relative error: * arithmetic domain # trials peak rms * IEEE 1,50 30000 5.5e-7 7.5e-8 * * */
/* zetaf.c * * Riemann zeta function of two arguments * * * * SYNOPSIS: * * float x, q, y, zetaf(); * * y = zetaf( x, q ); * * * * DESCRIPTION: * * * * 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. * * * * ACCURACY: * * Relative error: * arithmetic domain # trials peak rms * IEEE 0,25 10000 6.9e-7 1.0e-7 * * Large arguments may produce underflow in powf(), in which * case the results are inaccurate. * * REFERENCE: * * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, * Series, and Products, p. 1073; Academic Press, 1980. * */
To Cephes home page www.moshier.net:
Last update: 10 August 2000


Class: Aws::Redshift::Types::Snapshot
| Inherits: |
Struct
|
| Defined in: |
(unknown) |
Overview
Instance Attribute Summary collapse
-
#accounts_with_restore_access ⇒ Array<Types::AccountWithRestoreAccess>
A list of the AWS customer accounts authorized to restore the snapshot.
-
#actual_incremental_backup_size_in_mega_bytes ⇒ Float
The size of the incremental backup.
-
#availability_zone ⇒ String
The Availability Zone in which the cluster was created.
-
#backup_progress_in_mega_bytes ⇒ Float
The number of megabytes that have been transferred to the snapshot backup.
-
#cluster_create_time ⇒ Time
The time (UTC) when the cluster was originally created.
-
#cluster_identifier ⇒ String
The identifier of the cluster for which the snapshot was taken.
-
#cluster_version ⇒ String
The version ID of the Amazon Redshift engine that is running on the cluster.
-
#current_backup_rate_in_mega_bytes_per_second ⇒ Float
The number of megabytes per second being transferred to the snapshot backup.
-
#db_name ⇒ String
The name of the database that was created when the cluster was created.
-
#elapsed_time_in_seconds ⇒ Integer
The amount of time an in-progress snapshot backup has been running, or the amount of time it took a completed backup to finish.
-
#encrypted ⇒ Boolean
If
true
, the data in the snapshot is encrypted at rest. -
#encrypted_with_hsm ⇒ Boolean
A boolean that indicates whether the snapshot data is encrypted using the HSM keys of the source cluster.
-
#estimated_seconds_to_completion ⇒ Integer
The estimate of the time remaining before the snapshot backup will complete.
-
#kms_key_id ⇒ String
The AWS Key Management Service (KMS) key ID of the encryption key that was used to encrypt data in the cluster from which the snapshot was taken.
-
#master_username ⇒ String
The master user name for the cluster.
-
#node_type ⇒ String
The node type of the nodes in the cluster.
-
#number_of_nodes ⇒ Integer
The number of nodes in the cluster.
-
#owner_account ⇒ String
For manual snapshots, the AWS customer account used to create or copy the snapshot.
-
#port ⇒ Integer
The port that the cluster is listening on.
-
#restorable_node_types ⇒ Array<String>
The list of node types that this cluster snapshot is able to restore into.
-
#snapshot_create_time ⇒ Time
The time (UTC) when Amazon Redshift began the snapshot.
-
#snapshot_identifier ⇒ String
The snapshot identifier that is provided in the request.
-
#snapshot_type ⇒ String
The snapshot type.
-
#source_region ⇒ String
The source region from which the snapshot was copied.
-
#status ⇒ String
The snapshot status.
-
#tags ⇒ Array<Types::Tag>
The list of tags for the cluster snapshot.
-
#total_backup_size_in_mega_bytes ⇒ Float
The size of the complete set of backup data that would be used to restore the cluster.
-
#vpc_id ⇒ String
The VPC identifier of the cluster if the snapshot is from a cluster in a VPC.
Instance Attribute Details
#accounts_with_restore_access ⇒ Array<Types::AccountWithRestoreAccess>
A list of the AWS customer accounts authorized to restore the snapshot.
Returns null
if no accounts are authorized. Visible only to the
snapshot owner.
#actual_incremental_backup_size_in_mega_bytes ⇒ Float
The size of the incremental backup.
#availability_zone ⇒ String
The Availability Zone in which the cluster was created.
#backup_progress_in_mega_bytes ⇒ Float
The number of megabytes that have been transferred to the snapshot backup.
#cluster_create_time ⇒ Time
The time (UTC) when the cluster was originally created.
#cluster_identifier ⇒ String
The identifier of the cluster for which the snapshot was taken.
#cluster_version ⇒ String
The version ID of the Amazon Redshift engine that is running on the cluster.
#current_backup_rate_in_mega_bytes_per_second ⇒ Float
The number of megabytes per second being transferred to the snapshot
backup. Returns 0
for a completed backup.
#db_name ⇒ String
The name of the database that was created when the cluster was created.
#elapsed_time_in_seconds ⇒ Integer
The amount of time an in-progress snapshot backup has been running, or the amount of time it took a completed backup to finish.
#encrypted ⇒ Boolean
If true
, the data in the snapshot is encrypted at rest.
#encrypted_with_hsm ⇒ Boolean
A boolean that indicates whether the snapshot data is encrypted using
the HSM keys of the source cluster. true
indicates that the data is
encrypted using HSM keys.
#estimated_seconds_to_completion ⇒ Integer
The estimate of the time remaining before the snapshot backup will
complete. Returns 0
for a completed backup.
#kms_key_id ⇒ String
The AWS Key Management Service (KMS) key ID of the encryption key that was used to encrypt data in the cluster from which the snapshot was taken.
#master_username ⇒ String
The master user name for the cluster.
#node_type ⇒ String
The node type of the nodes in the cluster.
#number_of_nodes ⇒ Integer
The number of nodes in the cluster.
#owner_account ⇒ String
For manual snapshots, the AWS customer account used to create or copy the snapshot. For automatic snapshots, the owner of the cluster. The owner can perform all snapshot actions, such as sharing a manual snapshot.
#port ⇒ Integer
The port that the cluster is listening on.
#restorable_node_types ⇒ Array<String>
The list of node types that this cluster snapshot is able to restore into.
#snapshot_create_time ⇒ Time
The time (UTC) when Amazon Redshift began the snapshot. A snapshot contains a copy of the cluster data as of this exact time.
#snapshot_identifier ⇒ String
The snapshot identifier that is provided in the request.
#snapshot_type ⇒ String
The snapshot type. Snapshots created using CreateClusterSnapshot and CopyClusterSnapshot will be of type \"manual\".
#source_region ⇒ String
The source region from which the snapshot was copied.
#status ⇒ String
The snapshot status. The value of the status depends on the API operation used. * CreateClusterSnapshot and CopyClusterSnapshot returns status as \"creating\". * DescribeClusterSnapshots returns status as \"creating\", \"available\", \"final snapshot\", or \"failed\". * DeleteClusterSnapshot returns status as \"deleted\".
#tags ⇒ Array<Types::Tag>
The list of tags for the cluster snapshot.
#total_backup_size_in_mega_bytes ⇒ Float
The size of the complete set of backup data that would be used to restore the cluster.
#vpc_id ⇒ String
The VPC identifier of the cluster if the snapshot is from a cluster in a VPC. Otherwise, this field is not in the output.


-
None at all. This might come as a bit of a shocker, but I really don't want your data.
-
There are multiple reasons to not be collecting your data
1. I don't really care much for your personal information.
2. Collecting data requires setting up a secret server and sending data, which I have no idea how to do.
-
Let's not kid ourselves by acting like either of us understands the GDPR guidelines. Also, since no data is being collected, GDPR does not apply. Cool Question.
-
Currently, it doesn't make any money at all. You can support development through PayPal on the Support Development tab in the side navigation bar.
-
I neither own, nor claim to own any of the sounds used in the app. All credits to FOX Media and NBC.
-
There are no plans at the moment. Designing a new app takes a lot of time, which I currently can't spare. However, if you're willing to make one, I'll be glad to assist.
-
Hitchcock and Scully.
-
This is coming in a future update. At the moment, you can write a mail to me with the audio tune you want and I'll be glad to send it to you.
-
You currently cannot. Maybe in a future update.


MenuOptions - input masks and multi-column autocomplete combined
Example of multiple select lists with autocomplete
Change MenuOptions Configuration
Pizza Type
Toppings
Crust thickness
Cheese Amount
Bake Instructions
Delivery Instructions
Driver tip
Show serialize values


Privacy Policy
Schooner (henceforth referred to as website) is a service provided by University of Turku (henceforth referred to as UTU), Department of Computing (henceforth referred to as DFT). This privacy policy will explain how our organization uses the personal data we collect from you when you use our website.
Topics:
What data do we collect?
How do we collect your data?
How will we use your data?
How do we store your data?
Marketing
What are your data protection rights?
What are cookies?
How do we use cookies?
What types of cookies do we use?
How to manage your cookies
Privacy policies of other websites
Changes to our privacy policy
How to contact us
How to contact the appropriate authorities
What data do we collect?
Following data is imported to this service:
Students registered to a course managed by this system:
- UTU Single Sign-On User ID
- Name
- Student ID
- Email Address
This website collects the following data:
From students:
- GitHub account (used to submit course execises)
- Submission details (date and time, score, feedback from evaluating assistant)
From everyone else:
How do we collect your data?
Personal data (listed above) is not collected by this system, but imported from UTU Peppi system, with the exception of GitHub account. For anonymous and student access, this website does not collect nor store any other identifiable personal data. System will collect score and grading information from other systems and add them to user data.
Users identified as administrators or course assistants will agree to have his or her UTU Single Sign-On (henceforth referred to as SSO) ID stored in the website's list of administrators or course assistants permanently and to be used in security audit logging.
How will we use your data?
Student submissions are used to check that required assignments are passed and to calculate course grade when the course is completed. Email is used to send automated messages notifying of deadlines, provided feedback for submissions, and to warn of any problems found in submissions.
Administrators and course assistants have their activities logged for security reasons, in case students note irregularities in their data and/or submissions.
How do we store your data?
This website stores all data into a locally-only accessible RDBMS, in a virtual server instance, securely hosted and maintained by UTU IT-Services.
This website will keep student's personal data for as long as standing orders require to keep them (for audits or examining data for contested results). Student data is removed annually when no longer subject to univeristy archival requirements.
Marketing
This website neither engages in any type of marketing nor share any of the data with external systems or parties.
What are your data protection rights?
University of Turku would like to make sure you are fully aware of all of your data protection rights. Every user is entitled to the following:
The right to access – You have the right to request us for copies of your personal data. We may charge you a small fee for this service. Fees are subject to change without notification, on hourly basis.
The right to rectification – You have the right to request that we correct any information you believe is inaccurate. You also have the right to request us to complete the information you believe is incomplete.
The right to erasure – You have the right to request that this service erase your personal data, under certain conditions.
The right to restrict processing – You have the right to request that this website restrict the processing of your personal data, under certain conditions.
The right to object to processing – You have the right to object to this service processing of your personal data, under certain conditions.
The right to data portability – You have the right to request us to transfer the data that we have collected to another organization, or directly to you, under certain conditions.
If you make a request, we have one month to respond to you. If you would like to exercise any of these rights, please contact the lecturer.
What are cookies
Cookies are text files placed on your computer to collect standard Internet log information and visitor behavior information. When you visit our websites, we may collect information from you automatically through cookies or similar technology.
For further information, visit allaboutcookies.org.
How do we use cookies?
This website uses cookies in a range of ways to improve your experience on our website, including:
What types of cookies do we use?
There are a number of different types of cookies, however, our website uses:
- Functionality – This website uses these cookies so that we can validate your UTU SSO session validity. A mix of site specific and UTU SSO service cookies are used.
How to manage cookies
You can set your browser not to accept cookies, and the above website tells you how to remove cookies from your browser. However, this will render the features of this website non-functional and may prevent successful completion of the course.
Privacy policies of other websites
Our website contains links to other websites. Our privacy policy applies only to our website, so if you click on a link to another website, you should read their privacy policy.
Changes to our privacy policy
Our organization keeps its privacy policy under regular review and places any updates on this web page. This privacy policy was last updated on 17th of August, 2021.
If you have any questions about our privacy policy, the data we hold on you, or you would like to exercise one of your data protection rights, please do not hesitate to contact us.
Email your course instructor, at an address which you can find from the course Moodle page.
Should you wish to report a complaint or if you feel that our organization has not addressed your concern in a satisfactory manner, you may contact the Office of the Data Protection Ombudsman.
Email<EMAIL_ADDRESS>Address:
Ratapihantie 9
6th floor
00520 Helsinki