[
https://issues.apache.org/jira/browse/MATH-1660?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17717798#comment-17717798
]
Fran Lattanzio commented on MATH-1660:
--------------------------------------
I will try to fix this function and submit a PR. In terms of the documentation,
scalb/ldexp are not part of the IEEE 754 spec.
There is this:
[https://pubs.opengroup.org/onlinepubs/009696899/functions/scalb.html]
Basically I think it is OK to link to something like that, with the
understanding that the contract for scalb is that it behaves as if it computes
2^n*x exactly and then applying round-to-nearest/ties to even.
Or stated another way it should compute fl(2^n*x) (where * is real
multiplication) and not 2^n [mul] x (where [mul] is floating point
multiplication). The three enumerated cases above fall out naturally from that
algorithm when considering subnormal outputs.
It's probably also possible to write this function using only multiplications
rather than bit-level hacking, in which case the hardware would do everything
properly.
> 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)
