libm/float/cmplxf.c

   1 /*                                                      cmplxf.c
   2  *
   3  *      Complex number arithmetic
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * typedef struct {
  10  *      float r;     real part
  11  *      float i;     imaginary part
  12  *     }cmplxf;
  13  *
  14  * cmplxf *a, *b, *c;
  15  *
  16  * caddf( a, b, c );     c = b + a
  17  * csubf( a, b, c );     c = b - a
  18  * cmulf( a, b, c );     c = b * a
  19  * cdivf( a, b, c );     c = b / a
  20  * cnegf( c );           c = -c
  21  * cmovf( b, c );        c = b
  22  *
  23  *
  24  *
  25  * DESCRIPTION:
  26  *
  27  * Addition:
  28  *    c.r  =  b.r + a.r
  29  *    c.i  =  b.i + a.i
  30  *
  31  * Subtraction:
  32  *    c.r  =  b.r - a.r
  33  *    c.i  =  b.i - a.i
  34  *
  35  * Multiplication:
  36  *    c.r  =  b.r * a.r  -  b.i * a.i
  37  *    c.i  =  b.r * a.i  +  b.i * a.r
  38  *
  39  * Division:
  40  *    d    =  a.r * a.r  +  a.i * a.i
  41  *    c.r  = (b.r * a.r  + b.i * a.i)/d
  42  *    c.i  = (b.i * a.r  -  b.r * a.i)/d
  43 \f * ACCURACY:
  44  *
  45  * In DEC arithmetic, the test (1/z) * z = 1 had peak relative
  46  * error 3.1e-17, rms 1.2e-17.  The test (y/z) * (z/y) = 1 had
  47  * peak relative error 8.3e-17, rms 2.1e-17.
  48  *
  49  * Tests in the rectangle {-10,+10}:
  50  *                      Relative error:
  51  * arithmetic   function  # trials      peak         rms
  52  *    IEEE       cadd       30000       5.9e-8      2.6e-8
  53  *    IEEE       csub       30000       6.0e-8      2.6e-8
  54  *    IEEE       cmul       30000       1.1e-7      3.7e-8
  55  *    IEEE       cdiv       30000       2.1e-7      5.7e-8
  56  */
  57 \f/*                             cmplx.c
  58  * complex number arithmetic
  59  */
  60
  61
  62 /*
  63 Cephes Math Library Release 2.1:  December, 1988
  64 Copyright 1984, 1987, 1988 by Stephen L. Moshier
  65 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  66 */
  67
  68
  69 #include <math.h>
  70 extern float MAXNUMF, MACHEPF, PIF, PIO2F;
  71 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  72 #ifdef ANSIC
  73 float sqrtf(float), frexpf(float, int *);
  74 float ldexpf(float, int);
  75 float cabsf(cmplxf *), atan2f(float, float), cosf(float), sinf(float);
  76 #else
  77 float sqrtf(), frexpf(), ldexpf();
  78 float cabsf(), atan2f(), cosf(), sinf();
  79 #endif
  80 /*
  81 typedef struct
  82         {
  83         float r;
  84         float i;
  85         }cmplxf;
  86 */
  87 cmplxf czerof = {0.0, 0.0};
  88 extern cmplxf czerof;
  89 cmplxf conef = {1.0, 0.0};
  90 extern cmplxf conef;
  91
  92 /*      c = b + a       */
  93
  94 void caddf( a, b, c )
  95 register cmplxf *a, *b;
  96 cmplxf *c;
  97 {
  98
  99 c->r = b->r + a->r;
 100 c->i = b->i + a->i;
 101 }
 102
 103
 104 /*      c = b - a       */
 105
 106 void csubf( a, b, c )
 107 register cmplxf *a, *b;
 108 cmplxf *c;
 109 {
 110
 111 c->r = b->r - a->r;
 112 c->i = b->i - a->i;
 113 }
 114
 115 /*      c = b * a */
 116
 117 void cmulf( a, b, c )
 118 register cmplxf *a, *b;
 119 cmplxf *c;
 120 {
 121 register float y;
 122
 123 y    = b->r * a->r  -  b->i * a->i;
 124 c->i = b->r * a->i  +  b->i * a->r;
 125 c->r = y;
 126 }
 127
 128
 129
 130 /*      c = b / a */
 131
 132 void cdivf( a, b, c )
 133 register cmplxf *a, *b;
 134 cmplxf *c;
 135 {
 136 float y, p, q, w;
 137
 138
 139 y = a->r * a->r  +  a->i * a->i;
 140 p = b->r * a->r  +  b->i * a->i;
 141 q = b->i * a->r  -  b->r * a->i;
 142
 143 if( y < 1.0f )
 144         {
 145         w = MAXNUMF * y;
 146         if( (fabsf(p) > w) || (fabsf(q) > w) || (y == 0.0f) )
 147                 {
 148                 c->r = MAXNUMF;
 149                 c->i = MAXNUMF;
 150                 mtherr( "cdivf", OVERFLOW );
 151                 return;
 152                 }
 153         }
 154 c->r = p/y;
 155 c->i = q/y;
 156 }
 157
 158
 159 /*      b = a   */
 160
 161 void cmovf( a, b )
 162 register short *a, *b;
 163 {
 164 int i;
 165
 166
 167 i = 8;
 168 do
 169         *b++ = *a++;
 170 while( --i );
 171 }
 172
 173
 174 void cnegf( a )
 175 register cmplxf *a;
 176 {
 177
 178 a->r = -a->r;
 179 a->i = -a->i;
 180 }
 181
 182 /*                                                      cabsf()
 183  *
 184  *      Complex absolute value
 185  *
 186  *
 187  *
 188  * SYNOPSIS:
 189  *
 190  * float cabsf();
 191  * cmplxf z;
 192  * float a;
 193  *
 194  * a = cabsf( &z );
 195  *
 196  *
 197  *
 198  * DESCRIPTION:
 199  *
 200  *
 201  * If z = x + iy
 202  *
 203  * then
 204  *
 205  *       a = sqrt( x**2 + y**2 ).
 206  *
 207  * Overflow and underflow are avoided by testing the magnitudes
 208  * of x and y before squaring.  If either is outside half of
 209  * the floating point full scale range, both are rescaled.
 210  *
 211  *
 212  * ACCURACY:
 213  *
 214  *                      Relative error:
 215  * arithmetic   domain     # trials      peak         rms
 216  *    IEEE      -10,+10     30000       1.2e-7      3.4e-8
 217  */
 218 \f
 219
 220 /*
 221 Cephes Math Library Release 2.1:  January, 1989
 222 Copyright 1984, 1987, 1989 by Stephen L. Moshier
 223 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 224 */
 225
 226
 227 /*
 228 typedef struct
 229         {
 230         float r;
 231         float i;
 232         }cmplxf;
 233 */
 234 /* square root of max and min numbers */
 235 #define SMAX  1.3043817825332782216E+19
 236 #define SMIN  7.6664670834168704053E-20
 237 #define PREC 12
 238 #define MAXEXPF 128
 239
 240
 241 #define SMAXT (2.0f * SMAX)
 242 #define SMINT (0.5f * SMIN)
 243
 244 float cabsf( z )
 245 register cmplxf *z;
 246 {
 247 float x, y, b, re, im;
 248 int ex, ey, e;
 249
 250 re = fabsf( z->r );
 251 im = fabsf( z->i );
 252
 253 if( re == 0.0f )
 254         {
 255         return( im );
 256         }
 257 if( im == 0.0f )
 258         {
 259         return( re );
 260         }
 261
 262 /* Get the exponents of the numbers */
 263 x = frexpf( re, &ex );
 264 y = frexpf( im, &ey );
 265
 266 /* Check if one number is tiny compared to the other */
 267 e = ex - ey;
 268 if( e > PREC )
 269         return( re );
 270 if( e < -PREC )
 271         return( im );
 272
 273 /* Find approximate exponent e of the geometric mean. */
 274 e = (ex + ey) >> 1;
 275
 276 /* Rescale so mean is about 1 */
 277 x = ldexpf( re, -e );
 278 y = ldexpf( im, -e );
 279
 280 /* Hypotenuse of the right triangle */
 281 b = sqrtf( x * x  +  y * y );
 282
 283 /* Compute the exponent of the answer. */
 284 y = frexpf( b, &ey );
 285 ey = e + ey;
 286
 287 /* Check it for overflow and underflow. */
 288 if( ey > MAXEXPF )
 289         {
 290         mtherr( "cabsf", OVERFLOW );
 291         return( MAXNUMF );
 292         }
 293 if( ey < -MAXEXPF )
 294         return(0.0f);
 295
 296 /* Undo the scaling */
 297 b = ldexpf( b, e );
 298 return( b );
 299 }
 300 \f/*                                                     csqrtf()
 301  *
 302  *      Complex square root
 303  *
 304  *
 305  *
 306  * SYNOPSIS:
 307  *
 308  * void csqrtf();
 309  * cmplxf z, w;
 310  *
 311  * csqrtf( &z, &w );
 312  *
 313  *
 314  *
 315  * DESCRIPTION:
 316  *
 317  *
 318  * If z = x + iy,  r = |z|, then
 319  *
 320  *                       1/2
 321  * Im w  =  [ (r - x)/2 ]   ,
 322  *
 323  * Re w  =  y / 2 Im w.
 324  *
 325  *
 326  * Note that -w is also a square root of z.  The solution
 327  * reported is always in the upper half plane.
 328  *
 329  * Because of the potential for cancellation error in r - x,
 330  * the result is sharpened by doing a Heron iteration
 331  * (see sqrt.c) in complex arithmetic.
 332  *
 333  *
 334  *
 335  * ACCURACY:
 336  *
 337  *                      Relative error:
 338  * arithmetic   domain     # trials      peak         rms
 339  *    IEEE      -10,+10    100000       1.8e-7       4.2e-8
 340  *
 341  */
 342 \f
 343
 344 void csqrtf( z, w )
 345 cmplxf *z, *w;
 346 {
 347 cmplxf q, s;
 348 float x, y, r, t;
 349
 350 x = z->r;
 351 y = z->i;
 352
 353 if( y == 0.0f )
 354         {
 355         if( x < 0.0f )
 356                 {
 357                 w->r = 0.0f;
 358                 w->i = sqrtf(-x);
 359                 return;
 360                 }
 361         else
 362                 {
 363                 w->r = sqrtf(x);
 364                 w->i = 0.0f;
 365                 return;
 366                 }
 367         }
 368
 369 if( x == 0.0f )
 370         {
 371         r = fabsf(y);
 372         r = sqrtf(0.5f*r);
 373         if( y > 0 )
 374                 w->r = r;
 375         else
 376                 w->r = -r;
 377         w->i = r;
 378         return;
 379         }
 380
 381 /* Approximate  sqrt(x^2+y^2) - x  =  y^2/2x - y^4/24x^3 + ... .
 382  * The relative error in the first term is approximately y^2/12x^2 .
 383  */
 384 if( (fabsf(y) < fabsf(0.015f*x))
 385    && (x > 0) )
 386         {
 387         t = 0.25f*y*(y/x);
 388         }
 389 else
 390         {
 391         r = cabsf(z);
 392         t = 0.5f*(r - x);
 393         }
 394
 395 r = sqrtf(t);
 396 q.i = r;
 397 q.r = 0.5f*y/r;
 398
 399 /* Heron iteration in complex arithmetic:
 400  * q = (q + z/q)/2
 401  */
 402 cdivf( &q, z, &s );
 403 caddf( &q, &s, w );
 404 w->r *= 0.5f;
 405 w->i *= 0.5f;
 406 }
 407