libquadmath/math/complex.c

   1 #include "quadmath-imp.h"
   2
   3
   4 #define REALPART(z) (__real__(z))
   5 #define IMAGPART(z) (__imag__(z))
   6 #define COMPLEX_ASSIGN(z_, r_, i_) {__real__(z_) = (r_); __imag__(z_) = (i_);}
   7
   8
   9 // Horrible... GCC doesn't know how to multiply or divide these
  10 // __complex128 things. We have to do it on our own.
  11 // Protect it around macros so, some day, we can switch it on
  12
  13 #if 0
  14
  15 # define C128_MULT(x,y) ((x)*(y))
  16 # define C128_DIV(x,y) ((x)/(y))
  17
  18 #else
  19
  20 #define C128_MULT(x,y) mult_c128(x,y)
  21 #define C128_DIV(x,y) div_c128(x,y)
  22
  23 static inline __complex128 mult_c128 (__complex128 x, __complex128 y)
  24 {
  25   __float128 r1 = REALPART(x), i1 = IMAGPART(x);
  26   __float128 r2 = REALPART(y), i2 = IMAGPART(y);
  27   __complex128 res;
  28   COMPLEX_ASSIGN(res, r1*r2 - i1*i2, i2*r1 + i1*r2);
  29   return res;
  30 }
  31
  32
  33 // Careful: the algorithm for the division sucks. A lot.
  34 static inline __complex128 div_c128 (__complex128 x, __complex128 y)
  35 {
  36   __float128 n = hypotq (REALPART (y), IMAGPART (y));
  37   __float128 r1 = REALPART(x), i1 = IMAGPART(x);
  38   __float128 r2 = REALPART(y), i2 = IMAGPART(y);
  39   __complex128 res;
  40   COMPLEX_ASSIGN(res, r1*r2 + i1*i2, i1*r2 - i2*r1);
  41   return res / n;
  42 }
  43
  44 #endif
  45
  46
  47
  48 __float128
  49 cabsq (__complex128 z)
  50 {
  51   return hypotq (REALPART (z), IMAGPART (z));
  52 }
  53
  54
  55 __complex128
  56 cexpq (__complex128 z)
  57 {
  58   __float128 a, b;
  59   __complex128 v;
  60
  61   a = REALPART (z);
  62   b = IMAGPART (z);
  63   COMPLEX_ASSIGN (v, cosq (b), sinq (b));
  64   return expq (a) * v;
  65 }
  66
  67
  68 __complex128
  69 cexpiq (__float128 x)
  70 {
  71   __complex128 v;
  72   COMPLEX_ASSIGN (v, cosq (x), sinq (x));
  73   return v;
  74 }
  75
  76
  77 __float128
  78 cargq (__complex128 z)
  79 {
  80   return atan2q (IMAGPART (z), REALPART (z));
  81 }
  82
  83
  84 __complex128
  85 clogq (__complex128 z)
  86 {
  87   __complex128 v;
  88   COMPLEX_ASSIGN (v, logq (cabsq (z)), cargq (z));
  89   return v;
  90 }
  91
  92
  93 __complex128
  94 clog10q (__complex128 z)
  95 {
  96   __complex128 v;
  97   COMPLEX_ASSIGN (v, log10q (cabsq (z)), cargq (z));
  98   return v;
  99 }
 100
 101
 102 __complex128
 103 cpowq (__complex128 base, __complex128 power)
 104 {
 105   return cexpq (C128_MULT(power, clogq (base)));
 106 }
 107
 108
 109 __complex128
 110 csinq (__complex128 a)
 111 {
 112   __float128 r = REALPART (a), i = IMAGPART (a);
 113   __complex128 v;
 114   COMPLEX_ASSIGN (v, sinq (r) * coshq (i), cosq (r) * sinhq (i));
 115   return v;
 116 }
 117
 118
 119 __complex128
 120 csinhq (__complex128 a)
 121 {
 122   __float128 r = REALPART (a), i = IMAGPART (a);
 123   __complex128 v;
 124   COMPLEX_ASSIGN (v, sinhq (r) * cosq (i), coshq (r) * sinq (i));
 125   return v;
 126 }
 127
 128
 129 __complex128
 130 ccosq (__complex128 a)
 131 {
 132   __float128 r = REALPART (a), i = IMAGPART (a);
 133   __complex128 v;
 134   COMPLEX_ASSIGN (v, cosq (r) * coshq (i), - (sinq (r) * sinhq (i)));
 135   return v;
 136 }
 137
 138
 139 __complex128
 140 ccoshq (__complex128 a)
 141 {
 142   __float128 r = REALPART (a), i = IMAGPART (a);
 143   __complex128 v;
 144   COMPLEX_ASSIGN (v, coshq (r) * cosq (i),  sinhq (r) * sinq (i));
 145   return v;
 146 }
 147
 148
 149 __complex128
 150 ctanq (__complex128 a)
 151 {
 152   __float128 rt = tanq (REALPART (a)), it = tanhq (IMAGPART (a));
 153   __complex128 n, d;
 154   COMPLEX_ASSIGN (n, rt, it);
 155   COMPLEX_ASSIGN (d, 1, - (rt * it));
 156   return C128_DIV(n,d);
 157 }
 158
 159
 160 __complex128
 161 ctanhq (__complex128 a)
 162 {
 163   __float128 rt = tanhq (REALPART (a)), it = tanq (IMAGPART (a));
 164   __complex128 n, d;
 165   COMPLEX_ASSIGN (n, rt, it);
 166   COMPLEX_ASSIGN (d, 1, rt * it);
 167   return C128_DIV(n,d);
 168 }
 169
 170
 171 /* Square root algorithm from glibc.  */
 172 __complex128
 173 csqrtq (__complex128 z)
 174 {
 175   __float128 re = REALPART(z), im = IMAGPART(z);
 176   __complex128 v;
 177
 178   if (im == 0)
 179   {
 180     if (re < 0)
 181     {
 182       COMPLEX_ASSIGN (v, 0, copysignq (sqrtq (-re), im));
 183     }
 184     else
 185     {
 186       COMPLEX_ASSIGN (v, fabsq (sqrtq (re)), copysignq (0, im));
 187     }
 188   }
 189   else if (re == 0)
 190   {
 191     __float128 r = sqrtq (0.5 * fabsq (im));
 192     COMPLEX_ASSIGN (v, r, copysignq (r, im));
 193   }
 194   else
 195   {
 196     __float128 d = hypotq (re, im);
 197     __float128 r, s;
 198
 199     /* Use the identity   2  Re res  Im res = Im x
 200         to avoid cancellation error in  d +/- Re x.  */
 201     if (re > 0)
 202       r = sqrtq (0.5 * d + 0.5 * re), s = (0.5 * im) / r;
 203     else
 204       s = sqrtq (0.5 * d - 0.5 * re), r = fabsq ((0.5 * im) / s);
 205
 206     COMPLEX_ASSIGN (v, r, copysignq (s, im));
 207   }
 208   return v;
 209 }
 210