1 #include "quadmath-imp.h"
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_);}
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
15 # define C128_MULT(x,y) ((x)*(y))
16 # define C128_DIV(x,y) ((x)/(y))
20 #define C128_MULT(x,y) mult_c128(x,y)
21 #define C128_DIV(x,y) div_c128(x,y)
23 static inline __complex128 mult_c128 (__complex128 x, __complex128 y)
25 __float128 r1 = REALPART(x), i1 = IMAGPART(x);
26 __float128 r2 = REALPART(y), i2 = IMAGPART(y);
28 COMPLEX_ASSIGN(res, r1*r2 - i1*i2, i2*r1 + i1*r2);
33 // Careful: the algorithm for the division sucks. A lot.
34 static inline __complex128 div_c128 (__complex128 x, __complex128 y)
36 __float128 n = hypotq (REALPART (y), IMAGPART (y));
37 __float128 r1 = REALPART(x), i1 = IMAGPART(x);
38 __float128 r2 = REALPART(y), i2 = IMAGPART(y);
40 COMPLEX_ASSIGN(res, r1*r2 + i1*i2, i1*r2 - i2*r1);
49 cabsq (__complex128 z)
51 return hypotq (REALPART (z), IMAGPART (z));
56 cexpq (__complex128 z)
63 COMPLEX_ASSIGN (v, cosq (b), sinq (b));
72 COMPLEX_ASSIGN (v, cosq (x), sinq (x));
78 cargq (__complex128 z)
80 return atan2q (IMAGPART (z), REALPART (z));
85 clogq (__complex128 z)
88 COMPLEX_ASSIGN (v, logq (cabsq (z)), cargq (z));
94 clog10q (__complex128 z)
97 COMPLEX_ASSIGN (v, log10q (cabsq (z)), cargq (z));
103 cpowq (__complex128 base, __complex128 power)
105 return cexpq (C128_MULT(power, clogq (base)));
110 csinq (__complex128 a)
112 __float128 r = REALPART (a), i = IMAGPART (a);
114 COMPLEX_ASSIGN (v, sinq (r) * coshq (i), cosq (r) * sinhq (i));
120 csinhq (__complex128 a)
122 __float128 r = REALPART (a), i = IMAGPART (a);
124 COMPLEX_ASSIGN (v, sinhq (r) * cosq (i), coshq (r) * sinq (i));
130 ccosq (__complex128 a)
132 __float128 r = REALPART (a), i = IMAGPART (a);
134 COMPLEX_ASSIGN (v, cosq (r) * coshq (i), - (sinq (r) * sinhq (i)));
140 ccoshq (__complex128 a)
142 __float128 r = REALPART (a), i = IMAGPART (a);
144 COMPLEX_ASSIGN (v, coshq (r) * cosq (i), sinhq (r) * sinq (i));
150 ctanq (__complex128 a)
152 __float128 rt = tanq (REALPART (a)), it = tanhq (IMAGPART (a));
154 COMPLEX_ASSIGN (n, rt, it);
155 COMPLEX_ASSIGN (d, 1, - (rt * it));
156 return C128_DIV(n,d);
161 ctanhq (__complex128 a)
163 __float128 rt = tanhq (REALPART (a)), it = tanq (IMAGPART (a));
165 COMPLEX_ASSIGN (n, rt, it);
166 COMPLEX_ASSIGN (d, 1, rt * it);
167 return C128_DIV(n,d);
171 /* Square root algorithm from glibc. */
173 csqrtq (__complex128 z)
175 __float128 re = REALPART(z), im = IMAGPART(z);
182 COMPLEX_ASSIGN (v, 0, copysignq (sqrtq (-re), im));
186 COMPLEX_ASSIGN (v, fabsq (sqrtq (re)), copysignq (0, im));
191 __float128 r = sqrtq (0.5 * fabsq (im));
192 COMPLEX_ASSIGN (v, r, copysignq (r, im));
196 __float128 d = hypotq (re, im);
199 /* Use the identity 2 Re res Im res = Im x
200 to avoid cancellation error in d +/- Re x. */
202 r = sqrtq (0.5 * d + 0.5 * re), s = (0.5 * im) / r;
204 s = sqrtq (0.5 * d - 0.5 * re), r = fabsq ((0.5 * im) / s);
206 COMPLEX_ASSIGN (v, r, copysignq (s, im));