Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/complex.c @ 111:04ced10e8804

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc 7
author kono
date Fri, 27 Oct 2017 22:46:09 +0900
parents 561a7518be6b
children
comparison
equal deleted inserted replaced
68:561a7518be6b 111:04ced10e8804
1 /* GCC Quad-Precision Math Library
2 Copyright (C) 2010, 2011 Free Software Foundation, Inc.
3 Written by Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
4
5 This file is part of the libquadmath library.
6 Libquadmath is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Library General Public
8 License as published by the Free Software Foundation; either
9 version 2 of the License, or (at your option) any later version.
10
11 Libquadmath is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Library General Public License for more details.
15
16 You should have received a copy of the GNU Library General Public
17 License along with libquadmath; see the file COPYING.LIB. If
18 not, write to the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 Boston, MA 02110-1301, USA. */
20
1 #include "quadmath-imp.h" 21 #include "quadmath-imp.h"
22
23 #ifdef HAVE_FENV_H
24 # include <fenv.h>
25 #endif
2 26
3 27
4 #define REALPART(z) (__real__(z)) 28 #define REALPART(z) (__real__(z))
5 #define IMAGPART(z) (__imag__(z)) 29 #define IMAGPART(z) (__imag__(z))
6 #define COMPLEX_ASSIGN(z_, r_, i_) {__real__(z_) = (r_); __imag__(z_) = (i_);} 30 #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 31
47 32
48 __float128 33 __float128
49 cabsq (__complex128 z) 34 cabsq (__complex128 z)
50 { 35 {
51 return hypotq (REALPART (z), IMAGPART (z)); 36 return hypotq (REALPART (z), IMAGPART (z));
52 } 37 }
53 38
54 39
55 __complex128 40 __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) 41 cexpiq (__float128 x)
70 { 42 {
43 __float128 sinix, cosix;
71 __complex128 v; 44 __complex128 v;
72 COMPLEX_ASSIGN (v, cosq (x), sinq (x)); 45 sincosq (x, &sinix, &cosix);
46 COMPLEX_ASSIGN (v, cosix, sinix);
73 return v; 47 return v;
74 } 48 }
75 49
76 50
77 __float128 51 __float128
80 return atan2q (IMAGPART (z), REALPART (z)); 54 return atan2q (IMAGPART (z), REALPART (z));
81 } 55 }
82 56
83 57
84 __complex128 58 __complex128
85 clogq (__complex128 z) 59 cpowq (__complex128 base, __complex128 power)
86 { 60 {
87 __complex128 v; 61 return cexpq (power * clogq (base));
88 COMPLEX_ASSIGN (v, logq (cabsq (z)), cargq (z));
89 return v;
90 } 62 }
91 63
92 64
93 __complex128 65 __complex128
94 clog10q (__complex128 z) 66 ccosq (__complex128 x)
95 { 67 {
96 __complex128 v; 68 __complex128 y;
97 COMPLEX_ASSIGN (v, log10q (cabsq (z)), cargq (z)); 69
98 return v; 70 COMPLEX_ASSIGN (y, -IMAGPART (x), REALPART (x));
71 return ccoshq (y);
99 } 72 }
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