Fran Lattanzio created MATH-1660:
------------------------------------
Summary: FastMath/AccurateMath.scalb does not handle subnormal
results properly
Key: MATH-1660
URL: https://issues.apache.org/jira/browse/MATH-1660
Project: Commons Math
Issue Type: Bug
Components: core
Affects Versions: 4.0-beta1, 3.3
Reporter: Fran Lattanzio
FastMath.scalb does not compute subnormal values correctly. I have run this
against ACM 3.3, but I see the same code in AccurateMath.scalb in 4.0 A simple
example:
{code:java}
public class ScalbTest {
@Test
public void scalbSubnormal() {
double x = -0x1.8df4a353af3a8p-2;
int e = -1024;
double wrong = FastMath.scalb(x, e);
double right = StrictMath.scalb(x, e);
Assert.assertEquals(right, wrong, 0);
}
} {code}
The code that handles subnormal outputs is incorrect in this case: Rounding
causes the loss of exactly 1/2 an ulp, but the resulting scaled mantissa is
even. Thus, it should not be rounded up. Conceptually the code needs to check
for the following 3 cases.
# Less than 1/2 ulp was lost.
# Exactly 1/2 ulp lost.
# More than 1/2 ulp lost.
For case 1, there is nothing to do. For case 2, it needs to round up only if
the mantissa is odd. For case 3, always round up.
The code below is a rough guide to what the fix should be I believe.
{code:java}
// the input is a normal number and the result is a subnormal
number
// recover the hidden mantissa bit
mantissa |= 1L << 52;
// capture the lost bits before scaling down the mantissa.
final long lostBits = mantissa & ((1L << (-scaledExponent + 1))
- 1);
mantissa >>>= 1 - scaledExponent;
// there are 3 cases to consider:
// 1. We lost less than 1/2 an ulp -> nothing to do.
// 2. We lost exactly 1/2 ulp -> round up if mantissa is odd
because we are in ties-to-even.
// 3. We lost more than 1/2 ulp -> round up.
final long halfUlp = (1L << (-scaledExponent));
if((lostBits == halfUlp && (mantissa & 1L) == 1L) || lostBits >
halfUlp) {
mantissa++;
}
return Double.longBitsToDouble(sign | mantissa); {code}
This code is probably not perfect, and I have not tested it on a large range of
inputs, but it demonstrates the basic idea.
--
This message was sent by Atlassian Jira
(v8.20.10#820010)
