[ 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)