1 {-# LANGUAGE ScopedTypeVariables #-}
    2 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
    3 {-# LANGUAGE FlexibleInstances   #-}
    4 
    5 module Fixed (
    6  -- type classes
    7     Has_Radix(..),
    8  -- types
    9     Binary,
   10     Decimal,
   11     Fixed(..),
   12     Fixed_Frame(..),
   13     Result_Scale(..),
   14     Rounding_Mode(..),
   15     Scale(..),
   16  -- functions
   17     scale,
   18     shift,
   19     toFixed,
   20     align,
   21     lia_divides,    lia_gcd,    lia_lcm,
   22     lia_signum,
   23     lia_dim,
   24     convert,
   25     rescale,
   26     convertFixed,
   27     quotient,
   28     fit,
   29     divide,
   30     readFixedDecimal
   31 ) where
   32 
   33 import Numeric (showSigned)
   34 import Data.Char (isDigit)
   35 import Data.Ratio
   36 
   37 data Binary
   38 data Decimal
   39 
   40 class Has_Radix t
   41   where radix :: t -> Integer
   42 
   43 instance Has_Radix Binary
   44   where radix _ = 2
   45 
   46 instance Has_Radix Decimal
   47   where radix _ = 10
   48 
   49 newtype Scale t = Scale Integer
   50                   deriving (Eq, Ord, Show, Num)
   51 
   52 shift :: forall t . Has_Radix t => Integer -> Scale t -> Integer
   53 shift m (Scale p) =
   54   if p >= 0 then m * radix (undefined :: t) ^ p
   55   else m `quot` radix (undefined :: t) ^ negate p
   56 
   57 data Fixed t = FP !Integer !(Scale t)
   58 
   59 scale :: forall t . Has_Radix t => Fixed t -> Scale t
   60 scale (FP _ s) = s
   61 
   62 ulp :: Has_Radix t => Fixed t -> Fixed t
   63 ulp (FP _ s) = FP 1 s
   64 
   65 toFixed :: Has_Radix t => Integer -> Scale t -> Fixed t
   66 toFixed n s = FP (shift n s) s
   67 
   68 lift :: forall t a . Has_Radix t =>
   69     (Integer -> Integer) -> Fixed t -> Fixed t
   70 
   71 lift f (FP m s) = FP (f m) s
   72 
   73 align :: forall t a . Has_Radix t =>
   74     Fixed t -> Fixed t -> (Integer -> Integer -> a) -> a
   75 
   76 align (FP m1 s1) (FP m2 s2) f
   77   | s1 > s2 = f m1 (shift m2 (s1-s2))
   78   | s1 < s2 = f (shift m1 (s2-s1)) m2
   79   | True    = f m1 m2
   80 
   81 instance (Has_Radix t) => Eq (Fixed t)
   82   where x1 == x2 = align x1 x2 (==)
   83 
   84   -- automatically defines /=
   85 
   86 instance (Has_Radix t) => Ord (Fixed t)
   87   where compare x1 x2 = align x1 x2 compare
   88 
   89   -- automatically defines < <= > >= max min
   90 
   91 lia_divides :: Has_Radix t => Fixed t -> Fixed t -> Bool
   92 lia_divides x y = align x y (\m n -> m /= 0 && n`rem`m == 0)
   93 
   94 lia_gcd :: Has_Radix t => Fixed t -> Fixed t -> Fixed t
   95 lia_gcd x y = FP (align x y gcd) (scale x `max` scale y)
   96 
   97 lia_lcm :: Has_Radix t => Fixed t -> Fixed t -> Fixed t
   98 lia_lcm x y = FP (align x y lcm) (scale x `max` scale y)
   99 
  100 instance (Has_Radix t) => Show (Fixed t)
  101   where
  102     showsPrec p x = showSigned showPos p x
  103       where showPos x@(FP m s) rest =
  104               if s <= 0 then shows (shift m (negate s)) rest
  105               else shows i ('.' : show_fract s r)
  106               where p = shift 1 s
  107                     (i,r) = quotRem m p
  108                     show_fract 0 _ = rest
  109                     show_fract d f = shows h (show_fract (d-1) l)
  110                       where (h,l) = quotRem (10 * f) p
  111 
  112 
  113 instance (Has_Radix t) => Num (Fixed t)
  114   where
  115     x + y = FP (align x y (+)) (scale x `max` scale y)
  116     x - y = FP (align x y (-)) (scale x `max` scale y)
  117     negate = lift negate
  118     abs    = lift abs
  119 --  negate (FP m s) = FP (negate m) s
  120 --  abs    (FP m s) = FP (abs    m) s
  121     signum (FP m s) = FP (signum m) 0
  122     fromInteger n   = FP n 0
  123     (FP m1 s1) * (FP m2 s2)  = FP (m1*m2) (s1+s2)
  124 
  125 lia_signum :: (Ord t, Num t) => t -> t
  126 lia_signum x = if x < 0 then -1 else 1
  127 
  128 lia_dim :: (Ord t, Num t) => t -> t -> t
  129 lia_dim x y = if x < y then fromInteger 0 else x - y
  130 
  131 instance (Has_Radix t) => Real (Fixed t)
  132   where
  133     toRational (FP m s) =
  134       if s >= 0 then m % (shift 1 s) else (shift m (negate s)) % 1
  135 
  136 -- Integral has
  137   -- quotRem, quot, rem
  138   -- divMod, div, mod
  139   -- toInteger
  140 
  141 instance Has_Radix t => Enum (Fixed t)
  142   where
  143     succ                 = lift succ
  144     pred                 = lift pred
  145 --  succ (FP m s)        = FP (succ m) s
  146 --  pred (FP m s)        = FP (pred m) s
  147     toEnum n             = FP (fromIntegral n) 0
  148     fromEnum x           = error "range too small"
  149     enumFrom x           = iterate (+ ulp x) x
  150     enumFromThen x y     = iterate (+ (y-x)) x
  151     enumFromTo x z       = takeWhile (<= z) $ enumFrom x
  152     enumFromThenTo x y z
  153       | x > y            = takeWhile (>= z) $ enumFromThen x y
  154       | x <= y           = takeWhile (<= z) $ enumFromThen x y
  155 
  156 instance (Has_Radix t) => Integral (Fixed t)
  157   where
  158     quotRem x y = (q, x - q*y) where q = FP (align x y quot) 0
  159     divMod  x y = (q, x - q*y) where q = FP (align x y div ) 0
  160     toInteger x@(FP m s) =
  161       if s > 0 then m `quot` (shift 1 s) else shift m (negate s)
  162 
  163 -- Fractional has
  164   -- / recip fromRational
  165 
  166 -- RealFrac has
  167   -- properFraction   :: (Integral b) => a -> (b,a)
  168   -- properFraction x = i f   such that i integral, i+f = x
  169   -- truncate, round  :: (Integral b) => a -> b
  170   -- round is round to even
  171   -- ceiling, floor   :: (Integral b) => a -> b
  172 
  173 data Rounding_Mode
  174    = Down
  175    | Up
  176    | In
  177    | Out
  178    | Exact
  179    | Even
  180    | Nearest Rounding_Mode
  181 
  182 convert :: RealFrac t => Rounding_Mode -> t -> Integer
  183 
  184 convert round x =
  185   case round of
  186     Down      -> if d < 0 then i-1 else i
  187     Up        -> if d > 0 then i+1 else i
  188     In        -> i
  189     Out       -> i + d
  190     Exact     -> if d == 0 then i else error "inexact"
  191     Even      -> if even i then i else
  192                  if d /= 0 then i + d else
  193                  if mod (i + 1) 4 == 0 then i+1 else i-1
  194     Nearest h -> if d == 0 then i else
  195                  case compare (2*abs f) 1 of
  196                    LT -> i
  197                    GT -> i + d
  198                    EQ -> i + convert h f
  199   where (i,f) = properFraction x
  200         d     = if f < 0 then -1 else if f > 0 then 1 else 0
  201 
  202 rescale :: forall t r . (Has_Radix t, Real r) =>
  203            Rounding_Mode -> Integer -> r -> Fixed t
  204 
  205 rescale round digits x = FP (convert round y) (Scale digits)
  206   where y = toRational x * toRational (radix (undefined :: t) ^ digits)
  207 
  208 convertFixed :: Has_Radix t =>
  209     Rounding_Mode -> Fixed t -> (Integer, Fixed t)
  210 
  211 convertFixed round x = (q, x - fromInteger q)
  212   where q = convert round (toRational x)
  213 
  214 quotient :: Has_Radix t =>
  215     Rounding_Mode -> Fixed t -> Fixed t ->
  216     (Integer, Fixed t)
  217 quotient round x y = (q, x - fromInteger q * y)
  218   where q = convert round (toRational x / toRational y)
  219 
  220 data Result_Scale
  221    = Exact_Scale   !Integer
  222    | Maximum_Scale !Integer
  223 
  224 data -- Has_Radix t =>
  225      Fixed_Frame t = Fixed_Frame !Result_Scale !Rounding_Mode
  226                                  !(Maybe (Integer, Integer))
  227 
  228 fit :: Has_Radix t =>
  229     Fixed_Frame t -> Fixed t -> Fixed t
  230 fit (Fixed_Frame rs round bounds) x@(FP m (Scale s))
  231   = check_bounds bounds m' y
  232   where
  233     y@(FP m' _) =
  234       case rs of
  235         Exact_Scale e   -> if e == s then x else
  236                            rescale round e x
  237         Maximum_Scale e -> if e >= s then x else
  238                            rescale round e x
  239 
  240 check_bounds :: Ord t => Maybe (t,t) -> t -> r -> r
  241 
  242 check_bounds Nothing      _ y = y
  243 check_bounds (Just (l,u)) x y =
  244   if l <= x && x <= u then y else error "overflow"
  245 
  246 divide :: forall t u v . (Has_Radix t, Real u, Real v) =>
  247           Fixed_Frame t -> u -> v -> Fixed t
  248 divide (Fixed_Frame rs round bounds) x y
  249   = check_bounds bounds m (FP m (Scale s))
  250   where
  251     s = case rs of
  252           (Exact_Scale   z) -> z
  253           (Maximum_Scale z) -> z
  254     p = radix (undefined :: t) ^ s
  255     q = toRational x / toRational y
  256     m = convert round (q * fromInteger p)
  257 
  258 readFixedDecimal :: ReadS (Fixed Decimal)
  259 
  260 readFixedDecimal ('-':cs) = after_sign negate cs
  261 readFixedDecimal cs       = after_sign id     cs
  262 
  263 after_sign :: (Integer -> Integer) -> ReadS (Fixed Decimal)
  264 
  265 after_sign f (c:cs)
  266   | isDigit c = after_digit f (add_digit 0 c) cs
  267 after_sign _ _ = []
  268 
  269 after_digit :: (Integer -> Integer) -> Integer ->
  270                ReadS (Fixed Decimal)
  271 
  272 after_digit f n (c:cs)
  273   | isDigit c = after_digit f (add_digit n c) cs
  274 after_digit f n ('.':c:cs)
  275   | isDigit c = after_dot f (Scale 1) (add_digit n c) cs
  276 after_digit f n ('.':_) = []
  277 after_digit f n cs = [(FP (f n) (Scale 0), cs)]
  278 
  279 after_dot :: (Integer -> Integer) -> Scale Decimal -> Integer ->
  280              ReadS (Fixed Decimal)
  281 
  282 after_dot f s n (c:cs)
  283   | isDigit c = after_dot f (s+1) (add_digit n c) cs
  284 after_dot f s n cs = [(FP (f n) s, cs)]
  285 
  286 add_digit :: Integer -> Char -> Integer
  287 
  288 add_digit n c =
  289    n * 10 + fromIntegral (fromEnum c - fromEnum '0')
  290 
  291 instance Read (Fixed Decimal)
  292   where
  293     readsPrec _ = readFixedDecimal
  294 
  295 -- ReadS a = String -> [(a,String)]
  296