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 }