src/compiler/dmd/complex.c

   1 /*
   2  *  Placed into the public domain.
   3  *  Written by Walter Bright
   4  *  www.digitalmars.com
   5  */
   6
   7
   8 #include <math.h>
   9
  10 typedef struct Complex
  11 {
  12     long double re;
  13     long double im;
  14 } Complex;
  15
  16 Complex _complex_div(Complex x, Complex y)
  17 {
  18     Complex q;
  19     long double r;
  20     long double den;
  21
  22     if (fabs(y.re) < fabs(y.im))
  23     {
  24         r = y.re / y.im;
  25         den = y.im + r * y.re;
  26         q.re = (x.re * r + x.im) / den;
  27         q.im = (x.im * r - x.re) / den;
  28     }
  29     else
  30     {
  31         r = y.im / y.re;
  32         den = y.re + r * y.im;
  33         q.re = (x.re + r * x.im) / den;
  34         q.im = (x.im - r * x.re) / den;
  35     }
  36     return q;
  37 }
  38
  39 Complex _complex_mul(Complex x, Complex y)
  40 {
  41     Complex p;
  42
  43     p.re = x.re * y.re - x.im * y.im;
  44     p.im = x.im * y.re + x.re * y.im;
  45     return p;
  46 }
  47
  48 long double _complex_abs(Complex z)
  49 {
  50     long double x,y,ans,temp;
  51
  52     x = fabs(z.re);
  53     y = fabs(z.im);
  54     if (x == 0)
  55         ans = y;
  56     else if (y == 0)
  57         ans = x;
  58     else if (x > y)
  59     {
  60         temp = y / x;
  61         ans = x * sqrt(1 + temp * temp);
  62     }
  63     else
  64     {
  65         temp = x / y;
  66         ans = y * sqrt(1 + temp * temp);
  67     }
  68     return ans;
  69 }
  70
  71 Complex _complex_sqrt(Complex z)
  72 {
  73     Complex c;
  74     long double x,y,w,r;
  75
  76     if (z.re == 0 && z.im == 0)
  77     {
  78         c.re = 0;
  79         c.im = 0;
  80     }
  81     else
  82     {
  83         x = fabs(z.re);
  84         y = fabs(z.im);
  85         if (x >= y)
  86         {
  87             r = y / x;
  88             w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
  89         }
  90         else
  91         {
  92             r = x / y;
  93             w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
  94         }
  95         if (z.re >= 0)
  96         {
  97             c.re = w;
  98             c.im = z.im / (w + w);
  99         }
 100         else
 101         {
 102             c.im = (z.im >= 0) ? w : -w;
 103             c.re = z.im / (c.im + c.im);
 104         }
 105     }
 106     return c;
 107 }