summaryrefslogtreecommitdiff
path: root/libm/double/cmplx.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/cmplx.c')
-rw-r--r--libm/double/cmplx.c461
1 files changed, 461 insertions, 0 deletions
diff --git a/libm/double/cmplx.c b/libm/double/cmplx.c
new file mode 100644
index 000000000..dcd972bea
--- /dev/null
+++ b/libm/double/cmplx.c
@@ -0,0 +1,461 @@
+/* cmplx.c
+ *
+ * Complex number arithmetic
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * typedef struct {
+ * double r; real part
+ * double i; imaginary part
+ * }cmplx;
+ *
+ * cmplx *a, *b, *c;
+ *
+ * cadd( a, b, c ); c = b + a
+ * csub( a, b, c ); c = b - a
+ * cmul( a, b, c ); c = b * a
+ * cdiv( a, b, c ); c = b / a
+ * cneg( c ); c = -c
+ * cmov( 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.8: June, 2000
+Copyright 1984, 1995, 2000 by Stephen L. Moshier
+*/
+
+
+#include <math.h>
+
+#ifdef ANSIPROT
+extern double fabs ( double );
+extern double cabs ( cmplx * );
+extern double sqrt ( double );
+extern double atan2 ( double, double );
+extern double cos ( double );
+extern double sin ( double );
+extern double sqrt ( double );
+extern double frexp ( double, int * );
+extern double ldexp ( double, int );
+int isnan ( double );
+void cdiv ( cmplx *, cmplx *, cmplx * );
+void cadd ( cmplx *, cmplx *, cmplx * );
+#else
+double fabs(), cabs(), sqrt(), atan2(), cos(), sin();
+double sqrt(), frexp(), ldexp();
+int isnan();
+void cdiv(), cadd();
+#endif
+
+extern double MAXNUM, MACHEP, PI, PIO2, INFINITY, NAN;
+/*
+typedef struct
+ {
+ double r;
+ double i;
+ }cmplx;
+*/
+cmplx czero = {0.0, 0.0};
+extern cmplx czero;
+cmplx cone = {1.0, 0.0};
+extern cmplx cone;
+
+/* c = b + a */
+
+void cadd( a, b, c )
+register cmplx *a, *b;
+cmplx *c;
+{
+
+c->r = b->r + a->r;
+c->i = b->i + a->i;
+}
+
+
+/* c = b - a */
+
+void csub( a, b, c )
+register cmplx *a, *b;
+cmplx *c;
+{
+
+c->r = b->r - a->r;
+c->i = b->i - a->i;
+}
+
+/* c = b * a */
+
+void cmul( a, b, c )
+register cmplx *a, *b;
+cmplx *c;
+{
+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 cdiv( a, b, c )
+register cmplx *a, *b;
+cmplx *c;
+{
+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.0 )
+ {
+ w = MAXNUM * y;
+ if( (fabs(p) > w) || (fabs(q) > w) || (y == 0.0) )
+ {
+ c->r = MAXNUM;
+ c->i = MAXNUM;
+ mtherr( "cdiv", OVERFLOW );
+ return;
+ }
+ }
+c->r = p/y;
+c->i = q/y;
+}
+
+
+/* b = a
+ Caution, a `short' is assumed to be 16 bits wide. */
+
+void cmov( a, b )
+void *a, *b;
+{
+register short *pa, *pb;
+int i;
+
+pa = (short *) a;
+pb = (short *) b;
+i = 8;
+do
+ *pb++ = *pa++;
+while( --i );
+}
+
+
+void cneg( a )
+register cmplx *a;
+{
+
+a->r = -a->r;
+a->i = -a->i;
+}
+
+/* cabs()
+ *
+ * Complex absolute value
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double cabs();
+ * cmplx z;
+ * 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
+ {
+ double r;
+ double i;
+ }cmplx;
+*/
+
+#ifdef UNK
+#define PREC 27
+#define MAXEXP 1024
+#define MINEXP -1077
+#endif
+#ifdef DEC
+#define PREC 29
+#define MAXEXP 128
+#define MINEXP -128
+#endif
+#ifdef IBMPC
+#define PREC 27
+#define MAXEXP 1024
+#define MINEXP -1077
+#endif
+#ifdef MIEEE
+#define PREC 27
+#define MAXEXP 1024
+#define MINEXP -1077
+#endif
+
+
+double cabs( z )
+register cmplx *z;
+{
+double x, y, b, re, im;
+int ex, ey, e;
+
+#ifdef INFINITIES
+/* Note, cabs(INFINITY,NAN) = INFINITY. */
+if( z->r == INFINITY || z->i == INFINITY
+ || z->r == -INFINITY || z->i == -INFINITY )
+ return( INFINITY );
+#endif
+
+#ifdef NANS
+if( isnan(z->r) )
+ return(z->r);
+if( isnan(z->i) )
+ return(z->i);
+#endif
+
+re = fabs( z->r );
+im = fabs( z->i );
+
+if( re == 0.0 )
+ return( im );
+if( im == 0.0 )
+ return( re );
+
+/* Get the exponents of the numbers */
+x = frexp( re, &ex );
+y = frexp( im, &ey );
+
+/* Check if one number is tiny compared to the other */
+e = ex - ey;
+if( e > PREC )
+ return( re );
+if( e < -PREC )
+ return( im );
+
+/* Find approximate exponent e of the geometric mean. */
+e = (ex + ey) >> 1;
+
+/* Rescale so mean is about 1 */
+x = ldexp( re, -e );
+y = ldexp( im, -e );
+
+/* Hypotenuse of the right triangle */
+b = sqrt( x * x + y * y );
+
+/* Compute the exponent of the answer. */
+y = frexp( b, &ey );
+ey = e + ey;
+
+/* Check it for overflow and underflow. */
+if( ey > MAXEXP )
+ {
+ mtherr( "cabs", OVERFLOW );
+ return( INFINITY );
+ }
+if( ey < MINEXP )
+ return(0.0);
+
+/* Undo the scaling */
+b = ldexp( b, e );
+return( b );
+}
+ /* csqrt()
+ *
+ * Complex square root
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * void csqrt();
+ * cmplx z, w;
+ *
+ * csqrt( &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 csqrt( z, w )
+cmplx *z, *w;
+{
+cmplx q, s;
+double x, y, r, t;
+
+x = z->r;
+y = z->i;
+
+if( y == 0.0 )
+ {
+ if( x < 0.0 )
+ {
+ w->r = 0.0;
+ w->i = sqrt(-x);
+ return;
+ }
+ else
+ {
+ w->r = sqrt(x);
+ w->i = 0.0;
+ return;
+ }
+ }
+
+
+if( x == 0.0 )
+ {
+ r = fabs(y);
+ r = sqrt(0.5*r);
+ if( y > 0 )
+ 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( (fabs(y) < 2.e-4 * fabs(x))
+ && (x > 0) )
+ {
+ t = 0.25*y*(y/x);
+ }
+else
+ {
+ r = cabs(z);
+ t = 0.5*(r - x);
+ }
+
+r = sqrt(t);
+q.i = r;
+q.r = y/(2.0*r);
+/* Heron iteration in complex arithmetic */
+cdiv( &q, z, &s );
+cadd( &q, &s, w );
+w->r *= 0.5;
+w->i *= 0.5;
+}
+
+
+double hypot( x, y )
+double x, y;
+{
+cmplx z;
+
+z.r = x;
+z.i = y;
+return( cabs(&z) );
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-11 17:53:26 +0000