Mercurial > hg > CbC > CbC_gcc
comparison libquadmath/math/complex.c @ 68:561a7518be6b
update gcc-4.6
author | Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> |
---|---|
date | Sun, 21 Aug 2011 07:07:55 +0900 |
parents | |
children | 04ced10e8804 |
comparison
equal
deleted
inserted
replaced
67:f6334be47118 | 68:561a7518be6b |
---|---|
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 |