diff --git a/.ghci b/.ghci new file mode 100644 --- /dev/null +++ b/.ghci @@ -0,0 +1,2 @@ +:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h -optP-Iincludes +:set -v0 diff --git a/.gitignore b/.gitignore new file mode 100644 --- /dev/null +++ b/.gitignore @@ -0,0 +1,35 @@ +dist +docs +wiki +TAGS +tags +wip +contrib +.DS_Store +.*.swp +.*.swo +*.o +*.hi +.svn +*~ +*# +Makefile +Makefile.in +aclocal.m4 +analytics.buildinfo +autom4te.cache +autoscan.log +config.log +config.guess +config.status +config.sub +config.h +config.h.in +configure +configure.scan +depcomp +install-sh +missing +www +stamp-h1 +config.h.in diff --git a/.travis.yml b/.travis.yml new file mode 100644 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,26 @@ +language: haskell +before_install: + # Uncomment whenever hackage is down. + # - mkdir -p ~/.cabal && cp travis/config ~/.cabal/config && cabal update + - cabal update + + # Try installing some of the build-deps with apt-get for speed. + - travis/cabal-apt-install $mode + +install: + - cabal configure -flib-Werror $mode + - cabal build + +script: + - $script && hlint src --cpp-define HLINT + +notifications: + irc: + channels: + - "irc.freenode.org#haskell-lens" + skip_join: true + template: + - "\x0313compensated\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}" + +env: + - mode="--enable-tests" script="cabal test --show-details=always" diff --git a/.vim.custom b/.vim.custom new file mode 100644 --- /dev/null +++ b/.vim.custom @@ -0,0 +1,31 @@ +" Add the following to your .vimrc to automatically load this on startup + +" if filereadable(".vim.custom") +" so .vim.custom +" endif + +function StripTrailingWhitespace() + let myline=line(".") + let mycolumn = col(".") + silent %s/ *$// + call cursor(myline, mycolumn) +endfunction + +" enable syntax highlighting +syntax on + +" search for the tags file anywhere between here and / +set tags=TAGS;/ + +" highlight tabs and trailing spaces +set listchars=tab:‗‗,trail:‗ +set list + +" f2 runs hasktags +map :exec ":!hasktags -x -c --ignore src" + +" strip trailing whitespace before saving +" au BufWritePre *.hs,*.markdown silent! cal StripTrailingWhitespace() + +" rebuild hasktags after saving +au BufWritePost *.hs silent! :exec ":!hasktags -x -c --ignore src" diff --git a/AUTHORS.markdown b/AUTHORS.markdown new file mode 100644 --- /dev/null +++ b/AUTHORS.markdown @@ -0,0 +1,15 @@ +Analytics was started by [Edward Kmett](https://github.com/ekmett) in response to a question by [Alec Heller](https://github.com/deviant-logic) about if he should use `bound` to implement datalog. It has since somewhat expanded in scope. + +`compensated` was split out of the `analytics` repository. + +Improvements to the multiplication algorithm used in this package and many small tweaks like the ability to iterate the `Compensated` construction were made after prolonged discussions on the topic with Dan Doel. Any errors in the algorithm are my fault, however. + +* [Dan Doel](mailto:dan.doel@gmail.com) [@dolio](https://github.com/dolio) + +You can watch contributors carry on the quest for bragging rights in the [contributors graph](https://github.com/analytics/compensated/graphs/contributors). + +Omission from this list is by no means an attempt to discount your contribution. + +Thank you for all of your help! + +-Edward Kmett diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown new file mode 100644 --- /dev/null +++ b/CHANGELOG.markdown @@ -0,0 +1,3 @@ +0.1 +--- +* Split from [analytics](https://github.com/analytics/analytics) diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright 2013 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: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. 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. + +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS 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/README.markdown b/README.markdown new file mode 100644 --- /dev/null +++ b/README.markdown @@ -0,0 +1,23 @@ +compensated +=========== + +[![Build Status](https://secure.travis-ci.org/analytics/compensated.png)](http://travis-ci.org/analytics/compensated) + +This package provides a form of compensated arithmetic that lets you perform many operations on chains of doubles at higher precision. + +In particular addition, multiplication, etc. on a `Compensated Double` is done with twice as many bits of `significand` as on a raw `Double`. + +This construction is designed so that it can be iterated. + +This was forked from the [analytics](http://github.com/analytics/analytics.git) repository after repeated requests. + +Work still needs to be done to properly support most `Floating` operations, but everything below that level of the numeric hierarchy should work. + +Contact Information +------------------- + +Contributions and bug reports are welcome! + +Please feel free to contact me through github or on the #haskell IRC channel on irc.freenode.net. + +-Edward Kmett diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,55 @@ +#!/usr/bin/runhaskell +\begin{code} +{-# OPTIONS_GHC -Wall #-} +module Main (main) where + +import Data.List ( nub ) +import Data.Version ( showVersion ) +import Distribution.Package ( PackageName(PackageName), Package, PackageId, InstalledPackageId, packageVersion, packageName ) +import Distribution.PackageDescription ( PackageDescription(), TestSuite(..) ) +import Distribution.Simple ( defaultMainWithHooks, UserHooks(..), simpleUserHooks ) +import Distribution.Simple.Utils ( rewriteFile, createDirectoryIfMissingVerbose, copyFiles ) +import Distribution.Simple.BuildPaths ( autogenModulesDir ) +import Distribution.Simple.Setup ( BuildFlags(buildVerbosity), Flag(..), fromFlag, HaddockFlags(haddockDistPref)) +import Distribution.Simple.LocalBuildInfo ( withLibLBI, withTestLBI, LocalBuildInfo(), ComponentLocalBuildInfo(componentPackageDeps) ) +import Distribution.Text ( display ) +import Distribution.Verbosity ( Verbosity, normal ) +import System.FilePath ( () ) + +main :: IO () +main = defaultMainWithHooks simpleUserHooks + { buildHook = \pkg lbi hooks flags -> do + generateBuildModule (fromFlag (buildVerbosity flags)) pkg lbi + buildHook simpleUserHooks pkg lbi hooks flags + , postHaddock = \args flags pkg lbi -> do + copyFiles normal (haddockOutputDir flags pkg) [("images","Hierarchy.png")] + postHaddock simpleUserHooks args flags pkg lbi + } + +haddockOutputDir :: Package p => HaddockFlags -> p -> FilePath +haddockOutputDir flags pkg = destDir where + baseDir = case haddockDistPref flags of + NoFlag -> "." + Flag x -> x + destDir = baseDir "doc" "html" display (packageName pkg) + +generateBuildModule :: Verbosity -> PackageDescription -> LocalBuildInfo -> IO () +generateBuildModule verbosity pkg lbi = do + let dir = autogenModulesDir lbi + createDirectoryIfMissingVerbose verbosity True dir + withLibLBI pkg lbi $ \_ libcfg -> do + withTestLBI pkg lbi $ \suite suitecfg -> do + rewriteFile (dir "Build_" ++ testName suite ++ ".hs") $ unlines + [ "module Build_" ++ testName suite ++ " where" + , "deps :: [String]" + , "deps = " ++ (show $ formatdeps (testDeps libcfg suitecfg)) + ] + where + formatdeps = map (formatone . snd) + formatone p = case packageName p of + PackageName n -> n ++ "-" ++ showVersion (packageVersion p) + +testDeps :: ComponentLocalBuildInfo -> ComponentLocalBuildInfo -> [(InstalledPackageId, PackageId)] +testDeps xs ys = nub $ componentPackageDeps xs ++ componentPackageDeps ys + +\end{code} diff --git a/compensated.cabal b/compensated.cabal new file mode 100644 --- /dev/null +++ b/compensated.cabal @@ -0,0 +1,88 @@ +name: compensated +category: Numeric +version: 0.1 +license: BSD3 +cabal-version: >= 1.8 +license-file: LICENSE +author: Edward A. Kmett +maintainer: Edward A. Kmett +stability: provisional +homepage: http://github.com/analytics/compensated/ +bug-reports: http://github.com/analytics/compensated/issues +copyright: Copyright (C) 2013 Edward A. Kmett +build-type: Custom +tested-with: GHC == 7.4.1, GHC == 7.6.1 +synopsis: Compensated floating-point arithmetic +description: This package provides compensated floating point arithmetic. + +extra-source-files: + .travis.yml + .ghci + .gitignore + .vim.custom + travis/cabal-apt-install + travis/config + AUTHORS.markdown + README.markdown + CHANGELOG.markdown + +source-repository head + type: git + location: git://github.com/analytics/compensated.git + +-- You can disable the doctests test suite with -f-test-doctests +flag test-doctests + default: True + manual: True + +flag lib-Werror + default: False + manual: True + +library + build-depends: + base >= 4.3 && < 5, + bifunctors >= 3.2 && < 4, + binary >= 0.5 && < 0.6, + cereal >= 0.3.5 && < 0.4, + comonad >= 3 && < 4, + deepseq >= 1.3 && < 1.5, + distributive >= 0.3 && < 1, + generic-deriving == 1.4.*, + hashable >= 1.1.2.3 && < 1.3, + lens >= 3.8 && < 4, + log-domain >= 0.3 && < 1, + semigroupoids >= 3.0.2 && < 4, + semigroups >= 0.8.4 && < 1, + safecopy >= 0.8.1 && < 0.9, + vector >= 0.9 && < 0.11 + + exposed-modules: + Numeric.Compensated + + if flag(lib-Werror) + ghc-options: -Werror + + ghc-options: -Wall -fwarn-tabs -O2 + hs-source-dirs: src + +test-suite doctests + type: exitcode-stdio-1.0 + main-is: doctests.hs + ghc-options: -Wall -threaded + hs-source-dirs: tests + + if !flag(test-doctests) + buildable: False + else + build-depends: + base, + directory >= 1.0, + doctest >= 0.9.1, + filepath, + generic-deriving, + semigroups >= 0.9, + simple-reflect >= 0.3.1 + + if impl(ghc<7.6.1) + ghc-options: -Werror diff --git a/src/Numeric/Compensated.hs b/src/Numeric/Compensated.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Compensated.hs @@ -0,0 +1,662 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE PatternGuards #-} +-------------------------------------------------------------------- +-- | +-- Copyright : (c) Edward Kmett 2013 +-- License : BSD3 +-- Maintainer: Edward Kmett +-- Stability : experimental +-- Portability: non-portable +-- +-- This module provides a fairly extensive API for compensated +-- floating point arithmetic based on Knuth's error free +-- transformation, various algorithms by Ogita, Rump and Oishi, +-- Hida, Li and Bailey, Kahan summation, etc. with custom compensated +-- arithmetic circuits to do multiplication, division, etc. of compensated +-- numbers. +-- +-- In general if @a@ has x bits of significand, @Compensated a@ gives +-- you twice that. You can iterate this construction for arbitrary +-- precision. +-- +-- References: +-- +-- * +-- * +-- * Donald Knuth's \"The Art of Computer Programming, Volume 2: Seminumerical Algorithms\" +-- * +-------------------------------------------------------------------- +module Numeric.Compensated + ( Compensable(..) + , _Compensated + , Overcompensated + , primal + , residual + , uncompensated + , fadd + -- * lifting scalars + , add, times, squared, divide, split + , kahan, (+^), (*^) + -- * compensated operators + , square + ) where + +import Control.Applicative +import Control.Lens as L +import Control.Monad +import Data.Foldable as Foldable +import Data.Function (on) +import Data.Ratio +import Data.Binary as Binary +import Data.SafeCopy +import Data.Serialize as Serialize +import Data.Semigroup +import Data.Vector.Unboxed as U +import Data.Vector.Generic as G +import Data.Vector.Generic.Mutable as M +import Data.Data +import Foreign.Ptr +import Foreign.Storable +import Numeric.Log +import Text.Read as T +import Text.Show as T + +-- $setup +-- >>> :load Numeric.Compensated + +{-# ANN module "hlint: ignore Use -" #-} +{-# ANN module "hlint: ignore Use curry" #-} + +-- | @'add' a b k@ computes @k x y@ such that +-- +-- > x + y = a + b +-- > x = fl(a + b) +-- +-- Which is to say that @x@ is the floating point image of @(a + b)@ and +-- @y@ stores the residual error term. +add :: Num a => a -> a -> (a -> a -> r) -> r +add a b k = k x y where + x = a + b + z = x - a + y = (a - (x - z)) + (b - z) +{-# INLINE add #-} + +-- | @'fadd' a b k@ computes @k x y@ such that +-- +-- > x + y = a + b +-- > x = fl(a + b) +-- +-- but only under the assumption that @'abs' a '>=' 'abs' b@. If you +-- aren't sure, use 'add'. +-- +-- Which is to say that @x@ is the floating point image of @(a + b)@ and +-- @y@ stores the residual error term. +fadd :: Num a => a -> a -> (a -> a -> r) -> r +fadd a b k = k x (b - (x - a)) where + x = a + b +{-# INLINE fadd #-} + +-- | @'times' a b k@ computes @k x y@ such that +-- +-- > x + y = a * b +-- > x = fl(a * b) +-- +-- Which is to say that @x@ is the floating point image of @(a * b)@ and +-- @y@ stores the residual error term. +-- +-- This could be nicer if we had access to a hardware fused multiply-add. +times :: Compensable a => a -> a -> (a -> a -> r) -> r +times a b k = + split a $ \a1 a2 -> + split b $ \b1 b2 -> + let x = a * b in k x (a2*b2 - (((x - a1*b1) - a2*b1) - a1*b2)) + -- let x = a * b in k x (((a1*b1-x)+a1*b2+b2*b1)+a2*b2) +{-# INLINEABLE times #-} + +-- this is a variant on a division algorithm by Liddicoat and Flynn +divide :: Compensable a => a -> a -> (a -> a -> r) -> r +divide a b = with (aX * ms) where + x0 = recip b + aX = times a x0 compensated -- calculate aX + m = 1 <| negate (times b x0 compensated) + mm = m*m + ms = 1+((m+mm)+m*mm) +{-# INLINEABLE divide #-} + +-- | Priest's renormalization algorithm +-- +-- @renorm a b c@ generates a 'Compensated' number assuming @a '>=' b '>=' c@. +renorm :: Compensable a => a -> a -> a -> Compensated a +renorm a b c = + fadd b c $ \x1 y1 -> + fadd a x1 $ \x2 y2 -> + fadd x2 (y1 + y2) compensated +{-# INLINE renorm #-} + +-- | @'squared' a k@ computes @k x y@ such that +-- +-- > x + y = a * a +-- > x = fl(a * a) +-- +-- Which is to say that @x@ is the floating point image of @(a * a)@ and +-- @y@ stores the residual error term. +squared :: Compensable a => a -> (a -> a -> r) -> r +squared a k = + split a $ \a1 a2 -> + let x = a * a in k x (a2*a2 - ((x - a1*a1) - 2*(a2*a1))) +{-# INLINE squared #-} + +-- | Calculate a fast square of a compensated number. +square :: Compensable a => Compensated a -> Compensated a +square m = + with m $ \a b -> + squared a $ \x1 y1 -> + times a b $ \x2 y2 -> + add y1 (x2*2) $ \x3 y3 -> + renorm x1 x3 (b*b + 2*y2 + y3) +{-# INLINE square #-} + +-- | error-free split of a floating point number into two parts. +-- +-- Note: these parts do not satisfy the `compensated` contract +split :: Compensable a => a -> (a -> a -> r) -> r +split a k = k x y where + c = magic*a + x = c - (c - a) + y = a - x +{-# INLINEABLE split #-} + +-- | Calculate a scalar + compensated sum with Kahan summation. +(+^) :: Compensable a => a -> Compensated a -> Compensated a +a +^ m = with m $ \b c -> let y = a - c; t = b + y in compensated t ((t - b) - y) +{-# INLINE (+^) #-} + +-- | Compute @a * 'Compensated' a@ +(*^) :: Compensable a => a -> Compensated a -> Compensated a +c *^ m = + with m $ \ a b -> + times c a $ \x1 y1 -> + times c b $ \x2 y2 -> + fadd x1 x2 $ \x3 y3 -> + add y1 y3 $ \x4 y4 -> + renorm x3 x4 (y4 + y2) +{-# INLINE (*^) #-} + +class (RealFrac a, Precise a, Floating a) => Compensable a where + -- | This provides a numeric data type with effectively doubled precision by + -- using Knuth's error free transform and a number of custom compensated + -- arithmetic circuits. + -- + -- This construction can be iterated, doubling precision each time. + -- + -- >>> round (Prelude.product [2..100] :: Compensated (Compensated (Compensated Double))) + -- 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 + -- + -- >>> Prelude.product [2..100] + -- 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 + data Compensated a + + -- | This extracts both the 'primal' and 'residual' components of a 'Compensated' + -- number. + with :: Compensable a => Compensated a -> (a -> a -> r) -> r + + -- | Used internally to construct 'compensated' values that satisfy our residual contract. + -- + -- When in doubt, use @'add' a b 'compensated'@ instead of @'compensated' a b@ + compensated :: Compensable a => a -> a -> Compensated a + + -- | This 'magic' number is used to 'split' the significand in half, so we can multiply + -- them separately without losing precision in 'times'. + magic :: a + +instance Compensable Double where + data Compensated Double = CD {-# UNPACK #-} !Double {-# UNPACK #-} !Double + with (CD a b) k = k a b + {-# INLINE with #-} + compensated = CD + {-# INLINE compensated #-} + magic = 134217729 + {-# INLINE magic #-} + +instance Compensable Float where + data Compensated Float = CF {-# UNPACK #-} !Float {-# UNPACK #-} !Float + with (CF a b) k = k a b + {-# INLINE with #-} + compensated = CF + {-# INLINE compensated #-} + magic = 4097 + {-# INLINE magic #-} + +instance Compensable a => Compensable (Compensated a) where + data Compensated (Compensated a) = CC !(Compensated a) !(Compensated a) + with (CC a b) k = k a b + {-# INLINE with #-} + compensated = CC + {-# INLINE compensated #-} + magic = times (magic - 1) (magic - 1) $ \ x y -> compensated x (y + 1) + {-# INLINE magic #-} + +instance Typeable1 Compensated where + typeOf1 _ = mkTyConApp (mkTyCon3 "analytics" "Data.Analytics.Numeric.Compensated" "Compensated") [] + +instance (Compensable a, Data a) => Data (Compensated a) where + gfoldl f z m = with m $ \a b -> z compensated `f` a `f` b + toConstr _ = compensatedConstr + gunfold k z c = case constrIndex c of + 1 -> k (k (z compensated)) + _ -> error "gunfold" + dataTypeOf _ = compensatedDataType + dataCast1 f = gcast1 f + +compensatedConstr :: Constr +compensatedConstr = mkConstr compensatedDataType "compensated" [] Prefix +{-# NOINLINE compensatedConstr #-} + +compensatedDataType :: DataType +compensatedDataType = mkDataType "Data.Analytics.Numeric.Compensated" [compensatedConstr] +{-# NOINLINE compensatedDataType #-} + +instance (Compensable a, Show a) => Show (Compensated a) where + showsPrec d m = with m $ \a b -> showParen (d > 10) $ + showString "compensated " . T.showsPrec 11 a . showChar ' ' . T.showsPrec 11 b + +instance (Compensable a, Read a) => Read (Compensated a) where + readPrec = parens $ prec 10 $ do + Ident "compensated" <- lexP + a <- step T.readPrec + b <- step T.readPrec + return $ compensated a b + +type Overcompensated a = Compensated (Compensated a) + +-- | This provides the isomorphism between the compact representation we store these in internally +-- and the naive pair of the 'primal' and 'residual' components. +_Compensated :: Compensable a => Iso' (Compensated a) (a, a) +_Compensated = iso (`with` (,)) (uncurry compensated) +{-# INLINE _Compensated #-} + +-- | This 'Lens' lets us edit the 'primal' directly, leaving the 'residual' untouched. +primal :: Compensable a => Lens' (Compensated a) a +primal f c = with c $ \a b -> f a <&> \a' -> compensated a' b +{-# INLINE primal #-} + +-- | This 'Lens' lets us edit the 'residual' directly, leaving the 'primal' untouched. +residual :: Compensable a => Lens' (Compensated a) a +residual f c = with c $ \a b -> compensated a <$> f b +{-# INLINE residual #-} + +-- | Extract the 'primal' component of a 'compensated' value, when and if compensation +-- is no longer required. +uncompensated :: Compensable a => Compensated a -> a +uncompensated c = with c const +{-# INLINE uncompensated #-} + +type instance Index (Compensated a) = Int +instance (Applicative f, Compensable a, Compensable b) => Each f (Compensated a) (Compensated b) a b where + each f m = with m $ \a b -> compensated <$> L.indexed f (0 :: Int) a <*> L.indexed f (1 :: Int) b + {-# INLINE each #-} + +instance Compensable a => Eq (Compensated a) where + m == n = with m $ \a b -> with n $ \c d -> a == c && b == d + m /= n = with m $ \a b -> with n $ \c d -> a /= c && b /= d + {-# INLINE (==) #-} + +instance Compensable a => Ord (Compensated a) where + compare m n = with m $ \a b -> with n $ \c d -> compare a c <> compare b d + {-# INLINE compare #-} + m <= n = with m $ \a b -> with n $ \c d -> case compare a c of + LT -> True + EQ -> b <= d + GT -> False + {-# INLINE (<=) #-} + m >= n = with m $ \a b -> with n $ \c d -> case compare a c of + LT -> False + EQ -> b >= d + GT -> a >= c -- @compare x NaN@ and @compare naN x@ always return 'GT', but @m >= n@ should be 'False' + {-# INLINE (>=) #-} + m > n = with m $ \a b -> with n $ \c d -> case compare a c of + LT -> False + EQ -> b > d + GT -> a > c -- @compare x NaN@ and @compare naN x@ always return 'GT', but @m >= n@ should be 'False' + {-# INLINE (>) #-} + m < n = with m $ \a b -> with n $ \c d -> case compare a c of + LT -> True + EQ -> b < d + GT -> False + {-# INLINE (<) #-} + +instance Compensable a => Semigroup (Compensated a) where + (<>) = (+) + {-# INLINE (<>) #-} + +instance Compensable a => Monoid (Compensated a) where + mempty = compensated 0 0 + {-# INLINE mempty #-} + mappend = (+) + {-# INLINE mappend #-} + +-- | Perform Kahan summation over a list. +kahan :: (Foldable f, Compensable a) => f a -> Compensated a +kahan = Foldable.foldr (+^) mempty +{-# INLINE kahan #-} + +-- (<|) = (+^) +instance (Reviewable p, Functor f, Compensable a, a ~ b) => Cons p f (Compensated a) (Compensated b) a b where + _Cons = unto $ \(a, e) -> with e $ \b c -> let y = a - c; t = b + y in compensated t ((t - b) - y) + {-# INLINE _Cons #-} + +-- (|>) = (+^) +instance (Reviewable p, Functor f, Compensable a, a ~ b) => Snoc p f (Compensated a) (Compensated b) a b where + _Snoc = unto $ \(e, a) -> with e $ \b c -> let y = a - c; t = b + y in compensated t ((t - b) - y) + {-# INLINE _Snoc #-} + +instance Compensable a => Num (Compensated a) where + m + n = + with m $ \a b -> + with n $ \c d -> + add a c $ \x1 y1 -> + add y1 d $ \x2 y2 -> + add b x2 $ \x3 y3 -> + add x1 x3 $ \x4 y4 -> + add x4 (y2 + y3 + y4) compensated + {-# INLINE (+) #-} + +{- + m + n = + with m $ \a b -> + with n $ \c d -> + add a c $ \x1 y1 -> + add b d $ \x2 y2 -> + renorm x1 x2 (y1 + y2) + {-# INLINE (+) #-} +-} + + + m * n = + with m $ \a b -> + with n $ \c d -> + times a c $ \x1 y1 -> + times b c $ \x2 y2 -> + times a d $ \x3 y3 -> + add x1 x2 $ \x4 y4 -> + add x3 x4 $ \x5 y5 -> + add y1 y4 $ \x6 y6 -> + add y5 x6 $ \x7 y7 -> + add x5 x7 $ \x8 y8 -> + add x8 (b*d + y2 + y3 + y6 + y7 + y8) compensated + {-# INLINE (*) #-} + +{- + m * n = + with m $ \a b -> + with n $ \c d -> + times a c $ \x1 y1 -> + times b c $ \x2 y2 -> + times a d $ \x3 y3 -> + add y1 x2 $ \x4 y4 -> + add x3 x4 $ \x5 y5 -> + renorm x1 x5 (b * d + y2 + y4 + y3 + y5) + {-# INLINE (*) #-} +-} + + negate m = with m (on compensated negate) + -- {-# INLINE negate #-} + + x - y = x + negate y + {-# INLINE (-) #-} + + signum m = with m $ \a _ -> compensated (signum a) 0 + {-# INLINE signum #-} + + abs m = with m $ \a b -> + if a < 0 + then compensated (negate a) (negate b) + else compensated a b + {-# INLINE abs #-} + + fromInteger i = add x (fromInteger (i - round x)) compensated where + x = fromInteger i + {-# INLINE fromInteger #-} + +instance Compensable a => Enum (Compensated a) where + succ a = a + 1 + {-# INLINE succ #-} + pred a = a - 1 + {-# INLINE pred #-} + toEnum i = add x (fromIntegral (i - round x)) compensated where + x = fromIntegral i + {-# INLINE toEnum #-} + fromEnum = round + {-# INLINE fromEnum #-} + enumFrom a = a : Prelude.enumFrom (a + 1) + {-# INLINE enumFrom #-} + enumFromThen a b = a : Prelude.enumFromThen b (b - a + b) + {-# INLINE enumFromThen #-} + enumFromTo a b + | a <= b = a : Prelude.enumFromTo (a + 1) b + | otherwise = [] + {-# INLINE enumFromTo #-} + enumFromThenTo a b c + | a <= b = up a + | otherwise = down a + where + delta = b - a + up x | x <= c = x : up (x + delta) + | otherwise = [] + down x | c <= x = x : down (x + delta) + | otherwise = [] + {-# INLINE enumFromThenTo #-} + +instance Compensable a => Fractional (Compensated a) where + recip m = with m $ \a b -> add (recip a) (-b / (a * a)) compensated + {-# INLINE recip #-} + + -- | A variant on a hardware division algorithm by Liddicoat and Flynn + a / b = (a*x0) * (1+((m+mm)+m*mm)) where + x0 = recip b + m = 1 - b*x0 + mm = m*m + -- {-# INLINE (/) #-} + + fromRational r = fromInteger (numerator r) / fromInteger (denominator r) + -- {-# INLINE fromRational #-} + +instance Compensable a => Real (Compensated a) where + toRational m = with m (on (+) toRational) + -- {-# INLINE toRational #-} + +instance Compensable a => RealFrac (Compensated a) where + properFraction m = with m $ \a b -> case properFraction a of + (w, p) -> add p b $ \ x y -> case properFraction x of + (w',q) -> (w + w', add q y compensated) + -- {-# INLINE properFraction #-} + +instance (Compensable a, Binary a) => Binary (Compensated a) where + get = compensated <$> Binary.get <*> Binary.get + put m = with m $ \a b -> do + Binary.put a + Binary.put b + +instance (Compensable a, Serialize a) => Serialize (Compensated a) where + get = compensated <$> Serialize.get <*> Serialize.get + put m = with m $ \a b -> do + Serialize.put a + Serialize.put b + +-- ಠ_ಠ this unnecessarily expects that the format won't change, because I can't derive a better instance. +instance (Compensable a, Serialize a) => SafeCopy (Compensated a) + +instance (Compensable a, Storable a) => Storable (Compensated a) where + sizeOf _ = sizeOf (undefined :: a) * 2 + -- {-# INLINE sizeOf #-} + alignment _ = alignment (undefined :: a) + -- {-# INLINE alignment #-} + peekElemOff p o | q <- castPtr p, o2 <- o * 2 = + compensated <$> peekElemOff q o2 <*> peekElemOff q (o2+1) + -- {-# INLINE peekElemOff #-} + pokeElemOff p o m | q <- castPtr p, o2 <- o * 2 = with m $ \a b -> do + pokeElemOff q o2 a + pokeElemOff q (o2+1) b + -- {-# INLINE pokeElemOff #-} + peekByteOff p o | q <- castPtr p = + compensated <$> peekByteOff q o <*> peekByteOff q (o + sizeOf (undefined :: a)) + -- {-# INLINE peekByteOff #-} + pokeByteOff p o m | q <- castPtr p = with m $ \a b -> do + pokeByteOff q o a + pokeByteOff q (o+sizeOf (undefined :: a)) b + -- {-# INLINE pokeByteOff #-} + peek p | q <- castPtr p = compensated <$> peek q <*> peekElemOff q 1 + -- {-# INLINE peek #-} + poke p m | q <- castPtr p = with m $ \a b -> do + poke q a + pokeElemOff q 1 b + -- {-# INLINE poke #-} + +newtype instance U.MVector s (Compensated a) = MV_Compensated (U.MVector s (a,a)) +newtype instance U.Vector (Compensated a) = V_Compensated (U.Vector (a, a)) + +instance (Compensable a, Unbox a) => M.MVector U.MVector (Compensated a) where + basicLength (MV_Compensated v) = M.basicLength v + {-# INLINE basicLength #-} + basicUnsafeSlice i n (MV_Compensated v) = MV_Compensated $ M.basicUnsafeSlice i n v + {-# INLINE basicUnsafeSlice #-} + basicOverlaps (MV_Compensated v1) (MV_Compensated v2) = M.basicOverlaps v1 v2 + {-# INLINE basicOverlaps #-} + basicUnsafeNew n = MV_Compensated `liftM` M.basicUnsafeNew n + {-# INLINE basicUnsafeNew #-} + basicUnsafeReplicate n m = with m $ \x y -> MV_Compensated `liftM` M.basicUnsafeReplicate n (x,y) + {-# INLINE basicUnsafeReplicate #-} + basicUnsafeRead (MV_Compensated v) i = uncurry compensated `liftM` M.basicUnsafeRead v i + {-# INLINE basicUnsafeRead #-} + basicUnsafeWrite (MV_Compensated v) i m = with m $ \ x y -> M.basicUnsafeWrite v i (x,y) + {-# INLINE basicUnsafeWrite #-} + basicClear (MV_Compensated v) = M.basicClear v + {-# INLINE basicClear #-} + basicSet (MV_Compensated v) m = with m $ \ x y -> M.basicSet v (x,y) + {-# INLINE basicSet #-} + basicUnsafeCopy (MV_Compensated v1) (MV_Compensated v2) = M.basicUnsafeCopy v1 v2 + {-# INLINE basicUnsafeCopy #-} + basicUnsafeMove (MV_Compensated v1) (MV_Compensated v2) = M.basicUnsafeMove v1 v2 + {-# INLINE basicUnsafeMove #-} + basicUnsafeGrow (MV_Compensated v) n = MV_Compensated `liftM` M.basicUnsafeGrow v n + {-# INLINE basicUnsafeGrow #-} + +instance (Compensable a, Unbox a) => G.Vector U.Vector (Compensated a) where + basicUnsafeFreeze (MV_Compensated v) = V_Compensated `liftM` G.basicUnsafeFreeze v + {-# INLINE basicUnsafeFreeze #-} + basicUnsafeThaw (V_Compensated v) = MV_Compensated `liftM` G.basicUnsafeThaw v + {-# INLINE basicUnsafeThaw #-} + basicLength (V_Compensated v) = G.basicLength v + {-# INLINE basicLength #-} + basicUnsafeSlice i n (V_Compensated v) = V_Compensated $ G.basicUnsafeSlice i n v + {-# INLINE basicUnsafeSlice #-} + basicUnsafeIndexM (V_Compensated v) i + = uncurry compensated `liftM` G.basicUnsafeIndexM v i + {-# INLINE basicUnsafeIndexM #-} + basicUnsafeCopy (MV_Compensated mv) (V_Compensated v) + = G.basicUnsafeCopy mv v + {-# INLINE basicUnsafeCopy #-} + elemseq _ m z = with m $ \x y -> G.elemseq (undefined :: U.Vector a) x + $ G.elemseq (undefined :: U.Vector a) y z + {-# INLINE elemseq #-} + +-- | /NB:/ Experimental and partially implemented. +-- +-- Other than sqrt, the accuracy of these is basically uncalculated! In fact many of these are known to be wrong! Patches and improvements are welcome. +instance Compensable a => Floating (Compensated a) where +#ifdef SPECIALIZE_INSTANCES + {-# SPECIALIZE instance Floating (Compensated Double) #-} + {-# SPECIALIZE instance Floating (Compensated Float) #-} + {-# SPECIALIZE instance Compensable a => Floating (Compensated (Compensated a)) #-} +#endif + + exp m = + with m $ \a b -> + times (exp a) (exp b) compensated + + sin m = + with m $ \a b -> + times (sin a) (cos b) $ \x1 y1 -> + times (sin b) (cos a) $ \x2 y2 -> + add x1 x2 $ \x3 y3 -> + add y1 y2 $ \x4 y4 -> + add x4 y3 $ \x5 y5 -> + add x5 x3 $ \x6 y6 -> + add (y4 + y5 + y6) x6 compensated + + cos m = + with m $ \a b -> + times (cos a) (cos b) $ \x1 y1 -> + times (-sin b) (sin a) $ \x2 y2 -> + add x1 x2 $ \x3 y3 -> + add y1 y2 $ \x4 y4 -> + add x4 y3 $ \x5 y5 -> + add x5 x3 $ \x6 y6 -> + add (y4 + y5 + y6) x6 compensated + + tan m = + with m $ \a b -> + add (tan a) (tan b) compensated / + (1 +^ times (tan a) (tan b) compensated) + + sinh m = + with m $ \a b -> + times (sinh a) (cosh b) $ \x1 y1 -> + times (cosh a) (sinh b) $ \x2 y2 -> + add x1 x2 $ \x3 y3 -> + add y1 y2 $ \x4 y4 -> + add x4 y3 $ \x5 y5 -> + add x5 x3 $ \x6 y6 -> + add (y4 + y5 + y6) x6 compensated + + cosh m = + with m $ \a b -> + times (cosh a) (cosh b) $ \x1 y1 -> + times (sinh b) (sinh a) $ \x2 y2 -> + add x1 x2 $ \x3 y3 -> + add y1 y2 $ \x4 y4 -> + add x4 y3 $ \x5 y5 -> + add x5 x3 $ \x6 y6 -> + add (y4 + y5 + y6) x6 compensated + + tanh m = + with m $ \a b -> + fadd (tanh a) (tanh b) compensated / + (1 +^ times (tanh a) (tanh b) compensated) + + -- This requires an accurate 'exp', which we currently lack. + log m = + with m $ \ a b -> let + xy1 = add (log a) (b/a) compensated + xy2 = xy1 + m * exp (-xy1) - 1 -- Newton Raphson step 1 + in xy2 + m * exp (-xy2) - 1 -- Newton Raphson step 2 + + -- | Hardware sqrt improved by the Babylonian algorithm (Newton Raphson) + sqrt m = with (z4 + m/z4) $ on compensated (/2) where + z0 = sqrt (m^.primal) + z1 = with (z0 <| (m / compensated z0 0)) $ on compensated (/2) + z2 = with (z1 + m/z1) $ on compensated (/2) + z3 = with (z2 + m/z2) $ on compensated (/2) + z4 = with (z3 + m/z3) $ on compensated (/2) + + -- (**) = error "TODO" + pi = error "TODO" + asin = error "TODO" + atan = error "TODO" + acos = error "TODO" + asinh = error "TODO" + atanh = error "TODO" + acosh = error "TODO" + +-- | TODO: do this right! +instance (Compensable a, Precise a) => Precise (Compensated a) where + log1p a = log (1 + a) + {-# INLINE log1p #-} + expm1 a = exp a - 1 + {-# INLINE expm1 #-} diff --git a/tests/doctests.hsc b/tests/doctests.hsc new file mode 100644 --- /dev/null +++ b/tests/doctests.hsc @@ -0,0 +1,74 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +-- | +-- Module : Main (doctests) +-- Copyright : (C) 2012-13 Edward Kmett +-- License : BSD-style (see the file LICENSE) +-- Maintainer : Edward Kmett +-- Stability : provisional +-- Portability : portable +-- +-- This module provides doctests for a project based on the actual versions +-- of the packages it was built with. It requires a corresponding Setup.lhs +-- to be added to the project +----------------------------------------------------------------------------- +module Main where + +import Build_doctests (deps) +import Control.Applicative +import Control.Monad +import Data.List +import System.Directory +import System.FilePath +import Test.DocTest + +##if defined(mingw32_HOST_OS) +##if defined(i386_HOST_ARCH) +##define USE_CP +import Control.Applicative +import Control.Exception +import Foreign.C.Types +foreign import stdcall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool +foreign import stdcall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt +##elif defined(x86_64_HOST_ARCH) +##define USE_CP +import Control.Applicative +import Control.Exception +import Foreign.C.Types +foreign import ccall "windows.h SetConsoleCP" c_SetConsoleCP :: CUInt -> IO Bool +foreign import ccall "windows.h GetConsoleCP" c_GetConsoleCP :: IO CUInt +##endif +##endif + +-- | Run in a modified codepage where we can print UTF-8 values on Windows. +withUnicode :: IO a -> IO a +##ifdef USE_CP +withUnicode m = do + cp <- c_GetConsoleCP + (c_SetConsoleCP 65001 >> m) `finally` c_SetConsoleCP cp +##else +withUnicode m = m +##endif + +main :: IO () +main = withUnicode $ getSources >>= \sources -> doctest $ + "-isrc" + : "-idist/build/autogen" + : "-optP-include" + : "-optPdist/build/autogen/cabal_macros.h" + : "-hide-all-packages" + : "-Iincludes" + : map ("-package="++) deps ++ sources + +getSources :: IO [FilePath] +getSources = filter (isSuffixOf ".hs") <$> go "src" + where + go dir = do + (dirs, files) <- getFilesAndDirectories dir + (files ++) . concat <$> mapM go dirs + +getFilesAndDirectories :: FilePath -> IO ([FilePath], [FilePath]) +getFilesAndDirectories dir = do + c <- map (dir ) . filter (`notElem` ["..", "."]) <$> getDirectoryContents dir + (,) <$> filterM doesDirectoryExist c <*> filterM doesFileExist c diff --git a/travis/cabal-apt-install b/travis/cabal-apt-install new file mode 100644 --- /dev/null +++ b/travis/cabal-apt-install @@ -0,0 +1,27 @@ +#! /bin/bash +set -eu + +APT="sudo apt-get -q -y" +CABAL_INSTALL_DEPS="cabal install --only-dependencies --force-reinstall" + +$APT update +$APT install dctrl-tools + +# Find potential system packages to satisfy cabal dependencies +deps() +{ + local M='^\([^ ]\+\)-[0-9.]\+ (.*$' + local G=' -o ( -FPackage -X libghc-\L\1\E-dev )' + local E="$($CABAL_INSTALL_DEPS "$@" --dry-run -v 2> /dev/null \ + | sed -ne "s/$M/$G/p" | sort -u)" + grep-aptavail -n -sPackage \( -FNone -X None \) $E | sort -u +} + +$APT install $(deps "$@") libghc-quickcheck2-dev # QuickCheck is special +$CABAL_INSTALL_DEPS "$@" # Install the rest via Hackage + +if ! $APT install hlint ; then + $APT install $(deps hlint) + cabal install hlint +fi + diff --git a/travis/config b/travis/config new file mode 100644 --- /dev/null +++ b/travis/config @@ -0,0 +1,16 @@ +-- This provides a custom ~/.cabal/config file for use when hackage is down that should work on unix +-- +-- This is particularly useful for travis-ci to get it to stop complaining +-- about a broken build when everything is still correct on our end. +-- +-- This uses Luite Stegeman's mirror of hackage provided by his 'hdiff' site instead +-- +-- To enable this, uncomment the before_script in .travis.yml + +remote-repo: hdiff.luite.com:http://hdiff.luite.com/packages/archive +remote-repo-cache: ~/.cabal/packages +world-file: ~/.cabal/world +build-summary: ~/.cabal/logs/build.log +remote-build-reporting: anonymous +install-dirs user +install-dirs global