From 3f71fe3998c947cd80218406ebe9b78a4e27399e Mon Sep 17 00:00:00 2001
From: Nick Bowler <nbowler@draconx.ca>
Date: Fri, 19 Feb 2010 17:47:59 -0500
Subject: [PATCH] floating: Make fromRational work properly.

The current implementation of fromRational (for both floating types) is
a simple floating point division, which does not work correctly when
either the numerator or the denominator cannot be represented exactly.

We do a correct conversion by scaling the rational by a power of two to
get a large integer that makes maximum use of the available precision.
Since such an integer is exactly representable, we can use scalb to get
the actual value.
---
 Data/Floating/Helpers.hs      | 40 ++++++++++++++++++++++++++++++++++++++++
 Data/Floating/Types/Double.hs |  7 ++++---
 Data/Floating/Types/Float.hs  |  7 ++++---
 altfloat.cabal                |  3 ++-
 4 files changed, 50 insertions(+), 7 deletions(-)
 create mode 100644 Data/Floating/Helpers.hs

diff --git a/Data/Floating/Helpers.hs b/Data/Floating/Helpers.hs
new file mode 100644
index 0000000..760c84e
--- /dev/null
+++ b/Data/Floating/Helpers.hs
@@ -0,0 +1,40 @@
+{-
+ - Copyright (C) 2010 Nick Bowler.
+ -
+ - License BSD2:  2-clause BSD license.  See LICENSE for full terms.
+ - This is free software: you are free to change and redistribute it.
+ - There is NO WARRANTY, to the extent permitted by law.
+ -}
+
+-- | Internal helper functions needed by at least two modules.
+module Data.Floating.Helpers where
+
+import Prelude hiding (RealFloat(..), RealFrac(..))
+import Data.Floating.Classes
+import Data.Floating.Instances
+import Data.Ratio
+import Data.Maybe
+
+-- | @binarySearch p low high@ computes the least integer on the interval
+-- [low, high] satisfying the given predicate, by binary search.  The
+-- predicate must partition the interval into two contiguous regions.
+binarySearch :: Integral a => (a -> Bool) -> a -> a -> a
+binarySearch p l u
+    | l >  u    = error "empty interval"
+    | l == u    = m
+    | p m       = binarySearch p l m
+    | otherwise = binarySearch p (m+1) u
+    where
+        m = l + div (u-l) 2
+
+-- | Find a power of two such that the given rational number, when multiplied
+-- by that power and rounded to an integer, has exactly as many digits as the
+-- precision of the floating point type.  The search may stop for values with
+-- extremely large (or small) magnitude, in which case the result shall
+-- overflow (or underflow) when scaled to the floating type.
+scaleRational :: PrimFloat a => a -> Rational -> (Integer, Int)
+scaleRational t x = (fromJust . toIntegral . round . scale x $ e, e) where
+    e         = binarySearch ((lbound <) . scale x) l (u*2)
+    (l, u)    = floatRange t
+    lbound    = floatRadix t ^ (floatPrecision t - 1)
+    scale x y = x * 2^^y
diff --git a/Data/Floating/Types/Double.hs b/Data/Floating/Types/Double.hs
index a27b52b..22c9bd6 100644
--- a/Data/Floating/Types/Double.hs
+++ b/Data/Floating/Types/Double.hs
@@ -30,6 +30,7 @@ import Foreign.C
 import System.IO.Unsafe
 
 import Data.Floating.Types
+import Data.Floating.Helpers
 import Data.Floating.Classes
 import Data.Floating.CMath
 
@@ -96,9 +97,9 @@ instance Sortable Double where
 
 instance Fractional Double where
     (D# x) / (D# y) = D# (x /## y)
-    fromRational = liftM2 (/)
-        (fromInteger . numerator)
-        (fromInteger . denominator)
+    fromRational x = scalb (toFloating s) (negate e) where
+        scale  = scaleRational (undefined :: Double)
+        (s, e) = scale x
 
 -- | Internal function which discards the fractional component of a Double.
 -- The results are meaningful only for finite input.
diff --git a/Data/Floating/Types/Float.hs b/Data/Floating/Types/Float.hs
index 101d005..a18ee5c 100644
--- a/Data/Floating/Types/Float.hs
+++ b/Data/Floating/Types/Float.hs
@@ -29,6 +29,7 @@ import Foreign.C
 import System.IO.Unsafe
 
 import Data.Floating.Types
+import Data.Floating.Helpers
 import Data.Floating.Types.Double
 import Data.Floating.Classes
 import Data.Floating.CMath
@@ -91,9 +92,9 @@ instance Sortable Float where
 
 instance Fractional Float where
     (F# x) / (F# y) = F# (x `divideFloat#` y)
-    fromRational = liftM2 (/)
-        (fromInteger . numerator)
-        (fromInteger . denominator)
+    fromRational x = scalb (toFloating s) (negate e) where
+        scale  = scaleRational (undefined :: Float)
+        (s, e) = scale x
 
 -- | Internal function which discards the fractional component of a Float.
 -- The results are meaningful only for finite input.
diff --git a/altfloat.cabal b/altfloat.cabal
index a7ef510..125d4ed 100644
--- a/altfloat.cabal
+++ b/altfloat.cabal
@@ -84,4 +84,5 @@ Library
         Data.Floating,
         Data.Poset
     Other-Modules:
-        Data.Floating.Instances, Data.Poset.Internal, Data.Poset.Instances
+        Data.Floating.Instances, Data.Floating.Helpers,
+        Data.Poset.Internal, Data.Poset.Instances
-- 
2.11.4.GIT