/* 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 "mconf.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) );
}
syntax highlighted by Code2HTML, v. 0.9.1