From 1077fa4d772832f77a677ce7fb7c2d513b959e3f Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 10 May 2001 00:40:28 +0000 Subject: uClibc now has a math library. muahahahaha! -Erik --- libm/double/clog.c | 1043 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1043 insertions(+) create mode 100644 libm/double/clog.c (limited to 'libm/double/clog.c') diff --git a/libm/double/clog.c b/libm/double/clog.c new file mode 100644 index 000000000..70a318a50 --- /dev/null +++ b/libm/double/clog.c @@ -0,0 +1,1043 @@ +/* clog.c + * + * Complex natural logarithm + * + * + * + * SYNOPSIS: + * + * void clog(); + * cmplx z, w; + * + * clog( &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 + * DEC -10,+10 7000 8.5e-17 1.9e-17 + * IEEE -10,+10 30000 5.0e-15 1.1e-16 + * + * Larger relative error can be observed for z near 1 +i0. + * In IEEE arithmetic the peak absolute error is 5.2e-16, rms + * absolute error 1.0e-16. + */ + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1984, 1995, 2000 by Stephen L. Moshier +*/ +#include +#ifdef ANSIPROT +static void cchsh ( double x, double *c, double *s ); +static double redupi ( double x ); +static double ctans ( cmplx *z ); +/* These are supposed to be in some standard place. */ +double fabs (double); +double sqrt (double); +double pow (double, double); +double log (double); +double exp (double); +double atan2 (double, double); +double cosh (double); +double sinh (double); +double asin (double); +double sin (double); +double cos (double); +double cabs (cmplx *); +void cadd ( cmplx *, cmplx *, cmplx * ); +void cmul ( cmplx *, cmplx *, cmplx * ); +void csqrt ( cmplx *, cmplx * ); +static void cchsh ( double, double *, double * ); +static double redupi ( double ); +static double ctans ( cmplx * ); +void clog ( cmplx *, cmplx * ); +void casin ( cmplx *, cmplx * ); +void cacos ( cmplx *, cmplx * ); +void catan ( cmplx *, cmplx * ); +#else +static void cchsh(); +static double redupi(); +static double ctans(); +double cabs(), fabs(), sqrt(), pow(); +double log(), exp(), atan2(), cosh(), sinh(); +double asin(), sin(), cos(); +void cadd(), cmul(), csqrt(); +void clog(), casin(), cacos(), catan(); +#endif + + +extern double MAXNUM, MACHEP, PI, PIO2; + +void clog( z, w ) +register cmplx *z, *w; +{ +double p, rr; + +/*rr = sqrt( z->r * z->r + z->i * z->i );*/ +rr = cabs(z); +p = log(rr); +#if ANSIC +rr = atan2( z->i, z->r ); +#else +rr = atan2( z->r, z->i ); +if( rr > PI ) + rr -= PI + PI; +#endif +w->i = rr; +w->r = p; +} + /* cexp() + * + * Complex exponential function + * + * + * + * SYNOPSIS: + * + * void cexp(); + * cmplx z, w; + * + * cexp( &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 + * DEC -10,+10 8700 3.7e-17 1.1e-17 + * IEEE -10,+10 30000 3.0e-16 8.7e-17 + * + */ + +void cexp( z, w ) +register cmplx *z, *w; +{ +double r; + +r = exp( z->r ); +w->r = r * cos( z->i ); +w->i = r * sin( z->i ); +} + /* csin() + * + * Complex circular sine + * + * + * + * SYNOPSIS: + * + * void csin(); + * cmplx z, w; + * + * csin( &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 + * 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 csin( z, w ) +register cmplx *z, *w; +{ +double ch, sh; + +cchsh( z->i, &ch, &sh ); +w->r = sin( z->r ) * ch; +w->i = cos( z->r ) * sh; +} + + + +/* calculate cosh and sinh */ + +static void cchsh( x, c, s ) +double x, *c, *s; +{ +double e, ei; + +if( fabs(x) <= 0.5 ) + { + *c = cosh(x); + *s = sinh(x); + } +else + { + e = exp(x); + ei = 0.5/e; + e = 0.5 * e; + *s = e - ei; + *c = e + ei; + } +} + + /* ccos() + * + * Complex circular cosine + * + * + * + * SYNOPSIS: + * + * void ccos(); + * cmplx z, w; + * + * ccos( &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 + * DEC -10,+10 8400 4.5e-17 1.3e-17 + * IEEE -10,+10 30000 3.8e-16 1.0e-16 + */ + +void ccos( z, w ) +register cmplx *z, *w; +{ +double ch, sh; + +cchsh( z->i, &ch, &sh ); +w->r = cos( z->r ) * ch; +w->i = -sin( z->r ) * sh; +} + /* ctan() + * + * Complex circular tangent + * + * + * + * SYNOPSIS: + * + * void ctan(); + * cmplx z, w; + * + * ctan( &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 + * 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 ctan( z, w ) +register cmplx *z, *w; +{ +double d; + +d = cos( 2.0 * z->r ) + cosh( 2.0 * z->i ); + +if( fabs(d) < 0.25 ) + d = ctans(z); + +if( d == 0.0 ) + { + mtherr( "ctan", OVERFLOW ); + w->r = MAXNUM; + w->i = MAXNUM; + return; + } + +w->r = sin( 2.0 * z->r ) / d; +w->i = sinh( 2.0 * z->i ) / d; +} + /* ccot() + * + * Complex circular cotangent + * + * + * + * SYNOPSIS: + * + * void ccot(); + * cmplx z, w; + * + * ccot( &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 + * 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 ccot( z, w ) +register cmplx *z, *w; +{ +double d; + +d = cosh(2.0 * z->i) - cos(2.0 * z->r); + +if( fabs(d) < 0.25 ) + d = ctans(z); + +if( d == 0.0 ) + { + mtherr( "ccot", OVERFLOW ); + w->r = MAXNUM; + w->i = MAXNUM; + return; + } + +w->r = sin( 2.0 * z->r ) / d; +w->i = -sinh( 2.0 * z->i ) / d; +} + +/* Program to subtract nearest integer multiple of PI */ +/* extended precision value of PI: */ +#ifdef UNK +static double DP1 = 3.14159265160560607910E0; +static double DP2 = 1.98418714791870343106E-9; +static double DP3 = 1.14423774522196636802E-17; +#endif + +#ifdef DEC +static unsigned short P1[] = {0040511,0007732,0120000,0000000,}; +static unsigned short P2[] = {0031010,0055060,0100000,0000000,}; +static unsigned short P3[] = {0022123,0011431,0105056,0001560,}; +#define DP1 *(double *)P1 +#define DP2 *(double *)P2 +#define DP3 *(double *)P3 +#endif + +#ifdef IBMPC +static unsigned short P1[] = {0x0000,0x5400,0x21fb,0x4009}; +static unsigned short P2[] = {0x0000,0x1000,0x0b46,0x3e21}; +static unsigned short P3[] = {0xc06e,0x3145,0x6263,0x3c6a}; +#define DP1 *(double *)P1 +#define DP2 *(double *)P2 +#define DP3 *(double *)P3 +#endif + +#ifdef MIEEE +static unsigned short P1[] = { +0x4009,0x21fb,0x5400,0x0000 +}; +static unsigned short P2[] = { +0x3e21,0x0b46,0x1000,0x0000 +}; +static unsigned short P3[] = { +0x3c6a,0x6263,0x3145,0xc06e +}; +#define DP1 *(double *)P1 +#define DP2 *(double *)P2 +#define DP3 *(double *)P3 +#endif + +static double redupi(x) +double x; +{ +double t; +long i; + +t = x/PI; +if( t >= 0.0 ) + t += 0.5; +else + t -= 0.5; + +i = t; /* the multiple */ +t = i; +t = ((x - t * DP1) - t * DP2) - t * DP3; +return(t); +} + +/* Taylor series expansion for cosh(2y) - cos(2x) */ + +static double ctans(z) +cmplx *z; +{ +double f, x, x2, y, y2, rn, t; +double d; + +x = fabs( 2.0 * z->r ); +y = fabs( 2.0 * z->i ); + +x = redupi(x); + +x = x * x; +y = y * y; +x2 = 1.0; +y2 = 1.0; +f = 1.0; +rn = 0.0; +d = 0.0; +do + { + rn += 1.0; + f *= rn; + rn += 1.0; + f *= rn; + x2 *= x; + y2 *= y; + t = y2 + x2; + t /= f; + d += t; + + rn += 1.0; + f *= rn; + rn += 1.0; + f *= rn; + x2 *= x; + y2 *= y; + t = y2 - x2; + t /= f; + d += t; + } +while( fabs(t/d) > MACHEP ); +return(d); +} + /* casin() + * + * Complex circular arc sine + * + * + * + * SYNOPSIS: + * + * void casin(); + * cmplx z, w; + * + * casin( &z, &w ); + * + * + * + * DESCRIPTION: + * + * Inverse complex sine: + * + * 2 + * w = -i clog( iz + csqrt( 1 - z ) ). + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 10100 2.1e-15 3.4e-16 + * IEEE -10,+10 30000 2.2e-14 2.7e-15 + * Larger relative error can be observed for z near zero. + * Also tested by csin(casin(z)) = z. + */ + +void casin( z, w ) +cmplx *z, *w; +{ +static cmplx ca, ct, zz, z2; +double x, y; + +x = z->r; +y = z->i; + +if( y == 0.0 ) + { + if( fabs(x) > 1.0 ) + { + w->r = PIO2; + w->i = 0.0; + mtherr( "casin", DOMAIN ); + } + else + { + w->r = asin(x); + w->i = 0.0; + } + return; + } + +/* Power series expansion */ +/* +b = cabs(z); +if( b < 0.125 ) +{ +z2.r = (x - y) * (x + y); +z2.i = 2.0 * x * y; + +cn = 1.0; +n = 1.0; +ca.r = x; +ca.i = y; +sum.r = x; +sum.i = y; +do + { + ct.r = z2.r * ca.r - z2.i * ca.i; + ct.i = z2.r * ca.i + z2.i * ca.r; + ca.r = ct.r; + ca.i = ct.i; + + cn *= n; + n += 1.0; + cn /= n; + n += 1.0; + b = cn/n; + + ct.r *= b; + ct.i *= b; + sum.r += ct.r; + sum.i += ct.i; + b = fabs(ct.r) + fabs(ct.i); + } +while( b > MACHEP ); +w->r = sum.r; +w->i = sum.i; +return; +} +*/ + + +ca.r = x; +ca.i = y; + +ct.r = -ca.i; /* iz */ +ct.i = ca.r; + + /* sqrt( 1 - z*z) */ +/* cmul( &ca, &ca, &zz ) */ +zz.r = (ca.r - ca.i) * (ca.r + ca.i); /*x * x - y * y */ +zz.i = 2.0 * ca.r * ca.i; + +zz.r = 1.0 - zz.r; +zz.i = -zz.i; +csqrt( &zz, &z2 ); + +cadd( &z2, &ct, &zz ); +clog( &zz, &zz ); +w->r = zz.i; /* mult by 1/i = -i */ +w->i = -zz.r; +return; +} + /* cacos() + * + * Complex circular arc cosine + * + * + * + * SYNOPSIS: + * + * void cacos(); + * cmplx z, w; + * + * cacos( &z, &w ); + * + * + * + * DESCRIPTION: + * + * + * w = arccos z = PI/2 - arcsin z. + * + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,+10 5200 1.6e-15 2.8e-16 + * IEEE -10,+10 30000 1.8e-14 2.2e-15 + */ + +void cacos( z, w ) +cmplx *z, *w; +{ + +casin( z, w ); +w->r = PIO2 - w->r; +w->i = -w->i; +} + /* catan() + * + * Complex circular arc tangent + * + * + * + * SYNOPSIS: + * + * void catan(); + * cmplx 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 + * 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 catan( z, w ) +cmplx *z, *w; +{ +double a, t, x, x2, y; + +x = z->r; +y = z->i; + +if( (x == 0.0) && (y > 1.0) ) + goto ovrf; + +x2 = x * x; +a = 1.0 - x2 - (y * y); +if( a == 0.0 ) + goto ovrf; + +#if ANSIC +t = atan2( 2.0 * x, a )/2.0; +#else +t = atan2( a, 2.0 * x )/2.0; +#endif +w->r = redupi( t ); + +t = y - 1.0; +a = x2 + (t * t); +if( a == 0.0 ) + goto ovrf; + +t = y + 1.0; +a = (x2 + (t * t))/a; +w->i = log(a)/4.0; +return; + +ovrf: +mtherr( "catan", OVERFLOW ); +w->r = MAXNUM; +w->i = MAXNUM; +} + + +/* csinh + * + * Complex hyperbolic sine + * + * + * + * SYNOPSIS: + * + * void csinh(); + * cmplx z, w; + * + * csinh( &z, &w ); + * + * + * DESCRIPTION: + * + * csinh z = (cexp(z) - cexp(-z))/2 + * = sinh x * cos y + i cosh x * sin y . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 3.1e-16 8.2e-17 + * + */ + +void +csinh (z, w) + cmplx *z, *w; +{ + double x, y; + + x = z->r; + y = z->i; + w->r = sinh (x) * cos (y); + w->i = cosh (x) * sin (y); +} + + +/* casinh + * + * Complex inverse hyperbolic sine + * + * + * + * SYNOPSIS: + * + * void casinh(); + * cmplx z, w; + * + * casinh (&z, &w); + * + * + * + * DESCRIPTION: + * + * casinh z = -i casin iz . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.8e-14 2.6e-15 + * + */ + +void +casinh (z, w) + cmplx *z, *w; +{ + cmplx u; + + u.r = 0.0; + u.i = 1.0; + cmul( z, &u, &u ); + casin( &u, w ); + u.r = 0.0; + u.i = -1.0; + cmul( &u, w, w ); +} + +/* ccosh + * + * Complex hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * void ccosh(); + * cmplx z, w; + * + * ccosh (&z, &w); + * + * + * + * DESCRIPTION: + * + * ccosh(z) = cosh x cos y + i sinh x sin y . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.9e-16 8.1e-17 + * + */ + +void +ccosh (z, w) + cmplx *z, *w; +{ + double x, y; + + x = z->r; + y = z->i; + w->r = cosh (x) * cos (y); + w->i = sinh (x) * sin (y); +} + + +/* cacosh + * + * Complex inverse hyperbolic cosine + * + * + * + * SYNOPSIS: + * + * void cacosh(); + * cmplx z, w; + * + * cacosh (&z, &w); + * + * + * + * DESCRIPTION: + * + * acosh z = i acos z . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.6e-14 2.1e-15 + * + */ + +void +cacosh (z, w) + cmplx *z, *w; +{ + cmplx u; + + cacos( z, w ); + u.r = 0.0; + u.i = 1.0; + cmul( &u, w, w ); +} + + +/* ctanh + * + * Complex hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * void ctanh(); + * cmplx z, w; + * + * ctanh (&z, &w); + * + * + * + * DESCRIPTION: + * + * tanh z = (sinh 2x + i sin 2y) / (cosh 2x + cos 2y) . + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 1.7e-14 2.4e-16 + * + */ + +/* 5.253E-02,1.550E+00 1.643E+01,6.553E+00 1.729E-14 21355 */ + +void +ctanh (z, w) + cmplx *z, *w; +{ + double x, y, d; + + x = z->r; + y = z->i; + d = cosh (2.0 * x) + cos (2.0 * y); + w->r = sinh (2.0 * x) / d; + w->i = sin (2.0 * y) / d; + return; +} + + +/* catanh + * + * Complex inverse hyperbolic tangent + * + * + * + * SYNOPSIS: + * + * void catanh(); + * cmplx z, w; + * + * catanh (&z, &w); + * + * + * + * DESCRIPTION: + * + * Inverse tanh, equal to -i catan (iz); + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 2.3e-16 6.2e-17 + * + */ + +void +catanh (z, w) + cmplx *z, *w; +{ + cmplx u; + + u.r = 0.0; + u.i = 1.0; + cmul (z, &u, &u); /* i z */ + catan (&u, w); + u.r = 0.0; + u.i = -1.0; + cmul (&u, w, w); /* -i catan iz */ + return; +} + + +/* cpow + * + * Complex power function + * + * + * + * SYNOPSIS: + * + * void cpow(); + * cmplx a, z, w; + * + * cpow (&a, &z, &w); + * + * + * + * DESCRIPTION: + * + * 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)). + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -10,+10 30000 9.4e-15 1.5e-15 + * + */ + + +void +cpow (a, z, w) + cmplx *a, *z, *w; +{ + double x, y, r, theta, absa, arga; + + x = z->r; + y = z->i; + absa = cabs (a); + if (absa == 0.0) + { + w->r = 0.0; + w->i = 0.0; + return; + } + arga = atan2 (a->i, a->r); + r = pow (absa, x); + theta = x * arga; + if (y != 0.0) + { + r = r * exp (-y * arga); + theta = theta + y * log (absa); + } + w->r = r * cos (theta); + w->i = r * sin (theta); + return; +} -- cgit v1.2.3