uclibc-ng/libm/ldouble/cmplxl.c

/*							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, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/


#include <math.h>

/*
typedef struct
	{
	long double r;
	long double i;
	}cmplxl;
*/

#ifdef ANSIPROT
extern 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 * );
#else
long double fabsl(), cabsl(), sqrtl(), atan2l(), cosl(), sinl();
long double frexpl(), ldexpl();
int isnanl();
void cdivl(), caddl();
#endif


extern 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, 1989
Copyright 1984, 1987, 1989 by Stephen L. Moshier
Direct 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
#endif


long 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 NANS
if( isnanl(z->r) )
  return(z->r);
if( isnanl(z->i) )
  return(z->i);
#endif

re = 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) );
}