arithmoi 0.12.1.0 → 0.13.0.0

raw patch · 58 files changed

+600/−663 lines, 58 filesdep +ghc-bignumdep +infinite-listdep −integer-gmpdep ~QuickCheckdep ~basedep ~chimeraPVP ok

version bump matches the API change (PVP)

Dependencies added: ghc-bignum, infinite-list

Dependencies removed: integer-gmp

Dependency ranges changed: QuickCheck, base, chimera, constraints, containers, deepseq, exact-pi, integer-logarithms, integer-roots, mod, quickcheck-classes, random, semirings, tasty, tasty-bench, tasty-rerun, vector

API changes (from Hackage documentation)

Files

• Math/NumberTheory/ArithmeticFunctions/Class.hs +0/−8
• Math/NumberTheory/ArithmeticFunctions/Inverse.hs +0/−4
• Math/NumberTheory/ArithmeticFunctions/Moebius.hs +2/−6
• Math/NumberTheory/Diophantine.hs +10/−2
• Math/NumberTheory/DirichletCharacters.hs +4/−23
• Math/NumberTheory/Euclidean/Coprimes.hs +0/−4
• Math/NumberTheory/Moduli/Equations.hs +5/−2
• Math/NumberTheory/Moduli/Internal.hs +41/−9
• Math/NumberTheory/Moduli/JacobiSymbol.hs +0/−4
• Math/NumberTheory/Moduli/Singleton.hs +0/−6
• Math/NumberTheory/Moduli/SomeMod.hs +0/−4
• Math/NumberTheory/Moduli/Sqrt.hs +31/−17
• Math/NumberTheory/Powers/Modular.hs +0/−84
• Math/NumberTheory/Primes.hs +12/−6
• Math/NumberTheory/Primes/Factorisation/Montgomery.hs +12/−15
• Math/NumberTheory/Primes/IntSet.hs +0/−7
• Math/NumberTheory/Primes/Testing/Certified.hs +10/−7
• Math/NumberTheory/Primes/Testing/Probabilistic.hs +7/−6
• Math/NumberTheory/Quadratic/EisensteinIntegers.hs +20/−8
• Math/NumberTheory/Quadratic/GaussianIntegers.hs +16/−7
• Math/NumberTheory/Recurrences/Bilinear.hs +82/−77
• Math/NumberTheory/Recurrences/Linear.hs +10/−7
• Math/NumberTheory/Recurrences/Pentagonal.hs +17/−15
• Math/NumberTheory/Utils.hs +29/−33
• Math/NumberTheory/Zeta/Dirichlet.hs +20/−17
• Math/NumberTheory/Zeta/Hurwitz.hs +31/−27
• Math/NumberTheory/Zeta/Riemann.hs +14/−10
• Math/NumberTheory/Zeta/Utils.hs +7/−15
• arithmoi.cabal +23/−22
• benchmark/Math/NumberTheory/RecurrencesBench.hs +10/−8
• benchmark/Math/NumberTheory/SieveBlockBench.hs +0/−4
• benchmark/Math/NumberTheory/ZetaBench.hs +3/−2
• changelog.md +11/−0
• test-suite/Math/NumberTheory/ArithmeticFunctions/InverseTests.hs +7/−6
• test-suite/Math/NumberTheory/ArithmeticFunctions/MertensTests.hs +0/−4
• test-suite/Math/NumberTheory/ArithmeticFunctions/SieveBlockTests.hs +0/−4
• test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs +2/−3
• test-suite/Math/NumberTheory/CurvesTests.hs +0/−1
• test-suite/Math/NumberTheory/DiophantineTests.hs +0/−2
• test-suite/Math/NumberTheory/EisensteinIntegersTests.hs +4/−3
• test-suite/Math/NumberTheory/EuclideanTests.hs +0/−1
• test-suite/Math/NumberTheory/GaussianIntegersTests.hs +6/−5
• test-suite/Math/NumberTheory/Moduli/CbrtTests.hs +0/−5
• test-suite/Math/NumberTheory/Moduli/JacobiTests.hs +0/−5
• test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs +0/−1
• test-suite/Math/NumberTheory/Moduli/SqrtTests.hs +0/−1
• test-suite/Math/NumberTheory/PrefactoredTests.hs +0/−2
• test-suite/Math/NumberTheory/Primes/SieveTests.hs +0/−1
• test-suite/Math/NumberTheory/PrimesTests.hs +0/−4
• test-suite/Math/NumberTheory/Recurrences/BilinearTests.hs +94/−92
• test-suite/Math/NumberTheory/Recurrences/LinearTests.hs +14/−15
• test-suite/Math/NumberTheory/Recurrences/PentagonalTests.hs +6/−5
• test-suite/Math/NumberTheory/SmoothNumbersTests.hs +6/−4
• test-suite/Math/NumberTheory/TestUtils.hs +0/−1
• test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +0/−1
• test-suite/Math/NumberTheory/UniqueFactorisationTests.hs +0/−1
• test-suite/Math/NumberTheory/Zeta/DirichletTests.hs +18/−16
• test-suite/Math/NumberTheory/Zeta/RiemannTests.hs +16/−14

@@ -8,7 +8,6 @@ -- factorisation domains. -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE GADTs               #-}  module Math.NumberTheory.ArithmeticFunctions.Class@@ -18,9 +17,6 @@   ) where  import Control.Applicative-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import Prelude hiding (Applicative(..))  import Math.NumberTheory.Primes@@ -66,11 +62,7 @@  instance Monoid a => Monoid (ArithmeticFunction n a) where   mempty  = pure mempty-#if __GLASGOW_HASKELL__ < 803-  mappend = liftA2 mappend-#else   mappend = (<>)-#endif  -- | Factorisation is expensive, so it is better to avoid doing it twice. -- Write 'runFunction (f + g) n' instead of 'runFunction f n + runFunction g n'.

@@ -9,7 +9,6 @@ -- <https://www.emis.de/journals/JIS/VOL19/Alekseyev/alek5.pdf Computing the Inverses, their Power Sums, and Extrema for Euler’s Totient and Other Multiplicative Functions> -- by M. A. Alekseyev. -{-# LANGUAGE CPP                 #-} {-# LANGUAGE FlexibleContexts    #-} {-# LANGUAGE RankNTypes          #-} {-# LANGUAGE ScopedTypeVariables #-}@@ -37,9 +36,6 @@ import qualified Data.Map as M import Data.Maybe import Data.Ord (Down(..))-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import Data.Semiring (Semiring(..), Mul(..)) import Data.Set (Set) import qualified Data.Set as S

