packages feed

rounding (empty) → 0.0.1

raw patch · 7 files changed

+589/−0 lines, 7 filesdep +arraydep +basesetup-changed

Dependencies added: array, base

Files

+ LICENSE view
@@ -0,0 +1,31 @@+Copyright (c) 2010, Edward Kmett++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are+met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Edward Kmett nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Numeric/Rounding.hsc view
@@ -0,0 +1,429 @@+{-# LANGUAGE CPP, ForeignFunctionInterface, EmptyDataDecls, FlexibleInstances, FlexibleContexts #-}+module Numeric.Rounding +    ( Round(..)+    , Rounding(..)+    , Up, Down, ToNearest, TowardZero+    , up, down, toNearest, towardZero+    , runUp, runDown, runToNearest, runTowardZero+    , double+    , float+    ) where++import Control.Applicative+import GHC.Real+import Data.Foldable+import Data.Traversable+import Data.Array+import Foreign+import Foreign.C.Types++#include <math.h>+#include <fenv.h>++-- TODO: tweak the lsbs of pi+-- modify the enum instance to properly round+-- implement complex numbers++newtype Round dir a = Round a+    deriving (Show, Eq, Ord, Bounded)++instance Functor (Round dir) where+    fmap f (Round a) = Round (f a)++instance Foldable (Round dir) where+    foldMap f (Round a) = f a++instance Traversable (Round dir) where+    traverse f (Round a) = Round <$> f a++class Rounding dir where+    mode :: Round dir a -> CInt+    rounding :: (Integral b, RealFrac a) => Round dir c -> (a -> b)++data ToNearest+data TowardZero+data Up+data Down++up :: a -> Round Up a +up = Round+{-# INLINE up #-}++down :: a -> Round Down a+down = Round+{-# INLINE down #-}++towardZero :: a -> Round TowardZero a +towardZero = Round+{-# INLINE towardZero #-}++toNearest :: a -> Round ToNearest a+toNearest = Round+{-# INLINE toNearest #-}++runUp :: Round Up a -> a+runUp (Round a) = a+{-# INLINE runUp #-}++runDown :: Round Down a -> a+runDown (Round a) = a+{-# INLINE runDown #-}++runTowardZero :: Round TowardZero a -> a+runTowardZero (Round a) = a+{-# INLINE runTowardZero #-}++runToNearest :: Round ToNearest a -> a+runToNearest (Round a) = a+{-# INLINE runToNearest #-}++double :: Double -> Double+double = id+{-# INLINE double #-}++float :: Float -> Float+float = id+{-# INLINE float #-}+++instance Rounding ToNearest where +    mode _ = #const FE_TONEAREST+    rounding _ = round++instance Rounding TowardZero where +    mode _ = #const FE_TOWARDZERO+    rounding _ = truncate++instance Rounding Up where +    mode _ = #const FE_UPWARD+    rounding _ = ceiling++instance Rounding Down where +    mode _ = #const FE_DOWNWARD+    rounding _ = floor++-- * Rounded Doubles++lift1D :: Rounding m => +          (CInt -> CDouble -> CDouble) -> +          Round m Double -> Round m Double+lift1D f r@(Round x) = Round (realToFrac (f (mode r) (realToFrac x)))++lift2D :: Rounding m => +          (CInt -> CDouble -> CDouble -> CDouble) -> +          Round m Double -> Round m Double -> Round m Double+lift2D f r@(Round x) (Round y) = Round (realToFrac (f (mode r) (realToFrac x) (realToFrac y)))++foreign import ccall unsafe "rounding.h madd" +    madd :: CInt -> CDouble -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mminus" +    mminus :: CInt -> CDouble -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mtimes"+    mtimes :: CInt -> CDouble -> CDouble -> CDouble +foreign import ccall unsafe "rounding.h mdiv"+    mdiv :: CInt -> CDouble -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mexp"+    mexp :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mpow"+    mpow :: CInt -> CDouble -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mlog"+    mlog :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h msqrt"+    msqrt :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h msin"+    msin :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mcos"+    mcos :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mtan"+    mtan :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h msinh"+    msinh :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mcosh"+    mcosh :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h mtanh"+    mtanh :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h masin"+    masin :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h macos"+    macos :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h matan"+    matan :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h masinh"+    masinh :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h macosh"+    macosh :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h matanh"+    matanh :: CInt -> CDouble -> CDouble+foreign import ccall unsafe "rounding.h matan2"+    matan2 :: CInt -> CDouble -> CDouble -> CDouble++instance Rounding dir => Num (Round dir Double) where+    fromInteger n = Round (fromInteger n)+    (+) = lift2D madd+    (-) = lift2D mminus+    (*) = lift2D mtimes+    abs (Round a) = Round (abs a)+    signum (Round a) = Round (signum a)++instance Rounding dir => Fractional (Round dir Double) where+    (/) = lift2D mdiv+    recip = lift2D mdiv 1+    fromRational = fromRat+    +instance Rounding dir => Enum (Round dir Double) where+    succ = (+1)+    pred = subtract 1+    toEnum n = Round (toEnum n)+    fromEnum (Round a) = fromEnum a+    enumFrom = numericEnumFrom+    enumFromThen = numericEnumFromThen+    enumFromTo = numericEnumFromTo+    enumFromThenTo = numericEnumFromThenTo+    +instance Rounding dir => Floating (Round dir Double) where+    pi = Round pi -- TODO twiddle LSB+    exp = lift1D mexp+    (**) = lift2D mpow+    log = lift1D mlog+    sqrt = lift1D msqrt+    sin = lift1D msin+    cos = lift1D mcos+    tan = lift1D mtan+    asin = lift1D masin+    acos = lift1D macos+    atan = lift1D matan+    sinh = lift1D msinh+    cosh = lift1D mcosh+    tanh = lift1D mtanh+    asinh = lift1D masinh+    acosh = lift1D macosh+    atanh = lift1D matanh++instance Rounding dir => Real (Round dir Double) where+    toRational (Round a) = toRational a -- tweak?++instance Rounding dir => RealFrac (Round dir Double) where+    properFraction = properFrac+    truncate (Round a) = truncate a+    round (Round a) = round a+    ceiling (Round a) = ceiling a+    floor (Round a) = floor a++instance Rounding dir => RealFloat (Round dir Double) where+    floatRadix (Round a) = floatRadix a+    floatDigits (Round a) = floatDigits a+    floatRange (Round a) = floatRange a+    decodeFloat (Round a) = decodeFloat a+    encodeFloat m e = Round (encodeFloat m e)+    exponent (Round a) = exponent a+    significand (Round a) = Round (significand a)+    scaleFloat n (Round a) = Round (scaleFloat n a)+    isNaN (Round a) = isNaN a+    isInfinite (Round a) = isInfinite a+    isDenormalized (Round a) = isDenormalized a+    isNegativeZero (Round a) = isNegativeZero a+    isIEEE (Round a) = isIEEE a+    atan2 = lift2D matan2++-- * Rounded Floats++lift1F :: Rounding m => +          (CInt -> CFloat -> CFloat ) -> +          Round m Float -> Round m Float+lift1F f r@(Round x) = Round (realToFrac (f (mode r) (realToFrac x)))++lift2F :: Rounding m => +          (CInt -> CFloat -> CFloat -> CFloat) -> +          Round m Float -> Round m Float -> Round m Float+lift2F f r@(Round x) (Round y) = Round (realToFrac (f (mode r) (realToFrac x) (realToFrac y)))++foreign import ccall unsafe "rounding.h madd" +    maddf :: CInt -> CFloat -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mminus" +    mminusf :: CInt -> CFloat -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mtimesf"+    mtimesf :: CInt -> CFloat -> CFloat -> CFloat +foreign import ccall unsafe "rounding.h mdivf"+    mdivf :: CInt -> CFloat -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mexpf"+    mexpf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mpowf"+    mpowf :: CInt -> CFloat -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mlogf"+    mlogf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h msqrtf"+    msqrtf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h msinf"+    msinf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mcosf"+    mcosf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mtanf"+    mtanf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h msinhf"+    msinhf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mcoshf"+    mcoshf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h mtanhf"+    mtanhf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h masinf"+    masinf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h macosf"+    macosf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h matanf"+    matanf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h masinhf"+    masinhf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h macoshf"+    macoshf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h matanhf"+    matanhf :: CInt -> CFloat -> CFloat+foreign import ccall unsafe "rounding.h matan2f"+    matan2f :: CInt -> CFloat -> CFloat -> CFloat++instance Rounding dir => Num (Round dir Float) where+    fromInteger n = Round (fromInteger n)+    (+) = lift2F maddf+    (-) = lift2F mminusf+    (*) = lift2F mtimesf+    abs (Round a) = Round (abs a)+    signum (Round a) = Round (signum a)++instance Rounding dir => Fractional (Round dir Float) where+    (/) = lift2F mdivf+    recip = lift2F mdivf 1+    fromRational = fromRat+    +instance Rounding dir => Enum (Round dir Float) where+    succ = (+1)+    pred = subtract 1+    toEnum n = Round (toEnum n)+    fromEnum (Round a) = fromEnum a+    enumFrom = numericEnumFrom+    enumFromThen = numericEnumFromThen+    enumFromTo = numericEnumFromTo+    enumFromThenTo = numericEnumFromThenTo+    +instance Rounding dir => Floating (Round dir Float) where+    pi = Round pi -- TODO twiddle LSB+    exp = lift1F mexpf+    (**) = lift2F mpowf+    log = lift1F mlogf+    sqrt = lift1F msqrtf+    sin = lift1F msinf+    cos = lift1F mcosf+    tan = lift1F mtanf+    asin = lift1F masinf+    acos = lift1F macosf+    atan = lift1F matanf+    sinh = lift1F msinhf+    cosh = lift1F mcoshf+    tanh = lift1F mtanhf+    asinh = lift1F masinhf+    acosh = lift1F macoshf+    atanh = lift1F matanhf++instance Rounding dir => Real (Round dir Float) where+    toRational (Round a) = toRational a -- tweak?++instance Rounding dir => RealFrac (Round dir Float) where+    properFraction = properFrac+    truncate (Round a) = truncate a+    round (Round a) = round a+    ceiling (Round a) = ceiling a+    floor (Round a) = floor a++instance Rounding dir => RealFloat (Round dir Float) where+    floatRadix (Round a) = floatRadix a+    floatDigits (Round a) = floatDigits a+    floatRange (Round a) = floatRange a+    decodeFloat (Round a) = decodeFloat a+    encodeFloat m e = Round (encodeFloat m e)+    exponent (Round a) = exponent a+    significand (Round a) = Round (significand a)+    scaleFloat n (Round a) = Round (scaleFloat n a)+    isNaN (Round a) = isNaN a+    isInfinite (Round a) = isInfinite a+    isDenormalized (Round a) = isDenormalized a+    isNegativeZero (Round a) = isNegativeZero a+    isIEEE (Round a) = isIEEE a+    atan2 = lift2F matan2f++-- * Fractional ++properFrac +    :: (Rounding dir, RealFrac a, Integral b) => +       Round dir a -> (b, Round dir a)+properFrac (Round a) = (b, Round c) +    where+    (b, c) = properFraction a++-- * Rounding Rationals++fromRat :: (Rounding dir, RealFloat a, RealFloat (Round dir a)) => Rational -> Round dir a+fromRat (n :% 0) | n > 0  =  1/0 -- +Infinity+                 | n == 0 =  0/0 -- NaN+                 | n < 0  = -1/0 -- -Infinity+fromRat (n :% d) | n > 0  = fromRat' (n :% d)+                 | n == 0 = encodeFloat 0 0 -- Zero+                 | n < 0  = - fromRat' ((-n) :% d)+++-- Conversion process:+-- Scale the rational number by the RealFloat base until+-- it lies in the range of the mantissa (as used by decodeFloat/encodeFloat).+-- Then round the rational to an Integer and encode it with the exponent+-- that we got from the scaling.+-- To speed up the scaling process we compute the log2 of the number to get+-- a first guess of the exponent.++fromRat' :: (Rounding dir, RealFloat a, RealFloat (Round dir a)) => Rational -> Round dir a+-- Invariant: argument is strictly positive+fromRat' x = r+  where b = floatRadix r+        p = floatDigits r+        (minExp0, _) = floatRange r+        minExp = minExp0 - p            -- the real minimum exponent+        xMin   = toRational (expt b (p-1))+        xMax   = toRational (expt b p)+        p0 = (integerLogBase b (numerator x) - integerLogBase b (denominator x) - p) `max` minExp+        f = if p0 < 0 then 1 % expt b (-p0) else expt b p0 % 1+        (x', p') = scaleRat (toRational b) minExp xMin xMax p0 (x / f)+        r = encodeFloat (rounding r x') p'++-- Scale x until xMin <= x < xMax, or p (the exponent) <= minExp.+scaleRat :: Rational -> Int -> Rational -> Rational -> Int -> Rational -> (Rational, Int)+scaleRat b minExp xMin xMax p x + | p <= minExp = (x, p)+ | x >= xMax   = scaleRat b minExp xMin xMax (p+1) (x/b)+ | x < xMin    = scaleRat b minExp xMin xMax (p-1) (x*b)+ | otherwise   = (x, p)+++-- Exponentiation with a cache for the most common numbers.+minExpt, maxExpt :: Int+minExpt = 0+maxExpt = 1100++expt :: Integer -> Int -> Integer+expt base n | base == 2 && n >= minExpt && n <= maxExpt = expts ! n+            | otherwise = base^n++expts :: Array Int Integer+expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]++-- Compute the (floor of the) log of i in base b.+-- Simplest way would be just divide i by b until it's smaller then b, but that would+-- be very slow!  We are just slightly more clever.+integerLogBase :: Integer -> Integer -> Int+integerLogBase b i+   | i < b     = 0+   | otherwise = doDiv (i `div` (b^l)) l+       where+        -- Try squaring the base first to cut down the number of divisions.+         l = 2 * integerLogBase (b*b) i++         doDiv :: Integer -> Int -> Int+         doDiv x y+            | x < b     = y+            | otherwise = doDiv (x `div` b) (y+1)+
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMainWithHooks simpleUserHooks
+ rounding-inc.h view
@@ -0,0 +1,55 @@+// fool GHC that this really is a valid header+#ifdef UNARY_T++#define UNARY(fun) \+    UNARY_T(double, m ## fun, fun) \+    UNARY_T(float , m ## fun ## f, fun) \+    UNARY_T(long double , m ## fun ## l, fun)++#define BINARY(fun) \+    BINARY_T(double, m ## fun, fun) \+    BINARY_T(float , m ## fun ## f, fun) \+    BINARY_T(long double , m ## fun ## l, fun)++#define INFIX(fun, op) \+    INFIX_T(double, m ## fun, op) \+    INFIX_T(float , m ## fun ## f, op) \+    INFIX_T(long double , m ## fun ## l, op)++UNARY(fabs)+UNARY(exp)+// UNARY(expm1)+// UNARY(exp2)+UNARY(log)+// UNARY(logp1)+UNARY(sqrt)+UNARY(sin)+UNARY(cos)+UNARY(tan)+UNARY(asin)+UNARY(acos)+UNARY(atan)+UNARY(sinh)+UNARY(cosh)+UNARY(tanh)+UNARY(asinh)+UNARY(acosh)+UNARY(atanh)+UNARY(ceil)+UNARY(floor)+UNARY(round)+UNARY(trunc)+BINARY(atan2)+// BINARY(hypot)+// BINARY(cbrt)+BINARY(pow)+INFIX(add, (x + y))+INFIX(times, (x * y))+INFIX(div, (x / y))+INFIX(minus, (x - y))++#undef UNARY+#undef BINARY+#undef INFIX++#endif
+ rounding.c view
@@ -0,0 +1,28 @@+#include "rounding.h"++#define UNARY_T(t,f,g) \+    t f (int mode, t x) { \+        fesetround(mode); \+        t result = g(x); \+        fesetround(FE_DEFAULT); \+        return result; \+    }+        +#define BINARY_T(t,f,g) \+    t f (int mode, t x, t y) { \+        fesetround(mode); \+        t result = g(x, y); \+        fesetround(FE_DEFAULT); \+        return result; \+    }+        +#define INFIX_T(t,f,op) \+    t f (int mode, t x, t y) { \+        fesetround(mode); \+        t result = (op); \+        fesetround(FE_DEFAULT); \+        return result; \+    }+        +#include "rounding-inc.h"+
+ rounding.cabal view
@@ -0,0 +1,18 @@+Name:              rounding+Version:           0.0.1+Synopsis:          Explicit floating point rounding mode wrappers+Homepage:          http://patch-tag.com/r/ekmett/rounding+License:           BSD3+License-file:      LICENSE+Author:            Edward Kmett+Maintainer:        ekmett@gmail.com+Category:          Math+Build-type:        Simple+Cabal-version:     >=1.6+Extra-Source-Files: rounding.h, rounding-inc.h++Library+  Exposed-modules: Numeric.Rounding+  Build-depends:   base >= 4 && < 5,+                   array >= 0.3.0 && < 0.4+  C-Sources:       rounding.c
+ rounding.h view
@@ -0,0 +1,25 @@+#ifndef INCLUDED_HASKELL_ROUNDING_H+#define INCLUDED_HASKELL_ROUNDING_H++#include <math.h>+#include <fenv.h>++#define FE_DEFAULT FE_TONEAREST++#define UNARY_T(t,f,g) \+    extern t f (int mode, t x);++#define BINARY_T(t,f,g) \+    extern t f (int mode, t x, t y);++#define INFIX_T(t,f,op) \+    extern t f (int mode, t x, t y);++// invoke the rounding "X-Macro"+#include "rounding-inc.h"++#undef UNARY_T+#undef BINARY_T+#undef INFIX_T++#endif