[
https://issues.apache.org/jira/browse/MATH-1660?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17717746#comment-17717746
]
Alex Herbert commented on MATH-1660:
------------------------------------
I confirm this is a bug in the bit manipulation of the floating point
representation. An alternative implementation of scalb can use lossless
multiplication by powers of 2. The final multiplication must be performed on a
normal number so any sub-normal result has the correct rounding (i.e. the only
rounding operation is the final one).
{noformat}
jshell> double x = -0x1.8df4a353af3a8p-2
x ==> -0.3886285323931653
jshell> int e = -1024
e ==> -1024
jshell> Math.scalb(x, e)
$25 ==> -2.161817970245126E-309
jshell> x * 0x1.0p-512 * 0x1.0p-512
$26 ==> -2.161817970245126E-309
// Pasted the AccurateMath scalb function into the jshell
jshell> scalb(x, e)
$29 ==> -2.16181797024513E-309
{noformat}
Note that any change to the scalb(double) function would have to be replicated
in the scalb(float) function as the code is essentially the same. It would be
useful to have equivalent float failing test cases.
> 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: 3.3, 4.0-beta1
> Reporter: Fran Lattanzio
> Priority: Major
>
> 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)
