dmd.backend.bcomplex source code

1 /**
2  * Compiler implementation of the
3  * $(LINK2 http://www.dlang.org, D programming language).
4  *
5  * Copyright:   public domain
6  * License:     public domain
7  * Source:      $(DMDSRC backend/_bcomplex.d)
8  */
9 
10 module dmd.backend.bcomplex;
11 
12 public import dmd.root.longdouble : targ_ldouble = longdouble;
13 import core.stdc.math : fabs, fabsl, sqrt;
14 version(CRuntime_Microsoft)
15     private import dmd.root.longdouble : fabsl, sqrt; // needed if longdouble is longdouble_soft
16 
17 extern (C++):
18 @nogc:
19 @safe:
20 nothrow:
21 
22 // Roll our own for reliable bootstrapping
23 
24 
25 struct Complex_f
26 {
27 nothrow:
28     float re, im;
29 
30     static Complex_f div(ref Complex_f x, ref Complex_f y)
31     {
32         if (fabs(y.re) < fabs(y.im))
33         {
34             const r = y.re / y.im;
35             const den = y.im + r * y.re;
36             return Complex_f(cast(float)((x.re * r + x.im) / den),
37                              cast(float)((x.im * r - x.re) / den));
38         }
39         else
40         {
41             const r = y.im / y.re;
42             const den = y.re + r * y.im;
43             return Complex_f(cast(float)((x.re + r * x.im) / den),
44                              cast(float)((x.im - r * x.re) / den));
45         }
46     }
47 
48     static Complex_f mul(ref Complex_f x, ref Complex_f y) pure
49     {
50         return Complex_f(x.re * y.re - x.im * y.im,
51                          x.im * y.re + x.re * y.im);
52     }
53 
54     static targ_ldouble abs(ref Complex_f z)
55     {
56         const targ_ldouble x = fabs(z.re);
57         const targ_ldouble y = fabs(z.im);
58         if (x == 0)
59             return y;
60         else if (y == 0)
61             return x;
62         else if (x > y)
63         {
64             const targ_ldouble temp = y / x;
65             return x * sqrt(1 + temp * temp);
66         }
67         else
68         {
69             const targ_ldouble temp = x / y;
70             return y * sqrt(1 + temp * temp);
71         }
72     }
73 
74     static Complex_f sqrtc(ref Complex_f z)
75     {
76         if (z.re == 0 && z.im == 0)
77         {
78             return Complex_f(0, 0);
79         }
80         else
81         {
82             const targ_ldouble x = fabs(z.re);
83             const targ_ldouble y = fabs(z.im);
84             targ_ldouble r, w;
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 
96             if (z.re >= 0)
97             {
98                 return Complex_f(cast(float)w, (z.im / cast(float)(w + w)));
99             }
100             else
101             {
102                 const cim = (z.im >= 0) ? w : -w;
103                 return Complex_f((z.im / cast(float)(cim + cim)), cast(float)cim);
104             }
105         }
106     }
107 }
108 
109 struct Complex_d
110 {
111 nothrow:
112     double re, im;
113 
114     static Complex_d div(ref Complex_d x, ref Complex_d y)
115     {
116         if (fabs(y.re) < fabs(y.im))
117         {
118             const targ_ldouble r = y.re / y.im;
119             const targ_ldouble den = y.im + r * y.re;
120             return Complex_d(cast(double)((x.re * r + x.im) / den),
121                              cast(double)((x.im * r - x.re) / den));
122         }
123         else
124         {
125             const targ_ldouble r = y.im / y.re;
126             const targ_ldouble den = y.re + r * y.im;
127             return Complex_d(cast(double)((x.re + r * x.im) / den),
128                              cast(double)((x.im - r * x.re) / den));
129         }
130     }
131 
132     static Complex_d mul(ref Complex_d x, ref Complex_d y) pure
133     {
134         return Complex_d(x.re * y.re - x.im * y.im,
135                          x.im * y.re + x.re * y.im);
136     }
137 
138     static targ_ldouble abs(ref Complex_d z)
139     {
140         const targ_ldouble x = fabs(z.re);
141         const targ_ldouble y = fabs(z.im);
142         if (x == 0)
143             return y;
144         else if (y == 0)
145             return x;
146         else if (x > y)
147         {
148             const targ_ldouble temp = y / x;
149             return x * sqrt(1 + temp * temp);
150         }
151         else
152         {
153             const targ_ldouble temp = x / y;
154             return y * sqrt(1 + temp * temp);
155         }
156     }
157 
158     static Complex_d sqrtc(ref Complex_d z)
159     {
160         if (z.re == 0 && z.im == 0)
161         {
162             return Complex_d(0, 0);
163         }
164         else
165         {
166             const targ_ldouble x = fabs(z.re);
167             const targ_ldouble y = fabs(z.im);
168             targ_ldouble r, w;
169             if (x >= y)
170             {
171                 r = y / x;
172                 w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
173             }
174             else
175             {
176                 r = x / y;
177                 w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
178             }
179 
180             if (z.re >= 0)
181             {
182                 return Complex_d(cast(double)w, (z.im / cast(double)(w + w)));
183             }
184             else
185             {
186                 const cim = (z.im >= 0) ? w : -w;
187                 return Complex_d((z.im / cast(double)(cim + cim)), cast(double)cim);
188             }
189         }
190     }
191 }
192 
193 
194 struct Complex_ld
195 {
196 nothrow:
197     targ_ldouble re, im;
198 
199     static Complex_ld div(ref Complex_ld x, ref Complex_ld y)
200     {
201         if (fabsl(y.re) < fabsl(y.im))
202         {
203             const targ_ldouble r = y.re / y.im;
204             const targ_ldouble den = y.im + r * y.re;
205             return Complex_ld((x.re * r + x.im) / den,
206                               (x.im * r - x.re) / den);
207         }
208         else
209         {
210             const targ_ldouble r = y.im / y.re;
211             const targ_ldouble den = y.re + r * y.im;
212             return Complex_ld((x.re + r * x.im) / den,
213                               (x.im - r * x.re) / den);
214         }
215     }
216 
217     static Complex_ld mul(ref Complex_ld x, ref Complex_ld y) pure
218     {
219         return Complex_ld(x.re * y.re - x.im * y.im,
220                           x.im * y.re + x.re * y.im);
221     }
222 
223     static targ_ldouble abs(ref Complex_ld z)
224     {
225         const targ_ldouble x = fabsl(z.re);
226         const targ_ldouble y = fabsl(z.im);
227         if (x == 0)
228             return y;
229         else if (y == 0)
230             return x;
231         else if (x > y)
232         {
233             const targ_ldouble temp = y / x;
234             return x * sqrt(1 + temp * temp);
235         }
236         else
237         {
238             const targ_ldouble temp = x / y;
239             return y * sqrt(1 + temp * temp);
240         }
241     }
242 
243     static Complex_ld sqrtc(ref Complex_ld z)
244     {
245         if (z.re == 0 && z.im == 0)
246         {
247             return Complex_ld(targ_ldouble(0), targ_ldouble(0));
248         }
249         else
250         {
251             const targ_ldouble x = fabsl(z.re);
252             const targ_ldouble y = fabsl(z.im);
253             targ_ldouble r, w;
254             if (x >= y)
255             {
256                 r = y / x;
257                 w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
258             }
259             else
260             {
261                 r = x / y;
262                 w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
263             }
264 
265             if (z.re >= 0)
266             {
267                 return Complex_ld(w, z.im / (w + w));
268             }
269             else
270             {
271                 const cim = (z.im >= 0) ? w : -w;
272                 return Complex_ld(z.im / (cim + cim), cim);
273             }
274         }
275     }
276 }