@@ -8,7 +8,6 @@ --  {-# LANGUAGE BangPatterns          #-}-{-# LANGUAGE CPP                   #-} {-# LANGUAGE MagicHash             #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE TypeFamilies          #-}@@ -24,16 +23,13 @@ import Data.Bits import Data.Int import Data.Word-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import qualified Data.Vector.Generic         as G import qualified Data.Vector.Generic.Mutable as M import qualified Data.Vector.Primitive as P import qualified Data.Vector.Unboxed         as U import qualified Data.Vector.Unboxed.Mutable as MU import GHC.Exts-import GHC.Integer.GMP.Internals+import GHC.Num.Integer import Unsafe.Coerce  import Math.NumberTheory.Roots (integerSquareRoot)@@ -51,7 +47,7 @@  -- | Convert to any numeric type. runMoebius :: Num a => Moebius -> a-runMoebius m = fromInteger (S# (dataToTag# m -# 1#))+runMoebius m = fromInteger (IS (dataToTag# m -# 1#))  fromMoebius :: Moebius -> Int8 fromMoebius m = intToInt8 $ I# (dataToTag# m)

@@ -6,6 +6,9 @@   ) where +import Data.List.Infinite (Infinite(..))+import qualified Data.List.Infinite as Inf+ import           Math.NumberTheory.Moduli.Sqrt  ( sqrtsModFactorisation ) import           Math.NumberTheory.Primes       ( factorise                                                 , unPrime@@ -18,12 +21,17 @@ -- | as described at https://en.wikipedia.org/wiki/Cornacchia%27s_algorithm  cornacchiaPrimitive' :: Integer -> Integer -> [(Integer, Integer)] cornacchiaPrimitive' d m = concatMap-  (findSolution . head . dropWhile (\r -> r * r >= m) . gcdSeq m)+  (findSolution . Inf.head . Inf.dropWhile (\r -> r * r >= m) . gcdSeq m)   roots  where+  roots :: [Integer]   roots = filter (<= m `div` 2) $ sqrtsModFactorisation (m - d) (factorise m)-  gcdSeq a b = a : gcdSeq b (mod a b)++  gcdSeq :: Integer -> Integer -> Infinite Integer+  gcdSeq a b = a :< gcdSeq b (mod a b)+   -- If s = sqrt((m - r*r) / d) is an integer then (r, s) is a solution+  findSolution :: Integer -> [(Integer, Integer)]   findSolution r = [ (r, s) | rem1 == 0 && s * s == s2 ]    where     (s2, rem1) = divMod (m - r * r) d

@@ -7,10 +7,8 @@ -- Implementation and enumeration of Dirichlet characters. -- -{-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE GADTs                      #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE KindSignatures             #-} {-# LANGUAGE LambdaCase                 #-} {-# LANGUAGE PatternSynonyms            #-}@@ -59,9 +57,6 @@   , validChar   ) where -#if !MIN_VERSION_base(4,12,0)-import Control.Applicative                                 (liftA2)-#endif import Data.Bits                                           (Bits(..)) import Data.Constraint import Data.Foldable@@ -70,9 +65,7 @@ import Data.List                                           (sort, unfoldr) import Data.Maybe                                          (mapMaybe, fromJust, fromMaybe) import Data.Mod-#if MIN_VERSION_base(4,12,0) import Data.Monoid                                         (Ap(..))-#endif import Data.Proxy                                          (Proxy(..)) import Data.Ratio                                          ((%), numerator, denominator) import Data.Semigroup                                      (Semigroup(..),Product(..))@@ -143,10 +136,10 @@ generator :: Prime Natural -> Word -> Natural generator p k = case cyclicGroupFromFactors [(p, k)] of   Nothing -> error "illegal"-  Just (Some cg)-    | Sub Dict <- proofFromCyclicGroup cg ->-      unMod $ multElement $ unPrimitiveRoot $ head $-        mapMaybe (isPrimitiveRoot cg) [2..maxBound]+  Just (Some cg) -> case proofFromCyclicGroup cg of+    Sub Dict -> case mapMaybe (isPrimitiveRoot cg) [2..maxBound] of+      [] -> error "illegal"+      hd : _ -> unMod $ multElement $ unPrimitiveRoot hd  -- | Implement the function \(\lambda\) from page 5 of -- https://www2.eecs.berkeley.edu/Pubs/TechRpts/1984/CSD-84-186.pdf@@ -465,18 +458,6 @@                                   Just (i,_) -> Just (bit i :: Natural, TwoPower i a b)           where options = [(i, bit (i-2) :: Natural) | i <- [2..k]]                 worksb (_,phi) = phi `stimes` b == mempty--#if !MIN_VERSION_base(4,12,0)-newtype Ap f a = Ap { getAp :: f a }-  deriving (Eq, Functor, Applicative, Monad)--instance (Applicative f, Semigroup a) => Semigroup (Ap f a) where-  (<>) = liftA2 (<>)--instance (Applicative f, Semigroup a, Monoid a) => Monoid (Ap f a) where-  mempty = pure mempty-  mappend = (<>)-#endif  -- | Similar to Maybe, but with different Semigroup and Monoid instances. type OrZero a = Ap Maybe a

@@ -6,7 +6,6 @@ -- -- Container for pairwise coprime numbers. -{-# LANGUAGE CPP                 #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TupleSections       #-} @@ -23,9 +22,6 @@ import Data.Euclidean import Data.List (tails) import Data.Maybe-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import Data.Semiring (Semiring(..), isZero) import Data.Traversable 

@@ -7,8 +7,11 @@ -- Polynomial modular equations. -- +{-# LANGUAGE MagicHash           #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnboxedSums         #-} {-# LANGUAGE ViewPatterns        #-}+{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}  module Math.NumberTheory.Moduli.Equations   ( solveLinear@@ -18,7 +21,7 @@ import Data.Constraint import Data.Maybe import Data.Mod-import GHC.Integer.GMP.Internals+import GHC.Num.Integer import GHC.TypeNats (KnownNat, natVal)  import Math.NumberTheory.Moduli.Chinese@@ -130,7 +133,7 @@   | a `rem` p' == 0   = solveLinear' p' b c   | otherwise-  = map (\n -> (n - b) * recipModInteger (2 * a) p' `mod` p')+  = map (\n -> let (# t | #) = integerRecipMod# (2 * a) (fromInteger p') in (n - b) * toInteger t `mod` p')   $ sqrtsModPrime (b * b - 4 * a * c) p     where       p' :: Integer

@@ -8,7 +8,10 @@ --  {-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MagicHash #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnboxedTuples #-}+{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}  module Math.NumberTheory.Moduli.Internal   ( isPrimitiveRoot'@@ -20,7 +23,7 @@ import Data.Mod import Data.Proxy import GHC.TypeNats (SomeNat(..), someNatVal)-import GHC.Integer.GMP.Internals+import GHC.Num.Integer import Numeric.Natural  import Math.NumberTheory.Moduli.Chinese@@ -64,7 +67,9 @@     thetaA     = theta p pkMinusOne a     thetaB     = theta p pkMinusOne b     pkMinusOne = p^(k-1)-    c          = (recipModInteger thetaA pkMinusOne * thetaB) `rem` pkMinusOne+    c          = (toInteger t * thetaB) `rem` pkMinusOne+      where+        (# t | #) = integerRecipMod# thetaA (fromInteger pkMinusOne)     (result, pkMinusPk1) = fromJust $ chinese (baseSol, p-1) (c, pkMinusOne)  -- compute the homomorphism theta given in https://math.stackexchange.com/a/1864495/418148@@ -74,7 +79,9 @@   where     pk           = pkMinusOne * p     p2kMinusOne  = pkMinusOne * pk-    numerator    = (powModInteger a (pk - pkMinusOne) p2kMinusOne - 1) `rem` p2kMinusOne+    numerator    = (toInteger t  - 1) `rem` p2kMinusOne+      where+        (# t | #) = integerPowMod# a (pk - pkMinusOne) (fromInteger p2kMinusOne)  -- TODO: Use Pollig-Hellman to reduce the problem further into groups of prime order. -- While Bach reduction simplifies the problem into groups of the form (Z/pZ)*, these@@ -88,14 +95,34 @@   | otherwise     = discreteLogarithmPrimePollard p a b  discreteLogarithmPrimeBSGS :: Int -> Int -> Int -> Int-discreteLogarithmPrimeBSGS p a b = head [i*m + j | (v,i) <- zip giants [0..m-1], j <- maybeToList (M.lookup v table)]+discreteLogarithmPrimeBSGS p a b =+  case [i*m + j | (v,i) <- zip giants [0..m-1], j <- maybeToList (M.lookup v table)] of+    [] -> error ("discreteLogarithmPrimeBSGS: failed, please report this as a bug. Inputs: " ++ show [p,a,b])+    hd : _ -> hd   where+    m :: Int     m        = integerSquareRoot (p - 2) + 1 -- simple way of ceiling (sqrt (p-1))++    babies :: [Int]     babies   = iterate (.* a) 1++    table :: M.Map Int Int     table    = M.fromList (zip babies [0..m-1])-    aInv     = recipModInteger (toInteger a) (toInteger p)-    bigGiant = fromInteger $ powModInteger aInv (toInteger m) (toInteger p)++    aInv :: Integer+    aInv     = fromIntegral ap+      where+        (# ap | #) = integerRecipMod# (toInteger a) (fromIntegral p)++    bigGiant :: Int+    bigGiant = fromIntegral aInvmp+      where+        (# aInvmp | #) = integerPowMod# aInv (toInteger m) (fromIntegral p)++    giants :: [Int]     giants   = iterate (.* bigGiant) b++    (.*) :: Int -> Int -> Int     x .* y   = x * y `rem` p  -- TODO: Use more advanced walks, in order to reduce divisions, cf@@ -106,7 +133,7 @@ discreteLogarithmPrimePollard p a b =   case concatMap runPollard [(x,y) | x <- [0..n], y <- [0..n]] of     (t:_)  -> fromInteger t-    []     -> error ("discreteLogarithm: pollard's rho failed, please report this as a bug. inputs " ++ show [p,a,b])+    []     -> error ("discreteLogarithm: pollard's rho failed, please report this as a bug. Inputs: " ++ show [p,a,b])   where     n                 = p-1 -- order of the cyclic group     halfN             = n `quot` 2@@ -116,9 +143,14 @@                           0 -> (xi*xi `rem` p, mul2 ai, mul2 bi)                           1 -> ( a*xi `rem` p,    ai+1,      bi)                           _ -> ( b*xi `rem` p,      ai,    bi+1)-    initialise (x,y)  = (powModInteger a x n * powModInteger b y n `rem` n, x, y)+    initialise (x,y)  = (toInteger axn * toInteger byn `rem` n, x, y)+      where+        (# axn | #) = integerPowMod# a x (fromInteger n)+        (# byn | #) = integerPowMod# b y (fromInteger n)     begin t           = go (step t) (step (step t))-    check t           = powModInteger a t p == b+    check t           = case integerPowMod# a t (fromInteger p) of+      (# atp | #) -> toInteger atp == b+      (# | _ #) -> False     go tort@(xi,ai,bi) hare@(x2i,a2i,b2i)       | xi == x2i, gcd (bi - b2i) n < sqrtN = case someNatVal (fromInteger n) of         SomeNat (Proxy :: Proxy n) -> map (toInteger . unMod) $ solveLinear (fromInteger (bi - b2i) :: Mod n) (fromInteger (ai - a2i))

@@ -11,7 +11,6 @@ --  {-# LANGUAGE BangPatterns #-}-{-# LANGUAGE CPP          #-} {-# LANGUAGE LambdaCase   #-}  module Math.NumberTheory.Moduli.JacobiSymbol@@ -21,9 +20,6 @@   ) where  import Data.Bits-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import Numeric.Natural  import Math.NumberTheory.Utils

@@ -7,7 +7,6 @@ -- Singleton data types. -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE DataKinds           #-} {-# LANGUAGE DeriveGeneric       #-} {-# LANGUAGE FlexibleInstances   #-}@@ -51,9 +50,6 @@ import Data.List (sort) import qualified Data.Map as M import Data.Proxy-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import GHC.Generics import GHC.TypeNats (KnownNat, Nat, natVal) import Numeric.Natural@@ -315,6 +311,4 @@ pattern CGDoubleOddPrimePower :: Prime a -> Word -> CyclicGroup a m pattern CGDoubleOddPrimePower p k <- CGDoubleOddPrimePower' p k -#if __GLASGOW_HASKELL__ > 801 {-# COMPLETE CG2, CG4, CGOddPrimePower, CGDoubleOddPrimePower #-}-#endif

@@ -7,7 +7,6 @@ -- Safe modular arithmetic with modulo on type level. -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE GADTs               #-} {-# LANGUAGE LambdaCase          #-} {-# LANGUAGE RankNTypes          #-}@@ -23,9 +22,6 @@ import Data.Euclidean (GcdDomain(..), Euclidean(..), Field) import Data.Mod import Data.Proxy-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import Data.Semiring (Semiring(..), Ring(..)) import Data.Type.Equality import GHC.TypeNats (KnownNat, SomeNat(..), sameNat, natVal, someNatVal)

@@ -8,11 +8,15 @@ -- <https://en.wikipedia.org/wiki/Jacobi_symbol Jacobi symbol>. -- -{-# LANGUAGE BangPatterns  #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE PostfixOperators #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-}-{-# LANGUAGE ViewPatterns  #-}+{-# LANGUAGE ViewPatterns #-}  module Math.NumberTheory.Moduli.Sqrt   ( -- * Modular square roots@@ -29,10 +33,13 @@ import Control.Monad (liftM2) import Data.Bits import Data.Constraint+import Data.List.Infinite (Infinite(..), (...))+import qualified Data.List.Infinite as Inf import Data.Maybe import Data.Mod import Data.Proxy-import GHC.TypeNats (KnownNat, SomeNat(..), natVal, someNatVal)+import GHC.TypeNats (KnownNat, SomeNat(..), natVal, someNatVal, Nat)+import Numeric.Natural (Natural)  import Math.NumberTheory.Moduli.Chinese import Math.NumberTheory.Moduli.JacobiSymbol@@ -139,15 +146,21 @@     prime = natVal square  -- | @p@ must be of form @4k + 1@-sqrtOfMinusOne :: KnownNat p => Mod p-sqrtOfMinusOne = res+sqrtOfMinusOne :: forall (p :: Nat). KnownNat p => Mod p+sqrtOfMinusOne = case results of+  [] -> error "sqrtOfMinusOne: internal invariant violated"+  hd : _ -> hd   where-    p = natVal res+    p :: Natural+    p = natVal (Proxy :: Proxy p)++    k :: Natural     k = (p - 1) `quot` 4-    res = head-      $ dropWhile (\n -> n == 1 || n == maxBound)-      $ map (^ k) [2 .. maxBound - 1] +    results :: [Mod p]+    results = dropWhile (\n -> n == 1 || n == maxBound) $+      map (^ k) [2 .. maxBound - 1]+ -- | @tonelliShanks square prime@ calculates a square root of @square@ --   modulo @prime@, where @prime@ is a prime of the form @4*k + 1@ and --   @square@ is a positive quadratic residue modulo @prime@, using the@@ -243,17 +256,18 @@ rem8 :: Integral a => a -> Int rem8 n = fromIntegral n .&. 7 -findNonSquare :: KnownNat n => Mod n-findNonSquare = res+findNonSquare :: forall (n :: Nat). KnownNat n => Mod n+findNonSquare+  | rem8 n == 3 || rem8 n == 5 = 2+  | otherwise = fromIntegral $ Inf.head $+    Inf.dropWhile (\p -> jacobi p n /= MinusOne) candidates   where-    n = natVal res-    res-      | rem8 n == 3 || rem8 n == 5 = 2-      | otherwise = fromIntegral $ head $-        dropWhile (\p -> jacobi p n /= MinusOne) candidates+    n = natVal (Proxy :: Proxy n)      -- It is enough to consider only prime candidates, but     -- the probability that the smallest non-residue is > 67     -- is small and 'jacobi' test is fast,     -- so we use [71..n] instead of filter isPrime [71..n].-    candidates = 3:5:7:11:13:17:19:23:29:31:37:41:43:47:53:59:61:67:[71..n]+    candidates :: Infinite Natural+    candidates = 3 :< 5 :< 7 :< 11 :< 13 :< 17 :< 19 :< 23 :< 29 :< 31 :<+      37 :< 41 :< 43 :< 47 :< 53 :< 59 :< 61 :< 67 :< (71...)

@@ -1,84 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers.Modular--- Copyright:   (c) 2017 Andrew Lelechenko--- Licence:     MIT--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ Modular powers (a. k. a. modular exponentiation).-----module Math.NumberTheory.Powers.Modular-  {-# DEPRECATED "Use Data.Mod or Data.Mod.Word instead" #-}-  ( powMod-  , powModWord-  , powModInt-  ) where--import GHC.Natural (powModNatural)-import qualified GHC.Integer.GMP.Internals as GMP (powModInteger)-import Math.NumberTheory.Utils.FromIntegral---- | @powMod@ @b@ @e@ @m@ computes (@b^e@) \`mod\` @m@ in effective way.--- An exponent @e@ must be non-negative, a modulo @m@ must be positive.--- Otherwise the behaviour of @powMod@ is undefined.------ >>> powMod 2 3 5--- 3--- >>> powMod 3 12345678901234567890 1001--- 1------ See also 'Math.NumberTheory.Moduli.Class.powMod' and 'Math.NumberTheory.Moduli.Class.powSomeMod'.------ For finite numeric types ('Int', 'Word', etc.)--- modulo @m@ should be such that @m^2@ does not overflow,--- otherwise the behaviour is undefined. If you--- need both to fit into machine word and to handle large moduli,--- take a look at 'powModInt' and 'powModWord'.------ >>> powMod 3 101 (2^60-1 :: Integer) -- correct--- 1018105167100379328--- >>> powMod 3 101 (2^60-1 :: Int) -- incorrect due to overflow--- 1115647832265427613--- >>> powModInt 3 101 (2^60-1 :: Int) -- correct--- 1018105167100379328-powMod :: (Integral a, Integral b) => a -> b -> a -> a-powMod x y m-  | m <= 0    = error "powModInt: non-positive modulo"-  | y <  0    = error "powModInt: negative exponent"-  | otherwise = f (x `rem` m) y 1 `mod` m-  where-    f _ 0 acc = acc-    f b e acc = f (b * b `rem` m) (e `quot` 2)-      (if odd e then b * acc `rem` m else acc)--{-# INLINE [1] powMod #-}-{-# RULES-"powMod/Integer" powMod = powModInteger-  #-}---- Work around https://ghc.haskell.org/trac/ghc/ticket/14085-powModInteger :: Integer -> Integer -> Integer -> Integer-powModInteger b e m = GMP.powModInteger (b `mod` m) e m--{-# RULES-"powMod/Natural" powMod = powModNatural-"powMod/Word"    powMod = powModWord-"powMod/Int"     powMod = powModInt-  #-}---- | Specialised version of 'powMod', able to handle large moduli correctly.------ >>> powModWord 3 101 (2^60-1)--- 1018105167100379328-powModWord :: Word -> Word -> Word -> Word-powModWord b e m = fromInteger $ GMP.powModInteger (toInteger b) (toInteger e) (toInteger m)---- | Specialised version of 'powMod', able to handle large moduli correctly.------ >>> powModInt 3 101 (2^60-1)--- 1018105167100379328-powModInt :: Int -> Int -> Int -> Int-powModInt x y m-  | m <= 0 = error "powModInt: non-positive modulo"-  | y <  0 = error "powModInt: negative exponent"-  | otherwise = wordToInt $ powModWord (intToWord (x `mod` m)) (intToWord y) (intToWord m)

@@ -7,6 +7,7 @@  {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE LambdaCase        #-}+{-# LANGUAGE PostfixOperators  #-}  {-# OPTIONS_GHC -fno-warn-orphans #-} @@ -24,6 +25,8 @@  import Data.Bits import Data.Coerce+import Data.List.Infinite (Infinite(..), (...), (....))+import qualified Data.List.Infinite as Inf import Data.Maybe import Data.Word import Numeric.Natural@@ -112,8 +115,8 @@   | n <= 2    = Prime 2   | n <= 3    = Prime 3   | n <= 5    = Prime 5-  | otherwise = head $ mapMaybe isPrime $-                  dropWhile (< n) $ map fromWheel30 [toWheel30 n ..]+  | otherwise = Inf.head $ mapMaybeInf isPrime $+                  Inf.dropWhile (< n) $ fmap fromWheel30 (toWheel30 n ...)                   -- dropWhile is important, because fromWheel30 (toWheel30 n) may appear to be < n.                   -- E. g., fromWheel30 (toWheel30 94) == 97 @@ -127,11 +130,14 @@   | n < 3     = Prime 2   | n < 5     = Prime 3   | n < 7     = Prime 5-  | otherwise = head $ mapMaybe isPrime $-                  dropWhile (> n) $ map fromWheel30 [toWheel30 n, toWheel30 n - 1 ..]+  | otherwise = Inf.head $ mapMaybeInf isPrime $+                  Inf.dropWhile (> n) $ fmap fromWheel30 ((toWheel30 n, toWheel30 n - 1) ....)                   -- dropWhile is important, because fromWheel30 (toWheel30 n) may appear to be > n.                   -- E. g., fromWheel30 (toWheel30 100) == 101 +mapMaybeInf :: (a -> Maybe b) -> Infinite a -> Infinite b+mapMaybeInf = Inf.foldr . (maybe id (:<) .)+ ------------------------------------------------------------------------------- -- Prime sequences @@ -153,7 +159,7 @@   Prime 2 -> Prime 3   Prime 3 -> Prime 5   Prime 5 -> Prime 7-  Prime p -> head $ mapMaybe (isPrime . fromWheel30) [toWheel30 p + 1 ..]+  Prime p -> Inf.head $ mapMaybeInf (isPrime . fromWheel30) ((toWheel30 p + 1) ...)  succGenericBounded   :: (Bits a, Integral a, UniqueFactorisation a, Bounded a)@@ -173,7 +179,7 @@   Prime 3 -> Prime 2   Prime 5 -> Prime 3   Prime 7 -> Prime 5-  Prime p -> head $ mapMaybe (isPrime . fromWheel30) [toWheel30 p - 1, toWheel30 p - 2 ..]+  Prime p -> Inf.head $ mapMaybeInf (isPrime . fromWheel30) ((toWheel30 p - 1, toWheel30 p - 2) ....)  -- 'dropWhile' is important, because 'psieveFrom' can actually contain primes less than p. enumFromGeneric :: Integral a => Prime a -> [Prime a]

@@ -43,12 +43,9 @@ import Data.Maybe import Data.Mod import Data.Proxy-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import Data.Traversable import GHC.Exts-import GHC.Integer.GMP.Internals hiding (integerToInt, wordToInteger)+import GHC.Num.BigNat import GHC.Natural import GHC.TypeNats (KnownNat, SomeNat(..), natVal, someNatVal) import System.Random@@ -334,7 +331,7 @@     pThen = multiply thn  p     pStep = multiply step p -    progression = pFrom : pThen : zipWith (`add` pStep) progression (tail progression)+    progression = pFrom : pThen : zipWith (`add` pStep) progression (drop 1 progression)  -- primes, compactly stored as a bit sieve primeStore :: [PrimeSieve]@@ -353,32 +350,32 @@   NatS# n#  -> case shiftToOddCount# n# of     (# 0##, m# #) -> goWord m# 1     (# k#,  m# #) -> (2, W# k#) <: goWord m# 1-  NatJ# n -> case shiftToOddCountBigNat n of-    (0, m) -> goBigNat m 1-    (k, m) -> (2, k) <: goBigNat m 1+  NatJ# (BN# n) -> case shiftToOddCountBigNat n of+    (# 0, m #) -> goBigNat m 1+    (# k, m #) -> (2, k) <: goBigNat m 1   where     x <: ~(l,b) = (x:l,b)      !(Ptr smallPrimesAddr#) = smallPrimesPtr -    goBigNat :: BigNat -> Int -> ([(Natural, Word)], Maybe Natural)+    goBigNat :: BigNat# -> Int -> ([(Natural, Word)], Maybe Natural)     goBigNat !m i@(I# i#)-      | isTrue# (sizeofBigNat# m ==# 1#)-      = goWord (bigNatToWord m) i+      | isTrue# (bigNatSize# m ==# 1#)+      = goWord (bigNatToWord# m) i       | i >= smallPrimesLength-      = ([], Just (NatJ# m))+      = ([], Just (NatJ# (BN# m)))       | otherwise       = let p# = #if MIN_VERSION_base(4,16,0)               word16ToWord# #endif               (indexWord16OffAddr# smallPrimesAddr# i#) in-      case m `quotRemBigNatWord` p# of+      case m `bigNatQuotRemWord#` p# of         (# mp, 0## #) ->-          let (# k, r #) = splitOff 1 mp in+          let !(# k, r #) = splitOff 1 mp in             (NatS# p#, k) <: goBigNat r (i + 1)           where-            splitOff !k x = case x `quotRemBigNatWord` p# of+            splitOff !k x = case x `bigNatQuotRemWord#` p# of               (# xp, 0## #) -> splitOff (k + 1) xp               _             -> (# k, x #)         _ -> goBigNat m (i + 1)

@@ -13,7 +13,6 @@ --  {-# LANGUAGE BangPatterns               #-}-{-# LANGUAGE CPP                        #-} {-# LANGUAGE DeriveDataTypeable         #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE ScopedTypeVariables        #-}@@ -179,15 +178,9 @@ isProperSubsetOf :: PrimeIntSet -> PrimeIntSet -> Bool isProperSubsetOf = coerce IS.isProperSubsetOf -#if MIN_VERSION_containers(0,5,11) -- | Check whether two sets are disjoint. disjoint :: PrimeIntSet -> PrimeIntSet -> Bool disjoint = coerce IS.disjoint-#else--- | Check whether two sets are disjoint.-disjoint :: PrimeIntSet -> PrimeIntSet -> Bool-disjoint (PrimeIntSet x) (PrimeIntSet y) = IS.null (IS.intersection x y)-#endif  -- | Difference between a set of primes and a set of integers. difference :: PrimeIntSet -> IntSet -> PrimeIntSet

@@ -6,8 +6,11 @@ -- -- Deterministic primality testing. +{-# LANGUAGE MagicHash           #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnboxedTuples       #-} {-# LANGUAGE ViewPatterns        #-}+{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}  module Math.NumberTheory.Primes.Testing.Certified   ( isCertifiedPrime@@ -17,7 +20,7 @@ import Data.Bits ((.&.)) import Data.Mod import Data.Proxy-import GHC.Integer.GMP.Internals (powModInteger)+import GHC.Num.Integer import GHC.TypeNats (SomeNat(..), someNatVal)  import Math.NumberTheory.Roots (integerSquareRoot)@@ -67,10 +70,10 @@     pm1 = p-1     ppProd pps = product [pf^e | (pf,e,_,_) <- pps]     verify (pf,_,base,proof) = pf == cprime proof && crit pf base && checkPrimalityProof proof-    crit pf base = gcd p (x-1) == 1 && y == 1+    crit pf base = gcd p (toInteger x-1) == 1 && y == 1       where-        x = powModInteger base (pm1 `quot` pf) p-        y = powModInteger x pf p+        (# x | #) = integerPowMod# base (pm1 `quot` pf) (fromInteger p)+        (# y | #) = integerPowMod# (toInteger x) pf (fromInteger p)  -- | @'isTrivialPrime'@ checks whether its argument is a trivially --   known prime.@@ -144,9 +147,9 @@                             Just ppr ->(p,e,bs,ppr)               where                 q = nm1 `quot` p-                x = powModInteger bs q n-                y = powModInteger x p n-                g = gcd n (x-1)+                (# x | #) = integerPowMod# bs q (fromInteger n)+                (# y | #) = integerPowMod# (toInteger x) p (fromInteger n)+                g = gcd n (toInteger x-1)  -- | Find a decomposition of p-1 for the pocklington certificate. --   Usually bloody slow if p-1 has two (or more) /large/ prime divisors.

@@ -22,8 +22,9 @@ import Data.Bits import Data.Mod import Data.Proxy+import GHC.Num.BigNat+import GHC.Num.Integer import GHC.Exts (Word(..), Int(..), (-#), (<#), isTrue#)-import GHC.Integer.GMP.Internals import GHC.TypeNats (KnownNat, SomeNat(..), someNatVal)  import Math.NumberTheory.Moduli.JacobiSymbol@@ -174,7 +175,7 @@  -- n odd positive, n > abs q, index odd testLucas :: Integer -> Integer -> Integer -> (Integer, Integer, Integer)-testLucas n q (S# i#) = look (finiteBitSize (0 :: Word) - 2)+testLucas n q (IS i#) = look (finiteBitSize (0 :: Word) - 2)   where     j = I# i#     look k@@ -192,10 +193,10 @@           v2n1  = ((un1 - (2*q)*un)*vn-qn) `rem` n           q2n   = (qn*qn) `rem` n           q2n1  = (qn*qn*q) `rem` n-testLucas n q (Jp# bn#) = test (s# -# 1#)+testLucas n q (IP bn#) = test (s# -# 1#)   where-    s# = sizeofBigNat# bn#-    test j# = case indexBigNat# bn# j# of+    s# = bigNatSize# bn#+    test j# = case bigNatIndex# bn# j# of                 0## -> test (j# -# 1#)                 w# -> look (j# -# 1#) (W# w#) (finiteBitSize (0 :: Word) - 1)     look j# w i@@ -204,7 +205,7 @@     go k# w i un un1 vn qn       | i < 0       = if isTrue# (k# <# 0#)                          then (un,vn,qn)-                         else go (k# -# 1#) (W# (indexBigNat# bn# k#)) (finiteBitSize (0 :: Word) - 1) un un1 vn qn+                         else go (k# -# 1#) (W# (bigNatIndex# bn# k#)) (finiteBitSize (0 :: Word) - 1) un un1 vn qn       | testBit w i = go k# w (i-1) u2n1 u2n2 v2n1 q2n1       | otherwise   = go k# w (i-1) u2n u2n1 v2n q2n         where

@@ -10,6 +10,7 @@  {-# LANGUAGE BangPatterns   #-} {-# LANGUAGE DeriveGeneric  #-}+{-# LANGUAGE PostfixOperators #-} {-# LANGUAGE RankNTypes     #-} {-# LANGUAGE TypeFamilies   #-} @@ -30,7 +31,10 @@ import Control.DeepSeq import Data.Coerce import Data.Euclidean-import Data.List                                       (mapAccumL, partition)+import Data.List                                       (mapAccumL)+import Data.List.Infinite (Infinite(..), (...))+import qualified Data.List.Infinite as Inf+import Data.List.NonEmpty (NonEmpty(..)) import Data.Maybe import Data.Ord                                        (comparing) import qualified Data.Semiring as S@@ -191,14 +195,22 @@ -- -- >>> take 10 primes -- [Prime 2+ω,Prime 2,Prime 3+2*ω,Prime 3+ω,Prime 4+3*ω,Prime 4+ω,Prime 5+3*ω,Prime 5+2*ω,Prime 5,Prime 6+5*ω]-primes :: [Prime EisensteinInteger]-primes = coerce $ (2 :+ 1) : mergeBy (comparing norm) l r+primes :: Infinite (Prime EisensteinInteger)+primes = coerce $ (2 :+ 1) :< mergeBy (comparing norm) l r   where-    leftPrimes, rightPrimes :: [Prime Integer]-    (leftPrimes, rightPrimes) = partition (\p -> unPrime p `mod` 3 == 2) [U.nextPrime 2 ..]-    rightPrimes' = filter (\prime -> unPrime prime `mod` 3 == 1) $ tail rightPrimes-    l = [unPrime p :+ 0 | p <- leftPrimes]-    r = [g | p <- rightPrimes', let x :+ y = unPrime (findPrime p), g <- [x :+ y, x :+ (x - y)]]+    leftPrimes, rightPrimes :: Infinite (Prime Integer)+    (leftPrimes, rightPrimes) = Inf.partition (\p -> unPrime p `mod` 3 == 2) (U.nextPrime 2...)++    rightPrimes' :: Infinite (Prime Integer)+    rightPrimes' = Inf.filter (\prime -> unPrime prime `mod` 3 == 1) $ Inf.tail rightPrimes++    l :: Infinite EisensteinInteger+    l = fmap (\p -> unPrime p :+ 0) leftPrimes++    r :: Infinite EisensteinInteger+    r = Inf.concatMap+        (\p -> let x :+ y = unPrime (findPrime p) in (x :+ y) :| [x :+ (x - y)])+        rightPrimes'  -- | [Implementation notes for factorise function] --

@@ -9,6 +9,7 @@ --  {-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE PostfixOperators #-} {-# LANGUAGE TypeFamilies  #-}  module Math.NumberTheory.Quadratic.GaussianIntegers (@@ -24,7 +25,10 @@ import Control.DeepSeq (NFData) import Data.Coerce import Data.Euclidean-import Data.List (mapAccumL, partition)+import Data.List (mapAccumL)+import Data.List.Infinite (Infinite(..), (...))+import qualified Data.List.Infinite as Inf+import Data.List.NonEmpty (NonEmpty(..)) import Data.Maybe import Data.Ord (comparing) import qualified Data.Semiring as S@@ -131,14 +135,19 @@ -- -- >>> take 10 primes -- [Prime 1+ι,Prime 2+ι,Prime 1+2*ι,Prime 3,Prime 3+2*ι,Prime 2+3*ι,Prime 4+ι,Prime 1+4*ι,Prime 5+2*ι,Prime 2+5*ι]-primes :: [U.Prime GaussianInteger]-primes = coerce $ (1 :+ 1) : mergeBy (comparing norm) l r+primes :: Infinite (U.Prime GaussianInteger)+primes = coerce $ (1 :+ 1) :< mergeBy (comparing norm) l r   where-    leftPrimes, rightPrimes :: [Prime Integer]-    (leftPrimes, rightPrimes) = partition (\p -> unPrime p `mod` 4 == 3) [U.nextPrime 3 ..]-    l = [unPrime p :+ 0 | p <- leftPrimes]-    r = [g | p <- rightPrimes, let Prime (x :+ y) = findPrime p, g <- [x :+ y, y :+ x]]+    leftPrimes, rightPrimes :: Infinite (Prime Integer)+    (leftPrimes, rightPrimes) = Inf.partition (\p -> unPrime p `mod` 4 == 3) (U.nextPrime 3 ...) +    l :: Infinite (GaussianInteger)+    l = fmap (\p -> unPrime p :+ 0) leftPrimes++    r :: Infinite (GaussianInteger)+    r = Inf.concatMap+        (\p -> let x :+ y = unPrime (findPrime p) in (x :+ y) :| [y :+ x])+        rightPrimes  -- |Find a Gaussian integer whose norm is the given prime number -- of form 4k + 1 using

@@ -26,6 +26,8 @@ -- 1  {-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE PostfixOperators    #-}+{-# LANGUAGE RankNTypes          #-} {-# LANGUAGE ScopedTypeVariables #-}  module Math.NumberTheory.Recurrences.Bilinear@@ -49,6 +51,9 @@  import Data.Euclidean (GcdDomain(..)) import Data.List (scanl', zipWith4)+import Data.List.Infinite (Infinite(..), (...))+import qualified Data.List.Infinite as Inf+import Data.List.NonEmpty (NonEmpty(..)) import Data.Maybe import Data.Ratio import Data.Semiring (Semiring(..))@@ -67,12 +72,12 @@ -- -- >>> take 6 binomial :: [[Int]] -- [[1],[1,1],[1,2,1],[1,3,3,1],[1,4,6,4,1],[1,5,10,10,5,1]]-binomial :: Semiring a => [[a]]-binomial = iterate (\l -> zipWith plus (l ++ [zero]) (zero : l)) [one]-{-# SPECIALIZE binomial :: [[Int]]     #-}-{-# SPECIALIZE binomial :: [[Word]]    #-}-{-# SPECIALIZE binomial :: [[Integer]] #-}-{-# SPECIALIZE binomial :: [[Natural]] #-}+binomial :: Semiring a => Infinite [a]+binomial = Inf.iterate (\l -> zipWith plus (l ++ [zero]) (zero : l)) [one]+{-# SPECIALIZE binomial :: Infinite [Int]     #-}+{-# SPECIALIZE binomial :: Infinite [Word]    #-}+{-# SPECIALIZE binomial :: Infinite [Integer] #-}+{-# SPECIALIZE binomial :: Infinite [Natural] #-}  -- | Pascal triangle, rotated by 45 degrees. --@@ -84,12 +89,12 @@ -- -- >>> take 6 (map (take 6) binomialRotated) :: [[Int]] -- [[1,1,1,1,1,1],[1,2,3,4,5,6],[1,3,6,10,15,21],[1,4,10,20,35,56],[1,5,15,35,70,126],[1,6,21,56,126,252]]-binomialRotated :: Semiring a => [[a]]-binomialRotated = iterate (tail . scanl' plus zero) (repeat one)-{-# SPECIALIZE binomialRotated :: [[Int]]     #-}-{-# SPECIALIZE binomialRotated :: [[Word]]    #-}-{-# SPECIALIZE binomialRotated :: [[Integer]] #-}-{-# SPECIALIZE binomialRotated :: [[Natural]] #-}+binomialRotated :: Semiring a => Infinite (Infinite a)+binomialRotated = Inf.iterate (Inf.tail . Inf.scanl' plus zero) (Inf.repeat one)+{-# SPECIALIZE binomialRotated :: Infinite (Infinite Int)     #-}+{-# SPECIALIZE binomialRotated :: Infinite (Infinite Word)    #-}+{-# SPECIALIZE binomialRotated :: Infinite (Infinite Integer) #-}+{-# SPECIALIZE binomialRotated :: Infinite (Infinite Natural) #-}  -- | The n-th (zero-based) line of 'binomial' -- (and the n-th diagonal of 'binomialRotated').@@ -111,15 +116,15 @@ -- -- >>> take 6 (binomialDiagonal 5) -- [1,6,21,56,126,252]-binomialDiagonal :: (Enum a, GcdDomain a) => a -> [a]-binomialDiagonal n = scanl'+binomialDiagonal :: (Enum a, GcdDomain a) => a -> Infinite a+binomialDiagonal n = Inf.scanl'   (\x k -> fromJust (x `times` (n `plus` k) `divide` k))   one-  [one..]-{-# SPECIALIZE binomialDiagonal :: Int     -> [Int]     #-}-{-# SPECIALIZE binomialDiagonal :: Word    -> [Word]    #-}-{-# SPECIALIZE binomialDiagonal :: Integer -> [Integer] #-}-{-# SPECIALIZE binomialDiagonal :: Natural -> [Natural] #-}+  (one...)+{-# SPECIALIZE binomialDiagonal :: Int     -> Infinite Int     #-}+{-# SPECIALIZE binomialDiagonal :: Word    -> Infinite Word    #-}+{-# SPECIALIZE binomialDiagonal :: Integer -> Infinite Integer #-}+{-# SPECIALIZE binomialDiagonal :: Natural -> Infinite Natural #-}  -- | Prime factors of a binomial coefficient. --@@ -153,14 +158,14 @@ -- takes O(k n^2 ln n) time and forces thunks @stirling1 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. -- -- One could also consider 'Math.Combinat.Numbers.unsignedStirling1st' from <http://hackage.haskell.org/package/combinat combinat> package to compute stand-alone values.-stirling1 :: (Num a, Enum a) => [[a]]-stirling1 = scanl f [1] [0..]+stirling1 :: (Num a, Enum a) => Infinite [a]+stirling1 = Inf.scanl f [1] (0...)   where     f xs n = 0 : zipIndexedListWithTail (\_ x y -> x + n * y) 1 xs 0-{-# SPECIALIZE stirling1 :: [[Int]]     #-}-{-# SPECIALIZE stirling1 :: [[Word]]    #-}-{-# SPECIALIZE stirling1 :: [[Integer]] #-}-{-# SPECIALIZE stirling1 :: [[Natural]] #-}+{-# SPECIALIZE stirling1 :: Infinite [Int]     #-}+{-# SPECIALIZE stirling1 :: Infinite [Word]    #-}+{-# SPECIALIZE stirling1 :: Infinite [Integer] #-}+{-# SPECIALIZE stirling1 :: Infinite [Natural] #-}  -- | Infinite zero-based table of <https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind Stirling numbers of the second kind>. --@@ -171,14 +176,14 @@ -- takes O(k n^2 ln n) time and forces thunks @stirling2 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. -- -- One could also consider 'Math.Combinat.Numbers.stirling2nd' from <http://hackage.haskell.org/package/combinat combinat> package to compute stand-alone values.-stirling2 :: (Num a, Enum a) => [[a]]-stirling2 = iterate f [1]+stirling2 :: (Num a, Enum a) => Infinite [a]+stirling2 = Inf.iterate f [1]   where     f xs = 0 : zipIndexedListWithTail (\k x y -> x + k * y) 1 xs 0-{-# SPECIALIZE stirling2 :: [[Int]]     #-}-{-# SPECIALIZE stirling2 :: [[Word]]    #-}-{-# SPECIALIZE stirling2 :: [[Integer]] #-}-{-# SPECIALIZE stirling2 :: [[Natural]] #-}+{-# SPECIALIZE stirling2 :: Infinite [Int]     #-}+{-# SPECIALIZE stirling2 :: Infinite [Word]    #-}+{-# SPECIALIZE stirling2 :: Infinite [Integer] #-}+{-# SPECIALIZE stirling2 :: Infinite [Natural] #-}  -- | Infinite one-based table of <https://en.wikipedia.org/wiki/Lah_number Lah numbers>. -- @lah !! n !! k@ equals to lah(n + 1, k + 1).@@ -188,15 +193,15 @@ -- -- Complexity: @lah !! n !! k@ is O(n ln n) bits long, its computation -- takes O(k n ln n) time and forces thunks @lah !! n !! i@ for @0 <= i <= k@.-lah :: Integral a => [[a]]+lah :: Integral a => Infinite [a] -- Implementation was derived from code by https://github.com/grandpascorpion-lah = zipWith f (tail factorial) [1..]+lah = Inf.zipWith f (Inf.tail factorial) (1...)   where     f nf n = scanl (\x k -> x * (n - k) `div` (k * (k + 1))) nf [1..n-1]-{-# SPECIALIZE lah :: [[Int]]     #-}-{-# SPECIALIZE lah :: [[Word]]    #-}-{-# SPECIALIZE lah :: [[Integer]] #-}-{-# SPECIALIZE lah :: [[Natural]] #-}+{-# SPECIALIZE lah :: Infinite [Int]     #-}+{-# SPECIALIZE lah :: Infinite [Word]    #-}+{-# SPECIALIZE lah :: Infinite [Integer] #-}+{-# SPECIALIZE lah :: Infinite [Natural] #-}  -- | Infinite zero-based table of <https://en.wikipedia.org/wiki/Eulerian_number Eulerian numbers of the first kind>. --@@ -206,14 +211,14 @@ -- Complexity: @eulerian1 !! n !! k@ is O(n ln n) bits long, its computation -- takes O(k n^2 ln n) time and forces thunks @eulerian1 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. ---eulerian1 :: (Num a, Enum a) => [[a]]-eulerian1 = scanl f [] [1..]+eulerian1 :: (Num a, Enum a) => Infinite [a]+eulerian1 = Inf.scanl f [] (1...)   where     f xs n = 1 : zipIndexedListWithTail (\k x y -> (n - k) * x + (k + 1) * y) 1 xs 0-{-# SPECIALIZE eulerian1 :: [[Int]]     #-}-{-# SPECIALIZE eulerian1 :: [[Word]]    #-}-{-# SPECIALIZE eulerian1 :: [[Integer]] #-}-{-# SPECIALIZE eulerian1 :: [[Natural]] #-}+{-# SPECIALIZE eulerian1 :: Infinite [Int]     #-}+{-# SPECIALIZE eulerian1 :: Infinite [Word]    #-}+{-# SPECIALIZE eulerian1 :: Infinite [Integer] #-}+{-# SPECIALIZE eulerian1 :: Infinite [Natural] #-}  -- | Infinite zero-based table of <https://en.wikipedia.org/wiki/Eulerian_number#Eulerian_numbers_of_the_second_kind Eulerian numbers of the second kind>. --@@ -223,14 +228,14 @@ -- Complexity: @eulerian2 !! n !! k@ is O(n ln n) bits long, its computation -- takes O(k n^2 ln n) time and forces thunks @eulerian2 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. ---eulerian2 :: (Num a, Enum a) => [[a]]-eulerian2 = scanl f [] [1..]+eulerian2 :: (Num a, Enum a) => Infinite [a]+eulerian2 = Inf.scanl f [] (1...)   where     f xs n = 1 : zipIndexedListWithTail (\k x y -> (2 * n - k - 1) * x + (k + 1) * y) 1 xs 0-{-# SPECIALIZE eulerian2 :: [[Int]]     #-}-{-# SPECIALIZE eulerian2 :: [[Word]]    #-}-{-# SPECIALIZE eulerian2 :: [[Integer]] #-}-{-# SPECIALIZE eulerian2 :: [[Natural]] #-}+{-# SPECIALIZE eulerian2 :: Infinite [Int]     #-}+{-# SPECIALIZE eulerian2 :: Infinite [Word]    #-}+{-# SPECIALIZE eulerian2 :: Infinite [Integer] #-}+{-# SPECIALIZE eulerian2 :: Infinite [Natural] #-}  -- | Infinite zero-based sequence of <https://en.wikipedia.org/wiki/Bernoulli_number Bernoulli numbers>, -- computed via <https://en.wikipedia.org/wiki/Bernoulli_number#Connection_with_Stirling_numbers_of_the_second_kind connection>@@ -243,10 +248,10 @@ -- takes O(n^3 ln n) time and forces thunks @stirling2 !! i !! j@ for @0 <= i <= n@ and @0 <= j <= i@. -- -- One could also consider 'Math.Combinat.Numbers.bernoulli' from <http://hackage.haskell.org/package/combinat combinat> package to compute stand-alone values.-bernoulli :: Integral a => [Ratio a]-bernoulli = helperForBEEP id (map recip [1..])-{-# SPECIALIZE bernoulli :: [Ratio Int] #-}-{-# SPECIALIZE bernoulli :: [Rational] #-}+bernoulli :: Integral a => Infinite (Ratio a)+bernoulli = helperForBEEP id (Inf.map recip (1...))+{-# SPECIALIZE bernoulli :: Infinite (Ratio Int) #-}+{-# SPECIALIZE bernoulli :: Infinite (Rational) #-}  -- | <https://en.wikipedia.org/wiki/Faulhaber%27s_formula Faulhaber's formula>. --@@ -259,10 +264,10 @@ faulhaberPoly p   = zipWith (*) ((0:)   $ reverse-  $ take (p + 1) bernoulli)+  $ Inf.take (p + 1) bernoulli)   $ map (% (fromIntegral p+1))   $ zipWith (*) (iterate negate (if odd p then 1 else -1))-  $ binomial !! (p+1)+  $ binomial Inf.!! (fromIntegral (p+1))  -- | Infinite zero-based list of <https://en.wikipedia.org/wiki/Euler_number Euler numbers>. -- The algorithm used was derived from <http://www.emis.ams.org/journals/JIS/VOL4/CHEN/AlgBE2.pdf Algorithms for Bernoulli numbers and Euler numbers>@@ -271,28 +276,28 @@ -- -- >>> take 10 euler' :: [Rational] -- [1 % 1,0 % 1,(-1) % 1,0 % 1,5 % 1,0 % 1,(-61) % 1,0 % 1,1385 % 1,0 % 1]-euler' :: forall a . Integral a => [Ratio a]-euler' = tail $ helperForBEEP tail as+euler' :: forall a . Integral a => Infinite (Ratio a)+euler' = Inf.tail $ helperForBEEP (drop 1) as   where-    as :: [Ratio a]-    as = zipWith3+    as :: Infinite (Ratio a)+    as = Inf.zipWith3         (\sgn frac ones -> (sgn * ones) % frac)-        (cycle [1, 1, 1, 1, -1, -1, -1, -1])-        (dups (iterate (2 *) 1))-        (cycle [1, 1, 1, 0])+        (Inf.cycle (1 :| [1, 1, 1, -1, -1, -1, -1]))+        (dups (Inf.iterate (2 *) 1))+        (Inf.cycle (1 :| [1, 1, 0])) -    dups :: forall x . [x] -> [x]-    dups = foldr (\n list -> n : n : list) []-{-# SPECIALIZE euler' :: [Ratio Int]     #-}-{-# SPECIALIZE euler' :: [Rational]      #-}+    dups :: forall x . Infinite x -> Infinite x+    dups = Inf.foldr (\n list -> n :< n :< list)+{-# SPECIALIZE euler' :: Infinite (Ratio Int)     #-}+{-# SPECIALIZE euler' :: Infinite (Rational)      #-}  -- | The same sequence as @euler'@, but with type @[a]@ instead of @[Ratio a]@ -- as the denominators in @euler'@ are always @1@. -- -- >>> take 10 euler :: [Integer] -- [1,0,-1,0,5,0,-61,0,1385,0]-euler :: forall a . Integral a => [a]-euler = map numerator euler'+euler :: forall a . Integral a => Infinite a+euler = Inf.map numerator euler'  -- | Infinite zero-based list of the @n@-th order Euler polynomials evaluated at @1@. -- The algorithm used was derived from <http://www.emis.ams.org/journals/JIS/VOL4/CHEN/AlgBE2.pdf Algorithms for Bernoulli numbers and Euler numbers>@@ -302,10 +307,10 @@ -- -- >>> take 10 eulerPolyAt1 :: [Rational] -- [1 % 1,1 % 2,0 % 1,(-1) % 4,0 % 1,1 % 2,0 % 1,(-17) % 8,0 % 1,31 % 2]-eulerPolyAt1 :: forall a . Integral a => [Ratio a]-eulerPolyAt1 = tail $ helperForBEEP tail (map recip (iterate (2 *) 1))-{-# SPECIALIZE eulerPolyAt1 :: [Ratio Int]     #-}-{-# SPECIALIZE eulerPolyAt1 :: [Rational]      #-}+eulerPolyAt1 :: forall a . Integral a => Infinite (Ratio a)+eulerPolyAt1 = Inf.tail $ helperForBEEP (drop 1) (Inf.map recip (Inf.iterate (2 *) 1))+{-# SPECIALIZE eulerPolyAt1 :: Infinite (Ratio Int)     #-}+{-# SPECIALIZE eulerPolyAt1 :: Infinite (Rational)      #-}  ------------------------------------------------------------------------------- -- Utils@@ -334,9 +339,9 @@ -- @stirling2 !! n@, while @euler, eulerPolyAt1@ only use -- @tail $ stirling2 !! n@. As such, this argument serves to pass @id@ -- in the former case, and @tail@ in the latter.-helperForBEEP :: Integral a => ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]-helperForBEEP g xs = map (f . g) stirling2+helperForBEEP :: Integral a => (forall b. [b] -> [b]) -> Infinite (Ratio a) -> Infinite (Ratio a)+helperForBEEP g xs = Inf.map (f . g) stirling2   where-    f = sum . zipWith4 (\sgn fact x stir -> sgn * fact * x * stir) (cycle [1, -1]) factorial xs-{-# SPECIALIZE helperForBEEP :: ([Ratio Int] -> [Ratio Int]) -> [Ratio Int] -> [Ratio Int] #-}-{-# SPECIALIZE helperForBEEP :: ([Rational] -> [Rational]) -> [Rational] -> [Rational]     #-}+    f = sum . zipWith4 (\sgn fact x stir -> sgn * fact * x * stir) (cycle [1, -1]) (Inf.toList factorial) (Inf.toList xs)+{-# SPECIALIZE helperForBEEP :: (forall b. [b] -> [b]) -> Infinite (Ratio Int) -> Infinite (Ratio Int) #-}+{-# SPECIALIZE helperForBEEP :: (forall b. [b] -> [b]) -> Infinite (Rational) -> Infinite (Rational)     #-}

@@ -6,7 +6,8 @@ -- -- Efficient calculation of linear recurrent sequences, including Fibonacci and Lucas sequences. -{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE BangPatterns     #-}+{-# LANGUAGE PostfixOperators #-}  module Math.NumberTheory.Recurrences.Linear   ( factorial@@ -19,6 +20,8 @@   ) where  import Data.Bits+import Data.List.Infinite (Infinite(..), (...))+import qualified Data.List.Infinite as Inf import Numeric.Natural import Math.NumberTheory.Primes @@ -29,12 +32,12 @@ -- -- The time-and-space behaviour of 'factorial' is similar to described in -- "Math.NumberTheory.Recurrences.Bilinear#memory".-factorial :: (Num a, Enum a) => [a]-factorial = scanl (*) 1 [1..]-{-# SPECIALIZE factorial :: [Int]     #-}-{-# SPECIALIZE factorial :: [Word]    #-}-{-# SPECIALIZE factorial :: [Integer] #-}-{-# SPECIALIZE factorial :: [Natural] #-}+factorial :: (Num a, Enum a) => Infinite a+factorial = Inf.scanl (*) 1 (1...)+{-# SPECIALIZE factorial :: Infinite Int     #-}+{-# SPECIALIZE factorial :: Infinite Word    #-}+{-# SPECIALIZE factorial :: Infinite Integer #-}+{-# SPECIALIZE factorial :: Infinite Natural #-}  -- | Prime factors of a factorial. --

@@ -7,6 +7,7 @@ -- Values of <https://en.wikipedia.org/wiki/Partition_(number_theory)#Partition_function partition function>. -- +{-# LANGUAGE PostfixOperators #-} {-# LANGUAGE TypeApplications #-}  module Math.NumberTheory.Recurrences.Pentagonal@@ -14,6 +15,8 @@   ) where  import qualified Data.Chimera as Ch+import Data.List.Infinite (Infinite(..), (...))+import qualified Data.List.Infinite as Inf import Data.Vector (Vector) import Numeric.Natural (Natural) @@ -22,13 +25,10 @@ -- -- >>> take 10 pents -- [0,1,2,5,7,12,15,22,26,35]-pents :: (Enum a, Num a) => [a]-pents = interleave (scanl (\acc n -> acc + 3 * n - 1) 0 [1..])-                   (scanl (\acc n -> acc + 3 * n - 2) 1 [2..])-  where-    interleave :: [a] -> [a] -> [a]-    interleave (n : ns) (m : ms) = n : m : interleave ns ms-    interleave _ _ = []+pents :: (Enum a, Num a) => Infinite a+pents = Inf.interleave+  (Inf.scanl (\acc n -> acc + 3 * n - 1) 0 (1...))+  (Inf.scanl (\acc n -> acc + 3 * n - 2) 1 (2...))  -- | @p(n) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + p(n-11) + ...@, where @p(0) = 1@ -- and @p(k) = 0@ for a negative integer @k@. Uses a @Chimera@ from the@@ -39,8 +39,8 @@   = sum   $ zipWith (*) (cycle [1, 1, -1, -1])   $ map (f . (n -))-  $ takeWhile (<= n)-  $ tail pents+  $ Inf.takeWhile (<= n)+  $ Inf.tail pents  -- | Infinite zero-based table of <https://oeis.org/A000041 partition numbers>. --@@ -51,9 +51,11 @@ -- >>> import Data.Mod -- >>> partition !! 1000 :: Mod 1000 -- (991 `modulo` 1000)-partition :: Num a => [a]-partition = Ch.toList $ Ch.tabulateFix @Vector partitionF-{-# SPECIALIZE partition :: [Int]     #-}-{-# SPECIALIZE partition :: [Word]    #-}-{-# SPECIALIZE partition :: [Integer] #-}-{-# SPECIALIZE partition :: [Natural] #-}+partition :: Num a => Infinite a+partition = Inf.tabulate (Ch.index ch)+  where+    ch = Ch.tabulateFix @Vector partitionF+{-# SPECIALIZE partition :: Infinite Int     #-}+{-# SPECIALIZE partition :: Infinite Word    #-}+{-# SPECIALIZE partition :: Infinite Integer #-}+{-# SPECIALIZE partition :: Infinite Natural #-}

@@ -40,9 +40,12 @@  import Data.Bits import Data.Euclidean+import Data.List.Infinite (Infinite(..)) import Data.Semiring (Semiring(..), isZero) import GHC.Base-import GHC.Integer.GMP.Internals+import GHC.Num.BigNat+import GHC.Num.Integer+import GHC.Num.Natural import qualified Math.NumberTheory.Utils.FromIntegral as UT import GHC.Natural import GHC.TypeNats@@ -77,16 +80,16 @@ -- | Specialised version for @'Integer'@. --   Precondition: argument nonzero (not checked). shiftOCInteger :: Integer -> (Word, Integer)-shiftOCInteger n@(S# i#) =+shiftOCInteger n@(IS i#) =     case shiftToOddCount# (int2Word# i#) of       (# 0##, _ #) -> (0, n)-      (# z#, w# #) -> (W# z#, wordToInteger w#)-shiftOCInteger n@(Jp# bn#) = case bigNatZeroCount bn# of+      (# z#, w# #) -> (W# z#, integerFromWord# w#)+shiftOCInteger n@(IP bn#) = case bigNatZeroCount bn# of                                  0## -> (0, n)-                                 z#  -> (W# z#, bigNatToInteger (bn# `shiftRBigNat` word2Int# z#))-shiftOCInteger n@(Jn# bn#) = case bigNatZeroCount bn# of+                                 z#  -> (W# z#, integerFromBigNat# (bn# `bigNatShiftR#` z#))+shiftOCInteger n@(IN bn#) = case bigNatZeroCount bn# of                                  0## -> (0, n)-                                 z#  -> (W# z#, bigNatToNegInteger (bn# `shiftRBigNat` word2Int# z#))+                                 z#  -> (W# z#, integerFromBigNatNeg# (bn# `bigNatShiftR#` z#))  -- | Specialised version for @'Natural'@. --   Precondition: argument nonzero (not checked).@@ -95,23 +98,23 @@     case shiftToOddCount# i# of       (# 0##, _ #) -> (0, n)       (# z#, w# #) -> (W# z#, NatS# w#)-shiftOCNatural n@(NatJ# bn#) = case bigNatZeroCount bn# of+shiftOCNatural n@(NatJ# (BN# bn#)) = case bigNatZeroCount bn# of                                  0## -> (0, n)-                                 z#  -> (W# z#, bigNatToNatural (bn# `shiftRBigNat` word2Int# z#))+                                 z#  -> (W# z#, naturalFromBigNat# (bn# `bigNatShiftR#` z#)) -shiftToOddCountBigNat :: BigNat -> (Word, BigNat)+shiftToOddCountBigNat :: BigNat# -> (# Word, BigNat# #) shiftToOddCountBigNat bn# = case bigNatZeroCount bn# of-  0## -> (0, bn#)-  z#  -> (W# z#, bn# `shiftRBigNat` word2Int# z#)+  0## -> (# 0, bn# #)+  z#  -> (# W# z#, bn# `bigNatShiftR#` z# #)  -- | Count trailing zeros in a @'BigNat'@. --   Precondition: argument nonzero (not checked, Integer invariant).-bigNatZeroCount :: BigNat -> Word#+bigNatZeroCount :: BigNat# -> Word# bigNatZeroCount bn# = count 0## 0#   where     !(W# bitSize#) = intToWord (finiteBitSize (0 :: Word))     count a# i# =-          case indexBigNat# bn# i# of+          case bigNatIndex# bn# i# of             0## -> count (a# `plusWord#` bitSize#) (i# +# 1#)             w#  -> a# `plusWord#` ctz# w# @@ -139,13 +142,13 @@ -- | Specialised version for @'Int'@. --   Precondition: argument nonzero (not checked). shiftOInteger :: Integer -> Integer-shiftOInteger (S# i#) = wordToInteger (shiftToOdd# (int2Word# i#))-shiftOInteger n@(Jp# bn#) = case bigNatZeroCount bn# of+shiftOInteger (IS i#) = integerFromWord# (shiftToOdd# (int2Word# i#))+shiftOInteger n@(IP bn#) = case bigNatZeroCount bn# of                                  0## -> n-                                 z#  -> bigNatToInteger (bn# `shiftRBigNat` word2Int# z#)-shiftOInteger n@(Jn# bn#) = case bigNatZeroCount bn# of+                                 z#  -> integerFromBigNat# (bn# `bigNatShiftR#` z#)+shiftOInteger n@(IN bn#) = case bigNatZeroCount bn# of                                  0## -> n-                                 z#  -> bigNatToNegInteger (bn# `shiftRBigNat` word2Int# z#)+                                 z#  -> integerFromBigNatNeg# (bn# `bigNatShiftR#` z#)  -- | Shift argument right until the result is odd. --   Precondition: argument not @0@, not checked.@@ -180,26 +183,19 @@  -- | Merges two ordered lists into an ordered list. Checks for neither its -- precondition or postcondition.-mergeBy :: (a -> a -> Ordering) -> [a] -> [a] -> [a]+mergeBy :: (a -> a -> Ordering) -> Infinite a -> Infinite a -> Infinite a mergeBy cmp = loop   where-    loop [] ys  = ys-    loop xs []  = xs-    loop (x:xs) (y:ys)+    loop ( x:< xs) (y :< ys)       = case cmp x y of-         GT -> y : loop (x:xs) ys-         _  -> x : loop xs (y:ys)+         GT -> y :< loop (x :< xs) ys+         _  -> x :< loop xs (y :< ys)  -- | Work around https://ghc.haskell.org/trac/ghc/ticket/14085 recipMod :: Integer -> Integer -> Maybe Integer-recipMod x m = case recipModInteger (x `P.mod` m) m of-  0 -> Nothing-  y -> Just y--bigNatToNatural :: BigNat -> Natural-bigNatToNatural bn-  | isTrue# (sizeofBigNat# bn ==# 1#) = NatS# (bigNatToWord bn)-  | otherwise = NatJ# bn+recipMod x m = case integerRecipMod# (x `P.mod` m) (fromInteger m) of+  (# | _ #) -> Nothing+  (# y | #) -> Just (toInteger y)  ------------------------------------------------------------------------------- -- Helpers for mapping to rough numbers and back.

@@ -6,6 +6,7 @@ -- -- Dirichlet beta-function. +{-# LANGUAGE PostfixOperators    #-} {-# LANGUAGE ScopedTypeVariables #-}  module Math.NumberTheory.Zeta.Dirichlet@@ -15,12 +16,14 @@   ) where  import Data.ExactPi-import Data.List                      (zipWith4)+import Data.List.Infinite (Infinite(..), (....))+import Data.List.NonEmpty (NonEmpty(..))+import qualified Data.List.Infinite as Inf import Data.Ratio                     ((%))  import Math.NumberTheory.Recurrences  (euler, factorial) import Math.NumberTheory.Zeta.Hurwitz (zetaHurwitz)-import Math.NumberTheory.Zeta.Utils   (intertwine, skipOdds)+import Math.NumberTheory.Zeta.Utils   (skipOdds)  -- | Infinite sequence of exact values of Dirichlet beta-function at odd arguments, starting with @β(1)@. --@@ -30,37 +33,37 @@ -- >>> import Data.Number.Fixed -- >>> approximateValue (betasOdd !! 25) :: Fixed Prec50 -- 0.99999999999999999999999960726927497384196726751694-betasOdd :: [ExactPi]-betasOdd = zipWith Exact [1, 3 ..] $ zipWith4+betasOdd :: Infinite ExactPi+betasOdd = Inf.zipWith Exact ((1, 3)....) $ Inf.zipWith4                                      (\sgn denom eul twos -> sgn * (eul % (twos * denom)))-                                     (cycle [1, -1])+                                     (Inf.cycle (1 :| [-1]))                                      (skipOdds factorial)                                      (skipOdds euler)-                                     (iterate (4 *) 4)+                                     (Inf.iterate (4 *) 4)  -- | Infinite sequence of approximate values of the Dirichlet @β@ function at -- positive even integer arguments, starting with @β(0)@.-betasEven :: forall a. (Floating a, Ord a) => a -> [a]-betasEven eps = (1 / 2) : hurwitz+betasEven :: forall a. (Floating a, Ord a) => a -> Infinite a+betasEven eps = (1 / 2) :< hurwitz   where-    hurwitz :: [a]+    hurwitz :: Infinite a     hurwitz =-        zipWith3 (\quarter threeQuarters four ->+        Inf.zipWith3 (\quarter threeQuarters four ->             (quarter - threeQuarters) / four)-        (tail . skipOdds $ zetaHurwitz eps 0.25)-        (tail . skipOdds $ zetaHurwitz eps 0.75)-        (iterate (16 *) 16)+        (Inf.tail . skipOdds $ zetaHurwitz eps 0.25)+        (Inf.tail . skipOdds $ zetaHurwitz eps 0.75)+        (Inf.iterate (16 *) 16)  -- | Infinite sequence of approximate (up to given precision) -- values of Dirichlet beta-function at integer arguments, starting with @β(0)@. -- -- >>> take 5 (betas 1e-14) :: [Double] -- [0.5,0.7853981633974483,0.9159655941772189,0.9689461462593694,0.9889445517411051]-betas :: (Floating a, Ord a) => a -> [a]-betas eps = e : o : scanl1 f (intertwine es os)+betas :: (Floating a, Ord a) => a -> Infinite a+betas eps = e :< o :< Inf.scanl1 f (Inf.interleave es os)   where-    e : es = betasEven eps-    o : os = map (getRationalLimit (\a b -> abs (a - b) < eps) . rationalApproximations) betasOdd+    e :< es = betasEven eps+    o :< os = Inf.map (getRationalLimit (\a b -> abs (a - b) < eps) . rationalApproximations) betasOdd      -- Cap-and-floor to improve numerical stability:     -- 1 > beta(n + 1) - 1 > (beta(n) - 1) / 2

@@ -6,12 +6,16 @@ -- -- Hurwitz zeta function. +{-# LANGUAGE PostfixOperators    #-} {-# LANGUAGE ScopedTypeVariables #-}  module Math.NumberTheory.Zeta.Hurwitz   ( zetaHurwitz   ) where +import Data.List.Infinite (Infinite(..), (....))+import qualified Data.List.Infinite as Inf+ import Math.NumberTheory.Recurrences (bernoulli, factorial) import Math.NumberTheory.Zeta.Utils  (skipEvens, skipOdds) @@ -33,8 +37,8 @@ -- -- >>> zetaHurwitz 1e-15 0.25 !! 5 -- 1024.3489745265808-zetaHurwitz :: forall a . (Floating a, Ord a) => a -> a -> [a]-zetaHurwitz eps a = zipWith3 (\s i t -> s + i + t) ss is ts+zetaHurwitz :: forall a . (Floating a, Ord a) => a -> a -> Infinite a+zetaHurwitz eps a = Inf.zipWith3 (\s i t -> s + i + t) ss is ts   where     -- When given @1e-14@ as the @eps@ argument, this'll be     -- @div (33 * (length . takeWhile (>= 1) . iterate (/ 10) . recip) 1e-14) 10 == div (33 * 14) 10@@@ -53,70 +57,70 @@     aPlusN = a + fromInteger digitsOfPrecision      -- @[(a + n)^s | s <- [0, 1, 2 ..]]@-    powsOfAPlusN :: [a]-    powsOfAPlusN = iterate (aPlusN *) 1+    powsOfAPlusN :: Infinite a+    powsOfAPlusN = Inf.iterate (aPlusN *) 1      -- [                   [      1      ] |                   ]     -- | \sum_{k=0}^\(n-1) | ----------- | | s <- [0, 1, 2 ..] |     -- [                   [ (a + k) ^ s ] |                   ]     -- @S@ value in 4.8.5 formula.-    ss :: [a]+    ss :: Infinite a     ss = let numbers = map ((a +) . fromInteger) [0..digitsOfPrecision-1]-             denoms  = replicate (fromInteger digitsOfPrecision) 1 :-                       iterate (zipWith (*) numbers) numbers-         in map (sum . map recip) denoms+             denoms  = replicate (fromInteger digitsOfPrecision) 1 :<+                       Inf.iterate (zipWith (*) numbers) numbers+         in Inf.map (sum . map recip) denoms      -- [ (a + n) ^ (1 - s)            a + n         |                   ]     -- | ----------------- = ---------------------- | s <- [0, 1, 2 ..] |     -- [       s - 1          (a + n) ^ s * (s - 1) |                   ]     -- @I@ value in 4.8.5 formula.-    is :: [a]-    is = let denoms = zipWith+    is :: Infinite a+    is = let denoms = Inf.zipWith                       (\powOfA int -> powOfA * fromInteger int)                       powsOfAPlusN-                      [-1, 0..]-         in map (aPlusN /) denoms+                      ((-1, 0)....)+         in Inf.map (aPlusN /) denoms      -- [      1      |             ]     -- [ ----------- | s <- [0 ..] ]     -- [ (a + n) ^ s |             ]-    constants2 :: [a]-    constants2 = map recip powsOfAPlusN+    constants2 :: Infinite a+    constants2 = Inf.map recip powsOfAPlusN      -- [ [(s)_(2*k - 1) | k <- [1 ..]], s <- [0 ..]], i.e. odd indices of     -- infinite rising factorial sequences, each sequence starting at a     -- positive integer.-    pochhammers :: [[Integer]]+    pochhammers :: Infinite (Infinite Integer)     pochhammers = let -- [ [(s)_k | k <- [1 ..]], s <- [1 ..]]-                      pochhs :: [[Integer]]-                      pochhs = iterate (\(x : xs) -> map (`div` x) xs) (tail factorial)+                      pochhs :: Infinite (Infinite Integer)+                      pochhs = Inf.iterate (\(x :< xs) -> Inf.map (`div` x) xs) (Inf.tail factorial)                   in -- When @s@ is @0@, the infinite sequence of rising                      -- factorials starting at @s@ is @[0,0,0,0..]@.-                     repeat 0 : map skipOdds pochhs+                     Inf.repeat 0 :< Inf.map skipOdds pochhs      -- [            B_2k           |             ]     -- | ------------------------- | k <- [1 ..] |     -- [ (2k)! (a + n) ^ (2*k - 1) |             ]     second :: [a]     second =-        take (fromInteger digitsOfPrecision) $-        zipWith3+        Inf.take (fromInteger digitsOfPrecision) $+        Inf.zipWith3         (\bern evenFac denom -> fromRational bern / (denom * fromInteger evenFac))-        (tail $ skipOdds bernoulli)-        (tail $ skipOdds factorial)+        (Inf.tail $ skipOdds bernoulli)+        (Inf.tail $ skipOdds factorial)         -- Recall that @powsOfAPlusN = [(a + n) ^ s | s <- [0 ..]]@, so this         -- is @[(a + n) ^ (2 * s - 1) | s <- [1 ..]]@         (skipEvens powsOfAPlusN) -    fracs :: [a]-    fracs = map-            (sum . zipWith (\s p -> s * fromInteger p) second)+    fracs :: Infinite a+    fracs = Inf.map+            (sum . zipWith (\s p -> s * fromInteger p) second . Inf.toList)             pochhammers      -- Infinite list of @T@ values in 4.8.5 formula, for every @s@ in     -- @[0, 1, 2 ..]@.-    ts :: [a]-    ts = zipWith+    ts :: Infinite a+    ts = Inf.zipWith          (\constant2 frac -> constant2 * (0.5 + frac))          constants2          fracs

@@ -6,6 +6,7 @@ -- -- Riemann zeta-function. +{-# LANGUAGE PostfixOperators    #-} {-# LANGUAGE ScopedTypeVariables #-}  module Math.NumberTheory.Zeta.Riemann@@ -15,11 +16,13 @@   ) where  import Data.ExactPi+import Data.List.Infinite (Infinite(..), (...), (....))+import qualified Data.List.Infinite as Inf import Data.Ratio                     ((%))  import Math.NumberTheory.Recurrences  (bernoulli) import Math.NumberTheory.Zeta.Hurwitz (zetaHurwitz)-import Math.NumberTheory.Zeta.Utils   (intertwine, skipEvens, skipOdds)+import Math.NumberTheory.Zeta.Utils   (skipEvens, skipOdds)  -- | Infinite sequence of exact values of Riemann zeta-function at even arguments, starting with @ζ(0)@. -- Note that due to numerical errors conversion to 'Double' may return values below 1:@@ -33,15 +36,16 @@ -- >>> approximateValue (zetasEven !! 25) :: Fixed Prec50 -- 1.00000000000000088817842111574532859293035196051773 ---zetasEven :: [ExactPi]-zetasEven = zipWith Exact [0, 2 ..] $ zipWith (*) (skipOdds bernoulli) cs+zetasEven :: Infinite ExactPi+zetasEven = Inf.zipWith Exact ((0, 2)....) $ Inf.zipWith (*) (skipOdds bernoulli) cs   where-    cs = (- 1 % 2) : zipWith (\i f -> i * (-4) / fromInteger (2 * f * (2 * f - 1))) cs [1..]+    cs :: Infinite Rational+    cs = (- 1 % 2) :< Inf.zipWith (\i f -> i * (-4) / fromInteger (2 * f * (2 * f - 1))) cs (1...)  -- | Infinite sequence of approximate values of Riemann zeta-function -- at odd arguments, starting with @ζ(1)@.-zetasOdd :: forall a. (Floating a, Ord a) => a -> [a]-zetasOdd eps = (1 / 0) : tail (skipEvens $ zetaHurwitz eps 1)+zetasOdd :: forall a. (Floating a, Ord a) => a -> Infinite a+zetasOdd eps = (1 / 0) :< Inf.tail (skipEvens $ zetaHurwitz eps 1)  -- | Infinite sequence of approximate (up to given precision) -- values of Riemann zeta-function at integer arguments, starting with @ζ(0)@.@@ -52,11 +56,11 @@ -- Beware to force evaluation of @zetas !! 1@ if the type @a@ does not support infinite values -- (for instance, 'Data.Number.Fixed.Fixed'). ---zetas :: (Floating a, Ord a) => a -> [a]-zetas eps = e : o : scanl1 f (intertwine es os)+zetas :: (Floating a, Ord a) => a -> Infinite a+zetas eps = e :< o :< Inf.scanl1 f (Inf.interleave es os)   where-    e : es = map (getRationalLimit (\a b -> abs (a - b) < eps) . rationalApproximations) zetasEven-    o : os = zetasOdd eps+    e :< es = Inf.map (getRationalLimit (\a b -> abs (a - b) < eps) . rationalApproximations) zetasEven+    o :< os = zetasOdd eps      -- Cap-and-floor to improve numerical stability:     -- 0 < zeta(n + 1) - 1 < (zeta(n) - 1) / 2

@@ -7,33 +7,25 @@ -- Shared utilities used by functions from @Math.NumberTheory.Zeta@.  module Math.NumberTheory.Zeta.Utils-  ( intertwine-  , skipEvens+  ( skipEvens   , skipOdds   ) where --- | Joins two lists element-by-element together into one, starting with the--- first one provided as argument.------ >>> take 10 (intertwine [0, 2 ..] [1, 3 ..])--- [0,1,2,3,4,5,6,7,8,9]-intertwine :: [a] -> [a] -> [a]-intertwine [] ys = ys-intertwine (x : xs) ys = x : intertwine ys xs+import Data.List.Infinite (Infinite(..))+import qualified Data.List.Infinite as Inf  -- | Skips every odd-indexed element from an infinite list. -- Do NOT use with finite lists. -- -- >>> take 10 (skipOdds [0, 1 ..]) -- [0,2,4,6,8,10,12,14,16,18]-skipOdds :: [a] -> [a]-skipOdds (x : _ : xs) = x : skipOdds xs-skipOdds xs = xs+skipOdds :: Infinite a -> Infinite a+skipOdds (x :< _ :< xs) = x :< skipOdds xs  -- | Skips every even-indexed element from an infinite list. -- Do NOT use with finite lists. -- -- >>> take 10 (skipEvens [0, 1 ..]) -- [1,3,5,7,9,11,13,15,17,19]-skipEvens :: [a] -> [a]-skipEvens = skipOdds . tail+skipEvens :: Infinite a -> Infinite a+skipEvens = skipOdds . Inf.tail

@@ -1,5 +1,5 @@ name:          arithmoi-version:       0.12.1.0+version:       0.13.0.0 cabal-version: 2.0 build-type:    Simple license:       MIT@@ -17,8 +17,8 @@   powers (integer roots and tests, modular exponentiation). category:      Math, Algorithms, Number Theory author:        Andrew Lelechenko, Daniel Fischer-tested-with:   GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.4 GHC ==8.10.7 GHC ==9.0.2 GHC ==9.2.5 GHC ==9.4.4 GHC ==9.6.1-extra-source-files:+tested-with:   GHC ==9.0.2 GHC ==9.2.8 GHC ==9.4.5 GHC ==9.6.2+extra-doc-files:   changelog.md  source-repository head@@ -27,21 +27,22 @@  library   build-depends:-    base >=4.10 && <5,+    base >=4.15 && <5,     array >=0.5 && <0.6,-    containers >=0.5.8 && <0.7,-    chimera >=0.3,-    constraints,-    deepseq,-    exact-pi >=0.5,-    integer-gmp <1.2,-    integer-logarithms >=1.0,-    integer-roots >=1.0,-    mod,+    containers >=0.5.11 && <0.7,+    chimera >=0.3 && <0.4,+    constraints <0.14,+    deepseq <1.6,+    exact-pi >=0.5 && <0.6,+    ghc-bignum <1.4,+    infinite-list <0.2,+    integer-logarithms >=1.0 && <1.1,+    integer-roots >=1.0 && <1.1,+    mod <0.3,     random >=1.0 && <1.3,     transformers >=0.4 && <0.7,-    semirings >=0.5.2,-    vector >=0.12+    semirings >=0.5.2 && <0.7,+    vector >=0.12 && <0.14   exposed-modules:     Math.NumberTheory.ArithmeticFunctions     Math.NumberTheory.ArithmeticFunctions.Inverse@@ -62,7 +63,6 @@     Math.NumberTheory.Moduli.Singleton     Math.NumberTheory.Moduli.Sqrt     Math.NumberTheory.MoebiusInversion-    Math.NumberTheory.Powers.Modular     Math.NumberTheory.Prefactored     Math.NumberTheory.Primes     Math.NumberTheory.Primes.Counting@@ -110,18 +110,18 @@     arithmoi,     containers,     exact-pi >=0.4.1.1,-    integer-gmp <1.2,+    infinite-list,     integer-roots >=1.0,     mod,-    QuickCheck >=2.10,-    quickcheck-classes >=0.6.3,+    QuickCheck >=2.10 && <2.15,+    quickcheck-classes >=0.6.3 && <0.7,     random >=1.0 && <1.3,     semirings >=0.2,     smallcheck >=1.2 && <1.3,-    tasty >=0.10,+    tasty >=0.10 && <1.5,     tasty-hunit >=0.9 && <0.11,     tasty-quickcheck >=0.9 && <0.11,-    tasty-rerun >=1.1.17,+    tasty-rerun >=1.1.17 && <1.2,     tasty-smallcheck >=0.8 && <0.9,     transformers >=0.5,     vector@@ -180,11 +180,12 @@     constraints,     containers,     deepseq,+    infinite-list,     integer-logarithms,     mod,     random,     semirings,-    tasty-bench,+    tasty-bench < 0.4,     vector   other-modules:     Math.NumberTheory.ArithmeticFunctionsBench

@@ -7,27 +7,29 @@ import Test.Tasty.Bench  import Data.Euclidean (GcdDomain)+import Data.List.Infinite (Infinite(..))+import qualified Data.List.Infinite as Inf import Math.NumberTheory.Recurrences -benchTriangle :: String -> (forall a. (GcdDomain a, Integral a) => [[a]]) -> Int -> Benchmark+benchTriangle :: String -> (forall a. (GcdDomain a, Integral a) => Infinite [a]) -> Word -> Benchmark benchTriangle name triangle n = bgroup name-  [ benchAt (10 * n)  (1 * n)-  , benchAt (10 * n)  (2 * n)-  , benchAt (10 * n)  (5 * n)-  , benchAt (10 * n)  (9 * n)+  [ benchAt (10 * n)  (1 * fromIntegral n)+  , benchAt (10 * n)  (2 * fromIntegral n)+  , benchAt (10 * n)  (5 * fromIntegral n)+  , benchAt (10 * n)  (9 * fromIntegral n)   ]   where     benchAt i j = bench ("!! " ++ show i ++ " !! " ++ show j)-                $ nf (\(x, y) -> triangle !! x !! y :: Integer) (i, j)+                $ nf (\(x, y) -> triangle Inf.!! x !! y :: Integer) (i, j) -benchPartition :: Int -> Benchmark+benchPartition :: Word -> Benchmark benchPartition n = bgroup "partition"   [ benchAt n   , benchAt (n * 10)   , benchAt (n * 100)   ]   where-    benchAt m = bench ("!!" ++ show m) $  nf (\k -> partition !! k :: Integer) m+    benchAt m = bench ("!!" ++ show m) $  nf (\k -> partition Inf.!! k :: Integer) m  benchSuite :: Benchmark benchSuite = bgroup "Recurrences"

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP        #-} {-# LANGUAGE LambdaCase #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-}@@ -8,9 +7,6 @@   ) where  import Test.Tasty.Bench-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import qualified Data.Vector as V import qualified Data.Vector.Unboxed as U 

@@ -6,10 +6,11 @@  import Test.Tasty.Bench +import qualified Data.List.Infinite as Inf import Math.NumberTheory.Zeta  benchSuite :: Benchmark benchSuite = bgroup "Zeta"-  [ bench "riemann zeta"   $ nf (sum . take 20 . zetas) (1e-15 :: Double)-  , bench "dirichlet beta" $ nf (sum . take 20 . betas) (1e-15 :: Double)+  [ bench "riemann zeta"   $ nf (sum . Inf.take 20 . zetas) (1e-15 :: Double)+  , bench "dirichlet beta" $ nf (sum . Inf.take 20 . betas) (1e-15 :: Double)   ]

@@ -1,5 +1,16 @@ # Changelog +## 0.13.0.0++### Changed++* Migrate functions under `Math.NumberTheory.Recurrences` and `Math.NumberTheory.Zeta`, which operate on infinite lists, to use `Infinite` from `infinite-list` package.+* Migrate functions under `Math.NumberTheory.Quadratic` to return an `Infinite` list of quadratic primes.++### Removed++* Remove deprecated `Math.NumberTheory.Powers.Modular`.+ ## 0.12.1.0  ### Fixed

@@ -27,6 +27,7 @@  import Data.Bits (Bits) import Data.Euclidean+import qualified Data.List.Infinite as Inf import Data.Semiring (Semiring) import qualified Data.Set as S import Numeric.Natural (Natural)@@ -89,7 +90,7 @@   ]  totientSpecialCases1 :: [Assertion]-totientSpecialCases1 = zipWith mkAssert (tail factorial) totientCountFactorial+totientSpecialCases1 = zipWith mkAssert (drop 1 (Inf.toList factorial)) totientCountFactorial   where     mkAssert n m = assertEqual "should be equal" m (totientCount n) @@ -120,7 +121,7 @@   ]  totientSpecialCases2 :: [Assertion]-totientSpecialCases2 = zipWith mkAssert (tail factorial) totientMinFactorial+totientSpecialCases2 = zipWith mkAssert (drop 1 (Inf.toList factorial)) totientMinFactorial   where     mkAssert n m = assertEqual "should be equal" m (totientMin n) @@ -151,7 +152,7 @@   ]  totientSpecialCases3 :: [Assertion]-totientSpecialCases3 = zipWith mkAssert (tail factorial) totientMaxFactorial+totientSpecialCases3 = zipWith mkAssert (drop 1 (Inf.toList factorial)) totientMaxFactorial   where     mkAssert n m = assertEqual "should be equal" m (totientMax n) @@ -250,7 +251,7 @@   ]  sigmaSpecialCases1 :: [Assertion]-sigmaSpecialCases1 = zipWith mkAssert (tail factorial) sigmaCountFactorial+sigmaSpecialCases1 = zipWith mkAssert (drop 1 (Inf.toList factorial)) sigmaCountFactorial   where     mkAssert n m = assertEqual "should be equal" m (sigmaCount n) @@ -279,7 +280,7 @@   ]  sigmaSpecialCases2 :: [Assertion]-sigmaSpecialCases2 = zipWith mkAssert (drop 3 factorial) sigmaMinFactorial+sigmaSpecialCases2 = zipWith mkAssert (drop 3 (Inf.toList factorial)) sigmaMinFactorial   where     mkAssert n m = assertEqual "should be equal" m (sigmaMin n) @@ -308,7 +309,7 @@   ]  sigmaSpecialCases3 :: [Assertion]-sigmaSpecialCases3 = zipWith mkAssert (drop 3 factorial) sigmaMaxFactorial+sigmaSpecialCases3 = zipWith mkAssert (drop 3 (Inf.toList factorial)) sigmaMaxFactorial   where     mkAssert n m = assertEqual "should be equal" m (sigmaMax n) 

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.ArithmeticFunctions.Mertens -- -{-# LANGUAGE CPP        #-} {-# LANGUAGE LambdaCase #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-}@@ -18,9 +17,6 @@  import Test.Tasty -#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import qualified Data.Vector as V import qualified Data.Vector.Unboxed as U 

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.ArithmeticFunctions.SieveBlock -- -{-# LANGUAGE CPP        #-} {-# LANGUAGE LambdaCase #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-}@@ -19,9 +18,6 @@ import Test.Tasty import Test.Tasty.HUnit -#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import qualified Data.Vector as V  import Math.NumberTheory.ArithmeticFunctions

@@ -7,8 +7,6 @@ -- Tests for Math.NumberTheory.ArithmeticFunctions -- -{-# LANGUAGE CPP       #-}- {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.ArithmeticFunctionsTests@@ -19,6 +17,7 @@ import Test.Tasty.HUnit  import Data.List (sort)+import qualified Data.List.Infinite as Inf import qualified Data.Set as S import qualified Data.IntSet as IS @@ -292,7 +291,7 @@ nFreedomProperty3 (Power n) (Positive m) =     let n' | n == maxBound = n            | otherwise     = n + 1-        zet = 1 / zetas 1e-14 !! fromIntegral n' :: Double+        zet = 1 / zetas 1e-14 Inf.!! n' :: Double         m' = 100 * m         nfree = fromIntegral m' /                 fromIntegral (head (drop (m' - 1) $ nFrees n' :: [Integer]))

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.Curves -- -{-# LANGUAGE CPP                        #-} {-# LANGUAGE LambdaCase                 #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} 

@@ -1,7 +1,5 @@ -- Tests for Math.NumberTheory.Diophantine -{-# LANGUAGE CPP       #-}- {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.DiophantineTests

@@ -15,6 +15,7 @@  import Prelude hiding (gcd, rem, quot, quotRem) import Data.Euclidean+import qualified Data.List.Infinite as Inf import Data.Maybe (fromJust, isJust) import Data.Proxy import Test.Tasty.QuickCheck as QC hiding (Positive(..), NonNegative(..))@@ -99,7 +100,7 @@ -- | Checks that the numbers produced by @primes@ are actually Eisenstein -- primes. primesProperty1 :: Positive Int -> Bool-primesProperty1 (Positive index) = all (isJust . isPrime . (unPrime :: Prime E.EisensteinInteger -> E.EisensteinInteger)) $ take index E.primes+primesProperty1 (Positive index) = all (isJust . isPrime . (unPrime :: Prime E.EisensteinInteger -> E.EisensteinInteger)) $ Inf.take index E.primes  -- | Checks that the infinite list of Eisenstein primes @primes@ is ordered -- by the numbers' norm.@@ -107,13 +108,13 @@ primesProperty2 (Positive index) =     let isOrdered :: [Prime E.EisensteinInteger] -> Bool         isOrdered xs = all (\(x, y) -> E.norm (unPrime x) <= E.norm (unPrime y)) . zip xs $ tail xs-    in isOrdered $ take index E.primes+    in isOrdered $ Inf.take index E.primes  -- | Checks that the numbers produced by @primes@ are all in the first -- sextant. primesProperty3 :: Positive Int -> Bool primesProperty3 (Positive index) =-    all (\e -> abs (unPrime e) == (unPrime e :: E.EisensteinInteger)) $ take index E.primes+    all (\e -> abs (unPrime e) == (unPrime e :: E.EisensteinInteger)) $ Inf.take index E.primes  -- | An Eisenstein integer is either zero or associated (i.e. equal up to -- multiplication by a unit) to the product of its factors raised to their

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.Euclidean.Coprimes -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications    #-} 

@@ -18,6 +18,7 @@ import Control.Monad (zipWithM_) import Data.Euclidean import Data.List (groupBy, sort)+import qualified Data.List.Infinite as Inf import Data.Maybe (fromJust, mapMaybe) import Data.Proxy import Test.Tasty.QuickCheck as QC hiding (Positive(..), NonNegative(..))@@ -99,26 +100,26 @@ primesSpecialCase1 :: Assertion primesSpecialCase1 = assertEqual "primes"   (f $ mapMaybe isPrime [1+ι,2+ι,1+2*ι,3,3+2*ι,2+3*ι,4+ι,1+4*ι,5+2*ι,2+5*ι,6+ι,1+6*ι,5+4*ι,4+5*ι,7,7+2*ι,2+7*ι,6+5*ι,5+6*ι,8+3*ι,3+8*ι,8+5*ι,5+8*ι,9+4*ι,4+9*ι,10+ι,1+10*ι,10+3*ι,3+10*ι,8+7*ι,7+8*ι,11,11+4*ι,4+11*ι,10+7*ι,7+10*ι,11+6*ι,6+11*ι,13+2*ι,2+13*ι,10+9*ι,9+10*ι,12+7*ι,7+12*ι,14+ι,1+14*ι,15+2*ι,2+15*ι,13+8*ι,8+13*ι,15+4*ι])-  (f $ take 51 primes)+  (f $ Inf.take 51 primes)   where     f :: [Prime GaussianInteger] -> [[Prime GaussianInteger]]     f = map sort . groupBy (\g1 g2 -> norm (unPrime g1) == norm (unPrime g2))  -- | The list of primes should include only primes.-primesGeneratesPrimesProperty :: NonNegative Int -> Bool-primesGeneratesPrimesProperty (NonNegative i) = case isPrime (unPrime (primes !! i) :: GaussianInteger) of+primesGeneratesPrimesProperty :: Word -> Bool+primesGeneratesPrimesProperty i = case isPrime (unPrime (primes Inf.!! i) :: GaussianInteger) of   Nothing -> False   Just{}  -> True  -- | Check that primes generates the primes in order. orderingPrimes :: Assertion orderingPrimes = assertBool "primes are ordered" (and $ zipWith (<=) xs (tail xs))-  where xs = map (norm . unPrime) $ take 1000 primes+  where xs = map (norm . unPrime) $ Inf.take 1000 primes  numberOfPrimes :: Assertion numberOfPrimes = assertEqual "counting primes: OEIS A091100"   [16,100,668,4928,38404,313752]-  [4 * length (takeWhile ((<= 10^n) . norm . unPrime) primes) | n <- [1..6]]+  [4 * length (Inf.takeWhile ((<= 10^n) . norm . unPrime) primes) | n <- [1..6]]  -- | signum and abs should satisfy: z == signum z * abs z signumAbsProperty :: GaussianInteger -> Bool

@@ -7,8 +7,6 @@ -- Test for Math.NumberTheory.Moduli.Cbrt -- -{-# LANGUAGE CPP #-}- module Math.NumberTheory.Moduli.CbrtTests   ( testSuite   ) where@@ -18,9 +16,6 @@ import Math.NumberTheory.Primes import qualified Data.Euclidean as A import Data.List (genericReplicate)-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif import Test.Tasty import Math.NumberTheory.TestUtils 

@@ -7,8 +7,6 @@ -- Tests for Math.NumberTheory.Moduli.Jacobi -- -{-# LANGUAGE CPP                        #-}- {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.Moduli.JacobiTests@@ -18,9 +16,6 @@ import Test.Tasty  import Data.Bits-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif  import Math.NumberTheory.Moduli.Sqrt import Math.NumberTheory.TestUtils

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.Moduli.PrimitiveRoot -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE ScopedTypeVariables #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-}

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.Moduli.Sqrt -- -{-# LANGUAGE CPP             #-} {-# LANGUAGE ViewPatterns    #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-}

@@ -7,8 +7,6 @@ -- Tests for Math.NumberTheory.Prefactored -- -{-# LANGUAGE CPP #-}- {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.PrefactoredTests

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.Primes.Sieve -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE FlexibleContexts    #-} {-# LANGUAGE ScopedTypeVariables #-} 

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.Primes -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications    #-} @@ -22,9 +21,6 @@ import Data.Bits import Data.Int import Data.Proxy-#if __GLASGOW_HASKELL__ < 803-import Data.Semigroup-#endif  import Math.NumberTheory.Primes import qualified Math.NumberTheory.Primes.IntSet as PS

@@ -20,205 +20,207 @@  import Control.Arrow import Data.List (sort)+import qualified Data.List.Infinite as Inf import Data.Ratio  import Math.NumberTheory.Primes import Math.NumberTheory.Recurrences.Bilinear import Math.NumberTheory.TestUtils -binomialProperty1 :: NonNegative Int -> Bool-binomialProperty1 (NonNegative i) = length (binomial @Integer !! i) == i + 1+binomialProperty1 :: Word -> Bool+binomialProperty1 i = length (binomial @Integer Inf.!! i) == fromIntegral i + 1 -binomialProperty2 :: NonNegative Int -> Bool-binomialProperty2 (NonNegative i) = head (binomial @Integer !! i) == 1+binomialProperty2 :: Word -> Bool+binomialProperty2 i = head (binomial @Integer Inf.!! i) == 1 -binomialProperty3 :: NonNegative Int -> Bool-binomialProperty3 (NonNegative i) = binomial @Integer !! i !! i == 1+binomialProperty3 :: Word -> Bool+binomialProperty3 i = binomial @Integer Inf.!! i !! fromIntegral i == 1  binomialProperty4 :: Positive Int -> Positive Int -> Bool binomialProperty4 (Positive i) (Positive j)   =  j >= i-  || binomial @Integer !! i !! j-  == binomial !! (i - 1) !! (j - 1)-  +  binomial !! (i - 1) !! j+  || binomial @Integer Inf.!! fromIntegral i !! j+  == binomial Inf.!! (fromIntegral i - 1) !! (j - 1)+  +  binomial Inf.!! (fromIntegral i - 1) !! j  binomialProperty5 :: Word -> Word -> Bool binomialProperty5 n m' = n > 100000 ||-  sort (map (first unPrime) (factorise (binomial !! fromIntegral n !! fromIntegral m))) ==+  sort (map (first unPrime) (factorise (binomial Inf.!! n !! fromIntegral m))) ==     sort (map (first (toInteger . unPrime)) (binomialFactors n m))   where     m = m' `mod` (n + 1)  binomialProperty6 :: Word -> Word -> Bool binomialProperty6 n m' = n > 100000 ||-  binomial !! fromIntegral n !! fromIntegral m ==+  binomial Inf.!! n !! fromIntegral m ==     product (map (\(p, k) -> toInteger (unPrime p) ^ k) (binomialFactors n m))   where     m = m' `mod` (n + 1) -binomialRotatedProperty2 :: NonNegative Int -> Bool-binomialRotatedProperty2 (NonNegative i) = head (binomialRotated @Integer !! i) == 1+binomialRotatedProperty2 :: Word -> Bool+binomialRotatedProperty2 i = Inf.head (binomialRotated @Integer Inf.!! i) == 1 -binomialRotatedProperty3 :: NonNegative Int -> Bool-binomialRotatedProperty3 (NonNegative i) = head (binomialRotated @Integer) !! i == 1+binomialRotatedProperty3 :: Word -> Bool+binomialRotatedProperty3 i = Inf.head (binomialRotated @Integer) Inf.!! i == 1  binomialRotatedProperty4 :: Positive Int -> Positive Int -> Bool binomialRotatedProperty4 (Positive i) (Positive j)-  =  binomialRotated @Integer !! i !! j-  == binomialRotated !! i !! (j - 1)-  +  binomialRotated !! (i - 1) !! j+  =  binomialRotated @Integer Inf.!! fromIntegral i Inf.!! fromIntegral j+  == binomialRotated Inf.!! fromIntegral i Inf.!! (fromIntegral j - 1)+  +  binomialRotated Inf.!! (fromIntegral i - 1) Inf.!! fromIntegral j -binomialLineProperty1 :: NonNegative Int -> NonNegative Int -> Bool-binomialLineProperty1 (NonNegative i) (NonNegative j)+binomialLineProperty1 :: Word -> Word -> Bool+binomialLineProperty1 i j   =  j >= i-  || binomial @Integer !! i !! j == binomialLine (toInteger i) !! j+  || binomial @Integer Inf.!! i !! fromIntegral j == binomialLine (toInteger i) !! fromIntegral j -binomialLineProperty2 :: NonNegative Int -> NonNegative Int -> Bool-binomialLineProperty2 (NonNegative i) (NonNegative j)-  = binomialRotated @Integer !! i !! j == binomialLine (toInteger (i + j)) !! j+binomialLineProperty2 :: Word -> Word -> Bool+binomialLineProperty2 i j+  = binomialRotated @Integer Inf.!! i Inf.!! j == binomialLine (toInteger (i + j)) !! fromIntegral j -binomialDiagonalProperty1 :: NonNegative Int -> NonNegative Int -> Bool-binomialDiagonalProperty1 (NonNegative i) (NonNegative j)-  = binomialRotated @Integer !! i !! j == binomialDiagonal (toInteger i) !! j+binomialDiagonalProperty1 :: Word -> Word -> Bool+binomialDiagonalProperty1 i j+  = binomialRotated @Integer Inf.!! i Inf.!! j == binomialDiagonal (toInteger i) Inf.!! j -binomialDiagonalProperty2 :: NonNegative Int -> NonNegative Int -> Bool-binomialDiagonalProperty2 (NonNegative i) (NonNegative j)-  = binomial @Integer !! (i + j) !! j == binomialDiagonal (toInteger i) !! j+binomialDiagonalProperty2 :: Word -> Word -> Bool+binomialDiagonalProperty2 i j+  = binomial @Integer Inf.!! (i + j) !! fromIntegral j == binomialDiagonal (toInteger i) Inf.!! j -stirling1Property1 :: NonNegative Int -> Bool-stirling1Property1 (NonNegative i) = length (stirling1 !! i) == i + 1+stirling1Property1 :: Word -> Bool+stirling1Property1 i = length (stirling1 Inf.!! i) == fromIntegral i + 1 -stirling1Property2 :: NonNegative Int -> Bool-stirling1Property2 (NonNegative i)-  =  head (stirling1 !! i)+stirling1Property2 :: Word -> Bool+stirling1Property2 i+  =  head (stirling1 Inf.!! i)   == if i == 0 then 1 else 0 -stirling1Property3 :: NonNegative Int -> Bool-stirling1Property3 (NonNegative i) = stirling1 !! i !! i == 1+stirling1Property3 :: Word -> Bool+stirling1Property3 i = stirling1 Inf.!! i !! fromIntegral i == 1  stirling1Property4 :: Positive Int -> Positive Int -> Bool stirling1Property4 (Positive i) (Positive j)   =  j >= i-  || stirling1 !! i !! j-  == stirling1 !! (i - 1) !! (j - 1)-  +  (toInteger i - 1) * stirling1 !! (i - 1) !! j+  || stirling1 Inf.!! fromIntegral i !! j+  == stirling1 Inf.!! (fromIntegral i - 1) !! (j - 1)+  +  (toInteger i - 1) * stirling1 Inf.!! (fromIntegral i - 1) !! j -stirling2Property1 :: NonNegative Int -> Bool-stirling2Property1 (NonNegative i) = length (stirling2 !! i) == i + 1+stirling2Property1 :: Word -> Bool+stirling2Property1 i = length (stirling2 Inf.!! i) == fromIntegral i + 1 -stirling2Property2 :: NonNegative Int -> Bool-stirling2Property2 (NonNegative i)-  =  head (stirling2 !! i)+stirling2Property2 :: Word -> Bool+stirling2Property2 i+  =  head (stirling2 Inf.!! i)   == if i == 0 then 1 else 0 -stirling2Property3 :: NonNegative Int -> Bool-stirling2Property3 (NonNegative i) = stirling2 !! i !! i == 1+stirling2Property3 :: Word -> Bool+stirling2Property3 i = stirling2 Inf.!! i !! fromIntegral i == 1  stirling2Property4 :: Positive Int -> Positive Int -> Bool stirling2Property4 (Positive i) (Positive j)   =  j >= i-  || stirling2 !! i !! j-  == stirling2 !! (i - 1) !! (j - 1)-  +  toInteger j * stirling2 !! (i - 1) !! j+  || stirling2 Inf.!! fromIntegral i !! j+  == stirling2 Inf.!! (fromIntegral i - 1) !! (j - 1)+  +  toInteger j * stirling2 Inf.!! (fromIntegral i - 1) !! j -lahProperty1 :: NonNegative Int -> Bool-lahProperty1 (NonNegative i) = length (lah !! i) == i + 1+lahProperty1 :: Word -> Bool+lahProperty1 i = length (lah Inf.!! i) == fromIntegral i + 1 -lahProperty2 :: NonNegative Int -> Bool-lahProperty2 (NonNegative i)-  =  head (lah !! i)+lahProperty2 :: Word -> Bool+lahProperty2 i+  =  head (lah Inf.!! i)   == product [1 .. i+1] -lahProperty3 :: NonNegative Int -> Bool-lahProperty3 (NonNegative i) = lah !! i !! i == 1+lahProperty3 :: Word -> Bool+lahProperty3 i = lah Inf.!! i !! fromIntegral i == 1  lahProperty4 :: Positive Int -> Positive Int -> Bool lahProperty4 (Positive i) (Positive j)   =  j >= i-  || lah !! i !! j-  == sum [ stirling1 !! (i + 1) !! k * stirling2 !! k !! (j + 1) | k <- [j + 1 .. i + 1] ]+  || lah Inf.!! fromIntegral i !! j+  == sum [ stirling1 Inf.!! (fromIntegral i + 1) !! k * stirling2 Inf.!! fromIntegral k !! (j + 1) | k <- [j + 1 .. i + 1] ] -eulerian1Property1 :: NonNegative Int -> Bool-eulerian1Property1 (NonNegative i) = length (eulerian1 !! i) == i+eulerian1Property1 :: Word -> Bool+eulerian1Property1 i = length (eulerian1 Inf.!! i) == fromIntegral i  eulerian1Property2 :: Positive Int -> Bool-eulerian1Property2 (Positive i) = head (eulerian1 !! i) == 1+eulerian1Property2 (Positive i) = head (eulerian1 Inf.!! fromIntegral i) == 1  eulerian1Property3 :: Positive Int -> Bool-eulerian1Property3 (Positive i) = eulerian1 !! i !! (i - 1) == 1+eulerian1Property3 (Positive i) = eulerian1 Inf.!! fromIntegral i !! (i - 1) == 1  eulerian1Property4 :: Positive Int -> Positive Int -> Bool eulerian1Property4 (Positive i) (Positive j)   =  j >= i - 1-  || eulerian1 !! i !! j-  == toInteger (i - j) * eulerian1 !! (i - 1) !! (j - 1)-  +  (toInteger j + 1) * eulerian1 !! (i - 1) !! j+  || eulerian1 Inf.!! fromIntegral i !! j+  == toInteger (i - j) * eulerian1 Inf.!! (fromIntegral i - 1) !! (j - 1)+  +  (toInteger j + 1) * eulerian1 Inf.!! (fromIntegral i - 1) !! j -eulerian2Property1 :: NonNegative Int -> Bool-eulerian2Property1 (NonNegative i) = length (eulerian2 !! i) == i+eulerian2Property1 :: Word -> Bool+eulerian2Property1 i = length (eulerian2 Inf.!! i) == fromIntegral i  eulerian2Property2 :: Positive Int -> Bool eulerian2Property2 (Positive i)-  = head (eulerian2 !! i) == 1+  = head (eulerian2 Inf.!! fromIntegral i) == 1  eulerian2Property3 :: Positive Int -> Bool eulerian2Property3 (Positive i)-  =  eulerian2 !! i !! (i - 1)+  =  eulerian2 Inf.!! fromIntegral i !! (i - 1)   == product [1 .. toInteger i]  eulerian2Property4 :: Positive Int -> Positive Int -> Bool eulerian2Property4 (Positive i) (Positive j)   =  j >= i - 1-  || eulerian2 !! i !! j-  == toInteger (2 * i - j - 1) * eulerian2 !! (i - 1) !! (j - 1)-  +  (toInteger j + 1) * eulerian2 !! (i - 1) !! j+  || eulerian2 Inf.!! fromIntegral i !! j+  == toInteger (2 * i - j - 1) * eulerian2 Inf.!! (fromIntegral i - 1) !! (j - 1)+  +  (toInteger j + 1) * eulerian2 Inf.!! (fromIntegral i - 1) !! j  bernoulliSpecialCase1 :: Assertion-bernoulliSpecialCase1 = assertEqual "B_0 = 1" (head bernoulli) 1+bernoulliSpecialCase1 = assertEqual "B_0 = 1" (Inf.head bernoulli) 1  bernoulliSpecialCase2 :: Assertion-bernoulliSpecialCase2 = assertEqual "B_1 = -1/2" (bernoulli !! 1) (- 1 % 2)+bernoulliSpecialCase2 = assertEqual "B_1 = -1/2" (bernoulli Inf.!! 1) (- 1 % 2) -bernoulliProperty1 :: NonNegative Int -> Bool-bernoulliProperty1 (NonNegative m)-  = case signum (bernoulli !! m) of+bernoulliProperty1 :: Word -> Bool+bernoulliProperty1 m+  = case signum (bernoulli Inf.!! m) of     1  -> m == 0 || m `mod` 4 == 2     0  -> m /= 1 && odd m     -1 -> m == 1 || (m /= 0 && m `rem` 4 == 0)     _  -> False -bernoulliProperty2 :: NonNegative Int -> Bool-bernoulliProperty2 (NonNegative m)-  =  bernoulli !! m-  == (if m == 0 then 1 else 0)-  -  sum [ bernoulli !! k-         * (binomial !! m !! k % toInteger (m - k + 1))-         | k <- [0 .. m - 1]+bernoulliProperty2 :: Word -> (Rational, Rational)+bernoulliProperty2 m =+  ( bernoulli Inf.!! m+  , (if m == 0 then 1 else 0)+  -  sum [ bernoulli Inf.!! k+         * (binomial Inf.!! m !! fromIntegral k % toInteger (m - k + 1))+         | k <- map fromIntegral [(0 :: Int) .. fromIntegral m - 1]          ]+  )  -- | For every odd positive integer @n@, @E_n@ is @0@. eulerProperty1 :: Positive Int -> Bool-eulerProperty1 (Positive n) = euler !! (2 * n - 1) == 0+eulerProperty1 (Positive n) = euler Inf.!! (2 * fromIntegral n - 1) == 0  -- | Every positive even index produces a negative result.-eulerProperty2 :: NonNegative Int -> Bool-eulerProperty2 (NonNegative n) = euler !! (2 + 4 * n) < 0+eulerProperty2 :: Word -> Bool+eulerProperty2 n = euler Inf.!! (2 + 4 * n) < 0  -- | The Euler number sequence is https://oeis.org/A122045 eulerSpecialCase1 :: Assertion eulerSpecialCase1 = assertEqual "euler"-    (take 20 euler)+    (Inf.take 20 euler)     [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0, -50521, 0, 2702765, 0, -199360981, 0, 19391512145, 0, -2404879675441, 0]  -- | For any even positive integer @n@, @E_n(1)@ is @0@. eulerPAt1Property1 :: Positive Int -> Bool-eulerPAt1Property1 (Positive n) = (eulerPolyAt1 !! (2 * n)) == 0+eulerPAt1Property1 (Positive n) = (eulerPolyAt1 Inf.!! (2 * fromIntegral n)) == 0  -- | The numerators in this sequence are from https://oeis.org/A198631 while the -- denominators are from https://oeis.org/A006519. eulerPAt1SpecialCase1 :: Assertion eulerPAt1SpecialCase1 = assertEqual "eulerPolyAt1"-    (take 20 eulerPolyAt1)+    (Inf.take 20 eulerPolyAt1)     (zipWith (%) [1, 1, 0, -1, 0, 1, 0, -17, 0, 31, 0, -691, 0, 5461, 0, -929569, 0, 3202291, 0, -221930581]                  [1, 2, 1, 4, 1, 2, 1, 8, 1, 2, 1, 4, 1, 2, 1, 16, 1, 2, 1, 4]) @@ -275,7 +277,7 @@     [ testCase "B_0"                           bernoulliSpecialCase1     , testCase "B_1"                           bernoulliSpecialCase2     , testSmallAndQuick "sign"                 bernoulliProperty1-    , testSmallAndQuick "recursive definition" bernoulliProperty2+    , testEqualSmallAndQuick "recursive definition" bernoulliProperty2     ]     , testGroup "Euler numbers"     [ testCase "First 20 elements of E_n are correct"           eulerSpecialCase1

@@ -7,8 +7,6 @@ -- Tests for Math.NumberTheory.Recurrences.Linear -- -{-# LANGUAGE CPP       #-}- {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.Recurrences.LinearTests@@ -20,6 +18,7 @@  import Control.Arrow import Data.List (sort)+import qualified Data.List.Infinite as Inf  import Math.NumberTheory.Primes import Math.NumberTheory.Recurrences.Linear@@ -30,8 +29,8 @@ fibonacciProperty1 (AnySign n) = fibonacci n + fibonacci (n + 1) == fibonacci (n +2)  -- | Check that 'fibonacci' for negative indices is correctly defined.-fibonacciProperty2 :: NonNegative Int -> Bool-fibonacciProperty2 (NonNegative n) = fibonacci n == (if even n then negate else id) (fibonacci (- n))+fibonacciProperty2 :: Word -> Bool+fibonacciProperty2 n = fibonacci (fromIntegral n) == (if even n then negate else id) (fibonacci (- fromIntegral n))  -- | Check that 'fibonacciPair' is a pair of consequent 'fibonacci'. fibonacciPairProperty :: AnySign Int -> Bool@@ -51,8 +50,8 @@ lucasProperty1 (AnySign n) = lucas n + lucas (n + 1) == lucas (n +2)  -- | Check that 'lucas' for negative indices is correctly defined.-lucasProperty2 :: NonNegative Int -> Bool-lucasProperty2 (NonNegative n) = lucas n == (if odd n then negate else id) (lucas (- n))+lucasProperty2 :: Word -> Bool+lucasProperty2 n = lucas (fromIntegral n) == (if odd n then negate else id) (lucas (- fromIntegral n))  -- | Check that 'lucasPair' is a pair of consequent 'lucas'. lucasPairProperty :: AnySign Int -> Bool@@ -67,17 +66,17 @@ lucasSpecialCase1 = assertEqual "lucas" (lucas 1) 1  -- | Check that 'generalLucas' matches its definition.-generalLucasProperty1 :: AnySign Integer -> AnySign Integer -> NonNegative Int -> Bool-generalLucasProperty1 (AnySign p) (AnySign q) (NonNegative n) = un1 == un1' && vn1 == vn1' && un2 == p * un1 - q * un && vn2 == p * vn1 - q * vn+generalLucasProperty1 :: AnySign Integer -> AnySign Integer -> Word -> Bool+generalLucasProperty1 (AnySign p) (AnySign q) n = un1 == un1' && vn1 == vn1' && un2 == p * un1 - q * un && vn2 == p * vn1 - q * vn   where-    (un, un1, vn, vn1) = generalLucas p q n-    (un1', un2, vn1', vn2) = generalLucas p q (n + 1)+    (un, un1, vn, vn1) = generalLucas p q (fromIntegral n)+    (un1', un2, vn1', vn2) = generalLucas p q (fromIntegral n + 1)  -- | Check that 'generalLucas' 1 (-1) is 'fibonacciPair' plus 'lucasPair'.-generalLucasProperty2 :: NonNegative Int -> Bool-generalLucasProperty2 (NonNegative n) = (un, un1) == fibonacciPair n && (vn, vn1) == lucasPair n+generalLucasProperty2 :: Word -> Bool+generalLucasProperty2 n = (un, un1) == fibonacciPair (fromIntegral n) && (vn, vn1) == lucasPair (fromIntegral n)   where-    (un, un1, vn, vn1) = generalLucas 1 (-1) n+    (un, un1, vn, vn1) = generalLucas 1 (-1) (fromIntegral n)  -- | Check that 'generalLucas' p _ 0 is (0, 1, 2, p). generalLucasProperty3 :: AnySign Integer -> AnySign Integer -> Bool@@ -85,12 +84,12 @@  factorialProperty1 :: Word -> Bool factorialProperty1 n = n > 100000 ||-  sort (map (first unPrime) (factorise (factorial !! fromIntegral n))) ==+  sort (map (first unPrime) (factorise (factorial Inf.!! n))) ==     sort (map (first (toInteger . unPrime)) (factorialFactors n))  factorialProperty2 :: Word -> Bool factorialProperty2 n = n > 100000 ||-  factorial !! fromIntegral n ==+  factorial Inf.!! n ==     product (map (\(p, k) -> toInteger (unPrime p) ^ k) (factorialFactors n))  testSuite :: TestTree

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.Recurrences.Pentagonal -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE ViewPatterns        #-} @@ -17,6 +16,8 @@   ( testSuite   ) where +import Data.List.Infinite (Infinite(..))+import qualified Data.List.Infinite as Inf import Data.Proxy                    (Proxy (..)) import GHC.Natural                   (Natural) import GHC.TypeNats                  (SomeNat (..), someNatVal)@@ -30,7 +31,7 @@  -- | Helper to avoid writing @partition !!@ too many times. partition' :: Num a => Int -> a-partition' = (partition !!)+partition' = (partition Inf.!!) . fromIntegral  -- | Check that the @k@-th generalized pentagonal number is -- @div (3 * k² - k) 2@, where @k ∈ {0, 1, -1, 2, -2, 3, -3, 4, ...}@.@@ -48,7 +49,7 @@ -- https://oeis.org/A000041. partitionSpecialCase20 :: Assertion partitionSpecialCase20 = assertEqual "partition"-    (take 20 partition)+    (Inf.take 20 partition)     [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, 176, 231, 297, 385, 490]  -- | Copied from @Math.NumberTheory.Recurrences.Pentagonal@ to test the@@ -81,8 +82,8 @@ partitionProperty2 :: NonNegative Integer -> Positive Natural -> Bool partitionProperty2 (NonNegative m)                    n@(someNatVal . getPositive -> (SomeNat (Proxy :: Proxy n))) =-    (take m' . map getVal $ (partition :: [Mod n])) ==-    map helper (take m' partition :: [Integer])+    (Inf.take m' . Inf.map getVal $ (partition :: Infinite (Mod n))) ==+    map helper (Inf.take m' partition :: [Integer])   where     m' = fromIntegral m     n' = fromIntegral n

@@ -22,6 +22,8 @@ import Data.Coerce import Data.Euclidean import Data.List (nub)+import Data.List.Infinite (Infinite(..))+import qualified Data.List.Infinite as Inf import Numeric.Natural  import Math.NumberTheory.Primes (Prime (..))@@ -33,22 +35,22 @@ isSmoothPropertyHelper   :: (Eq a, Num a, Euclidean a)   => (a -> Integer)-  -> [a]+  -> Infinite a   -> Int   -> Int   -> Bool isSmoothPropertyHelper norm primes' i1 i2 =-    let primes = take i1 primes'+    let primes = Inf.take i1 primes'         basis  = fromList primes     in all (isSmooth basis) $ take i2 $ smoothOver' norm basis  isSmoothProperty1 :: Positive Int -> Positive Int -> Bool isSmoothProperty1 (Positive i1) (Positive i2) =-    isSmoothPropertyHelper G.norm (map unPrime G.primes) i1 i2+    isSmoothPropertyHelper G.norm (fmap unPrime G.primes) i1 i2  isSmoothProperty2 :: Positive Int -> Positive Int -> Bool isSmoothProperty2 (Positive i1) (Positive i2) =-    isSmoothPropertyHelper E.norm (map unPrime E.primes) i1 i2+    isSmoothPropertyHelper E.norm (fmap unPrime E.primes) i1 i2  smoothOverInRange :: Integral a => SmoothBasis a -> a -> a -> [a] smoothOverInRange s lo hi

@@ -8,7 +8,6 @@ --  {-# LANGUAGE ConstraintKinds            #-}-{-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}

@@ -7,7 +7,6 @@ -- Utils to test Math.NumberTheory -- -{-# LANGUAGE CPP                        #-} {-# LANGUAGE DeriveTraversable          #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}

@@ -7,7 +7,6 @@ -- Tests for Math.NumberTheory.ArithmeticFunctions -- -{-# LANGUAGE CPP                 #-} {-# LANGUAGE ScopedTypeVariables #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-}

@@ -14,6 +14,8 @@   ) where  import Data.ExactPi                (approximateValue)+import Data.List.Infinite (Infinite(..))+import qualified Data.List.Infinite as Inf  import Test.Tasty import Test.Tasty.HUnit            (Assertion, assertEqual, testCase)@@ -24,25 +26,25 @@ epsilon :: Double epsilon = 1e-14 -betas' :: [Double]+betas' :: Infinite Double betas' = betas epsilon  betasOddSpecialCase1 :: Assertion betasOddSpecialCase1   = assertEqualUpToEps "beta(1) = pi/4" epsilon-    (approximateValue $ head betasOdd)+    (approximateValue $ Inf.head betasOdd)     (pi / 4)  betasOddSpecialCase2 :: Assertion betasOddSpecialCase2   = assertEqualUpToEps "beta(3) = pi^3/32" epsilon-    (approximateValue $ betasOdd !! 1)+    (approximateValue $ betasOdd Inf.!! 1)     (pi^3 / 32)  betasOddSpecialCase3 :: Assertion betasOddSpecialCase3   = assertEqualUpToEps "beta(5) = 5*pi^5/1536" epsilon-    (approximateValue $ betasOdd !! 2)+    (approximateValue $ betasOdd Inf.!! 2)     ((5 * pi^5) / 1536)  betasOddProperty1 :: Positive Int -> Bool@@ -50,32 +52,32 @@   =  betaM < 1   || betaM < betaM1   where-    betaM  = approximateValue (betasOdd !! m)-    betaM1 = approximateValue (betasOdd !! (m + 1))+    betaM  = approximateValue (betasOdd Inf.!! fromIntegral m)+    betaM1 = approximateValue (betasOdd Inf.!! (fromIntegral m + 1)) -betasOddProperty2 :: NonNegative Int -> Bool-betasOddProperty2 (NonNegative m)+betasOddProperty2 :: Word -> Bool+betasOddProperty2 m   = abs (betaM - betaM') < epsilon   where-    betaM  = approximateValue (betasOdd !! m)-    betaM' = betas' !! ((2 * m) + 1)+    betaM  = approximateValue (betasOdd Inf.!! m)+    betaM' = betas' Inf.!! ((2 * m) + 1)  betasSpecialCase1 :: Assertion betasSpecialCase1   = assertEqual "beta(0) = 1/2"-    (head betas')+    (Inf.head betas')     (1 / 2)  betasSpecialCase2 :: Assertion betasSpecialCase2   = assertEqualUpToEps "beta(2) = 0.9159655" epsilon-    (betas' !! 2)+    (betas' Inf.!! 2)     0.9159655941772190150546035149323841107  betasSpecialCase3 :: Assertion betasSpecialCase3   = assertEqualUpToEps "beta(4) = 0.9889445" epsilon-    (betas' !! 4)+    (betas' Inf.!! 4)     0.9889445517411053361084226332283778213  betasProperty1 :: Positive Int -> Bool@@ -83,12 +85,12 @@   =  betaM <= betaM1   && betaM1 <= 1   where-    betaM  = betas' !! m-    betaM1 = betas' !! (m + 1)+    betaM  = betas' Inf.!! fromIntegral m+    betaM1 = betas' Inf.!! (fromIntegral m + 1)  betasProperty2 :: NonNegative Int -> NonNegative Int -> Bool betasProperty2 (NonNegative e1) (NonNegative e2)-  = maximum (take 35 $ drop 2 $ zipWith ((abs .) . (-)) (betas eps1) (betas eps2)) <= eps1 + eps2+  = maximum (Inf.take 35 $ Inf.drop 2 $ Inf.zipWith ((abs .) . (-)) (betas eps1) (betas eps2)) <= eps1 + eps2   where     eps1, eps2 :: Double     eps1 = max ((1.0 / 2) ^ e1) ((1.0 / 2) ^ 53)

@@ -14,6 +14,8 @@   ) where  import Data.ExactPi                (approximateValue)+import Data.List.Infinite (Infinite(..))+import qualified Data.List.Infinite as Inf  import Test.Tasty import Test.Tasty.HUnit            (Assertion, assertEqual, testCase)@@ -27,19 +29,19 @@ zetasEvenSpecialCase1 :: Assertion zetasEvenSpecialCase1   = assertEqual "zeta(0) = -1/2"-    (approximateValue $ head zetasEven)+    (approximateValue $ Inf.head zetasEven)     (-1 / 2)  zetasEvenSpecialCase2 :: Assertion zetasEvenSpecialCase2   = assertEqualUpToEps "zeta(2) = pi^2/6" epsilon-    (approximateValue $ zetasEven !! 1)+    (approximateValue $ zetasEven Inf.!! 1)     (pi * pi / 6)  zetasEvenSpecialCase3 :: Assertion zetasEvenSpecialCase3   = assertEqualUpToEps "zeta(4) = pi^4/90" epsilon-    (approximateValue $ zetasEven !! 2)+    (approximateValue $ zetasEven Inf.!! 2)     (pi ^ 4 / 90)  zetasEvenProperty1 :: Positive Int -> Bool@@ -47,35 +49,35 @@   =  zetaM < 1   || zetaM > zetaM1   where-    zetaM  = approximateValue (zetasEven !! m)-    zetaM1 = approximateValue (zetasEven !! (m + 1))+    zetaM  = approximateValue (zetasEven Inf.!! fromIntegral m)+    zetaM1 = approximateValue (zetasEven Inf.!! (fromIntegral m + 1))  zetasEvenProperty2 :: Positive Int -> Bool zetasEvenProperty2 (Positive m)   = abs (zetaM - zetaM') < epsilon   where-    zetaM  = approximateValue (zetasEven !! m)-    zetaM' = zetas' !! (2 * m)+    zetaM  = approximateValue (zetasEven Inf.!! fromIntegral m)+    zetaM' = zetas' Inf.!! (2 * fromIntegral m) -zetas' :: [Double]+zetas' :: Infinite Double zetas' = zetas epsilon  zetasSpecialCase1 :: Assertion zetasSpecialCase1   = assertEqual "zeta(1) = Infinity"-    (zetas' !! 1)+    (zetas' Inf.!! 1)     (1 / 0)  zetasSpecialCase2 :: Assertion zetasSpecialCase2   = assertEqualUpToEps "zeta(3) = 1.2020569" epsilon-    (zetas' !! 3)+    (zetas' Inf.!! 3)     1.2020569031595942853997381615114499908  zetasSpecialCase3 :: Assertion zetasSpecialCase3   = assertEqualUpToEps "zeta(5) = 1.0369277" epsilon-    (zetas' !! 5)+    (zetas' Inf.!! 5)     1.0369277551433699263313654864570341681  zetasProperty1 :: Positive Int -> Bool@@ -83,8 +85,8 @@   =  zetaM >= zetaM1   && zetaM1 >= 1   where-    zetaM  = zetas' !! m-    zetaM1 = zetas' !! (m + 1)+    zetaM  = zetas' Inf.!! fromIntegral m+    zetaM1 = zetas' Inf.!! (fromIntegral m + 1)  -- | Let z1 be an approximation of z with precision eps1, -- and z2 be an approximation of the same value with precision eps2.@@ -92,7 +94,7 @@ -- abs (z1 - z2) < eps1 + eps2. zetasProperty2 :: NonNegative Int -> NonNegative Int -> Bool zetasProperty2 (NonNegative e1) (NonNegative e2)-  = maximum (take 35 $ drop 2 $ zipWith ((abs .) . (-)) (zetas eps1) (zetas eps2)) < eps1 + eps2+  = maximum (Inf.take 35 $ Inf.drop 2 $ Inf.zipWith ((abs .) . (-)) (zetas eps1) (zetas eps2)) < eps1 + eps2   where     eps1, eps2 :: Double     eps1 = max ((1.0 / 2) ^ e1) ((1.0 / 2) ^ 53)