exact-pi 0.5.0.1 → 0.5.0.2
raw patch · 9 files changed
+699/−694 lines, 9 filessetup-changedPVP: major bump suggested
API removals or changes: PVP suggests a major version bump
API changes (from Hackage documentation)
- Data.ExactPi.TypeLevel: instance (Numeric.NumType.DK.Integers.KnownTypeInt z, GHC.TypeNats.KnownNat p, GHC.TypeNats.KnownNat q, 1 GHC.TypeNats.<= q) => Data.ExactPi.TypeLevel.KnownExactPi ('Data.ExactPi.TypeLevel.ExactPi' z p q)
+ Data.ExactPi.TypeLevel: instance (Numeric.NumType.DK.Integers.KnownTypeInt z, GHC.TypeNats.KnownNat p, GHC.TypeNats.KnownNat q, 1 Data.Type.Ord.<= q) => Data.ExactPi.TypeLevel.KnownExactPi ('Data.ExactPi.TypeLevel.ExactPi' z p q)
- Data.ExactPi.TypeLevel: injMin :: forall v a. MinCtxt v a => Proxy v -> a
+ Data.ExactPi.TypeLevel: injMin :: forall (v :: ExactPi') a. MinCtxt v a => Proxy v -> a
- Data.ExactPi.TypeLevel: type ExactNatural n = 'ExactPi' 'Zero n 1
+ Data.ExactPi.TypeLevel: type ExactNatural (n :: Nat) = 'ExactPi' 'Zero n 1
- Data.ExactPi.TypeLevel: type MinCtxt v a = (KnownExactPi v, MinCtxt' v a, KnownMinCtxt (MinCtxt' v))
+ Data.ExactPi.TypeLevel: type MinCtxt (v :: ExactPi') a = (KnownExactPi v, MinCtxt' v a, KnownMinCtxt MinCtxt' v)
- Data.ExactPi.TypeLevel: type Pi = 'ExactPi' 'Pos1 1 1
+ Data.ExactPi.TypeLevel: type Pi = 'ExactPi' 'Pos1 1 1
- Data.ExactPi.TypeLevel: type family MinCtxt' (v :: ExactPi')
+ Data.ExactPi.TypeLevel: type family MinCtxt' (v :: ExactPi') :: Type -> Constraint
Files
- LICENSE +20/−20
- README.md +6/−6
- Setup.hs +2/−2
- changelog.md +66/−62
- exact-pi.cabal +54/−54
- src/Data/ExactPi.hs +205/−205
- src/Data/ExactPi/TypeLevel.hs +137/−136
- test-suite/Test.hs +146/−146
- test-suite/TestUtils.hs +63/−63
LICENSE view
@@ -1,20 +1,20 @@-Copyright (c) 2015 Douglas McClean - -Permission is hereby granted, free of charge, to any person obtaining -a copy of this software and associated documentation files (the -"Software"), to deal in the Software without restriction, including -without limitation the rights to use, copy, modify, merge, publish, -distribute, sublicense, and/or sell copies of the Software, and to -permit persons to whom the Software is furnished to do so, subject to -the following conditions: - -The above copyright notice and this permission notice shall be included -in all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +Copyright (c) 2015 Douglas McClean++Permission is hereby granted, free of charge, to any person obtaining+a copy of this software and associated documentation files (the+"Software"), to deal in the Software without restriction, including+without limitation the rights to use, copy, modify, merge, publish,+distribute, sublicense, and/or sell copies of the Software, and to+permit persons to whom the Software is furnished to do so, subject to+the following conditions:++The above copyright notice and this permission notice shall be included+in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
README.md view
@@ -1,6 +1,6 @@-# exact-pi -Exact rational multiples of pi (and integer powers of pi) in Haskell - -[](https://travis-ci.org/dmcclean/exact-pi) -[](http://hackage.haskell.org/package/exact-pi) -[](https://www.stackage.org/package/exact-pi) +# exact-pi+Exact rational multiples of pi (and integer powers of pi) in Haskell++[](https://travis-ci.org/dmcclean/exact-pi)+[](http://hackage.haskell.org/package/exact-pi)+[](https://www.stackage.org/package/exact-pi)
Setup.hs view
@@ -1,2 +1,2 @@-import Distribution.Simple -main = defaultMain +import Distribution.Simple+main = defaultMain
changelog.md view
@@ -1,62 +1,66 @@-0.5.0.1 -------- -* Bump base dependency. -* Resolve compiler warnings. - -0.5.0.0 -------- -* Change implementation of 'rationalApproximations' to use Chudnovsky's approximations. - -0.4.1.4 -------- -* Comply with NoStarIsType pragma. - -0.4.1.3 -------- -* Add Semigroup ExactPi instance. - -0.4.1.2 -------- -* Bump base dependency. - -0.4.1.1 -------- -* Fixed infinite loop in definition of negate. - -0.4.1.0 -------- -* Added function for computing rational approximations of ExactPi values. - -0.4.0.0 -------- -* Added simpler constraints for converting ExactPi types to terms with the minimal context. - -0.3.1.0 -------- -* Added support for exactly comparing values. - -0.3.0.0 -------- -* Added a type-level representation of ExactPi values. - -0.2.1.2 -------- -* Fixed a bug in recip. -* Fixed approximation of exact values with a negative exponent for pi. - -0.2.1.1 -------- -* Fixed a missing case in isZero. - -0.2.1.0 -------- -* Added support for converting to exact integers or exact rationals. - -0.2.0.0 -------- -* Removed dependency on groups package, since it appears not to be widely used. -* Fixed a missing case alternative in recip. - -0.1.2.0 -------- -* Added support for GHC 7.8. +0.5.0.2+-------+* Support GHC 9.4.++0.5.0.1+-------+* Bump base dependency.+* Resolve compiler warnings.++0.5.0.0+-------+* Change implementation of 'rationalApproximations' to use Chudnovsky's approximations.++0.4.1.4+-------+* Comply with NoStarIsType pragma.++0.4.1.3+-------+* Add Semigroup ExactPi instance.++0.4.1.2+-------+* Bump base dependency.++0.4.1.1+-------+* Fixed infinite loop in definition of negate.++0.4.1.0+-------+* Added function for computing rational approximations of ExactPi values.++0.4.0.0+-------+* Added simpler constraints for converting ExactPi types to terms with the minimal context.++0.3.1.0+-------+* Added support for exactly comparing values.++0.3.0.0+-------+* Added a type-level representation of ExactPi values.++0.2.1.2+-------+* Fixed a bug in recip.+* Fixed approximation of exact values with a negative exponent for pi.++0.2.1.1+-------+* Fixed a missing case in isZero.++0.2.1.0+-------+* Added support for converting to exact integers or exact rationals.++0.2.0.0+-------+* Removed dependency on groups package, since it appears not to be widely used.+* Fixed a missing case alternative in recip.++0.1.2.0+-------+* Added support for GHC 7.8.
exact-pi.cabal view
@@ -1,54 +1,54 @@-name: exact-pi -version: 0.5.0.1 -synopsis: Exact rational multiples of pi (and integer powers of pi) -description: Provides an exact representation for rational multiples of pi alongside an approximate representation of all reals. - Useful for storing and computing with conversion factors between physical units. -homepage: https://github.com/dmcclean/exact-pi/ -bug-reports: https://github.com/dmcclean/exact-pi/issues/ -license: MIT -license-file: LICENSE -author: Douglas McClean -maintainer: douglas.mcclean@gmail.com -category: Data -build-type: Simple -extra-source-files: README.md, - changelog.md -cabal-version: >=1.10 -tested-with: GHC == 7.8.4, - GHC == 7.10.3, - GHC == 8.0.2, - GHC == 8.2.2, - GHC == 8.4.3, - GHC == 8.6.1 -library - exposed-modules: Data.ExactPi, - Data.ExactPi.TypeLevel - build-depends: base >=4.7 && <5, - numtype-dk >= 0.5 - if impl(ghc <8.0) - build-depends: - semigroups >=0.8 - ghc-options: -Wall - hs-source-dirs: src - default-language: Haskell2010 - -test-suite spec - main-is: Test.hs - build-depends: base >=4.7 && <4.13, - exact-pi, - numtype-dk >= 0.5, - QuickCheck >=2.10 && <2.12, - tasty >=0.10 && <1.2, - tasty-hunit >=0.9 && <0.11, - tasty-quickcheck >= 0.9 && <0.11 - if impl(ghc < 8.0) - build-depends: semigroups >=0.9 && < 1.0 - other-modules: TestUtils - type: exitcode-stdio-1.0 - ghc-options: -Wall - hs-source-dirs: test-suite - default-language: Haskell2010 - -source-repository head - type: git - location: https://github.com/dmcclean/exact-pi.git +name: exact-pi+version: 0.5.0.2+synopsis: Exact rational multiples of pi (and integer powers of pi)+description: Provides an exact representation for rational multiples of pi alongside an approximate representation of all reals.+ Useful for storing and computing with conversion factors between physical units.+homepage: https://github.com/dmcclean/exact-pi/+bug-reports: https://github.com/dmcclean/exact-pi/issues/+license: MIT+license-file: LICENSE+author: Douglas McClean+maintainer: douglas.mcclean@gmail.com+category: Data+build-type: Simple+extra-source-files: README.md,+ changelog.md+cabal-version: >=1.10+tested-with: GHC == 7.8.4,+ GHC == 7.10.3,+ GHC == 8.0.2,+ GHC == 8.2.2,+ GHC == 8.4.3,+ GHC == 8.6.1+library+ exposed-modules: Data.ExactPi,+ Data.ExactPi.TypeLevel+ build-depends: base >=4.7 && <5,+ numtype-dk >= 0.5+ if impl(ghc <8.0)+ build-depends:+ semigroups >=0.8+ ghc-options: -Wall+ hs-source-dirs: src+ default-language: Haskell2010++test-suite spec+ main-is: Test.hs+ build-depends: base >=4.7 && <5,+ exact-pi,+ numtype-dk >= 0.5,+ QuickCheck >=2.10,+ tasty >=0.10,+ tasty-hunit >=0.9 && <0.11,+ tasty-quickcheck >= 0.9 && <0.11+ if impl(ghc < 8.0)+ build-depends: semigroups >=0.9 && < 1.0+ other-modules: TestUtils+ type: exitcode-stdio-1.0+ ghc-options: -Wall+ hs-source-dirs: test-suite+ default-language: Haskell2010++source-repository head+ type: git+ location: https://github.com/dmcclean/exact-pi.git
src/Data/ExactPi.hs view
@@ -1,205 +1,205 @@-{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE ParallelListComp #-} - -{-# OPTIONS_HADDOCK show-extensions #-} - -{-| -Module : Data.ExactPi -Description : Exact rational multiples of powers of pi -License : MIT -Maintainer : douglas.mcclean@gmail.com -Stability : experimental - -This type is sufficient to exactly express the closure of Q ∪ {π} under multiplication and division. -As a result it is useful for representing conversion factors -between physical units. Approximate values are included both to close the remainder -of the arithmetic operations in the `Num` typeclass and to encode conversion -factors defined experimentally. --} -module Data.ExactPi -( - ExactPi(..), - approximateValue, - isZero, - isExact, - isExactZero, - isExactOne, - areExactlyEqual, - isExactInteger, - toExactInteger, - isExactRational, - toExactRational, - rationalApproximations, - -- * Utils - getRationalLimit -) -where - -import Data.Monoid -import Data.Ratio ((%), numerator, denominator) -import Data.Semigroup -import Prelude - --- | Represents an exact or approximate real value. --- The exactly representable values are rational multiples of an integer power of pi. -data ExactPi = Exact Integer Rational -- ^ @'Exact' z q@ = q * pi^z. Note that this means there are many representations of zero. - | Approximate (forall a.Floating a => a) -- ^ An approximate value. This representation was chosen because it allows conversion to floating types using their native definition of 'pi'. - --- | Approximates an exact or approximate value, converting it to a `Floating` type. --- This uses the value of `pi` supplied by the destination type, to provide the appropriate --- precision. -approximateValue :: Floating a => ExactPi -> a -approximateValue (Exact z q) = (pi ^^ z) * (fromRational q) -approximateValue (Approximate x) = x - --- | Identifies whether an 'ExactPi' is an exact or approximate representation of zero. -isZero :: ExactPi -> Bool -isZero (Exact _ 0) = True -isZero (Approximate x) = x == (0 :: Double) -isZero _ = False - --- | Identifies whether an 'ExactPi' is an exact value. -isExact :: ExactPi -> Bool -isExact (Exact _ _) = True -isExact _ = False - --- | Identifies whether an 'ExactPi' is an exact representation of zero. -isExactZero :: ExactPi -> Bool -isExactZero (Exact _ 0) = True -isExactZero _ = False - --- | Identifies whether an 'ExactPi' is an exact representation of one. -isExactOne :: ExactPi -> Bool -isExactOne (Exact 0 1) = True -isExactOne _ = False - --- | Identifies whether two 'ExactPi' values are exactly equal. -areExactlyEqual :: ExactPi -> ExactPi -> Bool -areExactlyEqual (Exact z1 q1) (Exact z2 q2) = (z1 == z2 && q1 == q2) || (q1 == 0 && q2 == 0) -areExactlyEqual _ _ = False - --- | Identifies whether an 'ExactPi' is an exact representation of an integer. -isExactInteger :: ExactPi -> Bool -isExactInteger (Exact 0 q) | denominator q == 1 = True -isExactInteger _ = False - --- | Converts an 'ExactPi' to an exact 'Integer' or 'Nothing'. -toExactInteger :: ExactPi -> Maybe Integer -toExactInteger (Exact 0 q) | denominator q == 1 = Just $ numerator q -toExactInteger _ = Nothing - --- | Identifies whether an 'ExactPi' is an exact representation of a rational. -isExactRational :: ExactPi -> Bool -isExactRational (Exact 0 _) = True -isExactRational _ = False - --- | Converts an 'ExactPi' to an exact 'Rational' or 'Nothing'. -toExactRational :: ExactPi -> Maybe Rational -toExactRational (Exact 0 q) = Just q -toExactRational _ = Nothing - --- | Converts an 'ExactPi' to a list of increasingly accurate rational approximations. Note --- that 'Approximate' values are converted using the 'Real' instance for 'Double' into a --- singleton list. Note that exact rationals are also converted into a singleton list. --- --- Implementation is based on Chudnovsky's algorithm. -rationalApproximations :: ExactPi -> [Rational] -rationalApproximations (Approximate x) = [toRational (x :: Double)] -rationalApproximations (Exact _ 0) = [0] -rationalApproximations (Exact 0 q) = [q] -rationalApproximations (Exact z q) - | even z = [q * 10005^^k * c^^z | c <- chudnovsky] - | otherwise = [q * 10005^^k * c^^z * r | c <- chudnovsky | r <- rootApproximation] - where k = z `div` 2 - -chudnovsky :: [Rational] -chudnovsky = [426880 / s | s <- partials] - where lk = iterate (+545140134) 13591409 - xk = iterate (*(-262537412640768000)) 1 - kk = iterate (+12) 6 - mk = 1: [m * ((k^(3::Int) - 16*k) % (n+1)^(3::Int)) | m <- mk | k <- kk | n <- [0..]] - values = [m * l / x | m <- mk | l <- lk | x <- xk] - partials = scanl1 (+) values - --- | Given an infinite converging sequence of rationals, find their limit. --- Takes a comparison function to determine when convergence is close enough. --- --- >>> getRationalLimit (==) (rationalApproximations (Exact 1 1)) :: Double --- 3.141592653589793 -getRationalLimit :: Fractional a => (a -> a -> Bool) -> [Rational] -> a -getRationalLimit cmp = go . map fromRational - where go (x:y:xs) - | cmp x y = y - | otherwise = go (y:xs) - go [x] = x - go _ = error "did not converge" - --- | A sequence of convergents approximating @sqrt 10005@, intended to be zipped --- with 'chudnovsky' in 'rationalApproximations'. Carefully chosen so that --- the denominator does not increase too rapidly but approximations are still --- appropriately precise. --- --- Chudnovsky's series provides no more than 15 digits --- per iteration, so the root approximation should not --- have a more rapid rate of convergence. -rootApproximation :: [Rational] -rootApproximation = map head . iterate (drop 4) $ go 1 0 100 1 40 - where - go pk' qk' pk qk a = (pk % qk): go pk qk (pk' + a*pk) (qk' + a*qk) (240-a) - -instance Show ExactPi where - show (Exact z q) | z == 0 = "Exactly " ++ show q - | z == 1 = "Exactly pi * " ++ show q - | otherwise = "Exactly pi^" ++ show z ++ " * " ++ show q - show (Approximate x) = "Approximately " ++ show (x :: Double) - -instance Num ExactPi where - fromInteger n = Exact 0 (fromInteger n) - (Exact z1 q1) * (Exact z2 q2) = Exact (z1 + z2) (q1 * q2) - (Exact _ 0) * _ = 0 - _ * (Exact _ 0) = 0 - x * y = Approximate $ approximateValue x * approximateValue y - (Exact z1 q1) + (Exact z2 q2) | z1 == z2 = Exact z1 (q1 + q2) -- by distributive property - x + y = Approximate $ approximateValue x + approximateValue y - abs (Exact z q) = Exact z (abs q) - abs (Approximate x) = Approximate $ abs x - signum (Exact _ q) = Exact 0 (signum q) - signum (Approximate x) = Approximate $ signum x -- we leave this tagged as approximate because we don't know "how" approximate the input was. a case could be made for exact answers here. - negate x = (Exact 0 (-1)) * x - -instance Fractional ExactPi where - fromRational = Exact 0 - recip (Exact z q) = Exact (negate z) (recip q) - recip (Approximate x) = Approximate (recip x) - -instance Floating ExactPi where - pi = Exact 1 1 - exp x | isExactZero x = 1 - | otherwise = approx1 exp x - log (Exact 0 1) = 0 - log x = approx1 log x - -- It would be possible to give tighter bounds to the trig functions, preserving exactness for arguments that have an exactly representable result. - sin = approx1 sin - cos = approx1 cos - tan = approx1 tan - asin = approx1 asin - atan = approx1 atan - acos = approx1 acos - sinh = approx1 sinh - cosh = approx1 cosh - tanh = approx1 tanh - asinh = approx1 asinh - acosh = approx1 acosh - atanh = approx1 atanh - -approx1 :: (forall a.Floating a => a -> a) -> ExactPi -> ExactPi -approx1 f x = Approximate (f (approximateValue x)) - --- | The multiplicative semigroup over 'Rational's augmented with multiples of 'pi'. -instance Semigroup ExactPi where - (<>) = mappend - --- | The multiplicative monoid over 'Rational's augmented with multiples of 'pi'. -instance Monoid ExactPi where - mempty = 1 - mappend = (*) +{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ParallelListComp #-}++{-# OPTIONS_HADDOCK show-extensions #-}++{-|+Module : Data.ExactPi+Description : Exact rational multiples of powers of pi+License : MIT+Maintainer : douglas.mcclean@gmail.com+Stability : experimental++This type is sufficient to exactly express the closure of Q ∪ {π} under multiplication and division.+As a result it is useful for representing conversion factors+between physical units. Approximate values are included both to close the remainder+of the arithmetic operations in the `Num` typeclass and to encode conversion+factors defined experimentally.+-}+module Data.ExactPi+(+ ExactPi(..),+ approximateValue,+ isZero,+ isExact,+ isExactZero,+ isExactOne,+ areExactlyEqual,+ isExactInteger,+ toExactInteger,+ isExactRational,+ toExactRational,+ rationalApproximations,+ -- * Utils+ getRationalLimit+)+where++import Data.Monoid+import Data.Ratio ((%), numerator, denominator)+import Data.Semigroup+import Prelude++-- | Represents an exact or approximate real value.+-- The exactly representable values are rational multiples of an integer power of pi.+data ExactPi = Exact Integer Rational -- ^ @'Exact' z q@ = q * pi^z. Note that this means there are many representations of zero.+ | Approximate (forall a.Floating a => a) -- ^ An approximate value. This representation was chosen because it allows conversion to floating types using their native definition of 'pi'.++-- | Approximates an exact or approximate value, converting it to a `Floating` type.+-- This uses the value of `pi` supplied by the destination type, to provide the appropriate+-- precision.+approximateValue :: Floating a => ExactPi -> a+approximateValue (Exact z q) = (pi ^^ z) * (fromRational q)+approximateValue (Approximate x) = x++-- | Identifies whether an 'ExactPi' is an exact or approximate representation of zero.+isZero :: ExactPi -> Bool+isZero (Exact _ 0) = True+isZero (Approximate x) = x == (0 :: Double)+isZero _ = False++-- | Identifies whether an 'ExactPi' is an exact value.+isExact :: ExactPi -> Bool+isExact (Exact _ _) = True+isExact _ = False++-- | Identifies whether an 'ExactPi' is an exact representation of zero.+isExactZero :: ExactPi -> Bool+isExactZero (Exact _ 0) = True+isExactZero _ = False++-- | Identifies whether an 'ExactPi' is an exact representation of one.+isExactOne :: ExactPi -> Bool+isExactOne (Exact 0 1) = True+isExactOne _ = False++-- | Identifies whether two 'ExactPi' values are exactly equal.+areExactlyEqual :: ExactPi -> ExactPi -> Bool+areExactlyEqual (Exact z1 q1) (Exact z2 q2) = (z1 == z2 && q1 == q2) || (q1 == 0 && q2 == 0)+areExactlyEqual _ _ = False++-- | Identifies whether an 'ExactPi' is an exact representation of an integer.+isExactInteger :: ExactPi -> Bool+isExactInteger (Exact 0 q) | denominator q == 1 = True+isExactInteger _ = False++-- | Converts an 'ExactPi' to an exact 'Integer' or 'Nothing'.+toExactInteger :: ExactPi -> Maybe Integer+toExactInteger (Exact 0 q) | denominator q == 1 = Just $ numerator q+toExactInteger _ = Nothing++-- | Identifies whether an 'ExactPi' is an exact representation of a rational.+isExactRational :: ExactPi -> Bool+isExactRational (Exact 0 _) = True+isExactRational _ = False++-- | Converts an 'ExactPi' to an exact 'Rational' or 'Nothing'.+toExactRational :: ExactPi -> Maybe Rational+toExactRational (Exact 0 q) = Just q+toExactRational _ = Nothing++-- | Converts an 'ExactPi' to a list of increasingly accurate rational approximations. Note+-- that 'Approximate' values are converted using the 'Real' instance for 'Double' into a+-- singleton list. Note that exact rationals are also converted into a singleton list.+--+-- Implementation is based on Chudnovsky's algorithm.+rationalApproximations :: ExactPi -> [Rational]+rationalApproximations (Approximate x) = [toRational (x :: Double)]+rationalApproximations (Exact _ 0) = [0]+rationalApproximations (Exact 0 q) = [q]+rationalApproximations (Exact z q)+ | even z = [q * 10005^^k * c^^z | c <- chudnovsky]+ | otherwise = [q * 10005^^k * c^^z * r | c <- chudnovsky | r <- rootApproximation]+ where k = z `div` 2++chudnovsky :: [Rational]+chudnovsky = [426880 / s | s <- partials]+ where lk = iterate (+545140134) 13591409+ xk = iterate (*(-262537412640768000)) 1+ kk = iterate (+12) 6+ mk = 1: [m * ((k^(3::Int) - 16*k) % (n+1)^(3::Int)) | m <- mk | k <- kk | n <- [0..]]+ values = [m * l / x | m <- mk | l <- lk | x <- xk]+ partials = scanl1 (+) values++-- | Given an infinite converging sequence of rationals, find their limit.+-- Takes a comparison function to determine when convergence is close enough.+--+-- >>> getRationalLimit (==) (rationalApproximations (Exact 1 1)) :: Double+-- 3.141592653589793+getRationalLimit :: Fractional a => (a -> a -> Bool) -> [Rational] -> a+getRationalLimit cmp = go . map fromRational+ where go (x:y:xs)+ | cmp x y = y+ | otherwise = go (y:xs)+ go [x] = x+ go _ = error "did not converge"++-- | A sequence of convergents approximating @sqrt 10005@, intended to be zipped+-- with 'chudnovsky' in 'rationalApproximations'. Carefully chosen so that+-- the denominator does not increase too rapidly but approximations are still+-- appropriately precise.+--+-- Chudnovsky's series provides no more than 15 digits+-- per iteration, so the root approximation should not+-- have a more rapid rate of convergence.+rootApproximation :: [Rational]+rootApproximation = map head . iterate (drop 4) $ go 1 0 100 1 40+ where+ go pk' qk' pk qk a = (pk % qk): go pk qk (pk' + a*pk) (qk' + a*qk) (240-a)++instance Show ExactPi where+ show (Exact z q) | z == 0 = "Exactly " ++ show q+ | z == 1 = "Exactly pi * " ++ show q+ | otherwise = "Exactly pi^" ++ show z ++ " * " ++ show q+ show (Approximate x) = "Approximately " ++ show (x :: Double)++instance Num ExactPi where+ fromInteger n = Exact 0 (fromInteger n)+ (Exact z1 q1) * (Exact z2 q2) = Exact (z1 + z2) (q1 * q2)+ (Exact _ 0) * _ = 0+ _ * (Exact _ 0) = 0+ x * y = Approximate $ approximateValue x * approximateValue y+ (Exact z1 q1) + (Exact z2 q2) | z1 == z2 = Exact z1 (q1 + q2) -- by distributive property+ x + y = Approximate $ approximateValue x + approximateValue y+ abs (Exact z q) = Exact z (abs q)+ abs (Approximate x) = Approximate $ abs x+ signum (Exact _ q) = Exact 0 (signum q)+ signum (Approximate x) = Approximate $ signum x -- we leave this tagged as approximate because we don't know "how" approximate the input was. a case could be made for exact answers here.+ negate x = (Exact 0 (-1)) * x++instance Fractional ExactPi where+ fromRational = Exact 0+ recip (Exact z q) = Exact (negate z) (recip q)+ recip (Approximate x) = Approximate (recip x)++instance Floating ExactPi where+ pi = Exact 1 1+ exp x | isExactZero x = 1+ | otherwise = approx1 exp x+ log (Exact 0 1) = 0+ log x = approx1 log x+ -- It would be possible to give tighter bounds to the trig functions, preserving exactness for arguments that have an exactly representable result.+ sin = approx1 sin+ cos = approx1 cos+ tan = approx1 tan+ asin = approx1 asin+ atan = approx1 atan+ acos = approx1 acos+ sinh = approx1 sinh+ cosh = approx1 cosh+ tanh = approx1 tanh+ asinh = approx1 asinh+ acosh = approx1 acosh+ atanh = approx1 atanh++approx1 :: (forall a.Floating a => a -> a) -> ExactPi -> ExactPi+approx1 f x = Approximate (f (approximateValue x))++-- | The multiplicative semigroup over 'Rational's augmented with multiples of 'pi'.+instance Semigroup ExactPi where+ (<>) = mappend++-- | The multiplicative monoid over 'Rational's augmented with multiples of 'pi'.+instance Monoid ExactPi where+ mempty = 1+ mappend = (*)
src/Data/ExactPi/TypeLevel.hs view
@@ -1,136 +1,137 @@-{-# OPTIONS_HADDOCK show-extensions #-} - -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE CPP #-} -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE TypeOperators #-} - -#if __GLASGOW_HASKELL__ > 805 -{-# LANGUAGE NoStarIsType #-} -#endif - -{-| -Module : Data.ExactPi.TypeLevel -Description : Exact non-negative rational multiples of powers of pi at the type level -License : MIT -Maintainer : douglas.mcclean@gmail.com -Stability : experimental - -This kind is sufficient to exactly express the closure of Q⁺ ∪ {π} under multiplication and division. -As a result it is useful for representing conversion factors between physical units. --} -module Data.ExactPi.TypeLevel -( - -- * Type Level ExactPi Values - type ExactPi'(..), - KnownExactPi(..), - -- * Arithmetic - type (*), type (/), type Recip, - type ExactNatural, - type One, type Pi, - -- * Conversion to Term Level - type MinCtxt, type MinCtxt', - injMin -) -where - -import Data.ExactPi -import Data.Maybe (fromJust) -import Data.Proxy -import Data.Ratio -import GHC.TypeLits hiding (type (*), type (^)) -import qualified GHC.TypeLits as N -import Numeric.NumType.DK.Integers hiding (type (*), type (/)) -import qualified Numeric.NumType.DK.Integers as Z - --- | A type-level representation of a non-negative rational multiple of an integer power of pi. --- --- Each type in this kind can be exactly represented at the term level by a value of type 'ExactPi', --- provided that its denominator is non-zero. --- --- Note that there are many representations of zero, and many representations of dividing by zero. --- These are not excluded because doing so introduces a lot of extra machinery. Play nice! Future --- versions may not include a representation for zero. --- --- Of course there are also many representations of every value, because the numerator need not be --- comprime to the denominator. For many purposes it is not necessary to maintain the types in reduced --- form, they will be appropriately reduced when converted to terms. -data ExactPi' = ExactPi' TypeInt -- Exponent of pi - Nat -- Numerator - Nat -- Denominator - --- | A KnownDimension is one for which we can construct a term-level representation. --- --- Each validly constructed type of kind 'ExactPi'' has a 'KnownExactPi' instance, provided that --- its denominator is non-zero. -class KnownExactPi (v :: ExactPi') where - -- | Converts an 'ExactPi'' type to an 'ExactPi' value. - exactPiVal :: Proxy v -> ExactPi - --- | Determines the minimum context required for a numeric type to hold the value --- associated with a specific 'ExactPi'' type. -type family MinCtxt' (v :: ExactPi') where - MinCtxt' ('ExactPi' 'Zero p 1) = Num - MinCtxt' ('ExactPi' 'Zero p q) = Fractional - MinCtxt' ('ExactPi' z p q) = Floating - -type MinCtxt v a = (KnownExactPi v, MinCtxt' v a, KnownMinCtxt (MinCtxt' v)) - --- | A KnownMinCtxt is a contraint on values sufficient to allow us to inject certain --- 'ExactPi' values into types that satisfy the constraint. -class KnownMinCtxt c where - -- | Injects an 'ExactPi' value into a specified type satisfying this constraint. - -- - -- The injection is permitted to fail if type constraint does not entail the 'MinCtxt' - -- required by the 'ExactPi'' representation of the supplied 'ExactPi' value. - inj :: c a => Proxy c -- ^ A proxy for identifying the required constraint. - -> ExactPi -- ^ The value to inject. - -> a -- ^ A value of the constrained type corresponding to the supplied 'ExactPi' value. - -instance KnownMinCtxt Num where - inj _ = fromInteger . fromJust . toExactInteger - -instance KnownMinCtxt Fractional where - inj _ = fromRational . fromJust . toExactRational - -instance KnownMinCtxt Floating where - inj _ = approximateValue - --- | Converts an 'ExactPi'' type to a numeric value with the minimum required context. --- --- When the value is known to be an integer, it can be returned as any instance of 'Num'. Similarly, --- rationals require 'Fractional', and values that involve 'pi' require 'Floating'. -injMin :: forall v a.(MinCtxt v a) => Proxy v -> a -injMin = inj (Proxy :: Proxy (MinCtxt' v)) . exactPiVal - -instance (KnownTypeInt z, KnownNat p, KnownNat q, 1 <= q) => KnownExactPi ('ExactPi' z p q) where - exactPiVal _ = Exact z' (p' % q') - where - z' = toNum (Proxy :: Proxy z) - p' = natVal (Proxy :: Proxy p) - q' = natVal (Proxy :: Proxy q) - --- | Forms the product of 'ExactPi'' types (in the arithmetic sense). -type family (a :: ExactPi') * (b :: ExactPi') :: ExactPi' where - ('ExactPi' z p q) * ('ExactPi' z' p' q') = 'ExactPi' (z Z.+ z') (p N.* p') (q N.* q') - --- | Forms the quotient of 'ExactPi'' types (in the arithmetic sense). -type family (a :: ExactPi') / (b :: ExactPi') :: ExactPi' where - ('ExactPi' z p q) / ('ExactPi' z' p' q') = 'ExactPi' (z Z.- z') (p N.* q') (q N.* p') - --- | Forms the reciprocal of an 'ExactPi'' type. -type family Recip (a :: ExactPi') :: ExactPi' where - Recip ('ExactPi' z p q) = 'ExactPi' (Negate z) q p - --- | Converts a type-level natural to an 'ExactPi'' type. -type ExactNatural n = 'ExactPi' 'Zero n 1 - --- | The 'ExactPi'' type representing the number 1. -type One = ExactNatural 1 - --- | The 'ExactPi'' type representing the number 'pi'. -type Pi = 'ExactPi' 'Pos1 1 1 +{-# OPTIONS_HADDOCK show-extensions #-}++{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE UndecidableInstances #-}++#if __GLASGOW_HASKELL__ > 805+{-# LANGUAGE NoStarIsType #-}+#endif++{-|+Module : Data.ExactPi.TypeLevel+Description : Exact non-negative rational multiples of powers of pi at the type level+License : MIT+Maintainer : douglas.mcclean@gmail.com+Stability : experimental++This kind is sufficient to exactly express the closure of Q⁺ ∪ {π} under multiplication and division.+As a result it is useful for representing conversion factors between physical units.+-}+module Data.ExactPi.TypeLevel+(+ -- * Type Level ExactPi Values+ type ExactPi'(..),+ KnownExactPi(..),+ -- * Arithmetic+ type (*), type (/), type Recip,+ type ExactNatural,+ type One, type Pi,+ -- * Conversion to Term Level+ type MinCtxt, type MinCtxt',+ injMin+)+where++import Data.ExactPi+import Data.Maybe (fromJust)+import Data.Proxy+import Data.Ratio+import GHC.TypeLits hiding (type (*), type (^))+import qualified GHC.TypeLits as N+import Numeric.NumType.DK.Integers hiding (type (*), type (/))+import qualified Numeric.NumType.DK.Integers as Z++-- | A type-level representation of a non-negative rational multiple of an integer power of pi.+--+-- Each type in this kind can be exactly represented at the term level by a value of type 'ExactPi',+-- provided that its denominator is non-zero.+--+-- Note that there are many representations of zero, and many representations of dividing by zero.+-- These are not excluded because doing so introduces a lot of extra machinery. Play nice! Future+-- versions may not include a representation for zero.+--+-- Of course there are also many representations of every value, because the numerator need not be+-- comprime to the denominator. For many purposes it is not necessary to maintain the types in reduced+-- form, they will be appropriately reduced when converted to terms.+data ExactPi' = ExactPi' TypeInt -- Exponent of pi+ Nat -- Numerator+ Nat -- Denominator++-- | A KnownDimension is one for which we can construct a term-level representation.+--+-- Each validly constructed type of kind 'ExactPi'' has a 'KnownExactPi' instance, provided that+-- its denominator is non-zero.+class KnownExactPi (v :: ExactPi') where+ -- | Converts an 'ExactPi'' type to an 'ExactPi' value.+ exactPiVal :: Proxy v -> ExactPi++-- | Determines the minimum context required for a numeric type to hold the value+-- associated with a specific 'ExactPi'' type.+type family MinCtxt' (v :: ExactPi') where+ MinCtxt' ('ExactPi' 'Zero p 1) = Num+ MinCtxt' ('ExactPi' 'Zero p q) = Fractional+ MinCtxt' ('ExactPi' z p q) = Floating++type MinCtxt v a = (KnownExactPi v, MinCtxt' v a, KnownMinCtxt (MinCtxt' v))++-- | A KnownMinCtxt is a contraint on values sufficient to allow us to inject certain+-- 'ExactPi' values into types that satisfy the constraint.+class KnownMinCtxt c where+ -- | Injects an 'ExactPi' value into a specified type satisfying this constraint.+ --+ -- The injection is permitted to fail if type constraint does not entail the 'MinCtxt'+ -- required by the 'ExactPi'' representation of the supplied 'ExactPi' value.+ inj :: c a => Proxy c -- ^ A proxy for identifying the required constraint.+ -> ExactPi -- ^ The value to inject.+ -> a -- ^ A value of the constrained type corresponding to the supplied 'ExactPi' value.++instance KnownMinCtxt Num where+ inj _ = fromInteger . fromJust . toExactInteger++instance KnownMinCtxt Fractional where+ inj _ = fromRational . fromJust . toExactRational++instance KnownMinCtxt Floating where+ inj _ = approximateValue++-- | Converts an 'ExactPi'' type to a numeric value with the minimum required context.+--+-- When the value is known to be an integer, it can be returned as any instance of 'Num'. Similarly,+-- rationals require 'Fractional', and values that involve 'pi' require 'Floating'.+injMin :: forall v a.(MinCtxt v a) => Proxy v -> a+injMin = inj (Proxy :: Proxy (MinCtxt' v)) . exactPiVal++instance (KnownTypeInt z, KnownNat p, KnownNat q, 1 <= q) => KnownExactPi ('ExactPi' z p q) where+ exactPiVal _ = Exact z' (p' % q')+ where+ z' = toNum (Proxy :: Proxy z)+ p' = natVal (Proxy :: Proxy p)+ q' = natVal (Proxy :: Proxy q)++-- | Forms the product of 'ExactPi'' types (in the arithmetic sense).+type family (a :: ExactPi') * (b :: ExactPi') :: ExactPi' where+ ('ExactPi' z p q) * ('ExactPi' z' p' q') = 'ExactPi' (z Z.+ z') (p N.* p') (q N.* q')++-- | Forms the quotient of 'ExactPi'' types (in the arithmetic sense).+type family (a :: ExactPi') / (b :: ExactPi') :: ExactPi' where+ ('ExactPi' z p q) / ('ExactPi' z' p' q') = 'ExactPi' (z Z.- z') (p N.* q') (q N.* p')++-- | Forms the reciprocal of an 'ExactPi'' type.+type family Recip (a :: ExactPi') :: ExactPi' where+ Recip ('ExactPi' z p q) = 'ExactPi' (Negate z) q p++-- | Converts a type-level natural to an 'ExactPi'' type.+type ExactNatural n = 'ExactPi' 'Zero n 1++-- | The 'ExactPi'' type representing the number 1.+type One = ExactNatural 1++-- | The 'ExactPi'' type representing the number 'pi'.+type Pi = 'ExactPi' 'Pos1 1 1
test-suite/Test.hs view
@@ -1,146 +1,146 @@-{-# LANGUAGE DataKinds #-} -{-# OPTIONS_GHC -fno-warn-type-defaults #-} -import Data.Fixed (Fixed(..)) -import Data.Ratio ((%)) -import Test.Tasty (TestTree, testGroup, defaultMain) -import Test.Tasty.HUnit ((@?=), Assertion, testCase) -import Test.Tasty.QuickCheck (testProperty) -import Test.QuickCheck (Positive(..)) - -import Data.ExactPi -import TestUtils (E, getValue, getDigit, getDigitBBP) - --- test pi^2 first since it does not rely on square roots -piSquaredDouble :: Assertion -piSquaredDouble = getValue (Exact 2 1) @?= (pi^2 :: Double) - --- first 57 digits of pi^2 --- http://www.wolframalpha.com/input/?i=pi%5E2 -piSquaredWAstart :: Assertion -piSquaredWAstart = getValue (Exact 2 1) @?= piSquared - -piSquared :: Fixed (E 57) -piSquared = 9.869604401089358618834490999876151135313699407240790626413 - --- last 21 digits of pi^2 on wolfram alpha http://www.wolframalpha.com/input/?i=pi%5E2 --- by asking for more digits as much as possible -piSquaredWAend :: Assertion -piSquaredWAend = x `mod` (10^21) @?= 643271910414561208753 - where - MkFixed x = getValue (Exact 2 1) :: Fixed (E 3647) - --- test first term matches formula of chudnovsky's algorithm -firstApproximation :: Assertion -firstApproximation = head (rationalApproximations (Exact 2 1)) @?= (426880 % 13591409)^2 * 10005 - --- pi tests -piDouble :: Assertion -piDouble = getValue (Exact 1 1) @?= (pi :: Double) - -piMatchesOeis :: Assertion -piMatchesOeis = getValue (Exact 1 1) @?= oeisValue - --- https://oeis.org/A000796 -oeisValue :: Fixed (E 104) -oeisValue = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214 - --- digits 762 to 767 of pi are 999999 -feynmanPoint :: Assertion -feynmanPoint = x `mod` 1000000 @?= 999999 - where - MkFixed x = getValue (Exact 1 1) :: Fixed (E 767) - --- last 21 digits of pi on wolfram alpha (http://www.wolframalpha.com/input/?i=pi) --- by asking for more digits as much as possible -piWAend :: Assertion -piWAend = x `mod` (10^21) @?= 706420467525907091548 - where - MkFixed x = getValue (Exact 1 1) :: Fixed (E 3647) - --- pi power tests --- http://www.wolframalpha.com/input/?i=1000th+digit+of+pi%5E3%2F10 -pi3 :: Assertion -pi3 = x `mod` 100 @?= 98 - where - MkFixed x = getValue (Exact 3 (1 % 10)) :: Fixed (E 1000) - --- http://www.wolframalpha.com/input/?i=1000th+digit+of+pi%5E-1+*+10 -piNegOne :: Assertion -piNegOne = x `mod` 100 @?= 87 - where - MkFixed x = getValue (Exact (-1) 10) :: Fixed (E 1000) - --- http://www.wolframalpha.com/input/?i=1000th+digit+of+pi%5E10+%2F+10%5E4 -pi10 :: Assertion -pi10 = x `mod` 100 @?= 58 - where - MkFixed x = getValue (Exact 10 (1 % 10^4)) :: Fixed (E 1000) - --- http://www.wolframalpha.com/input/?i=1000th+digit+of+pi%5E-10+*+100000 -piNeg10 :: Assertion -piNeg10 = x `mod` 100 @?= 01 - where - MkFixed x = getValue (Exact (-10) (10^5)) :: Fixed (E 1000) - --- http://www.wolframalpha.com/input/?i=400th+digit+of+pi%5E51+*+10%5E-25 -pi51 :: Assertion -pi51 = x `mod` 100 @?= 39 - where - MkFixed x = getValue (Exact 51 (1 % 10^25)) :: Fixed (E 400) - --- http://www.wolframalpha.com/input/?i=400th+digit+of+pi%5E-51+*+10%5E26 -piNeg51 :: Assertion -piNeg51 = x `mod` 100 @?= 93 - where - MkFixed x = getValue (Exact (-51) (10^26)) :: Fixed (E 400) - --- exact value of riemann zeta(50): should be very near 1 -zeta50 :: ExactPi -zeta50 = Exact 50 (39604576419286371856998202 % 285258771457546764463363635252374414183254365234375) - -zeta200 :: ExactPi -zeta200 = Exact 200 (996768098856666829529857264280799324216991774914413349936111645234527339243047375137731604604421998265202825395226558782117309054290681031680198580956052700765605768743424718675968548245722319600560038220395777111787342302 % 2682678748792657844957504192313280657551803049278355275671666881580642758576467817615493645217977237214155689404787155170845497733836863647685885197919191727452679238952541411298115541287013688972773507748859386210346035176197388875022427877722880764252312145081723341902733317236524547144682628641021437942981719970703125) - --- value of zeta(50) - 1 from wolfram alpha (up to a Double) --- http://www.wolframalpha.com/input/?i=zeta(50)-1 -zeta50MinusOne :: Assertion -zeta50MinusOne = t @?= 8.8817842109308159e-16 - where - t = getRationalLimit (==) . map (subtract 1) . rationalApproximations $ zeta50 :: Double - --- http://www.wolframalpha.com/input/?i=zeta(200)-1 -zeta200MinusOne :: Assertion -zeta200MinusOne = t @?= 6.2230152778611417071e-61 - where - t = getRationalLimit (==) . map (subtract 1) . rationalApproximations $ zeta200 :: Double - --- test against bbp formula -prop :: Positive Integer -> Bool -prop (Positive n) = getDigit n == getDigitBBP (n - 1) - -tests :: TestTree -tests = testGroup "Rational approximation tests" - [ testGroup "π² tests" [ testCase "matches double precision" piSquaredDouble - , testCase "matches start of wolfram alpha" piSquaredWAstart - , testCase "matches end of wolfram alpha" piSquaredWAend - , testCase "first term matches chudnovsky" firstApproximation - ] - , testGroup "π tests" [ testCase "matches double precision" piDouble - , testCase "matches oeis digits" piMatchesOeis - , testCase "has feynman point" feynmanPoint - , testCase "matches end of wolfram alpha" piWAend - ] - , testGroup "πᵏ tests" [ testCase "digits near 1000 of k=3" pi3 - , testCase "digits near 1000 of k=-1" piNegOne - , testCase "digits near 1000 of k=10" pi10 - , testCase "digits near 1000 of k=-10" piNeg10 - , testCase "digits near 400 of k=51" pi51 - , testCase "digits near 400 of k=-51" piNeg51 - , testCase "ζ(50)-1 double precision" zeta50MinusOne - , testCase "ζ(500)-1 double precision" zeta200MinusOne - ] - , testProperty "hex digits match BBP formula" prop - ] - -main :: IO () -main = defaultMain tests +{-# LANGUAGE DataKinds #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}+import Data.Fixed (Fixed(..))+import Data.Ratio ((%))+import Test.Tasty (TestTree, testGroup, defaultMain)+import Test.Tasty.HUnit ((@?=), Assertion, testCase)+import Test.Tasty.QuickCheck (testProperty)+import Test.QuickCheck (Positive(..))++import Data.ExactPi+import TestUtils (E, getValue, getDigit, getDigitBBP)++-- test pi^2 first since it does not rely on square roots+piSquaredDouble :: Assertion+piSquaredDouble = getValue (Exact 2 1) @?= (pi^2 :: Double)++-- first 57 digits of pi^2+-- http://www.wolframalpha.com/input/?i=pi%5E2+piSquaredWAstart :: Assertion+piSquaredWAstart = getValue (Exact 2 1) @?= piSquared++piSquared :: Fixed (E 57)+piSquared = 9.869604401089358618834490999876151135313699407240790626413++-- last 21 digits of pi^2 on wolfram alpha http://www.wolframalpha.com/input/?i=pi%5E2+-- by asking for more digits as much as possible+piSquaredWAend :: Assertion+piSquaredWAend = x `mod` (10^21) @?= 643271910414561208753+ where+ MkFixed x = getValue (Exact 2 1) :: Fixed (E 3647)++-- test first term matches formula of chudnovsky's algorithm+firstApproximation :: Assertion+firstApproximation = head (rationalApproximations (Exact 2 1)) @?= (426880 % 13591409)^2 * 10005++-- pi tests+piDouble :: Assertion+piDouble = getValue (Exact 1 1) @?= (pi :: Double)++piMatchesOeis :: Assertion+piMatchesOeis = getValue (Exact 1 1) @?= oeisValue++-- https://oeis.org/A000796+oeisValue :: Fixed (E 104)+oeisValue = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214++-- digits 762 to 767 of pi are 999999+feynmanPoint :: Assertion+feynmanPoint = x `mod` 1000000 @?= 999999+ where+ MkFixed x = getValue (Exact 1 1) :: Fixed (E 767)++-- last 21 digits of pi on wolfram alpha (http://www.wolframalpha.com/input/?i=pi)+-- by asking for more digits as much as possible+piWAend :: Assertion+piWAend = x `mod` (10^21) @?= 706420467525907091548+ where+ MkFixed x = getValue (Exact 1 1) :: Fixed (E 3647)++-- pi power tests+-- http://www.wolframalpha.com/input/?i=1000th+digit+of+pi%5E3%2F10+pi3 :: Assertion+pi3 = x `mod` 100 @?= 98+ where+ MkFixed x = getValue (Exact 3 (1 % 10)) :: Fixed (E 1000)++-- http://www.wolframalpha.com/input/?i=1000th+digit+of+pi%5E-1+*+10+piNegOne :: Assertion+piNegOne = x `mod` 100 @?= 87+ where+ MkFixed x = getValue (Exact (-1) 10) :: Fixed (E 1000)++-- http://www.wolframalpha.com/input/?i=1000th+digit+of+pi%5E10+%2F+10%5E4+pi10 :: Assertion+pi10 = x `mod` 100 @?= 58+ where+ MkFixed x = getValue (Exact 10 (1 % 10^4)) :: Fixed (E 1000)++-- http://www.wolframalpha.com/input/?i=1000th+digit+of+pi%5E-10+*+100000+piNeg10 :: Assertion+piNeg10 = x `mod` 100 @?= 01+ where+ MkFixed x = getValue (Exact (-10) (10^5)) :: Fixed (E 1000)++-- http://www.wolframalpha.com/input/?i=400th+digit+of+pi%5E51+*+10%5E-25+pi51 :: Assertion+pi51 = x `mod` 100 @?= 39+ where+ MkFixed x = getValue (Exact 51 (1 % 10^25)) :: Fixed (E 400)++-- http://www.wolframalpha.com/input/?i=400th+digit+of+pi%5E-51+*+10%5E26+piNeg51 :: Assertion+piNeg51 = x `mod` 100 @?= 93+ where+ MkFixed x = getValue (Exact (-51) (10^26)) :: Fixed (E 400)++-- exact value of riemann zeta(50): should be very near 1+zeta50 :: ExactPi+zeta50 = Exact 50 (39604576419286371856998202 % 285258771457546764463363635252374414183254365234375)++zeta200 :: ExactPi+zeta200 = Exact 200 (996768098856666829529857264280799324216991774914413349936111645234527339243047375137731604604421998265202825395226558782117309054290681031680198580956052700765605768743424718675968548245722319600560038220395777111787342302 % 2682678748792657844957504192313280657551803049278355275671666881580642758576467817615493645217977237214155689404787155170845497733836863647685885197919191727452679238952541411298115541287013688972773507748859386210346035176197388875022427877722880764252312145081723341902733317236524547144682628641021437942981719970703125)++-- value of zeta(50) - 1 from wolfram alpha (up to a Double)+-- http://www.wolframalpha.com/input/?i=zeta(50)-1+zeta50MinusOne :: Assertion+zeta50MinusOne = t @?= 8.8817842109308159e-16+ where+ t = getRationalLimit (==) . map (subtract 1) . rationalApproximations $ zeta50 :: Double++-- http://www.wolframalpha.com/input/?i=zeta(200)-1+zeta200MinusOne :: Assertion+zeta200MinusOne = t @?= 6.2230152778611417071e-61+ where+ t = getRationalLimit (==) . map (subtract 1) . rationalApproximations $ zeta200 :: Double++-- test against bbp formula+prop :: Positive Integer -> Bool+prop (Positive n) = getDigit n == getDigitBBP (n - 1)++tests :: TestTree+tests = testGroup "Rational approximation tests"+ [ testGroup "π² tests" [ testCase "matches double precision" piSquaredDouble+ , testCase "matches start of wolfram alpha" piSquaredWAstart+ , testCase "matches end of wolfram alpha" piSquaredWAend+ , testCase "first term matches chudnovsky" firstApproximation+ ]+ , testGroup "π tests" [ testCase "matches double precision" piDouble+ , testCase "matches oeis digits" piMatchesOeis+ , testCase "has feynman point" feynmanPoint+ , testCase "matches end of wolfram alpha" piWAend+ ]+ , testGroup "πᵏ tests" [ testCase "digits near 1000 of k=3" pi3+ , testCase "digits near 1000 of k=-1" piNegOne+ , testCase "digits near 1000 of k=10" pi10+ , testCase "digits near 1000 of k=-10" piNeg10+ , testCase "digits near 400 of k=51" pi51+ , testCase "digits near 400 of k=-51" piNeg51+ , testCase "ζ(50)-1 double precision" zeta50MinusOne+ , testCase "ζ(500)-1 double precision" zeta200MinusOne+ ]+ , testProperty "hex digits match BBP formula" prop+ ]++main :: IO ()+main = defaultMain tests
test-suite/TestUtils.hs view
@@ -1,63 +1,63 @@-{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE DataKinds #-} -module TestUtils - ( getValue - , getDigit - , getDigitBBP - , E - ) where - -import Data.Proxy (Proxy) -import Data.List (foldl') -import Data.Fixed (mod', HasResolution(..), Fixed) - -import GHC.TypeLits (Nat, KnownNat, SomeNat(..), natVal, someNatVal) - -import Data.ExactPi - --- E n generalises E2/E3/E6/E12 from Data.Fixed to give more precise --- fixed-precision arithmetic: Fixed (E 30) has 30 decimal places. -data E (n :: Nat) - -instance KnownNat n => HasResolution (E n) where - resolution _ = 10^natVal (undefined :: E n) - --- this function is not necessarily in general safe but is fine in the cases used here -getValue :: (Eq a, Fractional a) => ExactPi -> a -getValue = getRationalLimit (==) . rationalApproximations - -getDigit :: Integer -> Int -getDigit n = case someNatVal d of - Just (SomeNat (_ :: Proxy m)) -> (floor $ 16^n * (getValue (Exact 1 1) :: Fixed (E m))) `mod` 16 - Nothing -> error "negative digit requested" - where d = fromInteger $ 4 * n `div` 3 + 1 --------------------------------------------------------------------------------- -powModInteger :: Integer -> Integer -> Integer -> Integer -powModInteger a k n = a^k `mod` n - -infTerms :: Integer -> Int -> Integer -> Float -infTerms n j k = 16^^(n-k) / (fromIntegral $ 8*k + fromIntegral j) - -finiteTerms :: Integer -> Int -> Integer -> Float -finiteTerms n j k = (fromIntegral $ powModInteger 16 (n-k) (8*k + j')) / (fromIntegral $ 8*k + j') - where j' = fromIntegral j - -summation :: Integer -> Int -> Float -summation n j = stabilise $ scanl plus finitePart [infTerms n j k | k <- [n+1..]] - where finitePart = foldl' plus 0 [finiteTerms n j k | k <- [0..n]] - -mod1 :: Float -> Float -mod1 x = mod' x 1 - -plus :: Float -> Float -> Float -plus x y = mod1 (x + y) - -stabilise :: Eq a => [a] -> a -stabilise (x:y:xs) - | x == y = x - | otherwise = stabilise (y:xs) -stabilise _ = error "finite list" - -getDigitBBP :: Integer -> Int -getDigitBBP n = floor . (16 *) . mod1 $ 4 * summation n 1 - 2 * summation n 4 - summation n 5 - summation n 6 +{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DataKinds #-}+module TestUtils+ ( getValue+ , getDigit+ , getDigitBBP+ , E+ ) where++import Data.Proxy (Proxy)+import Data.List (foldl')+import Data.Fixed (mod', HasResolution(..), Fixed)++import GHC.TypeLits (Nat, KnownNat, SomeNat(..), natVal, someNatVal)++import Data.ExactPi++-- E n generalises E2/E3/E6/E12 from Data.Fixed to give more precise+-- fixed-precision arithmetic: Fixed (E 30) has 30 decimal places.+data E (n :: Nat)++instance KnownNat n => HasResolution (E n) where+ resolution _ = 10^natVal (undefined :: E n)++-- this function is not necessarily in general safe but is fine in the cases used here+getValue :: (Eq a, Fractional a) => ExactPi -> a+getValue = getRationalLimit (==) . rationalApproximations++getDigit :: Integer -> Int+getDigit n = case someNatVal d of+ Just (SomeNat (_ :: Proxy m)) -> (floor $ 16^n * (getValue (Exact 1 1) :: Fixed (E m))) `mod` 16+ Nothing -> error "negative digit requested"+ where d = fromInteger $ 4 * n `div` 3 + 1+--------------------------------------------------------------------------------+powModInteger :: Integer -> Integer -> Integer -> Integer+powModInteger a k n = a^k `mod` n++infTerms :: Integer -> Int -> Integer -> Float+infTerms n j k = 16^^(n-k) / (fromIntegral $ 8*k + fromIntegral j)++finiteTerms :: Integer -> Int -> Integer -> Float+finiteTerms n j k = (fromIntegral $ powModInteger 16 (n-k) (8*k + j')) / (fromIntegral $ 8*k + j')+ where j' = fromIntegral j++summation :: Integer -> Int -> Float+summation n j = stabilise $ scanl plus finitePart [infTerms n j k | k <- [n+1..]]+ where finitePart = foldl' plus 0 [finiteTerms n j k | k <- [0..n]]++mod1 :: Float -> Float+mod1 x = mod' x 1++plus :: Float -> Float -> Float+plus x y = mod1 (x + y)++stabilise :: Eq a => [a] -> a+stabilise (x:y:xs)+ | x == y = x+ | otherwise = stabilise (y:xs)+stabilise _ = error "finite list"++getDigitBBP :: Integer -> Int+getDigitBBP n = floor . (16 *) . mod1 $ 4 * summation n 1 - 2 * summation n 4 - summation n 5 - summation n 6