diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,5 @@
+# Revision history for float128
+
+## 0.1  -- 2020-06-15
+
+* First release.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2020, Claude Heiland-Allen
+
+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 Claude Heiland-Allen 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/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/c/float128.c b/c/float128.c
new file mode 100644
--- /dev/null
+++ b/c/float128.c
@@ -0,0 +1,58 @@
+#define _GNU_SOURCE
+#include <math.h>
+#include <HsFFI.h>
+
+HsDouble f128_get_d(const _Float128 *a) { return *a; }
+HsInt    f128_get_i(const _Float128 *a) { return *a; }
+void     f128_set_d(_Float128 *r, HsDouble b) { *r = b; }
+void     f128_set_i(_Float128 *r, HsInt b) { *r = b; }
+
+void f128_add(_Float128 *r, const _Float128 *a, const _Float128 *b) { *r = *a + *b; }
+void f128_sub(_Float128 *r, const _Float128 *a, const _Float128 *b) { *r = *a - *b; }
+void f128_mul(_Float128 *r, const _Float128 *a, const _Float128 *b) { *r = *a * *b; }
+void f128_div(_Float128 *r, const _Float128 *a, const _Float128 *b) { *r = *a / *b; }
+void f128_pow(_Float128 *r, const _Float128 *a, const _Float128 *b) { *r = powf128(*a, *b); }
+void f128_atan2(_Float128 *r, const _Float128 *a, const _Float128 *b) { *r = atan2f128(*a, *b); }
+void f128_min(_Float128 *r, const _Float128 *a, const _Float128 *b) { *r = fminf128(*a, *b); }
+void f128_max(_Float128 *r, const _Float128 *a, const _Float128 *b) { *r = fmaxf128(*a, *b); }
+
+int f128_lt(const _Float128 *a, const _Float128 *b) { return *a <  *b; }
+int f128_le(const _Float128 *a, const _Float128 *b) { return *a <= *b; }
+int f128_eq(const _Float128 *a, const _Float128 *b) { return *a == *b; }
+int f128_ne(const _Float128 *a, const _Float128 *b) { return *a != *b; }
+int f128_ge(const _Float128 *a, const _Float128 *b) { return *a >= *b; }
+int f128_gt(const _Float128 *a, const _Float128 *b) { return *a >  *b; }
+
+void f128_abs(_Float128 *r, const _Float128 *a) { *r = fabsf128(*a); }
+void f128_sgn(_Float128 *r, const _Float128 *a) { *r = (*a > 0) - (0 > *a); }
+void f128_neg(_Float128 *r, const _Float128 *a) { *r = - *a; }
+void f128_recip(_Float128 *r, const _Float128 *a) { *r = 1.0L / *a; }
+void f128_sqrt(_Float128 *r, const _Float128 *a) { *r = sqrtf128(*a); }
+void f128_exp(_Float128 *r, const _Float128 *a) { *r = expf128(*a); }
+void f128_log(_Float128 *r, const _Float128 *a) { *r = logf128(*a); }
+void f128_sin(_Float128 *r, const _Float128 *a) { *r = sinf128(*a); }
+void f128_cos(_Float128 *r, const _Float128 *a) { *r = cosf128(*a); }
+void f128_tan(_Float128 *r, const _Float128 *a) { *r = tanf128(*a); }
+void f128_sinh(_Float128 *r, const _Float128 *a) { *r = sinhf128(*a); }
+void f128_cosh(_Float128 *r, const _Float128 *a) { *r = coshf128(*a); }
+void f128_tanh(_Float128 *r, const _Float128 *a) { *r = tanhf128(*a); }
+void f128_asin(_Float128 *r, const _Float128 *a) { *r = asinf128(*a); }
+void f128_acos(_Float128 *r, const _Float128 *a) { *r = acosf128(*a); }
+void f128_atan(_Float128 *r, const _Float128 *a) { *r = atanf128(*a); }
+void f128_asinh(_Float128 *r, const _Float128 *a) { *r = asinhf128(*a); }
+void f128_acosh(_Float128 *r, const _Float128 *a) { *r = acoshf128(*a); }
+void f128_atanh(_Float128 *r, const _Float128 *a) { *r = atanhf128(*a); }
+void f128_floor(_Float128 *r, const _Float128 *a) { *r = floorf128(*a); }
+void f128_ceil(_Float128 *r, const _Float128 *a) { *r = ceilf128(*a); }
+void f128_round(_Float128 *r, const _Float128 *a) { *r = roundf128(*a); }
+void f128_trunc(_Float128 *r, const _Float128 *a) { *r = truncf128(*a); }
+
+int f128_isnan(const _Float128 *a) { return fpclassify(*a) == FP_NAN; }
+int f128_isinf(const _Float128 *a) { return fpclassify(*a) == FP_INFINITE; }
+int f128_isdenorm(const _Float128 *a) { return fpclassify(*a) == FP_SUBNORMAL; }
+int f128_isnegzero(const _Float128 *a) { return fpclassify(*a) == FP_ZERO && signbit(*a); }
+
+void f128_ldexp(_Float128 *r, const _Float128 *a, int b) { *r = ldexpf128(*a, b); }
+void f128_frexp(_Float128 *r, const _Float128 *a, int *b) { *r = frexpf128(*a, b); }
+
+void f128_pi(_Float128 *r) { *r = M_PIf128; }
diff --git a/float128.cabal b/float128.cabal
new file mode 100644
--- /dev/null
+++ b/float128.cabal
@@ -0,0 +1,43 @@
+name:                float128
+version:             0.1
+synopsis:            FFI bindings for C _Float128
+description:
+  This package provides a Float128 type, wrapping C's _Float128, which
+  is typically a IEEE quadruple precision binary floating point type.
+  It does not provide a CFloat128 type usable for FFI without wrapping in Ptr,
+  this needs to be done by the compiler.
+
+homepage:            https://code.mathr.co.uk/float128
+license:             BSD3
+license-file:        LICENSE
+author:              Claude Heiland-Allen
+maintainer:          claude@mathr.co.uk
+copyright:           (c) 2020 Claude Heiland-Allen
+category:            Math
+build-type:          Simple
+extra-source-files:
+  ChangeLog.md
+
+cabal-version:       >=1.10
+
+source-repository head
+  type: git
+  location: https://code.mathr.co.uk/float128.git
+
+source-repository this
+  type: git
+  location: https://code.mathr.co.uk/float128.git
+  tag: float128-0.1
+
+library
+  exposed-modules:
+    Numeric.Float128
+  build-depends:
+    base >=4.6 && <4.15,
+    integer-gmp >=0.5 && <1.1
+  c-sources: c/float128.c
+  hs-source-dirs: hs
+  default-language: Haskell2010
+  other-extensions:
+    ForeignFunctionInterface,
+    MagicHash
diff --git a/hs/Numeric/Float128.hs b/hs/Numeric/Float128.hs
new file mode 100644
--- /dev/null
+++ b/hs/Numeric/Float128.hs
@@ -0,0 +1,333 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+{-# LANGUAGE MagicHash #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Float128.Float128
+-- Copyright   :  (C) 2020 Claude Heiland-Allen
+-- License     :  BSD3
+-- Maintainer  :  Claude Heiland-Allen <claude@mathr.co.uk>
+-- Stability   :  experimental
+-- Portability :  non-portable (needs C _Float128)
+--
+-- This module contains a Float128 type that can be used if you need the
+-- extra precision or range from the IEEE _Float128 quadruple precision type.
+-- It has 15bit signed exponent (compared to 11bit signed exponent for Double)
+-- and 113bit mantissa (compared to 53bit mantissa for Double).
+--
+-- Performance is likely to be poor, as the instances are implemented using
+-- FFI with Ptr Float128, copying to and from memory around each operation.
+-- If you need to bind to functions taking/returning _Float128 you need to
+-- write wrapper functions expecting pointers to _Float128 instead, as GHC
+-- does not expose a CFloat128 FFI type.
+module Numeric.Float128
+  ( 
+  -- * _Float128 data type
+    Float128(..)
+  -- * RealFrac alternatives
+  , truncate'
+  , round'
+  , ceiling'
+  , floor'
+  -- * Conversions
+  , fromDouble
+  , toDouble
+  , fromInt
+  , toInt
+  ) where
+
+import Data.Bits (bit, testBit, (.&.), (.|.), shiftL, shiftR)
+import Data.Ratio ((%), numerator, denominator)
+import Data.Word (Word64)
+import Foreign (Ptr, castPtr, with, alloca)
+import Foreign.C.Types (CIntMax(..), CInt(..), CDouble(..))
+import Foreign.Storable (Storable(..))
+import Numeric (showFloat, readFloat, readSigned)
+import System.IO.Unsafe (unsafePerformIO)
+import GHC.Exts (Int(..))
+import GHC.Integer.Logarithms (integerLog2#)
+
+-- | The _Float128 type.
+data Float128 = F128 !Word64 !Word64 -- most significant, least significant
+
+bigEndian :: Bool
+bigEndian = False -- FIXME
+
+instance Storable Float128 where
+  sizeOf _ = 2 * sizeOf (0 :: Word64)
+  alignment _ = 2 * alignment (0 :: Word64)
+  peek p = do
+    let q :: Ptr Word64
+        q = castPtr p
+    a <- peekElemOff q 0
+    b <- peekElemOff q 1
+    return $ if bigEndian then F128 a b else F128 b a
+  poke p (F128 msw lsw) = do
+    let q :: Ptr Word64
+        q = castPtr p
+        (a, b) = if bigEndian then (msw, lsw) else (lsw, msw)
+    pokeElemOff q 0 a
+    pokeElemOff q 1 b
+
+instance Eq Float128 where
+  (==) = cmp f128_eq
+  (/=) = cmp f128_ne
+
+instance Ord Float128 where
+  (<=) = cmp f128_le
+  (< ) = cmp f128_lt
+  (>=) = cmp f128_ge
+  (> ) = cmp f128_gt
+  min  = f2 f128_min
+  max  = f2 f128_max
+
+instance Num Float128 where
+  fromInteger z = encodeFloat z 0
+  negate = f1 f128_neg
+  (+) = f2 f128_add
+  (-) = f2 f128_sub
+  (*) = f2 f128_mul
+  abs = f1 f128_abs
+  signum = f1 f128_sgn
+
+instance Real Float128 where
+  toRational l = case decodeFloat l of
+    (m, e)
+      | e >= 0 -> m `shiftL` e % 1
+      | otherwise -> m % bit (negate e)
+
+instance Fractional Float128 where
+  fromRational q = -- FIXME accuracy?
+    let a = fromInteger (numerator q) / fromInteger (denominator q)
+        r = q - toRational a
+        b = fromInteger (numerator r) / fromInteger (denominator r)
+    in  a + b
+  (/) = f2 f128_div
+  recip = f1 f128_recip
+
+instance RealFrac Float128 where
+  properFraction l
+    | l >= 0 = let n' = floor' l
+                   f = l - n'
+               in  (fromInteger . toInteger' $ n', f)
+    
+    | l <  0 = let n' = ceiling' l
+                   f = l - n'
+               in  (fromInteger . toInteger' $ n', f)
+    | otherwise = (0, l) -- NaN
+  truncate = fromInteger . toInteger' . truncate'
+  round    = fromInteger . toInteger' . round'
+  ceiling  = fromInteger . toInteger' . ceiling'
+  floor    = fromInteger . toInteger' . floor'
+
+toInteger' :: Float128 -> Integer
+toInteger' l = case decodeFloat l of
+  (m, e)
+    | e >= 0 -> m `shiftL` e
+    | otherwise -> m `shiftR` negate e
+
+-- | Alternate versions of RealFrac methods that
+--   keep the value as a Float128.
+truncate', round', ceiling', floor' :: Float128 -> Float128
+truncate' = f1 f128_trunc
+round'    = f1 f128_round
+ceiling'  = f1 f128_ceil
+floor'    = f1 f128_floor
+
+instance Floating Float128 where
+  pi = unsafePerformIO $ do
+    alloca $ \lp -> do
+      f128_pi lp
+      peek lp
+  exp = f1 f128_exp
+  log = f1 f128_log
+  sqrt = f1 f128_sqrt
+  (**) = f2 f128_pow
+  -- logBase
+  sin = f1 f128_sin
+  cos = f1 f128_cos
+  tan = f1 f128_tan
+  sinh = f1 f128_sinh
+  cosh = f1 f128_cosh
+  tanh = f1 f128_tanh
+  asin = f1 f128_asin
+  acos = f1 f128_acos
+  atan = f1 f128_atan
+  asinh = f1 f128_asinh
+  acosh = f1 f128_acosh
+  atanh = f1 f128_atanh
+
+instance RealFloat Float128 where
+  floatRadix _ = 2
+  floatDigits _ = 113
+  floatRange _ = (-16381,16384) -- FIXME verify?
+
+  decodeFloat l@(F128 msw lsw)
+    | isNaN l = (0, 0)
+    | isInfinite l = (0, 0)
+    | l == 0 = (0, 0)
+    | isDenormalized l = case decodeFloat (scaleFloat 128 l) of
+        (m, e) -> (m, e - 128)
+    | otherwise =
+        ( (if s then negate else id) (shiftL (0x1000000000000 .|. toInteger msw .&. 0xFFFFffffFFFF) 64 .|. toInteger lsw)
+        , fromIntegral e0 - 16383 - 112 -- FIXME verify
+        )
+    where
+      s = shiftR msw 48 `testBit` 15
+      e0 = shiftR msw 48 .&. (bit 15 - 1)
+
+  encodeFloat m e
+    | m == 0 = F128 0 0
+    | m <  0 = negate (encodeFloat (negate m) e)
+    | b >= bit 15 - 1 = F128 (shiftL (bit 15 - 1) 48) 0 -- infinity
+    | b <= 0 = scaleFloat (b - 128) (encodeFloat m (e - b + 128)) -- denormal
+    | otherwise = F128 msw lsw -- normal
+    where
+      bigEndian = False -- FIXME
+      l = I# (integerLog2# m)
+      t = l - 112 -- FIXME verify
+      m' | t >= 0    = m `shiftR`        t
+         | otherwise = m `shiftL` negate t
+      -- FIXME: verify that m' `testBit` 112 == True
+      lsw = fromInteger (m' .&. 0xFFFFffffFFFFffff)
+      msw = fromInteger (shiftR m' 64 .&. 0xFFFFffffFFFF) .|. shiftL (fromIntegral b) 48 
+      b = e + t + 16383 + 112 -- FIXME verify
+
+  exponent l@(F128 msw lsw)
+    | isNaN l = 0
+    | isInfinite l = 0
+    | l == 0 = 0
+    | isDenormalized l = snd (decodeFloat l) + 113
+    | otherwise = fromIntegral e0 - 16383 - 112 + 113
+    where
+      e0 = shiftR msw 48 .&. (bit 15 - 1)
+
+  significand l = unsafePerformIO $ do
+    with l $ \lp -> alloca $ \ep -> do
+      f128_frexp lp lp ep
+      peek lp
+
+  scaleFloat e l = unsafePerformIO $ do
+    with l $ \lp -> do
+      f128_ldexp lp lp (fromIntegral e)
+      peek lp
+
+  isNaN = tst f128_isnan
+  isInfinite = tst f128_isinf
+  isDenormalized = tst f128_isdenorm
+  isNegativeZero = tst f128_isnegzero
+  isIEEE _= True
+
+  atan2 = f2 f128_atan2
+
+instance Read Float128 where
+  readsPrec _ = readSigned readFloat
+
+instance Show Float128 where
+  showsPrec p x = showParen (p >= 7 && take 1 s == "-") (s ++) -- FIXME: precedence issues?
+    where s = showFloat x ""
+
+fromInt :: Int -> Float128
+fromInt i = unsafePerformIO $ do
+  alloca $ \lp -> do
+    f128_set_i lp (fromIntegral i)
+    peek lp
+
+toInt :: Float128 -> Int
+toInt l = unsafePerformIO $ with l f128_get_i
+
+fromDouble :: Double -> Float128
+fromDouble i = unsafePerformIO $ do
+  alloca $ \lp -> do
+    f128_set_d lp i
+    peek lp
+
+toDouble :: Float128 -> Double
+toDouble l = unsafePerformIO $ with l f128_get_d
+
+f2 :: F2 -> Float128 -> Float128 -> Float128
+f2 f a b = unsafePerformIO $ do
+  with a $ \ap -> with b $ \bp -> alloca $ \rp -> do
+    f rp ap bp
+    peek rp
+
+type F2 = Ptr Float128 -> Ptr Float128 -> Ptr Float128 -> IO ()
+
+foreign import ccall unsafe "f128_add"   f128_add   :: F2
+foreign import ccall unsafe "f128_sub"   f128_sub   :: F2
+foreign import ccall unsafe "f128_mul"   f128_mul   :: F2
+foreign import ccall unsafe "f128_div"   f128_div   :: F2
+foreign import ccall unsafe "f128_pow"   f128_pow   :: F2
+foreign import ccall unsafe "f128_min"   f128_min   :: F2
+foreign import ccall unsafe "f128_max"   f128_max   :: F2
+foreign import ccall unsafe "f128_atan2" f128_atan2 :: F2
+
+f1 :: F1 -> Float128 -> Float128
+f1 f a = unsafePerformIO $ do
+  with a $ \ap -> alloca $ \rp -> do
+    f rp ap
+    peek rp
+
+type F1 = Ptr Float128 -> Ptr Float128 -> IO ()
+
+foreign import ccall unsafe "f128_abs"   f128_abs   :: F1
+foreign import ccall unsafe "f128_sgn"   f128_sgn   :: F1
+foreign import ccall unsafe "f128_neg"   f128_neg   :: F1
+foreign import ccall unsafe "f128_sqrt"  f128_sqrt  :: F1
+foreign import ccall unsafe "f128_recip" f128_recip :: F1
+foreign import ccall unsafe "f128_exp"   f128_exp   :: F1
+foreign import ccall unsafe "f128_log"   f128_log   :: F1
+foreign import ccall unsafe "f128_sin"   f128_sin   :: F1
+foreign import ccall unsafe "f128_cos"   f128_cos   :: F1
+foreign import ccall unsafe "f128_tan"   f128_tan   :: F1
+foreign import ccall unsafe "f128_sinh"  f128_sinh  :: F1
+foreign import ccall unsafe "f128_cosh"  f128_cosh  :: F1
+foreign import ccall unsafe "f128_tanh"  f128_tanh  :: F1
+foreign import ccall unsafe "f128_asin"  f128_asin  :: F1
+foreign import ccall unsafe "f128_acos"  f128_acos  :: F1
+foreign import ccall unsafe "f128_atan"  f128_atan  :: F1
+foreign import ccall unsafe "f128_asinh" f128_asinh :: F1
+foreign import ccall unsafe "f128_acosh" f128_acosh :: F1
+foreign import ccall unsafe "f128_atanh" f128_atanh :: F1
+foreign import ccall unsafe "f128_floor" f128_floor :: F1
+foreign import ccall unsafe "f128_ceil"  f128_ceil  :: F1
+foreign import ccall unsafe "f128_round" f128_round :: F1
+foreign import ccall unsafe "f128_trunc" f128_trunc :: F1
+
+type CMP = Ptr Float128 -> Ptr Float128 -> IO CInt
+
+cmp :: CMP -> Float128 -> Float128 -> Bool
+cmp f a b = unsafePerformIO $ do
+  with a $ \ap -> with b $ \bp -> do
+    r <- f ap bp
+    return (r /= 0)
+
+foreign import ccall unsafe "f128_eq" f128_eq :: CMP
+foreign import ccall unsafe "f128_ne" f128_ne :: CMP
+foreign import ccall unsafe "f128_lt" f128_lt :: CMP
+foreign import ccall unsafe "f128_le" f128_le :: CMP
+foreign import ccall unsafe "f128_gt" f128_gt :: CMP
+foreign import ccall unsafe "f128_ge" f128_ge :: CMP
+
+type TST = Ptr Float128 -> IO CInt
+
+tst :: TST -> Float128 -> Bool
+tst f a = unsafePerformIO $ do
+  with a $ \ap -> do
+    r <- f ap
+    return (r /= 0)
+
+foreign import ccall unsafe "f128_isnan" f128_isnan :: TST
+foreign import ccall unsafe "f128_isinf" f128_isinf :: TST
+foreign import ccall unsafe "f128_isdenorm" f128_isdenorm :: TST
+foreign import ccall unsafe "f128_isnegzero" f128_isnegzero :: TST
+
+foreign import ccall unsafe "f128_get_d" f128_get_d :: Ptr Float128 -> IO Double
+foreign import ccall unsafe "f128_get_i" f128_get_i :: Ptr Float128 -> IO Int
+
+foreign import ccall unsafe "f128_set_d" f128_set_d :: Ptr Float128 -> Double -> IO ()
+foreign import ccall unsafe "f128_set_i" f128_set_i :: Ptr Float128 -> Int -> IO ()
+
+foreign import ccall unsafe "f128_ldexp" f128_ldexp :: Ptr Float128 -> Ptr Float128 -> CInt -> IO ()
+foreign import ccall unsafe "f128_frexp" f128_frexp :: Ptr Float128 -> Ptr Float128 -> Ptr CInt -> IO ()
+
+foreign import ccall unsafe "f128_pi"    f128_pi    :: Ptr Float128 -> IO ()