exact-pi 0.4.1.4 → 0.5.1.0
raw patch · 8 files changed
Files
- README.md +0/−1
- Setup.hs +0/−2
- changelog.md +17/−0
- exact-pi.cabal +52/−34
- src/Data/ExactPi.hs +56/−23
- src/Data/ExactPi/TypeLevel.hs +1/−1
- test-suite/Test.hs +146/−0
- test-suite/TestUtils.hs +63/−0
README.md view
@@ -1,6 +1,5 @@ # 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
@@ -1,2 +0,0 @@-import Distribution.Simple-main = defaultMain
changelog.md view
@@ -1,3 +1,20 @@+0.5.1.0+-------+* Fix warnings.++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.
exact-pi.cabal view
@@ -1,37 +1,55 @@-name: exact-pi-version: 0.4.1.4-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+cabal-version: >=1.10+name: exact-pi+version: 0.5.1.0+license: MIT+license-file: LICENSE+maintainer: douglas.mcclean@gmail.com+author: Douglas McClean+tested-with:+ ghc ==9.14.1 ghc ==9.12.2 ghc ==9.10.3 ghc ==9.8.4 ghc ==9.6.7 ghc ==9.4.8+ ghc ==9.2.8 ghc ==9.0.2 ghc ==8.10.7 ghc ==8.8.4 ghc ==8.6.5 ghc ==8.4.4++homepage: https://github.com/dmcclean/exact-pi/+bug-reports: https://github.com/dmcclean/exact-pi/issues/+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.++category: Data+build-type: Simple+extra-source-files:+ README.md+ changelog.md++source-repository head+ type: git+ location: https://github.com/dmcclean/exact-pi.git+ library- exposed-modules: Data.ExactPi,- Data.ExactPi.TypeLevel- build-depends: base >=4.7 && <5,- numtype-dk >= 0.5- if impl(ghc <8.0)+ exposed-modules:+ Data.ExactPi+ Data.ExactPi.TypeLevel++ hs-source-dirs: src+ default-language: Haskell2010+ ghc-options: -Wall build-depends:- semigroups >=0.8- ghc-options: -Wall- hs-source-dirs: src- default-language: Haskell2010+ base >=4.11 && <5,+ numtype-dk >=0.5 && <0.6,+ infinite-list <0.2 -source-repository head- type: git- location: https://github.com/dmcclean/exact-pi.git+test-suite spec+ type: exitcode-stdio-1.0+ main-is: Test.hs+ hs-source-dirs: test-suite+ other-modules: TestUtils+ default-language: Haskell2010+ ghc-options: -Wall+ build-depends:+ base >=4.11 && <5,+ exact-pi,+ QuickCheck >=2.10,+ tasty >=0.10,+ tasty-hunit >=0.9 && <0.11,+ tasty-quickcheck >=0.9 && <0.12
src/Data/ExactPi.hs view
@@ -1,4 +1,6 @@-{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ParallelListComp #-}+{-# LANGUAGE PostfixOperators #-} {-# OPTIONS_HADDOCK show-extensions #-} @@ -28,10 +30,14 @@ toExactInteger, isExactRational, toExactRational,- rationalApproximations+ rationalApproximations,+ -- * Utils+ getRationalLimit ) where +import Data.List.Infinite (Infinite(..), (...))+import qualified Data.List.Infinite as Inf import Data.Monoid import Data.Ratio ((%), numerator, denominator) import Data.Semigroup@@ -46,7 +52,7 @@ -- 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 (Exact z q) = (pi ^^ z) * fromRational q approximateValue (Approximate x) = x -- | Identifies whether an 'ExactPi' is an exact or approximate representation of zero.@@ -95,29 +101,57 @@ toExactRational (Exact 0 q) = Just q toExactRational _ = Nothing --- | Converts an 'ExactPi' to a list of increasingly accurate rational approximations, on alternating--- sides of the actual value. 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.+-- | 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 based on work by Anders Kaseorg shared at http://qr.ae/RbXl8M.+-- Implementation is based on Chudnovsky's algorithm. rationalApproximations :: ExactPi -> [Rational] rationalApproximations (Approximate x) = [toRational (x :: Double)]-rationalApproximations (Exact 0 q) = [q]-rationalApproximations (Exact z q) = fmap (\pi' -> q * (pi' ^^ z)) piConvergents+rationalApproximations (Exact _ 0) = [0]+rationalApproximations (Exact 0 q) = [q]+rationalApproximations (Exact z q)+ | even z = Inf.toList $ fmap (\c -> q * 10005^^k * c^^z) chudnovsky+ | otherwise = Inf.toList $ Inf.zipWith (\c r -> q * 10005^^k * c^^z * r) chudnovsky rootApproximation+ where k = z `div` 2++chudnovsky :: Infinite Rational+chudnovsky = fmap (426880 /) partials where- piConvergents :: [Rational]- piConvergents = go True 2 4 where- go s p' q' | ltPi m = [q' | not s] ++ go True m q'- | otherwise = [p' | s] ++ go False p' m where- m = (numerator p' + numerator q')%(denominator p' + denominator q')- ltPi :: Rational -> Bool- ltPi x = ok x 1 where- ok y i =- y <= (27*i - 12)%5 ||- (y < (675*i - 216)%125 &&- ok ((y - fromInteger (5*i - 2))*(3*(3*i + 1)*(3*i + 2)%(i*(2*i - 1))))- (i + 1))+ lk = Inf.iterate (+545140134) 13591409+ xk = Inf.iterate (*(-262537412640768000)) 1+ kk = Inf.iterate (+12) 6+ mk = 1 :< Inf.zipWith3 (\m k n -> m * ((k^(3::Int) - 16*k) % (n+1)^(3::Int))) mk kk (0...)+ values = Inf.zipWith3 (\m l x -> m * l / x) mk lk xk+ partials = Inf.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 :: Infinite Rational+rootApproximation = fmap Inf.head . Inf.iterate (Inf.drop 4) $ go 1 0 100 1 40+ where+ go :: Integer -> Integer -> Integer -> Integer -> Integer -> Infinite Rational+ 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@@ -168,9 +202,8 @@ -- | 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
@@ -4,10 +4,10 @@ {-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE KindSignatures #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-}+{-# LANGUAGE UndecidableInstances #-} #if __GLASGOW_HASKELL__ > 805 {-# LANGUAGE NoStarIsType #-}
+ test-suite/Test.hs view
@@ -0,0 +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 = take 1 (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
@@ -0,0 +1,63 @@+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DataKinds #-}+module TestUtils+ ( getValue+ , getDigit+ , getDigitBBP+ , E+ ) where++import Prelude hiding (Foldable(..))+import Data.Fixed (mod', HasResolution(..), Fixed)+import Data.Foldable+import Data.Proxy (Proxy)+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