[
https://issues.apache.org/jira/browse/NUMBERS-156?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17350436#comment-17350436
]
Alex Herbert edited comment on NUMBERS-156 at 5/24/21, 1:42 PM:
----------------------------------------------------------------
The method {{extLinearSinglePass}} is missing from the performance chart for
lengths 1 to 5. I am guessing it is between the {{enormModExt}} and the
{{extLinearMod}} in that chart.
The single pass extended precision algorithm is a nice idea and the algorithm
looks good. I noted on review that the NaN detection can be changed to handle
infinite inputs:
{code:java}
if (Double.isNaN(x)) {
// found an NaN; no use in continuing
return x;
} ...
{code}
to
{code:java}
if (!Double.isFinite(x)) {
// found a non-finite input; no use in continuing.
// even if this is infinite we must check if a NaN is present
// later in the input for the correct result.
return nanOrPositiveInfinity(v, i);
} ...
private static double nanOrPositiveInfinity(double[] v, int start) {
for (int i = start; i < v.length; i++) {
if (Double.isNaN(v[i])) {
return Double.NaN;
}
}
return Double.POSITIVE_INFINITY;
}
{code}
This will just skip out if a non-finite is found and return either positive
infinity or NaN. It is an edge case that may be worth including. The effect on
performance for normal input will be the difference between Double.isNaN and
!Double.isFinite. I would think this is negligible. In the rare case the input
is long and contains non-finite values then the change will be useful.
The extended precision algorithms in [Dekker's original
paper|https://link.springer.com/article/10.1007/BF01397083] (see the last page)
include an extended precision square root. If you have a split number with high
and low parts you compute the square root using:
{code:java}
/**
* Compute the extended precision square root from the split number
* {@code x, xx}.
* This is a modification of Dekker's sqrt2 algorithm to ignore the
* roundoff of the square root.
*
* @param x the high part
* @param xx the low part
* @return the double
*/
static double sqrt2(double x, double xx) {
if (x > 0) {
double c = Math.sqrt(x);
double u = c * c;
//double uu = ExtendedPrecision.productLow(c, c, u);
// Here we use the optimised version:
double uu = productLowUnscaled(c, u);
double cc = (x - u - uu + xx) * 0.5 / c;
// Extended precision sqrt (y, yy)
// y = c + cc
// yy = c - y + cc (ignored)
return c + cc;
}
return x;
}
{code}
This could be applied to the extended precision sum. E.g.
{code:java}
return Math.sqrt(sum + c1) * SCALE_UP
{code}
Becomes:
{code:java}
return sqrt2(sum, c1) * SCALE_UP
{code}
However the sum of the round off parts may not be an non-overlapping two part
split number (x,xx) where xx is so small such that x + xx == x. This may be
required for the algorithm to work. This can be done using:
{code:java}
double x = sum + c1;
double xx = ExtendedPrecision.twoSumLow(sum, c1, x);
return sqrt2(x, xx) * SCALE_UP;
{code}
When I looked at this algorithm for the 2D case for use in the Complex class
it matched the result of BigDecimal.sqrt in all cases I tried. But the
performance overhead was not worth it for 2D as it gains only 1 ULP precision.
For longer arrays the performance impact is relatively less and it may be worth
considering.
was (Author: alexherbert):
The method {{extLinearSinglePass}} is missing from the performance chart for
lengths 1 to 5. I am guessing it is between the {{enormModExt}} and the
{{extLinearMod}} in that chart.
The single pass extended precision algorithm is a nice idea and the algorithm
looks good. I noted on review that the NaN detection can be changed to handle
infinite inputs:
{code:java}
if (Double.isNaN(x)) {
// found an NaN; no use in continuing
return x;
} ...
{code}
to
{code:java}
if (!Double.isFinite(x)) {
// found a non-finite input; no use in continuing.
// even if this is infinite we must check if a NaN is present
// later in the input for the correct result.
return nanOrPositiveInifinty(v, i);
} ...
private static double nanOrPositiveInfinity(double[] v, int start) {
for (int i = start; i < v.length; i++) {
if (Double.isNaN(v[i])) {
return Double.NaN;
}
}
return Double.POSITIVE_INFINITY;
}
{code}
This will just skip out if a non-finite is found and return either positive
infinity or NaN. It is an edge case that may be worth including. The effect on
performance for normal input will be the difference between Double.isNaN and
!Double.isFinite. I would think this is negligible. In the rare case the input
is long and contains non-finite values then the change will be useful.
The extended precision algorithms in [Dekker's original
paper|https://link.springer.com/article/10.1007/BF01397083] (see the last page)
include an extended precision square root. If you have a split number with high
and low parts you compute the square root using:
{code:java}
/**
* Compute the extended precision square root from the split number
* {@code x, xx}.
* This is a modification of Dekker's sqrt2 algorithm to ignore the
* roundoff of the square root.
*
* @param x the high part
* @param xx the low part
* @return the double
*/
static double sqrt2(double x, double xx) {
if (x > 0) {
double c = Math.sqrt(x);
double u = c * c;
//double uu = ExtendedPrecision.productLow(c, c, u);
// Here we use the optimised version:
double uu = productLowUnscaled(c, u);
double cc = (x - u - uu + xx) * 0.5 / c;
// Extended precision sqrt (y, yy)
// y = c + cc
// yy = c - y + cc (ignored)
return c + cc;
}
return x;
}
{code}
This could be applied to the extended precision sum. E.g.
{code:java}
return Math.sqrt(sum + c1) * SCALE_UP
{code}
Becomes:
{code:java}
return sqrt2(sum, c1) * SCALE_UP
{code}
However the sum of the round off parts may not be an non-overlapping two part
split number (x,xx) where xx is so small such that x + xx == x. This may be
required for the algorithm to work. This can be done using:
{code:java}
double x = sum + c1;
double xx = ExtendedPrecision.twoSumLow(sum, c1, x);
return sqrt2(x, xx) * SCALE_UP;
{code}
When I looked at this algorithm for the 2D case for use in the Complex class
it matched the result of BigDecimal.sqrt in all cases I tried. But the
performance overhead was not worth it for 2D as it gains only 1 ULP precision.
For longer arrays the performance impact is relatively less and it may be worth
considering.
> 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, performance3-1-5.png,
> performance3-all.png, performance4-1-5.png, performance4-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)
