diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/Numeric/Rounding.hsc b/Numeric/Rounding.hsc
new file mode 100644
--- /dev/null
+++ b/Numeric/Rounding.hsc
@@ -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)
+
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMainWithHooks simpleUserHooks
diff --git a/rounding-inc.h b/rounding-inc.h
new file mode 100644
--- /dev/null
+++ b/rounding-inc.h
@@ -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
diff --git a/rounding.c b/rounding.c
new file mode 100644
--- /dev/null
+++ b/rounding.c
@@ -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"
+
diff --git a/rounding.cabal b/rounding.cabal
new file mode 100644
--- /dev/null
+++ b/rounding.cabal
@@ -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
diff --git a/rounding.h b/rounding.h
new file mode 100644
--- /dev/null
+++ b/rounding.h
@@ -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