variable-precision 0.1.1 → 0.2
raw patch · 14 files changed
+1110/−360 lines, 14 filesdep +integer-gmp
Dependencies added: integer-gmp
Files
- CHANGES +1/−0
- Numeric/VariablePrecision.hs +27/−0
- Numeric/VariablePrecision/Algorithms.hs +300/−0
- Numeric/VariablePrecision/Aliases.hs +82/−0
- Numeric/VariablePrecision/Complex.hs +97/−14
- Numeric/VariablePrecision/Float.hs +500/−188
- Numeric/VariablePrecision/Float/Aliases.hs +0/−95
- Numeric/VariablePrecision/Precision.hs +41/−30
- Numeric/VariablePrecision/Precision/Reify.hs +12/−10
- TODO +2/−15
- TypeLevel/NaturalNumber/ExtraNumbers.hs +0/−2
- fast/Numeric/VariablePrecision/Integer/Logarithm.hs +8/−0
- pure/Numeric/VariablePrecision/Integer/Logarithm.hs +13/−0
- variable-precision.cabal +27/−6
CHANGES view
@@ -1,2 +1,3 @@+v0.2 generic IEEE-ish v0.1.1 fixed for ghc-7.0.4 v0.1 initial release
+ Numeric/VariablePrecision.hs view
@@ -0,0 +1,27 @@+{- |+Module : Numeric.VariablePrecision+Copyright : (c) Claude Heiland-Allen 2012+License : BSD3++Maintainer : claudiusmaximus@goto10.org+Stability : unstable+Portability : portable++Convenience module.++-}+module Numeric.VariablePrecision+ ( module Numeric.VariablePrecision.Float+ , module Numeric.VariablePrecision.Complex+ , module Numeric.VariablePrecision.Precision+ , module Numeric.VariablePrecision.Precision.Reify+ , module Numeric.VariablePrecision.Aliases+ , module Numeric.VariablePrecision.Algorithms+ ) where++import Numeric.VariablePrecision.Float+import Numeric.VariablePrecision.Complex+import Numeric.VariablePrecision.Precision+import Numeric.VariablePrecision.Precision.Reify+import Numeric.VariablePrecision.Aliases+import Numeric.VariablePrecision.Algorithms
+ Numeric/VariablePrecision/Algorithms.hs view
@@ -0,0 +1,300 @@+{-# LANGUAGE BangPatterns #-}+{- |+Module : Numeric.VariablePrecision.Algorithms+Copyright : (c) Claude Heiland-Allen 2012+License : BSD3++Maintainer : claudiusmaximus@goto10.org+Stability : unstable+Portability : BangPatterns++Implementations of various floating point algorithms. Accuracy has not+been extensively verified, and termination has not been proven.++Everything assumes that 'floatRadix' is 2. This is *not* checked.++Functions taking an @accuracy@ parameter may fail to terminate if+@accuracy@ is too small. Accuracy is measured in least significant+bits, similarly to '(=~=)'.++In this documentation, /basic functionality/ denotes that methods used+are from classes:++ * 'Num', 'Eq', 'Ord'.++Further, /basic RealFloat functionality/ denotes /basic functionality/+with the addition of:++ * Anything in 'RealFloat' except for 'atan2'.++The intention behind the used functionality documentation is to help+users decide when it is appropriate to use these generic implementations+to implement instances.++-}+module Numeric.VariablePrecision.Algorithms+ ( recodeFloat+ , viaDouble+ , (=~=)+ , genericRecip+ , genericSqrt+ , genericExp+ , genericLog+ , genericLog'+ , genericLog2+ , genericLog''+ , genericPi+ , genericPositiveZero+ , genericNegativeZero+ , genericPositiveInfinity+ , genericNegativeInfinity+ , genericNotANumber+ , sameSign+ ) where++import Data.Bits (bit, shiftR)+++-- | Special values implemented using basic RealFloat functionality.+genericPositiveZero, genericNegativeZero, genericPositiveInfinity, genericNegativeInfinity, genericNotANumber :: RealFloat a => a++genericPositiveZero = 0++genericNegativeZero = -0++genericPositiveInfinity = result+ where+ result = encodeFloat m e+ m = bit (floatDigits (undefined `asTypeOf` result))+ e = snd (floatRange (undefined `asTypeOf` result))++genericNegativeInfinity = result+ where+ result = encodeFloat (negate m) e+ m = bit (floatDigits (undefined `asTypeOf` result))+ e = snd (floatRange (undefined `asTypeOf` result))++genericNotANumber = genericPositiveInfinity + genericNegativeInfinity+++-- | Convert between generic 'RealFloat' types more efficiently than+-- 'realToFrac'. Tries hard to preserve special values like+-- infinities and negative zero, but any NaN payload is lost.+--+-- Uses only basic RealFloat functionality.+--+recodeFloat :: (RealFloat a, RealFloat b) => a -> b+recodeFloat !x+ | isNaN x = genericNotANumber+ | isInfinite x && x > 0 = genericPositiveInfinity+ | isInfinite x && x < 0 = genericNegativeInfinity+ | isNegativeZero x = genericNegativeZero+ | x == 0 = genericPositiveZero+ | otherwise = uncurry encodeFloat (decodeFloat x)+++-- | Check if two numbers have the same sign.+-- May give a nonsense result if an argument is NaN.+sameSign :: (Ord a, Num a) => a -> a -> Bool+sameSign a b = compare 0 a == compare 0 b+++-- | Approximate equality.+-- @(a =~= b) c@ when adding the difference to the larger in magnitude+-- changes at most @c@ least significant mantissa bits.+--+-- Uses only basic RealFloat functionality.+--+(=~=) :: RealFloat a => a -> a -> Int -> Bool+(=~=) !x !y !s+ | x == y = True+ | isNaN x && isNaN y = True+ | isNaN x || isNaN y = False+ | isInfinite x || isInfinite y = False+ | not (sameSign a b) = False+ | otherwise = abs (e - f) <= s && abs (x - y) <= encodeFloat 1 (s + (e `max` f))+ where+ (a, e) = decodeFloat x+ (b, f) = decodeFloat y+++-- | Compute a reciprocal using the Newton-Raphson division algorithm,+-- as described in+-- <http://en.wikipedia.org/wiki/Division_%28digital%29#Newton.E2.80.93Raphson_division>.+--+-- Uses only basic RealFloat functionality.+--+genericRecip :: RealFloat a => Int {- ^ accuracy -} -> a -> a+genericRecip accuracy y = recip' y+ where+ recip' f0+ | isNaN f0 = f0+ | isInfinite f0 && f0 > 0 = genericPositiveZero+ | isInfinite f0 && f0 < 0 = genericNegativeZero+ | isNegativeZero f0 = genericNegativeInfinity+ | f0 == 0 = genericPositiveInfinity+ | f0 < 0 = negate . recip' . negate $ f0+ | otherwise = scaleFloat sh (go d s0 x0)+ where+ x0 = k48 - k32 * d+ d = significand f0 -- in [0.5,1)+ sh = exponent d - exponent f0+ go !d !s !x+ | (x =~= x') accuracy = x'+ | s == 0 = x'+ | otherwise = go d (s - 1) x'+ where+ x' = scaleFloat 1 x - d * x * x -- x * (2 - d * x)+ -- an attempt to avoid recomputing per-type constants+ p = floatDigits (undefined `asTypeOf` y)+ s0 = ceiling (logBase 2 (fromIntegral (p + 1) / logBase 2 17) :: Double) :: Int+ k48 = recodeFloat (48/17 :: Double)+ k32 = recodeFloat (32/17 :: Double)+++-- | Compute a square root using Newton's method.+--+-- Uses basic RealFloat functionality and '(/)'.+--+genericSqrt :: RealFloat a => Int {- ^ accuracy -} -> a -> a+genericSqrt accuracy f0+ | f0 < 0 = genericNotANumber+ | f0 == 0 = f0 -- preserves negative zero+ | isNaN f0 = f0+ | isInfinite f0 = f0+ | otherwise = go (viaDouble sqrt f)+ where+ e = exponent f0+ d = if even e then 2 else 1+ s = e - d -- even+ f = scaleFloat (negate s) f0 -- in [1,4)+ go !r =+ let r' = scaleFloat (-1) (r + f / r)+ in if (r =~= r') accuracy then scaleFloat (s `shiftR` 1) r' else go r'+++-- | Compute an exponential using power series.+--+-- Uses basic RealFloat functionality, '(/)' and 'recip'.+--+genericExp :: RealFloat a => Int {-^ accuracy -} -> a -> a+genericExp accuracy x+ | isNaN x = x+ | isInfinite x && x < 0 = 0+ | isInfinite x = x+ | x == 0 = 1+ | x < 0 = recip . genericExp accuracy . negate $ x+ | otherwise = go 0 1 1+ where+ go !s !xnnf{- x^n / n! -} !n+ | (s =~= s') accuracy = s'+ | otherwise = go s' (xnnf * x / fromIntegral n) (n + 1 :: Int)+ where+ s' = s + xnnf+++-- | Compute a logarithm.+--+-- See 'genericLog''' for algorithmic references.+--+-- Uses basic RealFloat functionality, 'sqrt' and 'recip'.+--+genericLog :: RealFloat a => Int {- ^ accuracy -} -> a -> a+genericLog accuracy = genericLog' accuracy (genericLog2 accuracy)+++-- | Compute log 2.+--+-- See 'genericLog''' for algorithmic references.+--+-- Uses basic RealFloat functionality, 'sqrt' and 'recip'.+--+genericLog2 :: RealFloat a => Int {- ^ accuracy -} -> a+genericLog2 accuracy = negate (genericLog'' accuracy 0.5)+++-- | Compute a logarithm using decomposition and a value for @log 2@.+--+-- See 'genericLog''' for algorithmic references.+--+-- Uses basic RealFloat functionality, 'sqrt', and 'recip'.+--+genericLog' :: RealFloat a => Int {- ^ accuracy -} -> a {- ^ log 2 -} -> a -> a+genericLog' accuracy ln2 x+ | isNaN x = x+ | x == 0 = genericNegativeInfinity+ | x < 0 = genericNotANumber+ | isInfinite x = x+ | otherwise = mln2 + genericLog'' accuracy s+ where+ m = exponent x+ s = significand x+ mln2 -- micro-optimisation+ | m == 0 = 0+ | otherwise = fromIntegral m * ln2+++-- | Compute a logarithm for a value in [0.5,1) using the AGM method+-- as described in section 7 of+-- /The Logarithmic Constant: log 2/+-- Xavier Gourdon and Pascal Sebah, May 18, 2010,+-- <http://numbers.computation.free.fr/Constants/Log2/log2.ps>.+--+-- The precondition is not checked.+--+-- Uses basic RealFloat functionality, 'sqrt', and 'recip'.+--+genericLog'' :: RealFloat a => Int {- ^ accuracy -} -> a {- ^ value in [0.5,1) -} -> a+genericLog'' accuracy x = result+ where+ result = go (-1) 1 (encodeFloat 1 m) 0 1 (scaleFloat m x) 0+ m2 = accuracy - floatDigits (undefined `asTypeOf` result)+ m = m2 `shiftR` 1+ small y = y == 0 || exponent y <= m2+ go !n !a !b !s !c !d !t+ | small ds && small dt = recip (1 - s') - recip (1 - t')+ | otherwise = go n' a' b' s' c' d' t'+ where+ a' = scaleFloat (-1) (a + b)+ c' = scaleFloat (-1) (c + d)+ b' = sqrt (a * b)+ d' = sqrt (c * d)+ ds = scaleFloat n (a * a - b * b)+ dt = scaleFloat n (c * c - d * d)+ t' = t + dt+ s' = s + ds+ n' = n + 1+++-- | Compute pi using the method described in section 8 of+-- /Multiple-precision zero-finding methods and the complexity of elementary function evaluation/+-- Richard P Brent, 1975 (revised May 30, 2010),+-- <http://arxiv.org/abs/1004.3412>.+--+-- Uses basic RealFloat functionality, '(/)', and 'sqrt'.+--+genericPi :: RealFloat a => Int {- ^ accuracy -} -> a+-- Works ok up to around 600,000 bits (178,000 decimal digits) but after+-- that further increase to mantissa precision leads to problems.+-- Output compared against /Pi/ by Scott Hemphill <http://www.gutenberg.org/ebooks/50>.+genericPi accuracy = result+ where+ sqr x = x * x+ result = go 1 (sqrt 0.5) 0.25 0 1+ go !a !b !t !k !p+ | (p =~= p') accuracy = p'+ | otherwise = go a' b' t' k' p'+ where+ a' = scaleFloat (-1) (a + b)+ b' = sqrt (a * b)+ t' = t - scaleFloat k (sqr (a' - a))+ k' = k + 1+ p' = scaleFloat (-2) (sqr (a + b) / t)+++-- | Lift a function from Double to generic 'RealFloat' types.+viaDouble :: (RealFloat a, RealFloat b) => (Double -> Double) -> a -> b+viaDouble f = recodeFloat . f . recodeFloat+++-- FIXME everything assumes that floatRadix is 2 without checking
+ Numeric/VariablePrecision/Aliases.hs view
@@ -0,0 +1,82 @@+{-|+Module : Numeric.VariablePrecision.Aliases+Copyright : (c) Claude Heiland-Allen 2012+License : BSD3++Maintainer : claudiusmaximus@goto10.org+Stability : unstable+Portability : portable++Aliases for 'recodeFloat' and 'recodeComplex' with specialized types.++Aliases for commonly desired types.++-}+module Numeric.VariablePrecision.Aliases+ ( toFloat, fromFloat, toDouble, fromDouble+ , toComplexFloat, fromComplexFloat, toComplexDouble, fromComplexDouble+ , F8, F16, F24, F32, F40, F48, F53+ , f8, f16, f24, f32, f40, f48, f53+ , C8, C16, C24, C32, C40, C48, C53+ , c8, c16, c24, c32, c40, c48, c53+ , module TypeLevel.NaturalNumber+ , module TypeLevel.NaturalNumber.ExtraNumbers+ ) where++import TypeLevel.NaturalNumber (N8, n8)+import TypeLevel.NaturalNumber.ExtraNumbers+ (N16, n16, N24, n24, N32, n32, N40, n40, N48, n48, N53, n53)++import Numeric.VariablePrecision.Float (VFloat)+import Numeric.VariablePrecision.Complex (VComplex, recodeComplex, toComplex, fromComplex)+import Numeric.VariablePrecision.Algorithms (recodeFloat)++import Data.Complex (Complex)++-- | Convert to a Float from the same precision.+toFloat :: F24 -> Float+toFloat = recodeFloat++-- | Convert from a Float to the same precision.+fromFloat :: Float -> F24+fromFloat = recodeFloat++-- | Convert to a Double from the same precision.+toDouble :: F53 -> Double+toDouble = recodeFloat++-- | Convert from a Double to the same precision.+fromDouble :: Double -> F53+fromDouble = recodeFloat++-- | Convert to a Float from the same precision.+toComplexFloat :: C24 -> Complex Float+toComplexFloat = recodeComplex . toComplex++-- | Convert from a Float to the same precision.+fromComplexFloat :: Complex Float -> C24+fromComplexFloat = fromComplex . recodeComplex++-- | Convert to a Double from the same precision.+toComplexDouble :: C53 -> Complex Double+toComplexDouble = recodeComplex . toComplex++-- | Convert from a Double to the same precision.+fromComplexDouble :: Complex Double -> C53+fromComplexDouble = fromComplex . recodeComplex++type F8 = VFloat N8 ; f8 :: F8 ; f8 = 0+type F16 = VFloat N16 ; f16 :: F16 ; f16 = 0+type F24 = VFloat N24 ; f24 :: F24 ; f24 = 0+type F32 = VFloat N32 ; f32 :: F32 ; f32 = 0+type F40 = VFloat N40 ; f40 :: F40 ; f40 = 0+type F48 = VFloat N48 ; f48 :: F48 ; f48 = 0+type F53 = VFloat N53 ; f53 :: F53 ; f53 = 0++type C8 = VComplex N8 ; c8 :: C8 ; c8 = 0+type C16 = VComplex N16 ; c16 :: C16 ; c16 = 0+type C24 = VComplex N24 ; c24 :: C24 ; c24 = 0+type C32 = VComplex N32 ; c32 :: C32 ; c32 = 0+type C40 = VComplex N40 ; c40 :: C40 ; c40 = 0+type C48 = VComplex N48 ; c48 :: C48 ; c48 = 0+type C53 = VComplex N53 ; c53 :: C53 ; c53 = 0
Numeric/VariablePrecision/Complex.hs view
@@ -1,32 +1,45 @@-{-# LANGUAGE DeriveDataTypeable, GeneralizedNewtypeDeriving #-}+{-# LANGUAGE DeriveDataTypeable, GeneralizedNewtypeDeriving, Rank2Types #-} {- | Module : Numeric.VariablePrecision.Complex Copyright : (c) Claude Heiland-Allen 2012 License : BSD3 Maintainer : claudiusmaximus@goto10.org-Stability : provisional-Portability : DeriveDataTypeable, GeneralizedNewtypeDeriving+Stability : unstable+Portability : DeriveDataTypeable, GeneralizedNewtypeDeriving, Rank2Types Newtype wrapper around 'Data.Complex'. When both of 'Data.Complex' and this module need to be imported, use qualified imports. -} module Numeric.VariablePrecision.Complex- ( VComplex(..)+ ( VComplex() , (.+)+ , (.*)+ , (*.)+ , toComplex , fromComplex , withComplex+ , mapComplex+ , recodeComplex+ , scaleComplex , realPart , imagPart , conjugate , magnitude , magnitude2+ , sqr , phase , polar , cis , mkPolar- , module Numeric.VariablePrecision.Float+ , scaleComplex'+ , magnitude2'+ , sqr'+ , DComplex(..)+ , toDComplex+ , fromDComplex+ , withDComplex ) where import Data.Data (Data())@@ -34,35 +47,47 @@ import qualified Data.Complex as X import Numeric.VariablePrecision.Float+import Numeric.VariablePrecision.Precision+import Numeric.VariablePrecision.Algorithms (recodeFloat) -- | Newtype wrapper around 'X.Complex' so that instances can be written--- for 'HasPrecision' and 'VariablePrecision'. -newtype VComplex p = C{ toComplex :: X.Complex (VFloat p) }+-- for 'HasPrecision' and 'VariablePrecision'.+newtype VComplex p = FromComplex+ { -- | Convert 'VComplex' to 'X.Complex'.+ toComplex :: X.Complex (VFloat p)+ } deriving (Eq, Num, Fractional, Floating, Data, Typeable) +-- | Convert 'X.Complex' to 'VComplex'.+fromComplex :: X.Complex (VFloat p) -> VComplex p+fromComplex = FromComplex+ -- | Alike to 'X.:+', constructs a complex number from a real part and -- an imaginary part. (.+) :: NaturalNumber p => VFloat p -> VFloat p -> VComplex p-x .+ y = C (x X.:+ y)+x .+ y = fromComplex (x X.:+ y) infix 6 .+ instance HasPrecision VComplex+ instance VariablePrecision VComplex where- adjustPrecision (C (x X.:+ y)) = C (adjustPrecision x X.:+ adjustPrecision y)+ adjustPrecision = withComplex (mapComplex adjustPrecision) --- | Convert 'X.Complex' to 'VComplex'.-fromComplex :: X.Complex (VFloat p) -> VComplex p-fromComplex = C+instance Normed VComplex where+ norm1 z = abs (realPart z) + abs (imagPart z)+ norm2 = magnitude+ norm2Squared = magnitude2+ normInfinity z = abs (realPart z) `max` abs (imagPart z) -- | Lift an operation on 'X.Complex' to one on 'VComplex'. withComplex :: (X.Complex (VFloat p) -> X.Complex (VFloat q)) -> (VComplex p -> VComplex q) withComplex f = fromComplex . f . toComplex instance NaturalNumber p => Show (VComplex p) where- showsPrec p (C c) = showsPrec p c+ showsPrec p = showsPrec p . toComplex instance NaturalNumber p => Read (VComplex p) where- readsPrec p = map (first C) . readsPrec p+ readsPrec p = map (first fromComplex) . readsPrec p where first f (a, b) = (f a, b) -- | Unit at phase.@@ -101,5 +126,63 @@ magnitude2 :: NaturalNumber p => VComplex p -> VFloat p magnitude2 = magnitude2' . toComplex +-- | Apply a function to both components of a complex number.+mapComplex :: (RealFloat a, RealFloat b) => (a -> b) -> X.Complex a -> X.Complex b+mapComplex f (x X.:+ y) = f x X.:+ f y++-- | Much like 'mapComplex' 'recodeFloat'.+recodeComplex :: (RealFloat a, RealFloat b) => X.Complex a -> X.Complex b+recodeComplex = mapComplex recodeFloat++-- | Magnitude squared. magnitude2' :: RealFloat r => X.Complex r -> r magnitude2' (x X.:+ y) = x * x + y * y++-- | Complex square.+sqr :: NaturalNumber p => VComplex p -> VComplex p+sqr = withComplex sqr'++-- | Complex square.+sqr' :: RealFloat r => X.Complex r -> X.Complex r+sqr' (x X.:+ y) = (x + y) * (x - y) X.:+ scaleFloat 1 (x * y)++-- | Much like 'withComplex' 'scaleComplex''.+scaleComplex :: NaturalNumber p => Int -> VComplex p -> VComplex p+scaleComplex = withComplex . scaleComplex'++-- | Much like 'mapComplex' 'scaleFloat'.+scaleComplex' :: RealFloat r => Int -> X.Complex r -> X.Complex r+scaleComplex' = mapComplex . scaleFloat++-- | Real-complex multiplication.+(.*) :: NaturalNumber p => VFloat p -> VComplex p -> VComplex p+x .* y = withComplex (mapComplex (x *)) y+infixl 7 .*++-- | Complex-real multiplication.+(*.) :: NaturalNumber p => VComplex p -> VFloat p -> VComplex p+x *. y = withComplex (mapComplex (* y)) x+infixl 7 *.+++-- | A concrete format suitable for storage or wire transmission.+data DComplex = DComplex{ dRealPart :: !DFloat, dImagPart :: !DFloat }+ deriving (Eq, Ord, Read, Show, Data, Typeable)++-- | Freeze a 'VComplex'.+toDComplex :: NaturalNumber p => VComplex p -> DComplex+toDComplex v = DComplex (toDFloat (realPart v)) (toDFloat (imagPart v))++-- | Thaw a 'DComplex'. Results in 'Nothing' on precision mismatch.+fromDComplex :: NaturalNumber p => DComplex -> Maybe (VComplex p)+fromDComplex d = do+ r <- fromDFloat (dRealPart d)+ i <- fromDFloat (dImagPart d)+ return (r .+ i)++-- | Thaw a 'DComplex' to its natural precision. 'Nothing' is passed on+-- precision mismatch between real and imaginary parts.+withDComplex :: DComplex -> (forall p . NaturalNumber p => Maybe (VComplex p) -> r) -> r+withDComplex d f = withDFloat (dRealPart d) $ \r -> f $ do+ i <- fromDFloat (dImagPart d)+ return (r .+ i)
Numeric/VariablePrecision/Float.hs view
@@ -1,282 +1,594 @@-{-# LANGUAGE BangPatterns, DeriveDataTypeable #-}+{-# LANGUAGE BangPatterns, DeriveDataTypeable, Rank2Types #-} {- | Module : Numeric.VariablePrecision.Float Copyright : (c) Claude Heiland-Allen 2012 License : BSD3 Maintainer : claudiusmaximus@goto10.org-Stability : provisional-Portability : BangPatterns, DeriveDataTypeable+Stability : unstable+Portability : BangPatterns, DeriveDataTypeable, Rank2Types Variable precision software floating point based on @(Integer, Int)@ as-used by 'decodeFloat'.+used by 'decodeFloat'. Supports infinities and NaN, but not negative+zero or denormalization. Accuracy has not been extensively verified, and termination of numerical algorithms has not been proven. -'floatRange' is arbitrarily limited to mitigate the problems that-occur when enormous integers might be needed during some number-type conversions (worst case consequence: program abort in gmp).--No support for infinities, NaNs, negative zero or denormalization:-- * exponent overflow throws an error instead of resulting in infinity,-- * exponent underflow traces a warning and results in zero instead of- resulting in a denormalized number.--Some operations throw errors instead of resulting in an infinity or NaN:-- * @'recip' 0@,-- * @x '/' 0@,-- * @'sqrt' x | x < 0@,-- * @'log' x | x <= 0@.--The 'Floating' instance so far only implements algorithms for:-- * 'pi',-- * 'sqrt',-- * 'exp',-- * 'log'--with other 'Floating' methods transitting via 'Double', also 'log'-precision is limited due to internal use of @log 2 :: Double@.- -} module Numeric.VariablePrecision.Float ( VFloat()- , recodeFloat- , module Numeric.VariablePrecision.Precision- , module TypeLevel.NaturalNumber.ExtraNumbers+ , Normed(norm1, norm2, norm2Squared, normInfinity)+ , effectivePrecisionWith+ , effectivePrecision+ , (-@?)+ , DFloat(..)+ , toDFloat+ , fromDFloat+ , withDFloat ) where import Data.Data (Data()) import Data.Typeable (Typeable()) import Data.Bits (bit, shiftL, shiftR)-import Data.Monoid (mappend) import Data.Ratio ((%), numerator, denominator) import GHC.Float (showSignedFloat) import Numeric (readSigned, readFloat) import Text.FShow.RealFloat (DispFloat(), FShow(fshowsPrec), fshowFloat) -import Debug.Trace (trace) -- FIXME-+import Numeric.VariablePrecision.Algorithms import Numeric.VariablePrecision.Precision-import TypeLevel.NaturalNumber.ExtraNumbers (N24, n24, N53, n53)+import Numeric.VariablePrecision.Precision.Reify+import Numeric.VariablePrecision.Integer.Logarithm + -- | A software implementation of floating point arithmetic, using a strict--- pair of 'Integer' and 'Int', scaled similarly to 'decodeFloat'.-data VFloat p = F !Integer !Int deriving (Data, Typeable)+-- pair of 'Integer' and 'Int', scaled similarly to 'decodeFloat', along+-- with additional values representing:+--+-- * positive infinity (@1/0@),+--+-- * negative infinity (@-1/0@),+--+-- * not a number (@0/0@).+--+-- The 'Floating' instance so far only implements algorithms for:+--+-- * 'pi',+--+-- * 'sqrt',+--+-- * 'exp',+--+-- * 'log'.+--+-- These 'Floating' methods transit via 'Double' and so have limited+-- precision:+--+-- * 'sin', 'cos', 'tan',+--+-- * 'asin', 'acos', 'atan',+--+-- * 'sinh', 'cosh', 'tanh',+--+-- * 'asinh', 'acosh', 'atanh'.+--+-- 'floatRange' is arbitrarily limited to mitigate the problems that+-- occur when enormous integers might be needed during some number+-- type conversions (worst case consequence: program abort in gmp).+--+data VFloat p+ = F !Integer !Int+ -- invariant: matches decodeFloat spec+ -- if unsure, use encodeVFloat which maintains the invariant+ -- if sure, use checkVFloat which checks the invariant+ -- only construct with bare F when absolutely sure+ | FZero -- FIXME add negative zero+ | FPosInf+ | FNegInf+ | FNaN -- FIXME add payload+ deriving (Data, Typeable) --- | Convert between generic 'RealFloat' types--- more efficiently than 'realToFrac'.-recodeFloat :: (RealFloat a, RealFloat b) => a -> b-recodeFloat = uncurry encodeFloat . decodeFloat+encodeVFloat :: NaturalNumber p => VFloat p -> Integer -> Int -> VFloat p+encodeVFloat witness = self+ where+ b = fromIntegral $ precision (undefined `asTypeOf` witness)+ b' = b - 1+ self 0 !_ = FZero+ self m e = checkVFloat "encodeFloat'" $ encodeFloat'' (m > 0) m' (e - sh) l+ where+ absm = abs m+ m' = absm `shift` sh+ e2 = integerLog2 absm+ sh = b - e2+ l = integerLog2 m'+ encodeFloat'' !s' !m' !e' !l+ | m' <= 0 = failed -- FIXME+ | b' == l = F (if s' then m' else negate m') e'+ | b' < l = {-# SCC "encodeFloat''.shiftR" #-} encodeFloat'' s' (m' `shiftR` 1) (e' + 1) (l - 1)+ | b' > l = {-# SCC "encodeFloat''.shiftL" #-} encodeFloat'' s' (m' `shiftL` 1) (e' - 1) (l + 1)+ | otherwise = failed -- FIXME+ where+ failed = error $ "Numeric.VariablePrecision.Float.encodeVFloat: internal error (please report this bug): "+ ++ show (b, b', l, s', m', e') + instance NaturalNumber p => DispFloat (VFloat p) where++ instance NaturalNumber p => FShow (VFloat p) where+ fshowsPrec p = showSignedFloat fshowFloat p++ instance NaturalNumber p => Show (VFloat p) where+ showsPrec = fshowsPrec + instance NaturalNumber p => Read (VFloat p) where- readsPrec _ = readSigned readFloat -- FIXME ignores precedence + readsPrec _ = readSigned readFloat -- FIXME ignores precedence, NaN/Inf fail?++ instance HasPrecision VFloat+++minimumExponent, maximumExponent :: Int+minimumExponent = negate (bit 20)+maximumExponent = bit 20++asTypeIn :: (a -> b) -> a+asTypeIn _ = undefined++asTypeOut :: (a -> b) -> b+asTypeOut _ = undefined++asTypeOut2 :: (a -> b -> c) -> c+asTypeOut2 _ = undefined++ instance VariablePrecision VFloat where- adjustPrecision (F 0 _) = F 0 0- adjustPrecision x@(F m e) = result++ adjustPrecision = self where- result- | n > 0 = checkVFloat (F (m `shiftL` n) (e - n))- | n == 0 = checkVFloat (F m e)- | n < 0 = checkVFloat (F (m `shiftR` negate n) (e + negate n))+ p = asTypeIn self+ q = asTypeOut self+ np = floatDigits p+ nq = floatDigits q n = nq - np- np = precision x- nq = precision result+ self FZero = FZero+ self FPosInf = FPosInf+ self FNegInf = FNegInf+ self FNaN = FNaN+ self (F m e)+ | n > 0 = encodeVFloat q (m `shiftL` n) (e - n)+ | n == 0 = encodeVFloat q m e+ | n < 0 = encodeVFloat q (m `shiftR` negate n) (e + negate n)+ | otherwise = unreachable + instance Eq (VFloat p) where- F 0 _ == F 0 _ = True- F a b == F x y = a == x && b == y- F 0 _ /= F 0 _ = False- F a b /= F x y = a /= x || b /= y + FZero == FZero = True+ F a b == F x y = a == x && b == y+ FPosInf == FPosInf = True+ FNegInf == FNegInf = True+ -- everything else including NaN+ _ == _ = False++ a /= x = not (a == x)++ instance Ord (VFloat p) where- F 0 _ `compare` F x _ = 0 `compare` x- F a _ `compare` F 0 _ = a `compare` 0- F a b `compare` F x y- | a > 0 && x > 0 = (b `compare` y) `mappend` (a `compare` x)- | a > 0 && x < 0 = GT- | a < 0 && x > 0 = LT- | a < 0 && x < 0 = (y `compare` b) `mappend` (a `compare` x) + FZero < FZero = False+ FZero < F x _ = 0 < x+ F a _ < FZero = a < 0+ F a b < F x y+ | a > 0 && x > 0 && b < y = True+ | a > 0 && x > 0 && b == y = a < x+ | a > 0 && x > 0 && b > y = False+ | a > 0 && x < 0 = False+ | a < 0 && x > 0 = True+ | a < 0 && x < 0 && b < y = False+ | a < 0 && x < 0 && b == y = a < x+ | a < 0 && x < 0 && b > y = True+ | otherwise = unreachable+ FNaN < _ = False+ _ < FNaN = False+ FPosInf < _ = False+ _ < FPosInf = True+ _ < FNegInf = False+ FNegInf < _ = True++ a > x = x < a++ a <= x = a < x || a == x++ a >= x = a > x || a == x++ min a@FNaN !_ = a+ min !_ x@FNaN = x+ min a x+ | a <= x = a+ | otherwise = x++ max a@FNaN !_ = a+ max !_ x@FNaN = x+ max a x+ | a >= x = a+ | otherwise = x++ -- 'compare' uses default implementation in Ord++ instance NaturalNumber p => Num (VFloat p) where- F 0 _ + xy = xy- ab + F 0 _ = ab- F a b + F x y- | b > y = checkVFloat $ encodeFloat (a + (x `shiftR` (b - y))) b- | b == y = checkVFloat $ encodeFloat (a + x) b- | b < y = checkVFloat $ encodeFloat ((a `shiftR` (y - b)) + x) y- F 0 _ - xy = checkVFloat $ negate xy- ab - F 0 _ = checkVFloat $ ab- F a b - F x y- | b > y = checkVFloat $ encodeFloat (a - (x `shiftR` (b - y))) b- | b == y = checkVFloat $ encodeFloat (a - x) b- | b < y = checkVFloat $ encodeFloat ((a `shiftR` (y - b)) - x) y- ab@(F 0 _) * _ = checkVFloat $ ab- _ * xy@(F 0 _) = checkVFloat $ xy- ab@(F a b) * F x y = checkVFloat $ encodeFloat ((a * x) `shiftR` (k - 2)) (b + y + k - 2)- where k = precision ab- negate (F a b) = checkVFloat $ F (negate a) b- abs (F a b) = checkVFloat $ F (abs a) b- signum (F a _) = fromInteger (signum a)- fromInteger i = checkVFloat $ encodeFloat i 0 + f@(F a b) + F x y+ | b > y = encodeVFloat f (a + (x `shiftR` (b - y))) b+ | b == y = encodeVFloat f (a + x) b+ | b < y = encodeVFloat f ((a `shiftR` (y - b)) + x) y+ | otherwise = unreachable+ a@FNaN + _ = a+ _ + x@FNaN = x+ FZero + x = x+ a + FZero = a+ FPosInf + FNegInf = FNaN+ FNegInf + FPosInf = FNaN+ FPosInf + _ = FPosInf+ _ + FPosInf = FPosInf+ FNegInf + _ = FNegInf+ _ + FNegInf = FNegInf++ f@(F a b) - F x y+ | b > y = encodeVFloat f (a - (x `shiftR` (b - y))) b+ | b == y = encodeVFloat f (a - x) b+ | b < y = encodeVFloat f ((a `shiftR` (y - b)) - x) y+ | otherwise = unreachable+ a@FNaN - _ = a+ _ - x@FNaN = x+ FZero - x = negate x+ a - FZero = a+ FPosInf - FPosInf = FNaN+ FNegInf - FNegInf = FNaN+ FPosInf - _ = FPosInf+ _ - FPosInf = FNegInf+ FNegInf - _ = FNegInf+ _ - FNegInf = FPosInf++ negate (F a b) = checkVFloat "negate" $ F (negate a) b+ negate FZero = FZero+ negate FPosInf = FNegInf+ negate FNegInf = FPosInf+ negate a@FNaN = a++ abs !a+ | a < 0 = negate a+ | otherwise = a++ signum !a+ | a < 0 = -1+ | a > 0 = 1+ | otherwise = a++ f@(F a b) * F x y = encodeVFloat f ((a * x) `shiftR` (k - 1)) (b + y + k - 1) where k = fromIntegral $ precision f+ a@FNaN * _ = a+ _ * x@FNaN = x+ FZero * FPosInf = FNaN+ FZero * FNegInf = FNaN+ FZero * _ = FZero+ FPosInf * FZero = FNaN+ FNegInf * FZero = FNaN+ _ * FZero = FZero+ a * x+ | sameSign a x = FPosInf+ | otherwise = FNegInf++ fromInteger !i = encodeFloat i 0++ instance NaturalNumber p => Real (VFloat p) where- toRational (F 0 _) = 0++ toRational FZero = 0 toRational (F m e) | e > 0 = fromInteger (m `shiftL` e) | e == 0 = fromInteger m | e < 0 = m % bit (negate e)+ | otherwise = unreachable+ toRational FPosInf = 1 % 0+ toRational FNegInf = (-1) % 0+ toRational FNaN = 0 % 0 + instance NaturalNumber p => Fractional (VFloat p) where- _ / (F 0 _) = error "Numeric.VFloat./0" -- FIMXE- ab@(F 0 _) / _ = checkVFloat $ ab- ab@(F a b) / (F x y) = checkVFloat $ encodeFloat ((a `shiftL` (k + 2)) `quot` x) (b - y - k - 2) -- FIXME accuracy- where k = precision ab- recip (F 0 _) = error "Numeric.VFloat.recip 0" -- FIXME- recip xy@(F x y) = checkVFloat $ encodeFloat (bit (2 * k + 2) `quot` x) (negate y - 2 * k - 2) -- FIXME accuracy- where k = precision xy- fromRational r = checkVFloat $ fromInteger (numerator r) / fromInteger (denominator r) -- FIXME accuracy + f@(F _ _) / g@(F _ _) = f * recip g+ a@FNaN / _ = a+ _ / x@FNaN = x+ FPosInf / FPosInf = FNaN+ FPosInf / FNegInf = FNaN+ FNegInf / FPosInf = FNaN+ FNegInf / FNegInf = FNaN+ _ / FPosInf = FZero+ _ / FNegInf = FZero+ a / FZero+ | a > 0 = FPosInf+ | a < 0 = FNegInf+ | otherwise = FNaN+ FZero / _ = FZero+ a / x+ | a `sameSign` x = FPosInf+ | otherwise = FNegInf++ recip a@FNaN = a+ recip FZero = FPosInf+ recip FPosInf = FZero+ recip FNegInf = FZero+ recip f@(F m e) = encodeVFloat f (bit k `quot` m) (negate (k + e)) where k = 2 * fromIntegral (precision f)++ fromRational r = fromInteger (numerator r) / fromInteger (denominator r) -- FIXME accuracy++ instance NaturalNumber p => RealFrac (VFloat p) where- properFraction (F 0 _) = (0, checkVFloat $ 0)- properFraction me@(F m e)- | e >= 0 = (fromInteger m, checkVFloat $ 0)- | e < negate (precision me) = (0, checkVFloat $ me)- | otherwise = (fromInteger n', checkVFloat $ f')++ properFraction = self where- n = m `shiftR` (negate e)- d = F (n `shiftL` (negate e)) e- f = me - d- (n', f')- | (m >= 0) == (f >= 0) = (n, f)- | otherwise = (n + 1, f - 1)+ p = fromIntegral $ precision (asTypeIn self)+ self FZero = (0, FZero)+ self me@(F m e)+ | e >= 0 = (fromInteger m, FZero)+ | e < negate p = (0, me)+ | otherwise = (fromInteger n', f')+ where+ n = m `shiftR` (negate e)+ d = checkVFloat "properFraction" $ F (n `shiftL` (negate e)) e+ f = me - d+ (n', f')+ | (m >= 0) == (f >= 0) = (n, f)+ | otherwise = (n + 1, f - 1)+ self f = (error $ "Numeric.VariablePrecision.Float.properFraction: not finite: " ++ show f, f) + -- 'truncate' uses default implementation in RealFrac++ -- 'floor' uses default implementation in RealFrac++ -- 'ceiling' uses default implementation in RealFrac++ -- 'round' uses default implementation in RealFrac++ instance NaturalNumber p => RealFloat (VFloat p) where+ floatRadix _ = 2- floatDigits = precision- floatRange _ = (negate (bit 20), bit 20) -- FIXME++ floatDigits = self+ where+ prec = fromIntegral $ precision (asTypeIn self)+ self = const prec++ floatRange = const (minimumExponent, maximumExponent) -- FIXME arbitrary -- this floatRange is somewhat arbitrary, but toInteger gives integers -- with up to around (precision + maxExponent) bits, the value here -- gives rise to potentially more than 300k decimal digits...- isNaN _ = False- isInfinite _ = False++ isNaN FNaN = True+ isNaN _ = False++ isInfinite FPosInf = True+ isInfinite FNegInf = True+ isInfinite _ = False+ isDenormalized _ = False+ isNegativeZero _ = False- isIEEE _ = False- decodeFloat (F 0 _) = (0, 0)- decodeFloat (F m e) = (m, e)- encodeFloat 0 _ = F 0 0- encodeFloat m e = result- where- result = checkVFloat $ encodeFloat' (signum m) (abs m) e- b = precision result- hi = bit (b + 1)- lo = bit b- encodeFloat' !s' !m' !e'- | m' <= 0 = failed -- FIXME- | lo <= m' && m' < hi = F (s' * (m' `shiftR` 1)) (e' + 1)- | m' < lo = encodeFloat' s' (m' `shiftL` 1) (e' - 1)- | hi <= m' = encodeFloat' s' (m' `shiftR` 1) (e' + 1)- | otherwise = failed -- FIXME- where- failed = error $ "Numeric.VariablePrecision.VFloat.encodeFloat\n"- ++ show (m, e, b, lo, hi, s', m', e')- ++ "\nplease report this as a bug." -instance NaturalNumber p => Floating (VFloat p) where -- FIXME+ isIEEE _ = False -- FIXME what does this mean? - -- <http://en.wikipedia.org/wiki/AGM_method>- pi = checkVFloat $ go 1 (sqrt 0.5) 1 2 0+ decodeFloat FZero = (0, 0)+ decodeFloat (F m e) = (m, e)+ decodeFloat f = error $ "Numeric.VariablePrecision.Float.decodeFloat: not finite: " ++ show f++ encodeFloat = self where- go a b s k p- | p == p' = p'- | otherwise = go a' b' s' k' p'- where- a' = (a + b) / 2- b' = sqrt (a * b)- c = (a - b) / 2- s' = s - k' * c * c- k' = 2 * k- p' = 4 * a' * a' / s+ self = encodeVFloat (undefined `asTypeOf` asTypeOut2 self) - -- Newton's method- sqrt f- | 0 == f = F 0 0- | 0 < f = checkVFloat $ go 1+ exponent = self where- go !r =- let r' = (r + f / r) / 2- in if r == r' then r else go r'+ prec = fromIntegral $ precision (asTypeIn self)+ self FZero = 0+ self (F _ e) = e + prec+ self f = error $ "Numeric.VariablePrecision.Float.exponent: not finite: " ++ show f - -- power series- exp f = checkVFloat $ go 0 1 1 1+ significand = self where- go !e !nf !fn !n =- let e' = e + fn / nf- in if e == e' then e else go e' (nf * n) (f * fn) (n + 1)+ prec = fromIntegral $ precision (asTypeIn self)+ e = negate prec+ self (F m _) = checkVFloat "significand" $ F m e+ self f = f - -- <http://en.wikipedia.org/wiki/Logarithm#Arithmetic-geometric_mean_approximation>- log f@(F _ e)- | f > 0 = checkVFloat $ pi / (2 * agm 1 (encodeFloat 1 (2 - m) / f)) - fromIntegral m * ln2+ scaleFloat n (F m e) = checkVFloat "scaleFloat" $ F m (e + n)+ scaleFloat _ f = f++ -- 'atan2' uses default implementation in RealFloat+++shift :: Integer -> Int -> Integer+shift !n !k+ | k > 0 = n `shiftL` k+ | k == 0 = n+ | k < 0 = n `shiftR` (negate k)+ | otherwise = unreachable+++instance NaturalNumber p => Floating (VFloat p) where -- FIXME++ pi = genericPi 2++ sqrt = genericSqrt 2++ exp = genericExp 2++ log = self where- p = precision f- -- f ~= sqrt 2 * 2^(p + e)- -- f * 2^m > (sqrt 2) ^ p- -- sqrt 2 * 2 ^ (p + e) * 2 ^ m > sqrt 2 ^ p- -- 1/2 + p + e + m > p / 2- -- 1 + p + 2 e + 2 m > 0- m = negate $ p `div` 2 + e- agm !a! b =- let a' = (a + b) / 2- b' = sqrt (a * b)- in if a' == b' || (a == a' && b == b') then a' else agm a' b'- ln2 = viaDouble log 2 -- FIXME+ log2 = genericLog2 2+ self = genericLog' 2 log2 + -- '(**)' uses default implementation in Floating++ -- 'logBase' uses default implementation in Floating+ sin = viaDouble sin -- FIXME+ cos = viaDouble cos -- FIXME+ tan = viaDouble tan -- FIXME+ sinh = viaDouble sinh -- FIXME+ cosh = viaDouble cosh -- FIXME+ tanh = viaDouble tanh -- FIXME+ asin = viaDouble asin -- FIXME+ acos = viaDouble acos -- FIXME+ atan = viaDouble atan -- FIXME+ asinh = viaDouble asinh -- FIXME+ acosh = viaDouble acosh -- FIXME+ atanh = viaDouble atanh -- FIXME -viaDouble :: NaturalNumber p => (Double -> Double) -> (VFloat p -> VFloat p)-viaDouble f = recodeFloat . checkDouble . f . recodeFloat -checkDouble :: Double -> Double-checkDouble f- | isNaN f = error "Numeric.VariablePrecision.Float: isNaN" -- FIXME- | isInfinite f = error "Numeric.VariablePrecision.Float: isInfinite" -- FIXME- | otherwise = f+-- despite the name, using this is vital for correct behaviour+-- because it properly handles underflow and overflow as well as+-- checking that the invariant for F holds+checkVFloat :: NaturalNumber p => String -> VFloat p -> VFloat p+checkVFloat = self+ where+ prec = fromIntegral $ precision (asTypeOut2 self)+ prec' = prec - 1+ elo = minimumExponent+ ehi = maximumExponent+ self s x@(F m e)+ | not mok = error $ "Numeric.VariablePrecision.Float.checkVFloat." ++ s ++ ": internal error (please report this bug): " ++ show ((m, am, lm, prec, prec', mok), (elo, e, ehi, eok))+ | eok = x+ | e < elo = FZero -- underflow+ | m > 0 = FPosInf -- overflow+ | m < 0 = FNegInf -- overflow+ | otherwise = unreachable+ where+ eok = elo <= e && e <= ehi+ mok = lm == prec'+ lm = integerLog2 am+ am = abs m+ self _ x = x -checkVFloat :: NaturalNumber p => VFloat p -> VFloat p-checkVFloat x@(F _ e)- | lo <= e && e <= hi = x- | e < lo = trace ("Numeric.VariablePrecision.Float underflow: " ++ show x) 0 -- FIXME- | otherwise = error ("Numeric.VariablePrecision.Float overflow: " ++ show x) -- FIXME- where (lo, hi) = floatRange x++-- | A selection of norms.+class HasPrecision t => Normed t where+ norm1 :: NaturalNumber p => t p -> VFloat p+ norm2 :: NaturalNumber p => t p -> VFloat p+ norm2Squared :: NaturalNumber p => t p -> VFloat p+ normInfinity :: NaturalNumber p => t p -> VFloat p+++instance Normed VFloat where+ norm1 = abs+ norm2 = abs+ norm2Squared x = x * x+ normInfinity = abs+++-- | A measure of meaningful precision in the difference of two+-- finite non-zero values.+--+-- Values of very different magnitude have little meaningful+-- difference, because @a + b `approxEq` a@ when @|a| >> |b|@.+--+-- Very close values have little meaningful difference,+-- because @a + (a - b) `approxEq` a@ as @|a| >> |a - b|@.+--+-- 'effectivePrecisionWith' attempts to quantify this.+--+effectivePrecisionWith :: (Num t, RealFloat r) => (t -> r) {- ^ norm -} -> t -> t -> Int+effectivePrecisionWith n i j+ | t a && t b && t c = p - (d `max` (e - d))+ | otherwise = 0+ where+ t k = k > 0 && not (isInfinite k)+ d = (x `max` y) - z+ e = abs (x - y) `min` p+ p = floatDigits a+ x = exponent a+ y = exponent b+ z = exponent c+ a = n i+ b = n j+ c = n (i - j)+++-- | Much like 'effectivePrecisionWith' combined with 'normInfinity'.+effectivePrecision :: (NaturalNumber p, HasPrecision t, Normed t, Num (t p)) => t p -> t p -> Int+effectivePrecision = effectivePrecisionWith normInfinity+infix 6 `effectivePrecision`+++-- | An alias for 'effectivePrecision'.+(-@?) :: (NaturalNumber p, HasPrecision t, Normed t, Num (t p)) => t p -> t p -> Int+(-@?) = effectivePrecision+infix 6 -@?+++unreachable :: a+unreachable = error "Numeric.VariablePrecision.Float: internal error (please report this bug): unreachable code was reached"+++-- | A concrete format suitable for storage or wire transmission.+data DFloat+ = DFloat { dPrecision :: !Word, dMantissa :: !Integer, dExponent :: !Int }+ | DZero { dPrecision :: !Word }+ | DPositiveInfinity { dPrecision :: !Word }+ | DNegativeInfinity { dPrecision :: !Word }+ | DNotANumber { dPrecision :: !Word }+ deriving (Eq, Ord, Read, Show, Data, Typeable)++-- | Freeze a 'VFloat'.+toDFloat :: NaturalNumber p => VFloat p -> DFloat+toDFloat f@(F m e) = DFloat (precision f) m e+toDFloat f@FZero = DZero (precision f)+toDFloat f@FPosInf = DPositiveInfinity (precision f)+toDFloat f@FNegInf = DNegativeInfinity (precision f)+toDFloat f@FNaN = DNotANumber (precision f)++-- | Thaw a 'DFloat'. Results in 'Nothing' on precision mismatch.+fromDFloat :: NaturalNumber p => DFloat -> Maybe (VFloat p)+fromDFloat d+ | dPrecision d == precision result = Just result+ | otherwise = Nothing+ where+ result = case d of+ DFloat _ m e -> encodeVFloat undefined m e+ DZero _ -> FZero+ DPositiveInfinity _ -> FPosInf+ DNegativeInfinity _ -> FNegInf+ DNotANumber _ -> FNaN++-- | Thaw a 'DFloat' to its natural precision.+withDFloat :: DFloat -> (forall p . NaturalNumber p => VFloat p -> r) -> r+withDFloat (DFloat p m e) f = reifyPrecision p $ \prec -> f (encodeVFloat undefined m e `atPrecision` prec)+withDFloat d f = unsafeWithDFloat d f++-- | Thaw a 'DFloat' without guaranteeing a well-formed 'VFloat' value.+-- Possibly slightly faster.+unsafeWithDFloat :: DFloat -> (forall p . NaturalNumber p => VFloat p -> r) -> r+unsafeWithDFloat (DFloat p m e) f = reifyPrecision p $ \prec -> f (F m e `atPrecision` prec)+unsafeWithDFloat (DZero p) f = reifyPrecision p $ \prec -> f (FZero `atPrecision` prec)+unsafeWithDFloat (DPositiveInfinity p) f = reifyPrecision p $ \prec -> f (FPosInf `atPrecision` prec)+unsafeWithDFloat (DNegativeInfinity p) f = reifyPrecision p $ \prec -> f (FNegInf `atPrecision` prec)+unsafeWithDFloat (DNotANumber p) f = reifyPrecision p $ \prec -> f (FNaN `atPrecision` prec)
− Numeric/VariablePrecision/Float/Aliases.hs
@@ -1,95 +0,0 @@-{-|-Module : Numeric.VariablePrecision.Float.Aliases-Copyright : (c) Claude Heiland-Allen 2012-License : BSD3--Maintainer : claudiusmaximus@goto10.org-Stability : stable-Portability : portable--Boilerplate definitions generated by:--> flip mapM_ [1..53] $ \p -> let s = show p in-> putStrLn $ "type F" ++ s ++ " = VFloat N" ++ s ++-> " ; f" ++ s ++ " :: F" ++ s ++ " ; f" ++ s ++ " = 0"--Along with aliases for 'recodeFloat' with specialized types.--Using this module in ghc-7.0.4 might require @-fcontext-stack=100@.---}-module Numeric.VariablePrecision.Float.Aliases where--import Numeric.VariablePrecision.Float (VFloat, recodeFloat)-import TypeLevel.NaturalNumber-import TypeLevel.NaturalNumber.ExtraNumbers---- | Convert to a Float from the same precision.-toFloat :: F24 -> Float-toFloat = recodeFloat---- | Convert from a Float to the same precision.-fromFloat :: Float -> F24-fromFloat = recodeFloat---- | Convert to a Double from the same precision.-toDouble :: F53 -> Double-toDouble = recodeFloat---- | Convert from a Double to the same precision.-fromDouble :: Double -> F53-fromDouble = recodeFloat--type F1 = VFloat N1 ; f1 :: F1 ; f1 = 0-type F2 = VFloat N2 ; f2 :: F2 ; f2 = 0-type F3 = VFloat N3 ; f3 :: F3 ; f3 = 0-type F4 = VFloat N4 ; f4 :: F4 ; f4 = 0-type F5 = VFloat N5 ; f5 :: F5 ; f5 = 0-type F6 = VFloat N6 ; f6 :: F6 ; f6 = 0-type F7 = VFloat N7 ; f7 :: F7 ; f7 = 0-type F8 = VFloat N8 ; f8 :: F8 ; f8 = 0-type F9 = VFloat N9 ; f9 :: F9 ; f9 = 0-type F10 = VFloat N10 ; f10 :: F10 ; f10 = 0-type F11 = VFloat N11 ; f11 :: F11 ; f11 = 0-type F12 = VFloat N12 ; f12 :: F12 ; f12 = 0-type F13 = VFloat N13 ; f13 :: F13 ; f13 = 0-type F14 = VFloat N14 ; f14 :: F14 ; f14 = 0-type F15 = VFloat N15 ; f15 :: F15 ; f15 = 0-type F16 = VFloat N16 ; f16 :: F16 ; f16 = 0-type F17 = VFloat N17 ; f17 :: F17 ; f17 = 0-type F18 = VFloat N18 ; f18 :: F18 ; f18 = 0-type F19 = VFloat N19 ; f19 :: F19 ; f19 = 0-type F20 = VFloat N20 ; f20 :: F20 ; f20 = 0-type F21 = VFloat N21 ; f21 :: F21 ; f21 = 0-type F22 = VFloat N22 ; f22 :: F22 ; f22 = 0-type F23 = VFloat N23 ; f23 :: F23 ; f23 = 0-type F24 = VFloat N24 ; f24 :: F24 ; f24 = 0-type F25 = VFloat N25 ; f25 :: F25 ; f25 = 0-type F26 = VFloat N26 ; f26 :: F26 ; f26 = 0-type F27 = VFloat N27 ; f27 :: F27 ; f27 = 0-type F28 = VFloat N28 ; f28 :: F28 ; f28 = 0-type F29 = VFloat N29 ; f29 :: F29 ; f29 = 0-type F30 = VFloat N30 ; f30 :: F30 ; f30 = 0-type F31 = VFloat N31 ; f31 :: F31 ; f31 = 0-type F32 = VFloat N32 ; f32 :: F32 ; f32 = 0-type F33 = VFloat N33 ; f33 :: F33 ; f33 = 0-type F34 = VFloat N34 ; f34 :: F34 ; f34 = 0-type F35 = VFloat N35 ; f35 :: F35 ; f35 = 0-type F36 = VFloat N36 ; f36 :: F36 ; f36 = 0-type F37 = VFloat N37 ; f37 :: F37 ; f37 = 0-type F38 = VFloat N38 ; f38 :: F38 ; f38 = 0-type F39 = VFloat N39 ; f39 :: F39 ; f39 = 0-type F40 = VFloat N40 ; f40 :: F40 ; f40 = 0-type F41 = VFloat N41 ; f41 :: F41 ; f41 = 0-type F42 = VFloat N42 ; f42 :: F42 ; f42 = 0-type F43 = VFloat N43 ; f43 :: F43 ; f43 = 0-type F44 = VFloat N44 ; f44 :: F44 ; f44 = 0-type F45 = VFloat N45 ; f45 :: F45 ; f45 = 0-type F46 = VFloat N46 ; f46 :: F46 ; f46 = 0-type F47 = VFloat N47 ; f47 :: F47 ; f47 = 0-type F48 = VFloat N48 ; f48 :: F48 ; f48 = 0-type F49 = VFloat N49 ; f49 :: F49 ; f49 = 0-type F50 = VFloat N50 ; f50 :: F50 ; f50 = 0-type F51 = VFloat N51 ; f51 :: F51 ; f51 = 0-type F52 = VFloat N52 ; f52 :: F52 ; f52 = 0-type F53 = VFloat N53 ; f53 :: F53 ; f53 = 0
Numeric/VariablePrecision/Precision.hs view
@@ -4,7 +4,7 @@ License : BSD3 Maintainer : claudiusmaximus@goto10.org-Stability : provisional+Stability : unstable Portability : portable Classes for types with precision represented by a type-level natural@@ -21,67 +21,78 @@ , atPrecisionOf , (.@) , VariablePrecision(adjustPrecision)+ , auto , withPrecision , withPrecisionOf , (.@~) , module TypeLevel.NaturalNumber+ , module Data.Word ) where import TypeLevel.NaturalNumber+ ( NaturalNumber(..), Zero, SuccessorTo, n0, successorTo )+import Data.Word (Word) -- | A class for types with precision.+-- The methods must not evaluate their arguments, and their results+-- must not be evaluated. -- Minimal complete definition: (none). class HasPrecision t where- -- | Get the precision of a value. 'precisionOf' must not evaluate- -- its argument, and its result must not be evaluated. precisionOf :: NaturalNumber p => t p -> p- precisionOf _ = error "Numeric.VariablePrecision.Precision.HasPrecision.precisionOf: result evaluated"+ precisionOf _ = undefined + -- | Much like 'naturalNumberAsInt' combined with 'precisionOf'.-precision :: (HasPrecision t, NaturalNumber p) => t p -> Int-precision = naturalNumberAsInt . precisionOf+precision :: (NaturalNumber p, HasPrecision t) => t p -> Word+precision = fromIntegral . naturalNumberAsInt . precisionOf + -- | Much like 'const' with a restricted type.-atPrecision :: (HasPrecision t, NaturalNumber p) => t p -> p -> t p+atPrecision :: (NaturalNumber p, HasPrecision t) => t p -> p -> t p atPrecision = const + -- | Much like 'const' with a restricted type.-atPrecisionOf- :: (HasPrecision t, HasPrecision s, NaturalNumber p)- => t p -> s p -> t p+-- Precedence between '<' and '+'.+atPrecisionOf :: (HasPrecision t, HasPrecision s) => t p -> s p -> t p atPrecisionOf = const-infixl 5 `atPrecisionOf` -- precedence between Prelude.< and Prelude.++-- where _ = precisionOf t `asTypeOf` precisionOf s+infixl 5 `atPrecisionOf` + -- | An alias for 'atPrecisionOf'.-(.@)- :: (HasPrecision t, HasPrecision s, NaturalNumber p)- => t p -> s p -> t p+-- Precedence between '<' and '+'.+(.@) :: (HasPrecision t , HasPrecision s) => t p -> s p -> t p (.@) = atPrecisionOf-infixl 5 .@ -- precedence between Prelude.< and Prelude.++infixl 5 .@ --- | A class for types with variable precision.--- Minimal complete definition: (all).++-- | A class for types with adjustable precision.+-- Minimal complete definition: 'adjustPrecision'. class HasPrecision t => VariablePrecision t where -- | Adjust the precision of a value preserving as much accuracy as -- possible. adjustPrecision :: (NaturalNumber p, NaturalNumber q) => t p -> t q ++-- | Synonym for 'adjustPrecision'.+auto :: (VariablePrecision t, NaturalNumber p, NaturalNumber q) => t p -> t q+auto = adjustPrecision++ -- | Much like 'adjustPrecision' combined with 'atPrecision'.-withPrecision- :: (VariablePrecision t, NaturalNumber p, NaturalNumber q)- => t p -> q -> t q-withPrecision x q = adjustPrecision x `atPrecision` q+withPrecision :: (NaturalNumber p, NaturalNumber q, VariablePrecision t) => t p -> q -> t q+withPrecision s q = adjustPrecision s `atPrecision` q -- | Much like 'withPrecision' combined with 'precisionOf'.-withPrecisionOf- :: (VariablePrecision t, HasPrecision s, NaturalNumber p, NaturalNumber q)- => t p -> s q -> t q-withPrecisionOf x w = x `withPrecision` precisionOf w-infixl 5 `withPrecisionOf` -- precedence between Prelude.< and Prelude.++-- Precedence between '<' and '+'.+withPrecisionOf :: (NaturalNumber p, NaturalNumber q, VariablePrecision t, HasPrecision s) => t p -> s q -> t q+withPrecisionOf s w = s `withPrecision` precisionOf w+infixl 5 `withPrecisionOf` + -- | An alias for 'withPrecisionOf'.-(.@~)- :: (VariablePrecision t, HasPrecision s, NaturalNumber p, NaturalNumber q)- => t p -> s q -> t q+-- Precedence between '<' and '+'.+(.@~) :: (NaturalNumber p, NaturalNumber q, VariablePrecision t, HasPrecision s) => t p -> s q -> t q (.@~) = withPrecisionOf-infixl 5 .@~ -- precedence between Prelude.< and Prelude.++infixl 5 .@~
Numeric/VariablePrecision/Precision/Reify.hs view
@@ -5,7 +5,7 @@ License : BSD3 Maintainer : claudiusmaximus@goto10.org-Stability : provisional+Stability : unstable Portability : Rank2Types Reify from value-level to type-level using Rank2Types.@@ -18,36 +18,38 @@ , (.@$) ) where -import Numeric.VariablePrecision.Precision (NaturalNumber, n0, successorTo, VariablePrecision, withPrecision)+import Numeric.VariablePrecision.Precision+ ( VariablePrecision, withPrecision, Word+ , NaturalNumber, n0, successorTo+ ) -- | Reify a precision from value-level to type-level.-reifyPrecision :: Int -> (forall p . NaturalNumber p => p -> a) -> a+reifyPrecision :: Word -> (forall p . NaturalNumber p => p -> a) -> a -- Implemented as described in an email from Gregory Grosswhite -- <http://markmail.org/message/55iuty6axeljj2do> reifyPrecision = go n0 where- go :: NaturalNumber q => q -> Int -> (forall p . NaturalNumber p => p -> a) -> a+ go :: NaturalNumber q => q -> Word -> (forall p . NaturalNumber p => p -> a) -> a go n i f- | i < 0 = error $ "Numeric.VariablePrecision.Precision.Reify.reifyPrecision: negative argument: " ++ show i | i == 0 = f n- | i > 0 = go (successorTo n) (i - 1) f+ | otherwise = go (successorTo n) (i - 1) f -- | Much like 'reifyPrecision' combined with 'withPrecision'. withReifiedPrecision :: (VariablePrecision t, NaturalNumber p) => t p {-^ original value -}- -> Int {- ^ new precision -}+ -> Word {- ^ new precision -} -> (forall q. NaturalNumber q => t q -> a) {-^ operation -} -> a withReifiedPrecision x i f = reifyPrecision i (f . withPrecision x)-infixl 1 `withReifiedPrecision` -- same fixity as Prelude.$+infixl 1 `withReifiedPrecision` -- | An alias for 'withReifiedPrecision'. (.@$) :: (VariablePrecision t, NaturalNumber p) => t p {-^ original value -}- -> Int {- ^ new precision -}+ -> Word {- ^ new precision -} -> (forall q. NaturalNumber q => t q -> a) {-^ operation -} -> a (.@$) = withReifiedPrecision-infix 1 .@$ -- same fixity as Prelude.$+infixl 1 .@$
TODO view
@@ -5,28 +5,15 @@ examples of usage in documentation Numeric.VariablePrecision.Float- improve accuracy of log (bad: log 2 :: Double) readsPrec ignores precedence- proper exception types for- division by zero- sqrt negative- log nonpositive- floatRange overflow- reconsider warning on floatRange underflow- check accuracy of (/), recip check accuracy of fromRational check accuracy of (+), (-), (*)+ check accuracy of (/), recip proper (a)(sin|cos|tan)(h) (bad: viaDouble) profile and optimise space and time consider rounding modes consider mixed-precision operations- consider IEEE inf/nan/-0 semantics--Numeric.VariablePrecision.Complex.Aliases- if it proves necessary for convenience--Numeric.VariablePrecision.Algorithms- effective precision of difference+ consider IEEE -0 semantics TypeLevel.NaturalNumber.ExtraNumbers submit upstream and remove if accepted
TypeLevel/NaturalNumber/ExtraNumbers.hs view
@@ -13,8 +13,6 @@ > putStrLn $ "type N" ++ s ++ " = SuccessorTo N" ++ show (p - 1) ++ > " ; n" ++ s ++ " :: N" ++ s ++ " ; n" ++ s ++ " = undefined" -Using this module in ghc-7.0.4 might require @-fcontext-stack=100@.- -} module TypeLevel.NaturalNumber.ExtraNumbers where
+ fast/Numeric/VariablePrecision/Integer/Logarithm.hs view
@@ -0,0 +1,8 @@+{-# LANGUAGE MagicHash #-}+module Numeric.VariablePrecision.Integer.Logarithm where++import GHC.Exts (Int(I#))+import GHC.Integer.Logarithms (integerLog2#)++integerLog2 :: Integer -> Int+integerLog2 n = I# (integerLog2# n)
+ pure/Numeric/VariablePrecision/Integer/Logarithm.hs view
@@ -0,0 +1,13 @@+{-# LANGUAGE BangPatterns #-}+module Numeric.VariablePrecision.Integer.Logarithm where++import Data.Bits (shiftL)++integerLog2 :: Integer -> Int+integerLog2 n+ | n > 0 = go (-1) 1+ | otherwise = error $ "integerLog2: non-positive argument: " ++ show n+ where+ go !l !b+ | n < b = l+ | otherwise = go (l + 1) (b `shiftL` 1)
variable-precision.cabal view
@@ -1,5 +1,5 @@ Name: variable-precision-Version: 0.1.1+Version: 0.2 Synopsis: variable-precision floating point Description: Software floating point with type-tagged variable mantissa precision,@@ -13,6 +13,8 @@ The intention with this library is to be relatively simple but still useful, refer to the documentation for caveats concerning accuracy and assorted ill-behaviour.+ .+ Usage with ghc(i)-7.0.4 might require @-fcontext-stack=100@. Homepage: https://gitorious.org/variable-precision License: BSD3@@ -25,21 +27,40 @@ Cabal-version: >=1.6 -Extra-source-files: CHANGES README THANKS TODO+Extra-source-files:+ CHANGES+ README+ THANKS+ TODO+ pure/Numeric/VariablePrecision/Integer/Logarithm.hs+ fast/Numeric/VariablePrecision/Integer/Logarithm.hs +Flag fast+ Description: Enable optimisations requiring recent integer-gmp+ Default: True+ Library Exposed-modules:- Numeric.VariablePrecision.Complex+ Numeric.VariablePrecision+ Numeric.VariablePrecision.Algorithms Numeric.VariablePrecision.Float- Numeric.VariablePrecision.Float.Aliases+ Numeric.VariablePrecision.Complex Numeric.VariablePrecision.Precision Numeric.VariablePrecision.Precision.Reify+ Numeric.VariablePrecision.Aliases TypeLevel.NaturalNumber.ExtraNumbers+ Other-modules:+ Numeric.VariablePrecision.Integer.Logarithm Build-depends: base >= 3 && < 6, floatshow >= 0.2 && < 0.3, type-level-natural-number >= 1 && < 2- GHC-Options: -Wall -fno-warn-incomplete-patterns -fcontext-stack=100+ if (!flag(fast))+ HS-source-dirs: . pure+ if ( flag(fast))+ HS-source-dirs: . fast+ Build-depends: integer-gmp >= 0.4+ GHC-Options: -Wall -fcontext-stack=100 GHC-Prof-Options: -prof -auto-all -caf-all source-repository head@@ -49,4 +70,4 @@ source-repository this type: git location: git://gitorious.org/variable-precision/variable-precision.git- tag: v0.1.1+ tag: v0.2