properFraction implemented with modf primitive?

[guest]

I need a fast 'fraction' function for Double. I can use 'properFraction', but its Double instance uses 'decodeFloat', which is rather slow. It's also clumsy to use, because I have to provide an Integral type, although I'm not interested in the integral part. I can use (x - int2Double (double2Int x)) but this fails for x beyond the Int number range, and it is hard to fix that for an unknown implementation of Double.

What about a 'modf' primitive which either calls 'modf' from standard C library or invokes an appropriate FPU command? If Double is in IEEE format the 'fraction' command could also be implemented quite efficiently by some bit masking without the FPU.

[dons]

Can you use modf from the cmath library?

 http://hackage.haskell.org/packages/archive/cmath/0.3/doc/html/Foreign-C-Math-Double.html#v%3Amodf

for your immediate use case.

[guest]

Thanks for the pointer. Currently I'm using '(x - int2Double (double2Int x))' since I expect it's the fastest.

[guest]

We could insert code like the following to GHC.Float (attention, untested):

instance  RealFrac Double  where
    {-# INLINE properFraction #-}
    properFraction = fastProperFraction double2Int int2Double

{-# INLINE fastProperFraction #-}
fastProperFraction :: (Integral b) =>
   (a -> Int) -> (Int -> a) -> a -> (b,a)
fastProperFraction trunc toFloat x =
   if toFloat minBound <= x && x <= toFloat maxBound
     then let xt = trunc x
          in  (fromIntegral xt, x - toFloat xt)
     else decodeProperFraction x

decodeProperFraction :: (Integral b, RealFloat a) =>
   a -> (b, a)
decodeProperFraction x
  = case (decodeFloat x)      of { (m,n) ->
    let  b = floatRadix x     in
    if n >= 0 then
        (fromInteger m * fromInteger b ^ n, 0.0)
    else
        case (quotRem m (b^(negate n))) of { (w,r) ->
        (fromInteger w, encodeFloat r n)
        }
    }

How about that?

[igloo]

I've confirmed performance is worse with this program:

{-# LANGUAGE ForeignFunctionInterface #-}

import Foreign
import Foreign.C

main :: IO ()
main = f

count :: Int
count = 100000000

f :: IO ()
f = with 0 $ \dptr ->
    do let loop :: Int -> Double -> Double
           loop 0 d = d
           loop n d = modf (realToFrac d) dptr `seq` loop (n - 1) (d + 0.01)
       print $ loop count 100.456

g :: IO ()
g = do let loop :: Int -> Double -> Double
           loop 0 d = d
           loop n d = case properFraction d :: (Int, Double) of
                       (_, frac) -> frac `seq` loop (n - 1) (d + 0.01)
       print $ loop count 100.456

foreign import ccall "modf" modf :: CDouble -> Ptr CDouble -> CDouble

./f  8.22s user 0.00s system 99% cpu 8.231 total
./g  38.35s user 0.03s system 100% cpu 38.382 total

but I'm a bit nervous about the Integral type used affecting the result:

    properFraction x
      = case (decodeFloat x)      of { (m,n) ->
        let  b = floatRadix x     in
        if n >= 0 then
            (fromInteger m * fromInteger b ^ n, 0.0)
        else
            case (quotRem m (b^(negate n))) of { (w,r) ->
            (fromInteger w, encodeFloat r n)
            }
        }

I'd suggest the way forward is for someone to make a patch and then turn this ticket into a  library submission.

[igloo]

I'm going to close this ticket, but I suggest that if anyone is interested in this then they propose a patch on the libraries list, and argue why it is the right thing to do.