From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: Totally rework the math library, this time based on the MacOs X math library (which is itself based on the math lib from FreeBSD). -Erik --- libm/ldouble/clogl.c | 720 --------------------------------------------------- 1 file changed, 720 deletions(-) delete mode 100644 libm/ldouble/clogl.c (limited to 'libm/ldouble/clogl.c') diff --git a/libm/ldouble/clogl.c b/libm/ldouble/clogl.c deleted file mode 100644 index b3e6b25fb..000000000 --- a/libm/ldouble/clogl.c +++ /dev/null @@ -1,720 +0,0 @@ -/* clogl.c - * - * Complex natural logarithm - * - * - * - * SYNOPSIS: - * - * void clogl(); - * cmplxl z, w; - * - * clogl( &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. - */ - -#include -#ifdef ANSIPROT -static void cchshl ( long double x, long double *c, long double *s ); -static long double redupil ( long double x ); -static long double ctansl ( cmplxl *z ); -long double cabsl ( cmplxl *x ); -void csqrtl ( cmplxl *x, cmplxl *y ); -void caddl ( cmplxl *x, cmplxl *y, cmplxl *z ); -extern long double fabsl ( long double ); -extern long double sqrtl ( long double ); -extern long double logl ( long double ); -extern long double expl ( long double ); -extern long double atan2l ( long double, long double ); -extern long double coshl ( long double ); -extern long double sinhl ( long double ); -extern long double asinl ( long double ); -extern long double sinl ( long double ); -extern long double cosl ( long double ); -void clogl ( cmplxl *, cmplxl *); -void casinl ( cmplxl *, cmplxl *); -#else -static void cchshl(); -static long double redupil(); -static long double ctansl(); -long double cabsl(), fabsl(), sqrtl(); -lnog double logl(), expl(), atan2l(), coshl(), sinhl(); -long double asinl(), sinl(), cosl(); -void caddl(), csqrtl(), clogl(), casinl(); -#endif - -extern long double MAXNUML, MACHEPL, PIL, PIO2L; - -void clogl( z, w ) -register cmplxl *z, *w; -{ -long double p, rr; - -/*rr = sqrt( z->r * z->r + z->i * z->i );*/ -rr = cabsl(z); -p = logl(rr); -#if ANSIC -rr = atan2l( z->i, z->r ); -#else -rr = atan2l( z->r, z->i ); -if( rr > PIL ) - rr -= PIL + PIL; -#endif -w->i = rr; -w->r = p; -} - /* cexpl() - * - * Complex exponential function - * - * - * - * SYNOPSIS: - * - * void cexpl(); - * cmplxl z, w; - * - * cexpl( &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 cexpl( z, w ) -register cmplxl *z, *w; -{ -long double r; - -r = expl( z->r ); -w->r = r * cosl( z->i ); -w->i = r * sinl( z->i ); -} - /* csinl() - * - * Complex circular sine - * - * - * - * SYNOPSIS: - * - * void csinl(); - * cmplxl z, w; - * - * csinl( &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 csinl( z, w ) -register cmplxl *z, *w; -{ -long double ch, sh; - -cchshl( z->i, &ch, &sh ); -w->r = sinl( z->r ) * ch; -w->i = cosl( z->r ) * sh; -} - - - -/* calculate cosh and sinh */ - -static void cchshl( x, c, s ) -long double x, *c, *s; -{ -long double e, ei; - -if( fabsl(x) <= 0.5L ) - { - *c = coshl(x); - *s = sinhl(x); - } -else - { - e = expl(x); - ei = 0.5L/e; - e = 0.5L * e; - *s = e - ei; - *c = e + ei; - } -} - - /* ccosl() - * - * Complex circular cosine - * - * - * - * SYNOPSIS: - * - * void ccosl(); - * cmplxl z, w; - * - * ccosl( &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 ccosl( z, w ) -register cmplxl *z, *w; -{ -long double ch, sh; - -cchshl( z->i, &ch, &sh ); -w->r = cosl( z->r ) * ch; -w->i = -sinl( z->r ) * sh; -} - /* ctanl() - * - * Complex circular tangent - * - * - * - * SYNOPSIS: - * - * void ctanl(); - * cmplxl z, w; - * - * ctanl( &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 ctanl( z, w ) -register cmplxl *z, *w; -{ -long double d; - -d = cosl( 2.0L * z->r ) + coshl( 2.0L * z->i ); - -if( fabsl(d) < 0.25L ) - d = ctansl(z); - -if( d == 0.0L ) - { - mtherr( "ctan", OVERFLOW ); - w->r = MAXNUML; - w->i = MAXNUML; - return; - } - -w->r = sinl( 2.0L * z->r ) / d; -w->i = sinhl( 2.0L * z->i ) / d; -} - /* ccotl() - * - * Complex circular cotangent - * - * - * - * SYNOPSIS: - * - * void ccotl(); - * cmplxl z, w; - * - * ccotl( &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 ccotl( z, w ) -register cmplxl *z, *w; -{ -long double d; - -d = coshl(2.0L * z->i) - cosl(2.0L * z->r); - -if( fabsl(d) < 0.25L ) - d = ctansl(z); - -if( d == 0.0L ) - { - mtherr( "ccot", OVERFLOW ); - w->r = MAXNUML; - w->i = MAXNUML; - return; - } - -w->r = sinl( 2.0L * z->r ) / d; -w->i = -sinhl( 2.0L * 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 long double redupil(x) -long double x; -{ -long double t; -long i; - -t = x/PIL; -if( t >= 0.0L ) - t += 0.5L; -else - t -= 0.5L; - -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 long double ctansl(z) -cmplxl *z; -{ -long double f, x, x2, y, y2, rn, t; -long double d; - -x = fabsl( 2.0L * z->r ); -y = fabsl( 2.0L * z->i ); - -x = redupil(x); - -x = x * x; -y = y * y; -x2 = 1.0L; -y2 = 1.0L; -f = 1.0L; -rn = 0.0; -d = 0.0; -do - { - rn += 1.0L; - f *= rn; - rn += 1.0L; - f *= rn; - x2 *= x; - y2 *= y; - t = y2 + x2; - t /= f; - d += t; - - rn += 1.0L; - f *= rn; - rn += 1.0L; - f *= rn; - x2 *= x; - y2 *= y; - t = y2 - x2; - t /= f; - d += t; - } -while( fabsl(t/d) > MACHEPL ); -return(d); -} - /* casinl() - * - * Complex circular arc sine - * - * - * - * SYNOPSIS: - * - * void casinl(); - * cmplxl z, w; - * - * casinl( &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 casinl( z, w ) -cmplxl *z, *w; -{ -static cmplxl ca, ct, zz, z2; -long double x, y; - -x = z->r; -y = z->i; - -if( y == 0.0L ) - { - if( fabsl(x) > 1.0L ) - { - w->r = PIO2L; - w->i = 0.0L; - mtherr( "casinl", DOMAIN ); - } - else - { - w->r = asinl(x); - w->i = 0.0L; - } - return; - } - -/* Power series expansion */ -/* -b = cabsl(z); -if( b < 0.125L ) -{ -z2.r = (x - y) * (x + y); -z2.i = 2.0L * x * y; - -cn = 1.0L; -n = 1.0L; -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.0L; - cn /= n; - n += 1.0L; - b = cn/n; - - ct.r *= b; - ct.i *= b; - sum.r += ct.r; - sum.i += ct.i; - b = fabsl(ct.r) + fabs(ct.i); - } -while( b > MACHEPL ); -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.0L * ca.r * ca.i; - -zz.r = 1.0L - zz.r; -zz.i = -zz.i; -csqrtl( &zz, &z2 ); - -caddl( &z2, &ct, &zz ); -clogl( &zz, &zz ); -w->r = zz.i; /* mult by 1/i = -i */ -w->i = -zz.r; -return; -} - /* cacosl() - * - * Complex circular arc cosine - * - * - * - * SYNOPSIS: - * - * void cacosl(); - * cmplxl z, w; - * - * cacosl( &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 cacosl( z, w ) -cmplxl *z, *w; -{ - -casinl( z, w ); -w->r = PIO2L - w->r; -w->i = -w->i; -} - /* catanl() - * - * Complex circular arc tangent - * - * - * - * SYNOPSIS: - * - * void catanl(); - * cmplxl z, w; - * - * catanl( &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 catanl( z, w ) -cmplxl *z, *w; -{ -long double a, t, x, x2, y; - -x = z->r; -y = z->i; - -if( (x == 0.0L) && (y > 1.0L) ) - goto ovrf; - -x2 = x * x; -a = 1.0L - x2 - (y * y); -if( a == 0.0L ) - goto ovrf; - -#if ANSIC -t = atan2l( 2.0L * x, a ) * 0.5L; -#else -t = atan2l( a, 2.0 * x ) * 0.5L; -#endif -w->r = redupil( t ); - -t = y - 1.0L; -a = x2 + (t * t); -if( a == 0.0L ) - goto ovrf; - -t = y + 1.0L; -a = (x2 + (t * t))/a; -w->i = logl(a)/4.0; -return; - -ovrf: -mtherr( "catanl", OVERFLOW ); -w->r = MAXNUML; -w->i = MAXNUML; -} -- cgit v1.2.3