exact-pi 0.5.0.0 → 0.5.0.1
raw patch · 9 files changed
+694/−689 lines, 9 filesdep ~QuickCheckdep ~basedep ~tastysetup-changedPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: QuickCheck, base, tasty
API changes (from Hackage documentation)
Files
- LICENSE +20/−20
- README.md +6/−6
- Setup.hs +2/−2
- changelog.md +62/−57
- exact-pi.cabal +54/−54
- src/Data/ExactPi.hs +205/−205
- src/Data/ExactPi/TypeLevel.hs +136/−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,57 +1,62 @@-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.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.0-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.12,- 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.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
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 'chudnovksy' 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,136 @@-{-# 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 #-} + +#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