Opened 12 years ago

Closed 8 years ago

#2281 closed bug (wontfix)

properFraction implemented with modf primitive?

Reported by: guest Owned by:
Priority: normal Milestone: 7.2.1
Component: libraries/base Version: 6.8.2
Keywords: Cc: ghc@…
Operating System: Unknown/Multiple Architecture: Unknown/Multiple
Type of failure: None/Unknown Test Case:
Blocked By: Blocking:
Related Tickets: Differential Rev(s):
Wiki Page:

Description

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.

Change History (14)

comment:1 Changed 12 years ago by 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.

comment:2 Changed 12 years ago by guest

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

comment:3 Changed 12 years ago by guest

Component: Compilerlibraries/base

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?

comment:4 Changed 12 years ago by igloo

difficulty: Unknown
Milestone: 6.10 branch

comment:5 Changed 11 years ago by simonmar

Architecture: MultipleUnknown/Multiple

comment:6 Changed 11 years ago by simonmar

Operating System: MultipleUnknown/Multiple

comment:7 Changed 11 years ago by igloo

Milestone: 6.10 branch6.12.1
Owner: set to igloo

comment:8 Changed 11 years ago by igloo

Type: proposalbug

comment:9 Changed 10 years ago by igloo

Milestone: 6.12.16.12 branch
Type of failure: None/Unknown

comment:10 Changed 10 years ago by igloo

Milestone: 6.12 branch6.12.3

comment:11 Changed 10 years ago by igloo

Milestone: 6.12.36.14.1
Owner: igloo deleted

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.

comment:12 Changed 9 years ago by igloo

Milestone: 7.0.17.0.2

comment:13 Changed 9 years ago by igloo

Milestone: 7.0.27.2.1

comment:14 Changed 8 years ago by igloo

Resolution: wontfix
Status: newclosed

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.

Note: See TracTickets for help on using tickets.