[
https://issues.apache.org/jira/browse/NUMBERS-156?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17348809#comment-17348809
]
Alex Herbert commented on NUMBERS-156:
--------------------------------------
I noted some bugs in previous versions of {{enormMod}} that use scaling.
The loop must check if the absolute value is NaN. Previously this failed by
returning zero because the maximum absolute value detected was zero.
{code:java}
@Test
void testNaN() {
final double[] v = new double[] {0, Double.NaN, 0, 0, 0};
Assertions.assertEquals(Double.NaN, SafeNorm.value(v));
}
{code}
The methods that compute a compensation term fail this additional test:
{code:java}
@Test
void testInf() {
final double[] v = new double[] {0, Double.NEGATIVE_INFINITY, 0, 0, 0};
Assertions.assertEquals(Double.POSITIVE_INFINITY, SafeNorm.value(v));
}
{code}
The problem with enormModExt and enormModKahan is that computing the round-off
involves a {{Double.POSITIVE_INFINITY - Double.POSITIVE_INFINITY}} term which
produces a NaN.
So here are the fixed version of the methods.
enormMod
{code:java}
public static double value(double[] v) {
// Sum of big, normal and small numbers with 2-fold extended precision
summation
double s1 = 0;
double s2 = 0;
double s3 = 0;
for (int i = 0; i < v.length; i++) {
final double x = Math.abs(v[i]);
if (!(x <= Double.MAX_VALUE)) {
return x;
} else if (x > 0x1.0p500) {
// Scale down big numbers
s1 += square(x * 0x1.0p-600);
} else if (x < 0x1.0p-500) {
// Scale up small numbers
s3 += square(x * 0x1.0p600);
} else {
// Unscaled
s2 += square(x);
}
}
// The highest sum is the significant component. Add the next significant.
if (s1 != 0) {
return Math.sqrt(s1 + s2 * 0x1.0p-600 * 0x1.0p-600) * 0x1.0p600;
} else if (s2 != 0) {
return Math.sqrt(s2 + s3 * 0x1.0p-600 * 0x1.0p-600);
}
return Math.sqrt(s3) * 0x1.0p-600;
}
{code}
enormModExt:
{code:java}
public static double value(double[] v) {
// Sum of big, normal and small numbers with 2-fold extended precision
summation
double s1 = 0;
double s2 = 0;
double s3 = 0;
double c1 = 0;
double c2 = 0;
double c3 = 0;
for (int i = 0; i < v.length; i++) {
final double x = Math.abs(v[i]);
if (!(x <= Double.MAX_VALUE)) {
return x;
} else if (x > 0x1.0p500) {
// Scale down big numbers
final double y = square(x * 0x1.0p-600);
final double t = s1 + y;
c1 = ExtendedPrecision.twoSumLow(s1, y, t);
s1 = t;
} else if (x < 0x1.0p-500) {
// Scale up small numbers
final double y = square(x * 0x1.0p600);
final double t = s3 + y;
c3 = ExtendedPrecision.twoSumLow(s3, y, t);
s3 = t;
} else {
// Unscaled
final double y = square(x);
final double t = s2 + y;
c2 = ExtendedPrecision.twoSumLow(s2, y, t);
s2 = t;
}
}
// The highest sum is the significant component. Add the next significant.
// Adapted from LinearCombination dot2s summation.
if (s1 != 0) {
s2 = s2 * 0x1.0p-600 * 0x1.0p-600;
c2 = c2 * 0x1.0p-600 * 0x1.0p-600;
double sum = s1 + s2;
c1 += ExtendedPrecision.twoSumLow(s1, s2, sum) + c2;
return Math.sqrt(sum + c1) * 0x1.0p600;
} else if (s2 != 0) {
s3 = s3 * 0x1.0p-600 * 0x1.0p-600;
c3 = c3 * 0x1.0p-600 * 0x1.0p-600;
double sum = s2 + s3;
c2 += ExtendedPrecision.twoSumLow(s2, s3, sum) + c3;
return Math.sqrt(sum + c2);
}
return Math.sqrt(s3) * 0x1.0p-600;
}
{code}
extLinear:
{code:java}
public static double value(double[] v) {
// Find the magnitude limits ignoring zero
double max = 0;
double min = Double.POSITIVE_INFINITY;
for (int i = 0; i < v.length; i++) {
final double x = Math.abs(v[i]);
if (Double.isNaN(x)) {
return x;
} else if (x > max) {
max = x;
} else if (x < min && x != 0) {
min = x;
}
}
// Edge cases
if (max == 0 || max == Double.POSITIVE_INFINITY) {
return max;
}
// Use scaling if required
double[] x = v;
double rescale = 1;
if (max > 0x1.0p500) {
// Too big so scale down
x = x.clone();
for (int i = 0; i < x.length; i++) {
x[i] *= 0x1.0p-600;
}
rescale = 0x1.0p600;
} else if (min < 0x1.0p-500) {
// Too small so scale up
x = x.clone();
for (int i = 0; i < x.length; i++) {
x[i] *= 0x1.0p600;
}
rescale = 0x1.0p-600;
}
return Math.sqrt(LinearCombination.value(x, x)) * rescale;
}
{code}
Here is the adapted version of LinearCombination to remove some redundant
computations:
{code:java}
public static double value(double[] v) {
// Find the magnitude limits ignoring zero
double max = 0;
double min = Double.POSITIVE_INFINITY;
for (int i = 0; i < v.length; i++) {
final double x = Math.abs(v[i]);
if (Double.isNaN(x)) {
return x;
} else if (x > max) {
max = x;
} else if (x < min && x != 0) {
min = x;
}
}
// Edge cases
if (max == 0 || max == Double.POSITIVE_INFINITY) {
return max;
}
// Below here no value is infinite or NaN.
// Use scaling if required
double scale = 1;
double rescale = 1;
if (max > 0x1.0p500) {
// Too big so scale down
scale = 0x1.0p-600;
rescale = 0x1.0p600;
} else if (min < 0x1.0p-500) {
// Too small so scale up
scale = 0x1.0p600;
rescale = 0x1.0p-600;
}
// Same as LinearCombination but with scaling.
// Splitting is safe due to scaling and only one term requires splitting.
// Implement dot2s (Algorithm 5.4) from Ogita et al (2005).
final int len = v.length;
// p is the standard scalar product sum.
// s is the sum of round-off parts.
double a = v[0] * scale;
double p = a * a;
double s = productLowUnscaled(a, p);
// Remaining split products added to the current sum and round-off sum.
for (int i = 1; i < len; i++) {
a = v[i] * scale;
final double h = a * a;
final double r = productLowUnscaled(a, h);
final double x = p + h;
// s_i = s_(i-1) + (q_i + r_i)
s += ExtendedPrecision.twoSumLow(p, h, x) + r;
p = x;
}
p += s;
return Math.sqrt(p) * rescale;
}
/**
* Compute the low part of the double length number {@code (z,zz)} for the exact
* square of {@code x} using Dekker's mult12 algorithm. The standard
* precision product {@code x*x} must be provided. The number {@code x}
* is split into high and low parts using Dekker's algorithm.
*
* <p>Warning: This method does not perform scaling in Dekker's split and large
* finite numbers can create NaN results.
*
* @param x The factor.
* @param xx Square of the factor (x * x).
* @return the low part of the product double length number
*/
private static double productLowUnscaled(double x, double xx) {
// Split the numbers using Dekker's algorithm without scaling
final double hx = ExtendedPrecision.highPartUnscaled(x);
final double lx = x - hx;
return ExtendedPrecision.productLow(hx, lx, hx, lx, xx);
}
{code}
For this to work the private method {{ExtendedPrecision.productLow}} has to be
made package private. If in a different package then just copy the method from
ExtendedPrecision.
> SafeNorm 3D overload
> --------------------
>
> Key: NUMBERS-156
> URL: https://issues.apache.org/jira/browse/NUMBERS-156
> Project: Commons Numbers
> Issue Type: Improvement
> Reporter: Matt Juntunen
> Priority: Major
> Attachments: performance-all.png, performance-len-1-5.png,
> performance2-1-5.png, performance2-all.png
>
>
> We should create an overload of {{SafeNorm.value}} that accepts 3 arguments
> to potentially improve performance for 3D vectors.
--
This message was sent by Atlassian Jira
(v8.3.4#803005)
