Migrate Math.NumberTheory.Quadratic to infinite-list

Math/NumberTheory/Diophantine.hs

@@ -6,6 +6,9 @@ module Math.NumberTheory.Diophantine
  )
 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 @@ import Math.NumberTheory.Utils.FromIntegral
 -- | 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

Math/NumberTheory/DirichletCharacters.hs

@@ -136,10 +136,10 @@ instance Eq DirichletFactor where
 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

Math/NumberTheory/Moduli/Internal.hs

@@ -95,18 +95,34 @@ discreteLogarithmPrime p a b
  | 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 :: 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
@@ -117,7 +133,7 @@ discreteLogarithmPrimePollard :: Integer -> Integer -> Integer -> Natural
 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

Math/NumberTheory/Moduli/Sqrt.hs

@@ -9,6 +9,9 @@
 --
 
 {-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE PolyKinds #-}
+{-# LANGUAGE PostfixOperators #-}
 {-# LANGUAGE TupleSections #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
@@ -29,10 +32,13 @@ module Math.NumberTheory.Moduli.Sqrt
 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,14 +145,20 @@ sqrtModP' square
  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
@@ -243,17 +255,18 @@ rem4 n = fromIntegral n .&. 3
 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...)

Math/NumberTheory/Primes/Factorisation/Montgomery.hs

@@ -331,7 +331,7 @@ enumAndMultiplyFromThenTo p from thn to = zip [from, thn .. to] progression
  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]

Math/NumberTheory/Quadratic/EisensteinIntegers.hs

@@ -10,6 +10,7 @@
 
 {-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE PostfixOperators #-}
 {-# LANGUAGE RankNTypes #-}
 {-# LANGUAGE TypeFamilies #-}
 
@@ -30,7 +31,10 @@ import Prelude hiding (quot, quotRem, gcd)
 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 @@ findPrime p = case (r, sqrtsModPrime (9 * q * q - 1) p) of
 --
 -- >>> 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]
 --

Math/NumberTheory/Quadratic/GaussianIntegers.hs

@@ -9,6 +9,7 @@
 --
 
 {-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE PostfixOperators #-}
 {-# LANGUAGE TypeFamilies #-}
 
 module Math.NumberTheory.Quadratic.GaussianIntegers (
@@ -24,7 +25,10 @@ import Prelude hiding (quot, quotRem)
 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 @@ isPrime g@(x :+ y)
 --
 -- >>> 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

Math/NumberTheory/Recurrences/Bilinear.hs

@@ -277,7 +277,7 @@ faulhaberPoly p
 -- >>> 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 => Infinite (Ratio a)
-euler' = Inf.tail $ helperForBEEP tail as
+euler' = Inf.tail $ helperForBEEP (drop 1) as
  where
  as :: Infinite (Ratio a)
  as = Inf.zipWith3
@@ -308,7 +308,7 @@ euler = Inf.map numerator euler'
 -- >>> 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 => Infinite (Ratio a)
-eulerPolyAt1 = Inf.tail $ helperForBEEP tail (Inf.map recip (Inf.iterate (2 *) 1))
+eulerPolyAt1 = Inf.tail $ helperForBEEP (drop 1) (Inf.map recip (Inf.iterate (2 *) 1))
 {-# SPECIALIZE eulerPolyAt1 :: Infinite (Ratio Int) #-}
 {-# SPECIALIZE eulerPolyAt1 :: Infinite (Rational) #-}
 

Math/NumberTheory/Utils.hs

@@ -40,6 +40,7 @@ import qualified Prelude as P
 
 import Data.Bits
 import Data.Euclidean
+import Data.List.Infinite (Infinite(..))
 import Data.Semiring (Semiring(..), isZero)
 import GHC.Base
 import GHC.Num.BigNat
@@ -182,15 +183,13 @@ splitOff# p n = go 0## n
 
 -- | 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

changelog.md

@@ -5,6 +5,7 @@
 ### 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
 

test-suite/Math/NumberTheory/EisensteinIntegersTests.hs

@@ -15,6 +15,7 @@ module Math.NumberTheory.EisensteinIntegersTests
 
 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,21 +100,21 @@ euclideanDomainProperty1 e1 e2 = E.norm (e1 * e2) == E.norm e1 * E.norm e2
 -- | 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.
 primesProperty2 :: Positive Int -> Bool
 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