dmd.backend.bcomplex source code

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