Skip to content

Commit 030199f

Browse files
jeff5jeff5
committed
Re-work cmath.log, and cmath.log10 for accuracy and corner cases.
1 parent fe59267 commit 030199f

2 files changed

Lines changed: 116 additions & 36 deletions

File tree

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

Lines changed: 115 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,9 @@ public class cmath {
2323
/** ln({@link Double#MAX_VALUE}) or a little less */
2424
private static final double NEARLY_LN_DBL_MAX = 709.4361393;
2525

26+
/** log<sub>10</sub>e (Ref: Abramowitz &amp; Stegun [1972], p3). */
27+
private static final double LOG10E = 0.43429448190325182765;
28+
2629
private static PyComplex c_prodi(PyComplex x) {
2730
return (PyComplex)x.__mul__(i);
2831
}
@@ -204,18 +207,6 @@ public static PyComplex exp(PyObject z) {
204207
return exceptNaN(new PyComplex(u, v), zz);
205208
}
206209

207-
public static PyComplex log(PyObject in) {
208-
PyComplex x = complexFromPyObject(in);
209-
if (isNaN(x)) {
210-
if (Double.isInfinite(x.real) || Double.isInfinite(x.imag)) {
211-
return new PyComplex(Double.POSITIVE_INFINITY, Double.NaN);
212-
} else {
213-
return PyComplex.NaN;
214-
}
215-
}
216-
return new PyComplex(Math.log(abs(x)), Math.atan2(x.imag, x.real));
217-
}
218-
219210
public static double phase(PyObject in) {
220211
PyComplex x = complexFromPyObject(in);
221212
return Math.atan2(x.imag, x.real);
@@ -269,36 +260,125 @@ public static boolean isnan(PyObject in) {
269260
return Double.isNaN(x.real) || Double.isNaN(x.imag);
270261
}
271262

272-
public static PyComplex log10(PyObject in) {
273-
PyComplex x = complexFromPyObject(in);
274-
if (isNaN(x)) {
275-
if (Double.isInfinite(x.real) || Double.isInfinite(x.imag)) {
276-
return new PyComplex(Double.POSITIVE_INFINITY, Double.NaN);
277-
} else {
278-
return PyComplex.NaN;
279-
}
280-
}
281-
double l = abs(x);
282-
return new PyComplex(math.log10(new PyFloat(l)), Math.atan2(x.imag, x.real)
283-
/ Math.log(10.0));
263+
/**
264+
* Returns the natural logarithm of <i>w</i>.
265+
*
266+
* @param w
267+
* @return ln <i>w</i>
268+
*/
269+
public static PyComplex log(PyObject w) {
270+
PyComplex ww = complexFromPyObject(w);
271+
double u = ww.real, v = ww.imag;
272+
// The real part of the result is the log of the magnitude.
273+
double lnr = logHypot(u, v);
274+
// The imaginary part of the result is the arg. This may result in a nan.
275+
double theta = Math.atan2(v, u);
276+
PyComplex z = new PyComplex(lnr, theta);
277+
return exceptNaN(z, ww);
278+
}
279+
280+
/**
281+
* Returns the common logarithm of <i>w</i> (base 10 logarithm).
282+
*
283+
* @param w
284+
* @return log<sub>10</sub><i>w</i>
285+
*/
286+
public static PyComplex log10(PyObject w) {
287+
PyComplex ww = complexFromPyObject(w);
288+
double u = ww.real, v = ww.imag;
289+
// The expression is the same as for base e, scaled in magnitude.
290+
double logr = logHypot(u, v) * LOG10E;
291+
double theta = Math.atan2(v, u) * LOG10E;
292+
PyComplex z = new PyComplex(logr, theta);
293+
return exceptNaN(z, ww);
284294
}
285295

286-
public static PyComplex log(PyObject in, PyObject base) {
287-
return log(complexFromPyObject(in), complexFromPyObject(base));
296+
/**
297+
* Returns the logarithm of <i>w</i> to the given base. If the base is not specified, returns
298+
* the natural logarithm of <i>w</i>. There is one branch cut, from 0 along the negative real
299+
* axis to -&infin;, continuous from above.
300+
*
301+
* @param w
302+
* @param b
303+
* @return log<sub>b</sub><i>w</i>
304+
*/
305+
public static PyComplex log(PyObject w, PyObject b) {
306+
PyComplex ww = complexFromPyObject(w), bb = complexFromPyObject(b), z;
307+
double u = ww.real, v = ww.imag, br = bb.real, bi = bb.imag, x, y;
308+
// Natural log of w is (x,y)
309+
x = logHypot(u, v);
310+
y = Math.atan2(v, u);
311+
312+
if (bi != 0. || br <= 0.) {
313+
// Complex or negative real base requires complex log: general case.
314+
PyComplex lnb = log(bb);
315+
z = (PyComplex)(new PyComplex(x, y)).__div__(lnb);
316+
317+
} else {
318+
// Real positive base: frequent case. (b = inf or nan ends up here too.)
319+
double lnb = Math.log(br);
320+
z = new PyComplex(x / lnb, y / lnb);
321+
}
322+
323+
return exceptNaN(z, ww);
288324
}
289325

290-
public static PyComplex log(PyComplex x, PyComplex base) {
291-
if (isNaN(x)) {
292-
if (Double.isInfinite(x.real) || Double.isInfinite(x.imag)) {
293-
return new PyComplex(Double.POSITIVE_INFINITY, Double.NaN);
326+
/**
327+
* Helper function for the log of a complex number, dealing with the log magnitude, and without
328+
* intermediate overflow or underflow. It returns ln <i>r</i>, where <i>r<sup>2</sup> =
329+
* u<sup>2</sup>+v<sup>2</sup></i>. To do this it computes
330+
* &frac12;ln(u<sup>2</sup>+v<sup>2</sup>). Special cases are handled as follows:
331+
* <ul>
332+
* <li>if u or v is NaN, it returns NaN</li>
333+
* <li>if u or v is infinite, it returns positive infinity</li>
334+
* <li>if u and v are both zero, it raises a ValueError</li>
335+
* </ul>
336+
* We have this function instead of <code>Math.log(Math.hypot(u,v))</code> because a valid
337+
* result is still possible even when <code>hypot(u,v)</code> overflows, and because there's no
338+
* point in taking a square root when a log is to follow.
339+
*
340+
* @param u
341+
* @param v
342+
* @return &frac12;ln(u<sup>2</sup>+v<sup>2</sup>)
343+
*/
344+
private static double logHypot(double u, double v) {
345+
346+
if (Double.isInfinite(u) || Double.isInfinite(v)) {
347+
return Double.POSITIVE_INFINITY;
348+
349+
} else {
350+
// Cannot overflow, but if u=v=0 will return -inf.
351+
int scale = 0, ue = Math.getExponent(u), ve = Math.getExponent(v);
352+
double lnr;
353+
354+
if (ue < -511 && ve < -511) {
355+
// Both u and v are too small to square, or zero. (Just one would be ok.)
356+
scale = 600;
357+
} else if (ue > 510 || ve > 510) {
358+
// One of these is too big to square and double (or is nan or inf).
359+
scale = -600;
360+
}
361+
362+
if (scale == 0) {
363+
// Normal case: there is no risk of overflow or log of zero.
364+
lnr = 0.5 * Math.log(u * u + v * v);
365+
} else {
366+
// We must work with scaled values, us = u * 2**n etc..
367+
double us = Math.scalb(u, scale);
368+
double vs = Math.scalb(v, scale);
369+
// rs**2 = r**2 * 2**2n
370+
double rs2 = us * us + vs * vs;
371+
// So ln(r) = ln(u**2+v**2)/2 = ln(us**2+vs**2)/2 - n ln(2)
372+
lnr = 0.5 * Math.log(rs2) - scale * math.LN2;
373+
}
374+
375+
// (u,v) = 0 leads to ln(r) = -inf, but that's a domain error
376+
if (lnr == Double.NEGATIVE_INFINITY) {
377+
throw math.mathDomainError();
294378
} else {
295-
return PyComplex.NaN;
379+
return lnr;
296380
}
297381
}
298-
double l = abs(x);
299-
PyComplex log_base = log(base);
300-
return (PyComplex)new PyComplex(math.log(new PyFloat(l)), Math.atan2(x.imag, x.real))
301-
.__div__(log_base);
302382
}
303383

304384
public static PyComplex sin(PyObject in) {

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@ public class math implements ClassDictInit {
2424
private static final double MINUS_ONE = -1.0;
2525
private static final double TWO = 2.0;
2626
private static final double EIGHT = 8.0;
27-
private static final double LN2 = 0.693147180559945309417232121458; // Ref OEIS A002162
27+
static final double LN2 = 0.693147180559945309417232121458; // Ref OEIS A002162
2828

2929
private static final double INF = Double.POSITIVE_INFINITY;
3030
private static final double NINF = Double.NEGATIVE_INFINITY;

0 commit comments

Comments
 (0)