/*							cmplxl.c * *	Complex number arithmetic * * * * SYNOPSIS: * * typedef struct { *      long double r;     real part *      long double i;     imaginary part *     }cmplxl; * * cmplxl *a, *b, *c; * * caddl( a, b, c );     c = b + a * csubl( a, b, c );     c = b - a * cmull( a, b, c );     c = b * a * cdivl( a, b, c );     c = b / a * cnegl( c );           c = -c * cmovl( 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 *    DEC        cadd       10000       1.4e-17     3.4e-18 *    IEEE       cadd      100000       1.1e-16     2.7e-17 *    DEC        csub       10000       1.4e-17     4.5e-18 *    IEEE       csub      100000       1.1e-16     3.4e-17 *    DEC        cmul        3000       2.3e-17     8.7e-18 *    IEEE       cmul      100000       2.1e-16     6.9e-17 *    DEC        cdiv       18000       4.9e-17     1.3e-17 *    IEEE       cdiv      100000       3.7e-16     1.1e-16 *//*				cmplx.c * complex number arithmetic *//*Cephes Math Library Release 2.3:  March, 1995Copyright 1984, 1995 by Stephen L. Moshier*/#include <math.h>/*typedef struct	{	long double r;	long double i;	}cmplxl;*/#ifdef ANSIPROTextern long double fabsl ( long double );extern long double cabsl ( cmplxl * );extern long double sqrtl ( long double );extern long double atan2l ( long double, long double );extern long double cosl ( long double );extern long double sinl ( long double );extern long double frexpl ( long double, int * );extern long double ldexpl ( long double, int );extern int isnanl ( long double );void cdivl ( cmplxl *, cmplxl *, cmplxl * );void caddl ( cmplxl *, cmplxl *, cmplxl * );#elselong double fabsl(), cabsl(), sqrtl(), atan2l(), cosl(), sinl();long double frexpl(), ldexpl();int isnanl();void cdivl(), caddl();#endifextern double MAXNUML, MACHEPL, PIL, PIO2L, INFINITYL, NANL;cmplx czerol = {0.0L, 0.0L};cmplx conel = {1.0L, 0.0L};/*	c = b + a	*/void caddl( a, b, c )register cmplxl *a, *b;cmplxl *c;{c->r = b->r + a->r;c->i = b->i + a->i;}/*	c = b - a	*/void csubl( a, b, c )register cmplxl *a, *b;cmplxl *c;{c->r = b->r - a->r;c->i = b->i - a->i;}/*	c = b * a */void cmull( a, b, c )register cmplxl *a, *b;cmplxl *c;{long double y;y    = b->r * a->r  -  b->i * a->i;c->i = b->r * a->i  +  b->i * a->r;c->r = y;}/*	c = b / a */void cdivl( a, b, c )register cmplxl *a, *b;cmplxl *c;{long double y, p, q, w;y = a->r * a->r  +  a->i * a->i;p = b->r * a->r  +  b->i * a->i;q = b->i * a->r  -  b->r * a->i;if( y < 1.0L )	{	w = MAXNUML * y;	if( (fabsl(p) > w) || (fabsl(q) > w) || (y == 0.0L) )		{		c->r = INFINITYL;		c->i = INFINITYL;		mtherr( "cdivl", OVERFLOW );		return;		}	}c->r = p/y;c->i = q/y;}/*	b = a   Caution, a `short' is assumed to be 16 bits wide.  */void cmovl( a, b )void *a, *b;{register short *pa, *pb;int i;pa = (short *) a;pb = (short *) b;i = 16;do	*pb++ = *pa++;while( --i );}void cnegl( a )register cmplxl *a;{a->r = -a->r;a->i = -a->i;}/*							cabsl() * *	Complex absolute value * * * * SYNOPSIS: * * long double cabsl(); * cmplxl z; * long double a; * * a = cabs( &z ); * * * * DESCRIPTION: * * * If z = x + iy * * then * *       a = sqrt( x**2 + y**2 ). *  * Overflow and underflow are avoided by testing the magnitudes * of x and y before squaring.  If either is outside half of * the floating point full scale range, both are rescaled. * * * ACCURACY: * *                      Relative error: * arithmetic   domain     # trials      peak         rms *    DEC       -30,+30     30000       3.2e-17     9.2e-18 *    IEEE      -10,+10    100000       2.7e-16     6.9e-17 *//*Cephes Math Library Release 2.1:  January, 1989Copyright 1984, 1987, 1989 by Stephen L. MoshierDirect inquiries to 30 Frost Street, Cambridge, MA 02140*//*typedef struct	{	long double r;	long double i;	}cmplxl;*/#ifdef UNK#define PRECL 32#define MAXEXPL 16384#define MINEXPL -16384#endif#ifdef IBMPC#define PRECL 32#define MAXEXPL 16384#define MINEXPL -16384#endif#ifdef MIEEE#define PRECL 32#define MAXEXPL 16384#define MINEXPL -16384#endiflong double cabsl( z )register cmplxl *z;{long double x, y, b, re, im;int ex, ey, e;#ifdef INFINITIES/* Note, cabs(INFINITY,NAN) = INFINITY. */if( z->r == INFINITYL || z->i == INFINITYL   || z->r == -INFINITYL || z->i == -INFINITYL )  return( INFINITYL );#endif#ifdef NANSif( isnanl(z->r) )  return(z->r);if( isnanl(z->i) )  return(z->i);#endifre = fabsl( z->r );im = fabsl( z->i );if( re == 0.0 )	return( im );if( im == 0.0 )	return( re );/* Get the exponents of the numbers */x = frexpl( re, &ex );y = frexpl( im, &ey );/* Check if one number is tiny compared to the other */e = ex - ey;if( e > PRECL )	return( re );if( e < -PRECL )	return( im );/* Find approximate exponent e of the geometric mean. */e = (ex + ey) >> 1;/* Rescale so mean is about 1 */x = ldexpl( re, -e );y = ldexpl( im, -e );		/* Hypotenuse of the right triangle */b = sqrtl( x * x  +  y * y );/* Compute the exponent of the answer. */y = frexpl( b, &ey );ey = e + ey;/* Check it for overflow and underflow. */if( ey > MAXEXPL )	{	mtherr( "cabsl", OVERFLOW );	return( INFINITYL );	}if( ey < MINEXPL )	return(0.0L);/* Undo the scaling */b = ldexpl( b, e );return( b );}/*							csqrtl() * *	Complex square root * * * * SYNOPSIS: * * void csqrtl(); * cmplxl z, w; * * csqrtl( &z, &w ); * * * * DESCRIPTION: * * * If z = x + iy,  r = |z|, then * *                       1/2 * Im w  =  [ (r - x)/2 ]   , * * Re w  =  y / 2 Im w. * * * Note that -w is also a square root of z.  The root chosen * is always in the upper half plane. * * Because of the potential for cancellation error in r - x, * the result is sharpened by doing a Heron iteration * (see sqrt.c) in complex arithmetic. * * * * ACCURACY: * *                      Relative error: * arithmetic   domain     # trials      peak         rms *    DEC       -10,+10     25000       3.2e-17     9.6e-18 *    IEEE      -10,+10    100000       3.2e-16     7.7e-17 * *                        2 * Also tested by csqrt( z ) = z, and tested by arguments * close to the real axis. */void csqrtl( z, w )cmplxl *z, *w;{cmplxl q, s;long double x, y, r, t;x = z->r;y = z->i;if( y == 0.0L )	{	if( x < 0.0L )		{		w->r = 0.0L;		w->i = sqrtl(-x);		return;		}	else		{		w->r = sqrtl(x);		w->i = 0.0L;		return;		}	}if( x == 0.0L )	{	r = fabsl(y);	r = sqrtl(0.5L*r);	if( y > 0.0L )		w->r = r;	else		w->r = -r;	w->i = r;	return;	}/* Approximate  sqrt(x^2+y^2) - x  =  y^2/2x - y^4/24x^3 + ... . * The relative error in the first term is approximately y^2/12x^2 . */if( (fabsl(y) < 2.e-4L * fabsl(x))   && (x > 0) )	{	t = 0.25L*y*(y/x);	}else	{	r = cabsl(z);	t = 0.5L*(r - x);	}r = sqrtl(t);q.i = r;q.r = y/(2.0L*r);/* Heron iteration in complex arithmetic */cdivl( &q, z, &s );caddl( &q, &s, w );w->r *= 0.5L;w->i *= 0.5L;cdivl( &q, z, &s );caddl( &q, &s, w );w->r *= 0.5L;w->i *= 0.5L;}long double hypotl( x, y )long double x, y;{cmplxl z;z.r = x;z.i = y;return( cabsl(&z) );}