[jira] [Commented] (MATH-1660) FastMath/AccurateMath.scalb does not handle subnormal results properly

["Gilles Sadowski (Jira)" <ji...@apache.org>]

    [ https://issues.apache.org/jira/browse/MATH-1660?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17717698#comment-17717698 ] 

Gilles Sadowski commented on MATH-1660:
---------------------------------------

Thanks for the report!
Could you please provide the required modifications (code fix and unit test) in the form of a patch, or a PR on GitHub?
The comments are quite nice but, for ease of future maintenance, it would be great to add a link to a reference document that describes the correct behaviour (as a explanation for why the changes must be those proposed).

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