diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,6 +1,5 @@ # exact-pi Exact rational multiples of pi (and integer powers of pi) in Haskell -[![Build Status](https://travis-ci.org/dmcclean/exact-pi.svg?branch=master)](https://travis-ci.org/dmcclean/exact-pi) [![Hackage Version](https://img.shields.io/hackage/v/exact-pi.svg)](http://hackage.haskell.org/package/exact-pi) [![Stackage version](https://www.stackage.org/package/exact-pi/badge/lts?label=Stackage)](https://www.stackage.org/package/exact-pi) diff --git a/Setup.hs b/Setup.hs deleted file mode 100644 --- a/Setup.hs +++ /dev/null @@ -1,2 +0,0 @@ -import Distribution.Simple -main = defaultMain diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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. diff --git a/exact-pi.cabal b/exact-pi.cabal --- a/exact-pi.cabal +++ b/exact-pi.cabal @@ -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 diff --git a/src/Data/ExactPi.hs b/src/Data/ExactPi.hs --- a/src/Data/ExactPi.hs +++ b/src/Data/ExactPi.hs @@ -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 = (*) diff --git a/src/Data/ExactPi/TypeLevel.hs b/src/Data/ExactPi/TypeLevel.hs --- a/src/Data/ExactPi/TypeLevel.hs +++ b/src/Data/ExactPi/TypeLevel.hs @@ -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 #-} diff --git a/test-suite/Test.hs b/test-suite/Test.hs new file mode 100644 --- /dev/null +++ b/test-suite/Test.hs @@ -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 diff --git a/test-suite/TestUtils.hs b/test-suite/TestUtils.hs new file mode 100644 --- /dev/null +++ b/test-suite/TestUtils.hs @@ -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