diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -1,2 +1,3 @@ +v0.2 generic IEEE-ish v0.1.1 fixed for ghc-7.0.4 v0.1 initial release diff --git a/Numeric/VariablePrecision.hs b/Numeric/VariablePrecision.hs new file mode 100644 --- /dev/null +++ b/Numeric/VariablePrecision.hs @@ -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 diff --git a/Numeric/VariablePrecision/Algorithms.hs b/Numeric/VariablePrecision/Algorithms.hs new file mode 100644 --- /dev/null +++ b/Numeric/VariablePrecision/Algorithms.hs @@ -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 +-- . +-- +-- 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, +-- . +-- +-- 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), +-- . +-- +-- 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 . +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 diff --git a/Numeric/VariablePrecision/Aliases.hs b/Numeric/VariablePrecision/Aliases.hs new file mode 100644 --- /dev/null +++ b/Numeric/VariablePrecision/Aliases.hs @@ -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 diff --git a/Numeric/VariablePrecision/Complex.hs b/Numeric/VariablePrecision/Complex.hs --- a/Numeric/VariablePrecision/Complex.hs +++ b/Numeric/VariablePrecision/Complex.hs @@ -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) diff --git a/Numeric/VariablePrecision/Float.hs b/Numeric/VariablePrecision/Float.hs --- a/Numeric/VariablePrecision/Float.hs +++ b/Numeric/VariablePrecision/Float.hs @@ -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? - -- - 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 - -- - 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) diff --git a/Numeric/VariablePrecision/Float/Aliases.hs b/Numeric/VariablePrecision/Float/Aliases.hs deleted file mode 100644 --- a/Numeric/VariablePrecision/Float/Aliases.hs +++ /dev/null @@ -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 diff --git a/Numeric/VariablePrecision/Precision.hs b/Numeric/VariablePrecision/Precision.hs --- a/Numeric/VariablePrecision/Precision.hs +++ b/Numeric/VariablePrecision/Precision.hs @@ -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 .@~ diff --git a/Numeric/VariablePrecision/Precision/Reify.hs b/Numeric/VariablePrecision/Precision/Reify.hs --- a/Numeric/VariablePrecision/Precision/Reify.hs +++ b/Numeric/VariablePrecision/Precision/Reify.hs @@ -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 -- 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 .@$ diff --git a/TODO b/TODO --- a/TODO +++ b/TODO @@ -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 diff --git a/TypeLevel/NaturalNumber/ExtraNumbers.hs b/TypeLevel/NaturalNumber/ExtraNumbers.hs --- a/TypeLevel/NaturalNumber/ExtraNumbers.hs +++ b/TypeLevel/NaturalNumber/ExtraNumbers.hs @@ -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 diff --git a/fast/Numeric/VariablePrecision/Integer/Logarithm.hs b/fast/Numeric/VariablePrecision/Integer/Logarithm.hs new file mode 100644 --- /dev/null +++ b/fast/Numeric/VariablePrecision/Integer/Logarithm.hs @@ -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) diff --git a/pure/Numeric/VariablePrecision/Integer/Logarithm.hs b/pure/Numeric/VariablePrecision/Integer/Logarithm.hs new file mode 100644 --- /dev/null +++ b/pure/Numeric/VariablePrecision/Integer/Logarithm.hs @@ -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) diff --git a/variable-precision.cabal b/variable-precision.cabal --- a/variable-precision.cabal +++ b/variable-precision.cabal @@ -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