Performs a complex division operation. The standard complex division performs a set of operations that is suseptible to both overflow and underflow. This method is more numerically stable while still being relatively fast to execute. @param a the real part of the first number @param b the imag
(double a, double b, double c, double d, double[] results)
| 196 | * index is the real, 2nd is the imaginary. |
| 197 | */ |
| 198 | public static void cDiv(double a, double b, double c, double d, double[] results) |
| 199 | { |
| 200 | /** |
| 201 | * Douglas M. Priest. Efficient scaling for complex division. ACM Trans. |
| 202 | * Math. Softw., 30(4):389–401, 2004 |
| 203 | */ |
| 204 | long aa, bb, cc, dd, ss; |
| 205 | double t; |
| 206 | int ha, hb, hc, hd, hz, hw, hs; |
| 207 | |
| 208 | /*extract high-order 32 bits to estimate |z| and |w| */ |
| 209 | aa = Double.doubleToRawLongBits(a); |
| 210 | bb = Double.doubleToRawLongBits(b); |
| 211 | |
| 212 | ha = (int) ((aa >> 32) & 0x7fffffff); |
| 213 | hb = (int) ((bb >> 32) & 0x7fffffff); |
| 214 | hz = (ha > hb)? ha : hb; |
| 215 | |
| 216 | cc = Double.doubleToRawLongBits(c); |
| 217 | dd = Double.doubleToRawLongBits(d); |
| 218 | |
| 219 | hc = (int) ((cc >> 32) & 0x7fffffff); |
| 220 | hd = (int) ((dd >> 32) & 0x7fffffff); |
| 221 | hw = (hc > hd)? hc : hd; |
| 222 | |
| 223 | /* compute the scale factor */ |
| 224 | if (hz < 0x07200000 && hw >= 0x32800000 && hw < 0x47100000) |
| 225 | { |
| 226 | /* |z| < 2^-909 and 2^-215 <= |w| < 2^114 */ |
| 227 | hs = (((0x47100000 - hw) >> 1) & 0xfff00000) + 0x3ff00000; |
| 228 | } |
| 229 | else |
| 230 | hs = (((hw >> 2) - hw) + 0x6fd7ffff) & 0xfff00000; |
| 231 | ss = ((long) hs) << 32; |
| 232 | |
| 233 | /* scale c and d, and compute the quotient */ |
| 234 | double ssd = Double.longBitsToDouble(ss); |
| 235 | c *= ssd; |
| 236 | d *= ssd; |
| 237 | t = 1.0 / (c * c + d * d); |
| 238 | c *= ssd; |
| 239 | d *= ssd; |
| 240 | results[0] = (a * c + b * d) * t; |
| 241 | results[1] = (b * c - a * d) * t; |
| 242 | } |
| 243 | |
| 244 | /** |
| 245 | * Alters this complex number as if a division by another complex number was performed. |
no outgoing calls