tcomplex.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tcomplex.c (1592B)
       ---
            1 #include <u.h>
            2 #include <libc.h>
            3 #include "map.h"
            4 
            5 /*complex divide, defensive against overflow from
            6  *        * and /, but not from + and -
            7  *        assumes underflow yields 0.0
            8  *        uses identities:
            9  *        (a + bi)/(c + di) = ((a + bd/c) + (b - ad/c)i)/(c + dd/c)
           10  *        (a + bi)/(c + di) = (b - ai)/(d - ci)
           11 */
           12 void
           13 cdiv(double a, double b, double c, double d, double *u, double *v)
           14 {
           15         double r,t;
           16         if(fabs(c)<fabs(d)) {
           17                 t = -c; c = d; d = t;
           18                 t = -a; a = b; b = t;
           19         }
           20         r = d/c;
           21         t = c + r*d;
           22         *u = (a + r*b)/t;
           23         *v = (b - r*a)/t;
           24 }
           25 
           26 void
           27 cmul(double c1, double c2, double d1, double d2, double *e1, double *e2)
           28 {
           29         *e1 = c1*d1 - c2*d2;
           30         *e2 = c1*d2 + c2*d1;
           31 }
           32 
           33 void
           34 csq(double c1, double c2, double *e1, double *e2)
           35 {
           36         *e1 = c1*c1 - c2*c2;
           37         *e2 = c1*c2*2;
           38 }
           39 
           40 /* complex square root
           41  *        assumes underflow yields 0.0
           42  *        uses these identities:
           43  *        sqrt(x+_iy) = sqrt(r(cos(t)+_isin(t))
           44  *                   = sqrt(r)(cos(t/2)+_isin(t/2))
           45  *        cos(t/2) = sin(t)/2sin(t/2) = sqrt((1+cos(t)/2)
           46  *        sin(t/2) = sin(t)/2cos(t/2) = sqrt((1-cos(t)/2)
           47 */
           48 void
           49 csqrt(double c1, double c2, double *e1, double *e2)
           50 {
           51         double r,s;
           52         double x,y;
           53         x = fabs(c1);
           54         y = fabs(c2);
           55         if(x>=y) {
           56                 if(x==0) {
           57                         *e1 = *e2 = 0;
           58                         return;
           59                 }
           60                 r = x;
           61                 s = y/x;
           62         } else {
           63                 r = y;
           64                 s = x/y;
           65         }
           66         r *= sqrt(1+ s*s);
           67         if(c1>0) {
           68                 *e1 = sqrt((r+c1)/2);
           69                 *e2 = c2/(2* *e1);
           70         } else {
           71                 *e2 = sqrt((r-c1)/2);
           72                 if(c2<0)
           73                         *e2 = -*e2;
           74                 *e1 = c2/(2* *e2);
           75         }
           76 }
           77 
           78 
           79 void cpow(double c1, double c2, double *d1, double *d2, double pwr)
           80 {
           81         double theta = pwr*atan2(c2,c1);
           82         double r = pow(hypot(c1,c2), pwr);
           83         *d1 = r*cos(theta);
           84         *d2 = r*sin(theta);
           85 }