You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2020/01/28 17:19:56 UTC
[commons-numbers] 03/03: Remove methods for (a * c + b * d) and (b
* c - a * d) in Complex divide
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 980d80a40e67d4513f809a74e433fa9d58abfc64
Author: aherbert <ah...@apache.org>
AuthorDate: Tue Jan 28 17:19:48 2020 +0000
Remove methods for (a * c + b * d) and (b * c - a * d) in Complex divide
The scaling now reduces a and b before computing the division if there
is a chance of overflow. Thus the specialised methods to detect overflow
are redundant and cannot be hit by any input arguments, e.g.
divide(Double.MAX_VALUE, Double.MAX_VALUE, 1, 1).
This returns the listing for edge case handling to match ISO C99 G.5.1
(8).
---
.../apache/commons/numbers/complex/Complex.java | 70 ++++------------------
1 file changed, 13 insertions(+), 57 deletions(-)
diff --git a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java
index a78fe8a..04c7a6e 100644
--- a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java
+++ b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java
@@ -658,29 +658,29 @@ public final class Complex implements Serializable {
d = Math.scalb(d, -ilogbw);
}
final double denom = c * c + d * d;
- // Divisor is in the range [1, 2).
- // Avoid overflow if a or b are very big. Here we do not need to
- // handle sub-normal exponents so just get the maximum exponent.
+
+ // Note: Modification from the listing in ISO C99 G.5.1 (8):
+ // Avoid overflow if a or b are very big.
+ // Since (c, d) in the range [1, 2) the sum (ac + bd) could overflow
+ // when (a, b) are both above (Double.MAX_VALUE / 4). The same applies to
+ // (bc - ad) with large negative values.
+ // Use the maximum exponent as an approximation to the magnitude.
if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
ilogbw -= 2;
a /= 4;
b /= 4;
}
+
double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
// Recover infinities and zeros that computed as NaN+iNaN
// the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
- // --------------
- // Modification from the listing in ISO C99 G.5.1 (8):
- // Prevent overflow in (a * c + b * d) and (b * c - a * d).
- // It is only the sign that is important. not the magnitude.
- // --------------
if (Double.isNaN(x) && Double.isNaN(y)) {
if ((denom == 0.0) &&
(!Double.isNaN(a) || !Double.isNaN(b))) {
// nonzero/zero
// This case produces the same result as divide by a real-only zero
- // using divide(+/-0.0).
+ // using Complex.divide(+/-0.0)
x = Math.copySign(Double.POSITIVE_INFINITY, c) * a;
y = Math.copySign(Double.POSITIVE_INFINITY, c) * b;
} else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
@@ -688,65 +688,21 @@ public final class Complex implements Serializable {
// infinite/finite
a = boxInfinity(a);
b = boxInfinity(b);
- x = Double.POSITIVE_INFINITY * computeACplusBD(a, b, c, d);
- y = Double.POSITIVE_INFINITY * computeBCminusAD(a, b, c, d);
+ x = Double.POSITIVE_INFINITY * (a * c + b * d);
+ y = Double.POSITIVE_INFINITY * (b * c - a * d);
} else if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
Double.isFinite(a) && Double.isFinite(b)) {
// finite/infinite
c = boxInfinity(c);
d = boxInfinity(d);
- x = 0.0 * computeACplusBD(a, b, c, d);
- y = 0.0 * computeBCminusAD(a, b, c, d);
+ x = 0.0 * (a * c + b * d);
+ y = 0.0 * (b * c - a * d);
}
}
return new Complex(x, y);
}
/**
- * Compute {@code a*c + b*d} without overflow.
- * It is assumed: either {@code a} an\( b \)b} or {@code c} and {@code d} are
- * either zero or one (i.e. a boxed infinity); and the sign of the result is important,
- * not the value.
- *
- * @param a the a
- * @param b the b
- * @param c the c
- * @param d the d
- * @return The result
- */
- private static double computeACplusBD(double a, double b, double c, double d) {
- final double ac = a * c;
- final double bd = b * d;
- final double result = ac + bd;
- return Double.isFinite(result) ?
- result :
- // Overflow. Just divide by 2 as it is the sign of the result that matters.
- ac * 0.5 + bd * 0.5;
- }
-
- /**
- * Compute {@code b*c - a*d} without overflow.
- * It is assumed: either {@code a} and {@code b} or {@code c} and {@code d} are
- * either zero or one (i.e. a boxed infinity); and the sign of the result is important,
- * not the value.
- *
- * @param a the a
- * @param b the b
- * @param c the c
- * @param d the d
- * @return The result
- */
- private static double computeBCminusAD(double a, double b, double c, double d) {
- final double bc = b * c;
- final double ad = a * d;
- final double result = bc - ad;
- return Double.isFinite(result) ?
- result :
- // Overflow. Just divide by 2 as it is the sign of the result that matters.
- bc * 0.5 - ad * 0.5;
- }
-
- /**
* Returns a {@code Complex} whose value is {@code (this / divisor)},
* with {@code divisor} interpreted as a real number.
* Implements the formula: