compensated (empty) → 0.1
raw patch · 14 files changed
+1087/−0 lines, 14 filesdep +basedep +bifunctorsdep +binarybuild-type:Customsetup-changed
Dependencies added: base, bifunctors, binary, cereal, comonad, deepseq, directory, distributive, doctest, filepath, generic-deriving, hashable, lens, log-domain, safecopy, semigroupoids, semigroups, simple-reflect, vector
Files
- .ghci +2/−0
- .gitignore +35/−0
- .travis.yml +26/−0
- .vim.custom +31/−0
- AUTHORS.markdown +15/−0
- CHANGELOG.markdown +3/−0
- LICENSE +30/−0
- README.markdown +23/−0
- Setup.lhs +55/−0
- compensated.cabal +88/−0
- src/Numeric/Compensated.hs +662/−0
- tests/doctests.hsc +74/−0
- travis/cabal-apt-install +27/−0
- travis/config +16/−0
+ .ghci view
@@ -0,0 +1,2 @@+:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h -optP-Iincludes+:set -v0
+ .gitignore view
@@ -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
+ .travis.yml view
@@ -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"
+ .vim.custom view
@@ -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 <F2> :exec ":!hasktags -x -c --ignore src"<CR><CR>++" 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"
+ AUTHORS.markdown view
@@ -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
+ CHANGELOG.markdown view
@@ -0,0 +1,3 @@+0.1+---+* Split from [analytics](https://github.com/analytics/analytics)
+ LICENSE view
@@ -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.
+ README.markdown view
@@ -0,0 +1,23 @@+compensated+===========++[](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
+ Setup.lhs view
@@ -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}
+ compensated.cabal view
@@ -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 <ekmett@gmail.com>+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
+ src/Numeric/Compensated.hs view
@@ -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 <ekmett@gmail.com>+-- 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:+--+-- * <http://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4-old/docs/qd.pdf>+-- * <http://www.ti3.tuhh.de/paper/rump/OgRuOi05.pdf>+-- * Donald Knuth's \"The Art of Computer Programming, Volume 2: Seminumerical Algorithms\"+-- * <http://en.wikipedia.org/wiki/Kahan_summation_algorithm>+--------------------------------------------------------------------+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 #-}
+ tests/doctests.hsc view
@@ -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 <ekmett@gmail.com>+-- 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
+ travis/cabal-apt-install view
@@ -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+
+ travis/config view
@@ -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