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