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

["Fran Lattanzio (Jira)" <ji...@apache.org>]

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)