Skip to content

Commit a6bb110

Browse files
jeff5jeff5
committed
Re-work cmath.sqrt.
Addresses certain test failures in test_cmath.CMathTests.test_specific_values, part of issue #2237.
1 parent 9f0c136 commit a6bb110

File tree

1 file changed

+104
-20
lines changed

1 file changed

+104
-20
lines changed

‎src/org/python/modules/cmath.java‎

Lines changed: 104 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,9 @@ public class cmath {
1717
private static final PyComplex i = new PyComplex(0.0, 1.0);
1818
private static final PyComplex half_i = new PyComplex(0.0, 0.5);
1919

20+
/** 2<sup>-&#189;</sup> (Ref: Abramowitz &amp; Stegun [1972], p2). */
21+
private static final double ROOT_HALF = 0.70710678118654752440;
22+
2023
private static PyComplex c_prodi(PyComplex x) {
2124
return (PyComplex)x.__mul__(i);
2225
}
@@ -256,32 +259,113 @@ public static PyComplex sinh(PyObject in) {
256259
* math.cosh(x.real));
257260
}
258261

259-
public static PyComplex sqrt(PyObject in) {
260-
PyComplex x = complexFromPyObject(in);
262+
/**
263+
* Calculate <i>z = x+iy</i>, such that <i>z<sup>2</sup> = w</i>. In taking the square roots to
264+
* get <i>x</i> and <i>y</i>, we choose to have <i>x&ge;0</i> always, and <i>y</i> the same sign
265+
* as <i>v</i>.
266+
*
267+
* @param w to square-root
268+
* @return <i>w<sup>&#189;</sup></i>
269+
*/
270+
public static PyComplex sqrt(PyObject w) {
271+
/*
272+
* All the difficult parts are written for the first quadrant only (+,+), then the true sign
273+
* of the parts of w are factored in at the end, by flipping the result around the
274+
* diagonals.
275+
*/
276+
PyComplex ww = complexFromPyObject(w);
277+
double u = Math.abs(ww.real), v = Math.abs(ww.imag), x, y;
278+
279+
if (Double.isInfinite(u)) {
280+
// Special cases: u = inf
281+
x = Double.POSITIVE_INFINITY;
282+
y = (Double.isNaN(v) || Double.isInfinite(v)) ? v : 0.;
283+
284+
} else if (Double.isInfinite(v)) {
285+
// Special cases: v = inf, u != inf
286+
x = y = Double.POSITIVE_INFINITY;
287+
288+
} else if (Double.isNaN(u)) {
289+
// In the remaining cases, u == nan infects all.
290+
x = y = u;
261291

262-
if (Double.isInfinite(x.real) && Double.isNaN(x.imag)) {
263-
if (x.real == Double.NEGATIVE_INFINITY) {
264-
return new PyComplex(Double.NaN, Double.POSITIVE_INFINITY);
265-
} else {
266-
return new PyComplex(Double.POSITIVE_INFINITY, Double.NaN);
267-
}
268-
}
292+
} else {
269293

270-
if (x.imag == Double.POSITIVE_INFINITY) {
271-
return new PyComplex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
272-
} else if (x.imag == Double.NEGATIVE_INFINITY) {
273-
return new PyComplex(Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY);
274-
}
294+
if (v == 0.) {
295+
// Pure real (and positive since in first quadrant).
296+
x = (u == 0.) ? 0. : Math.sqrt(u);
297+
y = 0.;
275298

276-
if (x.real == 0.0 && x.imag == 0.0) {
277-
return new PyComplex(0.0, Math.copySign(0.0, x.imag));
299+
} else if (u == 0.) {
300+
// Pure imaginary, and v is positive.
301+
x = y = ROOT_HALF * Math.sqrt(v);
302+
303+
} else {
304+
/*
305+
* Let w = u + iv = 2a + 2ib, and define s**2 = a**2 + b**2. Then z = x + iy is
306+
* computed as x**2 = s + a, and y = b /x. Most of the logic here is about managing
307+
* the scaling.
308+
*/
309+
int ue = Math.getExponent(u), ve = Math.getExponent(v);
310+
int diff = ue - ve;
311+
312+
if (diff > 27) {
313+
// u is so much bigger than v we can ignore v in the square: s = u/2.
314+
x = Math.sqrt(u);
315+
316+
} else if (diff < -27) {
317+
// v is so much bigger than u we can ignore u in the square: s = v/2.
318+
if (ve >= Double.MAX_EXPONENT) {
319+
x = Math.sqrt(0.5 * u + 0.5 * v); // Avoid overflow in u+v
320+
} else {
321+
x = Math.sqrt(0.5 * (u + v));
322+
}
323+
324+
} else {
325+
/*
326+
* Use the full-fat formula: s = Math.sqrt(a * a + b * b). During calculation,
327+
* we will be squaring the components, so we scale by 2**n (small values up and
328+
* large values down).
329+
*/
330+
double s, a, b;
331+
final int LARGE = 510; // 1.999... * 2**LARGE is safe to square and double
332+
final int SMALL = -510; // 1.0 * 2**(SMALL-1) may squared with full precision
333+
final int SCALE = 600; // EVEN and > (52+SMALL-Double.MIN_EXPONENT)
334+
int n = 0;
335+
if (ue > LARGE || ve > LARGE) {
336+
// One of these is too big to square without overflow.
337+
a = Math.scalb(u, -(SCALE + 1)); // a = (u/2) * 2**n
338+
b = Math.scalb(v, -(SCALE + 1));
339+
n = -SCALE;
340+
} else if (ue < SMALL && ve < SMALL) {
341+
// Both of these are too small to square without loss of bits.
342+
a = Math.scalb(u, SCALE - 1); // a = (u/2) * 2**n
343+
b = Math.scalb(v, SCALE - 1);
344+
n = SCALE;
345+
} else {
346+
a = 0.5 * u; // a = u/2
347+
b = 0.5 * v;
348+
}
349+
350+
s = Math.sqrt(a * a + b * b);
351+
x = Math.sqrt(s + a);
352+
353+
// Restore x through the square root of the scale 2**(-n/2)
354+
if (n != 0) {
355+
x = Math.scalb(x, -n / 2);
356+
}
357+
}
358+
359+
// Finally, use y = v/2x
360+
y = v / (x + x);
361+
}
278362
}
279363

280-
double t = Math.sqrt((Math.abs(x.real) + abs(x)) / 2.0);
281-
if (x.real >= 0.0) {
282-
return new PyComplex(t, x.imag / (2.0 * t));
364+
// Flip according to the signs of the components of w.
365+
if (ww.real < 0.) {
366+
return new PyComplex(y, Math.copySign(x, ww.imag));
283367
} else {
284-
return new PyComplex(Math.abs(x.imag) / (2.0 * t), Math.copySign(1d, x.imag) * t);
368+
return new PyComplex(x, Math.copySign(y, ww.imag));
285369
}
286370
}
287371

0 commit comments

Comments
 (0)