-- Hoogle documentation, generated by Haddock
-- See Hoogle, http://www.haskell.org/hoogle/


-- | Efficient basic number-theoretic functions.
--   
--   A library of basic functionality needed for number-theoretic
--   calculations. The aim of this library is to provide efficient
--   implementations of the functions. Primes and related things (totients,
--   factorisation), powers (integer roots and tests, modular
--   exponentiation).
@package arithmoi
@version 0.8.0.0


-- | This module exports a class to represent Euclidean domains.
module Math.NumberTheory.Euclidean

-- | A class to represent a Euclidean domain, which is basically an
--   <a>Integral</a> without <a>toInteger</a>.
class (Eq a, Num a) => Euclidean a

-- | When restriced to a subring of the Euclidean domain <tt>a</tt>
--   isomorphic to <tt>Integer</tt>, this function should match
--   <tt>quotRem</tt> for Integers.
quotRem :: Euclidean a => a -> a -> (a, a)

-- | When restriced to a subring of the Euclidean domain <tt>a</tt>
--   isomorphic to <tt>Integer</tt>, this function should match
--   <tt>divMod</tt> for Integers.
divMod :: Euclidean a => a -> a -> (a, a)
quot :: Euclidean a => a -> a -> a
rem :: Euclidean a => a -> a -> a
div :: Euclidean a => a -> a -> a
mod :: Euclidean a => a -> a -> a

-- | <tt><a>gcd</a> x y</tt> is the greatest number that divides both
--   <tt>x</tt> and <tt>y</tt>.
gcd :: Euclidean a => a -> a -> a

-- | <tt><a>lcm</a> x y</tt> is the smallest number that both <tt>x</tt>
--   and <tt>y</tt> divide.
lcm :: Euclidean a => a -> a -> a

-- | Test whether two numbers are coprime .
coprime :: Euclidean a => a -> a -> Bool

-- | Calculate the greatest common divisor of two numbers and coefficients
--   for the linear combination.
--   
--   For signed types satisfies:
--   
--   <pre>
--   case extendedGCD a b of
--     (d, u, v) -&gt; u*a + v*b == d
--                  &amp;&amp; d == gcd a b
--   </pre>
--   
--   For unsigned and bounded types the property above holds, but since
--   <tt>u</tt> and <tt>v</tt> must also be unsigned, the result may look
--   weird. E. g., on 64-bit architecture
--   
--   <pre>
--   extendedGCD (2 :: Word) (3 :: Word) == (1, 2^64-1, 1)
--   </pre>
--   
--   For unsigned and unbounded types (like <a>Natural</a>) the result is
--   undefined.
--   
--   For signed types we also have
--   
--   <pre>
--   abs u &lt; abs b || abs b &lt;= 1
--   
--   abs v &lt; abs a || abs a &lt;= 1
--   </pre>
--   
--   (except if one of <tt>a</tt> and <tt>b</tt> is <a>minBound</a> of a
--   signed type).
extendedGCD :: Euclidean a => a -> a -> (a, a, a)

-- | Wrapper around <a>Integral</a>, which has an <tt>Eucledian</tt>
--   instance.
newtype WrappedIntegral a
WrappedIntegral :: a -> WrappedIntegral a
[unWrappedIntegral] :: WrappedIntegral a -> a
instance GHC.Enum.Enum a => GHC.Enum.Enum (Math.NumberTheory.Euclidean.WrappedIntegral a)
instance GHC.Real.Real a => GHC.Real.Real (Math.NumberTheory.Euclidean.WrappedIntegral a)
instance GHC.Real.Integral a => GHC.Real.Integral (Math.NumberTheory.Euclidean.WrappedIntegral a)
instance GHC.Num.Num a => GHC.Num.Num (Math.NumberTheory.Euclidean.WrappedIntegral a)
instance GHC.Show.Show a => GHC.Show.Show (Math.NumberTheory.Euclidean.WrappedIntegral a)
instance GHC.Classes.Ord a => GHC.Classes.Ord (Math.NumberTheory.Euclidean.WrappedIntegral a)
instance GHC.Classes.Eq a => GHC.Classes.Eq (Math.NumberTheory.Euclidean.WrappedIntegral a)
instance GHC.Real.Integral a => Math.NumberTheory.Euclidean.Euclidean (Math.NumberTheory.Euclidean.WrappedIntegral a)
instance Math.NumberTheory.Euclidean.Euclidean GHC.Types.Int
instance Math.NumberTheory.Euclidean.Euclidean GHC.Types.Word
instance Math.NumberTheory.Euclidean.Euclidean GHC.Integer.Type.Integer
instance Math.NumberTheory.Euclidean.Euclidean GHC.Natural.Natural


-- | Container for pairwise coprime numbers.
module Math.NumberTheory.Euclidean.Coprimes

-- | The input list is assumed to be a factorisation of some number into a
--   list of powers of (possibly, composite) non-zero factors. The output
--   list is a factorisation of the same number such that all factors are
--   coprime. Such transformation is crucial to continue factorisation
--   (lazily, in parallel or concurrent fashion) without having to merge
--   multiplicities of primes, which occurs more than in one composite
--   factor.
--   
--   <pre>
--   &gt;&gt;&gt; splitIntoCoprimes [(140, 1), (165, 1)]
--   Coprimes {unCoprimes = [(28,1),(33,1),(5,2)]}
--   
--   &gt;&gt;&gt; splitIntoCoprimes [(360, 1), (210, 1)]
--   Coprimes {unCoprimes = [(7,1),(5,2),(3,3),(2,4)]}
--   </pre>
splitIntoCoprimes :: (Euclidean a, Eq b, Num b) => [(a, b)] -> Coprimes a b

-- | A list of pairwise coprime numbers with their multiplicities.
data Coprimes a b

-- | Unwrap.
unCoprimes :: Coprimes a b -> [(a, b)]

-- | Wrap a non-zero number with its multiplicity into <a>Coprimes</a>.
--   
--   <pre>
--   &gt;&gt;&gt; singleton 210 1
--   Coprimes {unCoprimes = [(210,1)]}
--   </pre>
singleton :: (Eq a, Num a, Eq b, Num b) => a -> b -> Coprimes a b

-- | Add a non-zero number with its multiplicity to <a>Coprimes</a>.
--   
--   <pre>
--   &gt;&gt;&gt; insert 360 1 (singleton 210 1)
--   Coprimes {unCoprimes = [(7,1),(5,2),(3,3),(2,4)]}
--   
--   &gt;&gt;&gt; insert 2 4 (insert 7 1 (insert 5 2 (singleton 4 3)))
--   Coprimes {unCoprimes = [(7,1),(5,2),(2,10)]}
--   </pre>
insert :: (Euclidean a, Eq b, Num b) => a -> b -> Coprimes a b -> Coprimes a b
instance (GHC.Show.Show a, GHC.Show.Show b) => GHC.Show.Show (Math.NumberTheory.Euclidean.Coprimes.Coprimes a b)
instance (GHC.Classes.Eq a, GHC.Classes.Eq b) => GHC.Classes.Eq (Math.NumberTheory.Euclidean.Coprimes.Coprimes a b)
instance (Math.NumberTheory.Euclidean.Euclidean a, GHC.Classes.Eq b, GHC.Num.Num b) => GHC.Base.Semigroup (Math.NumberTheory.Euclidean.Coprimes.Coprimes a b)
instance (Math.NumberTheory.Euclidean.Euclidean a, GHC.Classes.Eq b, GHC.Num.Num b) => GHC.Base.Monoid (Math.NumberTheory.Euclidean.Coprimes.Coprimes a b)


-- | Safe modular arithmetic with modulo on type level.
module Math.NumberTheory.Moduli.Class

-- | Wrapper for residues modulo <tt>m</tt>.
--   
--   <tt>Mod 3 :: Mod 10</tt> stands for the class of integers, congruent
--   to 3 modulo 10 (…−17, −7, 3, 13, 23…). The modulo is stored on type
--   level, so it is impossible, for example, to add up by mistake residues
--   with different moduli.
--   
--   <pre>
--   &gt;&gt;&gt; :set -XDataKinds
--   
--   &gt;&gt;&gt; (3 :: Mod 10) + (4 :: Mod 12)
--   error: Couldn't match type ‘12’ with ‘10’...
--   
--   &gt;&gt;&gt; (3 :: Mod 10) + 8
--   (1 `modulo` 10)
--   </pre>
--   
--   Note that modulo cannot be negative.
data Mod (m :: Nat)

-- | The canonical representative of the residue class, always between 0
--   and <tt>m-1</tt> inclusively.
getVal :: Mod m -> Integer

-- | The canonical representative of the residue class, always between 0
--   and <tt>m-1</tt> inclusively.
getNatVal :: Mod m -> Natural

-- | Linking type and value levels: extract modulo <tt>m</tt> as a value.
getMod :: KnownNat m => Mod m -> Integer

-- | Linking type and value levels: extract modulo <tt>m</tt> as a value.
getNatMod :: KnownNat m => Mod m -> Natural

-- | Computes the modular inverse, if the residue is coprime with the
--   modulo.
--   
--   <pre>
--   &gt;&gt;&gt; :set -XDataKinds
--   
--   &gt;&gt;&gt; invertMod (3 :: Mod 10)
--   Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10
--   
--   &gt;&gt;&gt; invertMod (4 :: Mod 10)
--   Nothing
--   </pre>
invertMod :: KnownNat m => Mod m -> Maybe (Mod m)

-- | Drop-in replacement for <a>^</a>, with much better performance.
--   
--   <pre>
--   &gt;&gt;&gt; :set -XDataKinds
--   
--   &gt;&gt;&gt; powMod (3 :: Mod 10) 4
--   &gt; (1 `modulo` 10)
--   </pre>
powMod :: (KnownNat m, Integral a) => Mod m -> a -> Mod m

-- | Infix synonym of <a>powMod</a>.
(^%) :: (KnownNat m, Integral a) => Mod m -> a -> Mod m
infixr 8 ^%

-- | This type represents elements of the multiplicative group mod m, i.e.
--   those elements which are coprime to m. Use <tt>toMultElement</tt> to
--   construct.
data MultMod m
multElement :: MultMod m -> Mod m

-- | Attempt to construct a multiplicative group element.
isMultElement :: KnownNat m => Mod m -> Maybe (MultMod m)

-- | For elements of the multiplicative group, we can safely perform the
--   inverse without needing to worry about failure.
invertGroup :: KnownNat m => MultMod m -> MultMod m

-- | This type represents residues with unknown modulo and rational
--   numbers. One can freely combine them in arithmetic expressions, but
--   each operation will spend time on modulo's recalculation:
--   
--   <pre>
--   &gt;&gt;&gt; 2 `modulo` 10 + 4 `modulo` 15
--   (1 `modulo` 5)
--   
--   &gt;&gt;&gt; 2 `modulo` 10 * 4 `modulo` 15
--   (3 `modulo` 5)
--   
--   &gt;&gt;&gt; 2 `modulo` 10 + fromRational (3 % 7)
--   (1 `modulo` 10)
--   
--   &gt;&gt;&gt; 2 `modulo` 10 * fromRational (3 % 7)
--   (8 `modulo` 10)
--   </pre>
--   
--   If performance is crucial, it is recommended to extract <tt>Mod m</tt>
--   for further processing by pattern matching. E. g.,
--   
--   <pre>
--   case modulo n m of
--     SomeMod k -&gt; process k -- Here k has type Mod m
--     InfMod{}  -&gt; error "impossible"
--   </pre>
data SomeMod
[SomeMod] :: KnownNat m => Mod m -> SomeMod
[InfMod] :: Rational -> SomeMod

-- | Create modular value by representative of residue class and modulo.
--   One can use the result either directly (via functions from <a>Num</a>
--   and <a>Fractional</a>), or deconstruct it by pattern matching. Note
--   that <a>modulo</a> never returns <a>InfMod</a>.
modulo :: Integer -> Natural -> SomeMod
infixl 7 `modulo`

-- | Computes the inverse value, if it exists.
--   
--   <pre>
--   &gt;&gt;&gt; invertSomeMod (3 `modulo` 10)
--   Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10
--   
--   &gt;&gt;&gt; invertMod (4 `modulo` 10)
--   Nothing
--   
--   &gt;&gt;&gt; invertSomeMod (fromRational (2 % 5))
--   Just 5 % 2
--   </pre>
invertSomeMod :: SomeMod -> Maybe SomeMod

-- | Drop-in replacement for <a>^</a>, with much better performance. When
--   -O is enabled, there is a rewrite rule, which specialises <a>^</a> to
--   <a>powSomeMod</a>.
--   
--   <pre>
--   &gt;&gt;&gt; powSomeMod (3 `modulo` 10) 4
--   (1 `modulo` 10)
--   </pre>
powSomeMod :: Integral a => SomeMod -> a -> SomeMod

-- | This class gives the integer associated with a type-level natural.
--   There are instances of the class for every concrete literal: 0, 1, 2,
--   etc.
class KnownNat (n :: Nat)
instance GHC.TypeNats.KnownNat m => GHC.Show.Show (Math.NumberTheory.Moduli.Class.MultMod m)
instance GHC.Classes.Ord (Math.NumberTheory.Moduli.Class.MultMod m)
instance GHC.Classes.Eq (Math.NumberTheory.Moduli.Class.MultMod m)
instance GHC.Enum.Enum (Math.NumberTheory.Moduli.Class.Mod m)
instance GHC.Classes.Ord (Math.NumberTheory.Moduli.Class.Mod m)
instance GHC.Classes.Eq (Math.NumberTheory.Moduli.Class.Mod m)
instance GHC.Classes.Eq Math.NumberTheory.Moduli.Class.SomeMod
instance GHC.Show.Show Math.NumberTheory.Moduli.Class.SomeMod
instance GHC.Num.Num Math.NumberTheory.Moduli.Class.SomeMod
instance GHC.Real.Fractional Math.NumberTheory.Moduli.Class.SomeMod
instance GHC.TypeNats.KnownNat m => GHC.Base.Semigroup (Math.NumberTheory.Moduli.Class.MultMod m)
instance GHC.TypeNats.KnownNat m => GHC.Base.Monoid (Math.NumberTheory.Moduli.Class.MultMod m)
instance GHC.TypeNats.KnownNat m => GHC.Enum.Bounded (Math.NumberTheory.Moduli.Class.MultMod m)
instance GHC.TypeNats.KnownNat m => GHC.Show.Show (Math.NumberTheory.Moduli.Class.Mod m)
instance GHC.TypeNats.KnownNat m => GHC.Enum.Bounded (Math.NumberTheory.Moduli.Class.Mod m)
instance GHC.TypeNats.KnownNat m => GHC.Num.Num (Math.NumberTheory.Moduli.Class.Mod m)
instance GHC.TypeNats.KnownNat m => GHC.Real.Fractional (Math.NumberTheory.Moduli.Class.Mod m)


-- | Efficient calculation of linear recurrent sequences, including
--   Fibonacci and Lucas sequences.
module Math.NumberTheory.Recurrencies.Linear

-- | Infinite zero-based table of factorials.
--   
--   <pre>
--   &gt;&gt;&gt; take 5 factorial
--   [1,1,2,6,24]
--   </pre>
--   
--   The time-and-space behaviour of <a>factorial</a> is similar to
--   described in <a>Math.NumberTheory.Recurrencies.Bilinear#memory</a>.
factorial :: (Num a, Enum a) => [a]

-- | <tt><a>fibonacci</a> k</tt> calculates the <tt>k</tt>-th Fibonacci
--   number in <i>O</i>(<tt>log (abs k)</tt>) steps. The index may be
--   negative. This is efficient for calculating single Fibonacci numbers
--   (with large index), but for computing many Fibonacci numbers in close
--   proximity, it is better to use the simple addition formula starting
--   from an appropriate pair of successive Fibonacci numbers.
fibonacci :: Num a => Int -> a

-- | <tt><a>fibonacciPair</a> k</tt> returns the pair <tt>(F(k),
--   F(k+1))</tt> of the <tt>k</tt>-th Fibonacci number and its successor,
--   thus it can be used to calculate the Fibonacci numbers from some index
--   on without needing to compute the previous. The pair is efficiently
--   calculated in <i>O</i>(<tt>log (abs k)</tt>) steps. The index may be
--   negative.
fibonacciPair :: Num a => Int -> (a, a)

-- | <tt><a>lucas</a> k</tt> computes the <tt>k</tt>-th Lucas number. Very
--   similar to <tt><a>fibonacci</a></tt>.
lucas :: Num a => Int -> a

-- | <tt><a>lucasPair</a> k</tt> computes the pair <tt>(L(k), L(k+1))</tt>
--   of the <tt>k</tt>-th Lucas number and its successor. Very similar to
--   <tt><a>fibonacciPair</a></tt>.
lucasPair :: Num a => Int -> (a, a)

-- | <tt><a>generalLucas</a> p q k</tt> calculates the quadruple <tt>(U(k),
--   U(k+1), V(k), V(k+1))</tt> where <tt>U(i)</tt> is the Lucas sequence
--   of the first kind and <tt>V(i)</tt> the Lucas sequence of the second
--   kind for the parameters <tt>p</tt> and <tt>q</tt>, where <tt>p^2-4q /=
--   0</tt>. Both sequences satisfy the recurrence relation <tt>A(j+2) =
--   p*A(j+1) - q*A(j)</tt>, the starting values are <tt>U(0) = 0, U(1) =
--   1</tt> and <tt>V(0) = 2, V(1) = p</tt>. The Fibonacci numbers form the
--   Lucas sequence of the first kind for the parameters <tt>p = 1, q =
--   -1</tt> and the Lucas numbers form the Lucas sequence of the second
--   kind for these parameters. Here, the index must be non-negative, since
--   the terms of the sequence for negative indices are in general not
--   integers.
generalLucas :: Num a => a -> a -> Int -> (a, a, a, a)


-- | Bilinear recurrent sequences and Bernoulli numbers, roughly covering
--   Ch. 5-6 of <i>Concrete Mathematics</i> by R. L. Graham, D. E. Knuth
--   and O. Patashnik.
--   
--   <b>Note on memory leaks and memoization.</b> Top-level definitions in
--   this module are polymorphic, so the results of computations are not
--   retained in memory. Make them monomorphic to take advantages of
--   memoization. Compare
--   
--   <pre>
--   &gt;&gt;&gt; :set +s
--   
--   &gt;&gt;&gt; binomial !! 1000 !! 1000 :: Integer
--   1
--   (0.01 secs, 1,385,512 bytes)
--   
--   &gt;&gt;&gt; binomial !! 1000 !! 1000 :: Integer
--   1
--   (0.01 secs, 1,381,616 bytes)
--   </pre>
--   
--   against
--   
--   <pre>
--   &gt;&gt;&gt; let binomial' = binomial :: [[Integer]]
--   
--   &gt;&gt;&gt; binomial' !! 1000 !! 1000 :: Integer
--   1
--   (0.01 secs, 1,381,696 bytes)
--   
--   &gt;&gt;&gt; binomial' !! 1000 !! 1000 :: Integer
--   1
--   (0.01 secs, 391,152 bytes)
--   </pre>
module Math.NumberTheory.Recurrencies.Bilinear

-- | Infinite zero-based table of binomial coefficients (also known as
--   Pascal triangle): <tt>binomial !! n !! k == n! / k! / (n - k)!</tt>.
--   
--   <pre>
--   &gt;&gt;&gt; take 5 (map (take 5) binomial)
--   [[1],[1,1],[1,2,1],[1,3,3,1],[1,4,6,4,1]]
--   </pre>
--   
--   Complexity: <tt>binomial !! n !! k</tt> is O(n) bits long, its
--   computation takes O(k n) time and forces thunks <tt>binomial !! n !!
--   i</tt> for <tt>0 &lt;= i &lt;= k</tt>. Use the symmetry of Pascal
--   triangle <tt>binomial !! n !! k == binomial !! n !! (n - k)</tt> to
--   speed up computations.
--   
--   One could also consider <a>binomial</a> to compute stand-alone values.
binomial :: Integral a => [[a]]

-- | Infinite zero-based table of <a>Stirling numbers of the first
--   kind</a>.
--   
--   <pre>
--   &gt;&gt;&gt; take 5 (map (take 5) stirling1)
--   [[1],[0,1],[0,1,1],[0,2,3,1],[0,6,11,6,1]]
--   </pre>
--   
--   Complexity: <tt>stirling1 !! n !! k</tt> is O(n ln n) bits long, its
--   computation takes O(k n^2 ln n) time and forces thunks <tt>stirling1
--   !! i !! j</tt> for <tt>0 &lt;= i &lt;= n</tt> and <tt>max(0, k - n +
--   i) &lt;= j &lt;= k</tt>.
--   
--   One could also consider <a>unsignedStirling1st</a> to compute
--   stand-alone values.
stirling1 :: (Num a, Enum a) => [[a]]

-- | Infinite zero-based table of <a>Stirling numbers of the second
--   kind</a>.
--   
--   <pre>
--   &gt;&gt;&gt; take 5 (map (take 5) stirling2)
--   [[1],[0,1],[0,1,1],[0,1,3,1],[0,1,7,6,1]]
--   </pre>
--   
--   Complexity: <tt>stirling2 !! n !! k</tt> is O(n ln n) bits long, its
--   computation takes O(k n^2 ln n) time and forces thunks <tt>stirling2
--   !! i !! j</tt> for <tt>0 &lt;= i &lt;= n</tt> and <tt>max(0, k - n +
--   i) &lt;= j &lt;= k</tt>.
--   
--   One could also consider <a>stirling2nd</a> to compute stand-alone
--   values.
stirling2 :: (Num a, Enum a) => [[a]]

-- | Infinite one-based table of <a>Lah numbers</a>. <tt>lah !! n !! k</tt>
--   equals to lah(n + 1, k + 1).
--   
--   <pre>
--   &gt;&gt;&gt; take 5 (map (take 5) lah)
--   [[1],[2,1],[6,6,1],[24,36,12,1],[120,240,120,20,1]]
--   </pre>
--   
--   Complexity: <tt>lah !! n !! k</tt> is O(n ln n) bits long, its
--   computation takes O(k n ln n) time and forces thunks <tt>lah !! n !!
--   i</tt> for <tt>0 &lt;= i &lt;= k</tt>.
lah :: Integral a => [[a]]

-- | Infinite zero-based table of <a>Eulerian numbers of the first
--   kind</a>.
--   
--   <pre>
--   &gt;&gt;&gt; take 5 (map (take 5) eulerian1)
--   [[],[1],[1,1],[1,4,1],[1,11,11,1]]
--   </pre>
--   
--   Complexity: <tt>eulerian1 !! n !! k</tt> is O(n ln n) bits long, its
--   computation takes O(k n^2 ln n) time and forces thunks <tt>eulerian1
--   !! i !! j</tt> for <tt>0 &lt;= i &lt;= n</tt> and <tt>max(0, k - n +
--   i) &lt;= j &lt;= k</tt>.
eulerian1 :: (Num a, Enum a) => [[a]]

-- | Infinite zero-based table of <a>Eulerian numbers of the second
--   kind</a>.
--   
--   <pre>
--   &gt;&gt;&gt; take 5 (map (take 5) eulerian2)
--   [[],[1],[1,2],[1,8,6],[1,22,58,24]]
--   </pre>
--   
--   Complexity: <tt>eulerian2 !! n !! k</tt> is O(n ln n) bits long, its
--   computation takes O(k n^2 ln n) time and forces thunks <tt>eulerian2
--   !! i !! j</tt> for <tt>0 &lt;= i &lt;= n</tt> and <tt>max(0, k - n +
--   i) &lt;= j &lt;= k</tt>.
eulerian2 :: (Num a, Enum a) => [[a]]

-- | Infinite zero-based sequence of <a>Bernoulli numbers</a>, computed via
--   <a>connection</a> with <a>stirling2</a>.
--   
--   <pre>
--   &gt;&gt;&gt; take 5 bernoulli
--   [1 % 1,(-1) % 2,1 % 6,0 % 1,(-1) % 30]
--   </pre>
--   
--   Complexity: <tt>bernoulli !! n</tt> is O(n ln n) bits long, its
--   computation takes O(n^3 ln n) time and forces thunks <tt>stirling2 !!
--   i !! j</tt> for <tt>0 &lt;= i &lt;= n</tt> and <tt>0 &lt;= j &lt;=
--   i</tt>.
--   
--   One could also consider <a>bernoulli</a> to compute stand-alone
--   values.
bernoulli :: Integral a => [Ratio a]

-- | The same sequence as <tt>euler'</tt>, but with type <tt>[a]</tt>
--   instead of <tt>[Ratio a]</tt> as the denominators in <tt>euler'</tt>
--   are always <tt>1</tt>.
--   
--   <pre>
--   &gt;&gt;&gt; take 10 euler :: [Integer]
--   [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]
--   </pre>
euler :: forall a. Integral a => [a]

-- | Infinite zero-based list of the <tt>n</tt>-th order Euler polynomials
--   evaluated at <tt>1</tt>. The algorithm used was derived from
--   <a>Algorithms for Bernoulli numbers and Euler numbers</a> by Kwang-Wu
--   Chen, third formula of the Corollary in page 7. Element-by-element
--   division of sequences <a>A1986631</a> and <a>A006519</a> in OEIS.
--   
--   <pre>
--   &gt;&gt;&gt; take 10 eulerPolyAt1 :: [Rational]
--   [1 % 1,1 % 2,0 % 1,(-1) % 4,0 % 1,1 % 2,0 % 1,(-17) % 8,0 % 1,31 % 2]
--   </pre>
eulerPolyAt1 :: forall a. Integral a => [Ratio a]


module Math.NumberTheory.Recurrencies

-- | Infinite zero-based table of <a>partition numbers</a>.
--   
--   <pre>
--   &gt;&gt;&gt; take 10 partition
--   [1, 1, 2, 3, 5, 7, 11, 15, 22, 30]
--   </pre>
--   
--   <pre>
--   &gt;&gt;&gt; :set -XDataKinds
--   
--   &gt;&gt;&gt; partition !! 1000 :: Mod 1000
--   (991 `modulo` 1000)
--   </pre>
partition :: Num a => [a]


-- | Functions dealing with squares. Efficient calculation of integer
--   square roots and efficient testing for squareness.
module Math.NumberTheory.Powers.Squares

-- | Calculate the integer square root of a nonnegative number <tt>n</tt>,
--   that is, the largest integer <tt>r</tt> with <tt>r*r &lt;= n</tt>.
--   Throws an error on negative input.
integerSquareRoot :: Integral a => a -> a

-- | Calculate the integer square root of a nonnegative number <tt>n</tt>,
--   that is, the largest integer <tt>r</tt> with <tt>r*r &lt;= n</tt>. The
--   precondition <tt>n &gt;= 0</tt> is not checked.
integerSquareRoot' :: Integral a => a -> a

-- | Calculate the integer square root of a nonnegative number as well as
--   the difference of that number with the square of that root, that is if
--   <tt>(s,r) = integerSquareRootRem n</tt> then <tt>s^2 &lt;= n == s^2+r
--   &lt; (s+1)^2</tt>.
integerSquareRootRem :: Integral a => a -> (a, a)

-- | Calculate the integer square root of a nonnegative number as well as
--   the difference of that number with the square of that root, that is if
--   <tt>(s,r) = integerSquareRootRem' n</tt> then <tt>s^2 &lt;= n == s^2+r
--   &lt; (s+1)^2</tt>. The precondition <tt>n &gt;= 0</tt> is not checked.
integerSquareRootRem' :: Integral a => a -> (a, a)

-- | Returns <a>Nothing</a> if the argument is not a square,
--   <tt><a>Just</a> r</tt> if <tt>r*r == n</tt> and <tt>r &gt;= 0</tt>.
--   Avoids the expensive calculation of the square root if <tt>n</tt> is
--   recognized as a non-square before, prevents repeated calculation of
--   the square root if only the roots of perfect squares are needed.
--   Checks for negativity and <a>isPossibleSquare</a>.
exactSquareRoot :: Integral a => a -> Maybe a

-- | Test whether the argument is a square. After a number is found to be
--   positive, first <a>isPossibleSquare</a> is checked, if it is, the
--   integer square root is calculated.
isSquare :: Integral a => a -> Bool

-- | Test whether the input (a nonnegative number) <tt>n</tt> is a square.
--   The same as <a>isSquare</a>, but without the negativity test. Faster
--   if many known positive numbers are tested.
--   
--   The precondition <tt>n &gt;= 0</tt> is not tested, passing negative
--   arguments may cause any kind of havoc.
isSquare' :: Integral a => a -> Bool

-- | Test whether a non-negative number may be a square. Non-negativity is
--   not checked, passing negative arguments may cause any kind of havoc.
--   
--   First the remainder modulo 256 is checked (that can be calculated
--   easily without division and eliminates about 82% of all numbers).
--   After that, the remainders modulo 9, 25, 7, 11 and 13 are tested to
--   eliminate altogether about 99.436% of all numbers.
--   
--   This is the test used by <a>exactSquareRoot</a>. For large numbers,
--   the slower but more discriminating test <tt>isPossibleSqure2</tt> is
--   faster.
isPossibleSquare :: Integral a => a -> Bool

-- | Test whether a non-negative number may be a square. Non-negativity is
--   not checked, passing negative arguments may cause any kind of havoc.
--   
--   First the remainder modulo 256 is checked (that can be calculated
--   easily without division and eliminates about 82% of all numbers).
--   After that, the remainders modulo several small primes are tested to
--   eliminate altogether about 99.98954% of all numbers.
--   
--   For smallish to medium sized numbers, this hardly performs better than
--   <a>isPossibleSquare</a>, which uses smaller arrays, but for large
--   numbers, where calculating the square root becomes more expensive, it
--   is much faster (if the vast majority of tested numbers aren't
--   squares).
isPossibleSquare2 :: Integral a => a -> Bool


-- | Functions dealing with fourth powers. Efficient calculation of integer
--   fourth roots and efficient testing for being a square's square.
module Math.NumberTheory.Powers.Fourth

-- | Calculate the integer fourth root of a nonnegative number, that is,
--   the largest integer <tt>r</tt> with <tt>r^4 &lt;= n</tt>. Throws an
--   error on negaitve input.
integerFourthRoot :: Integral a => a -> a

-- | Calculate the integer fourth root of a nonnegative number, that is,
--   the largest integer <tt>r</tt> with <tt>r^4 &lt;= n</tt>. The
--   condition is <i>not</i> checked.
integerFourthRoot' :: Integral a => a -> a

-- | Returns <tt>Nothing</tt> if <tt>n</tt> is not a fourth power, <tt>Just
--   r</tt> if <tt>n == r^4</tt> and <tt>r &gt;= 0</tt>.
exactFourthRoot :: Integral a => a -> Maybe a

-- | Test whether an integer is a fourth power. First nonnegativity is
--   checked, then the unchecked test is called.
isFourthPower :: Integral a => a -> Bool

-- | Test whether a nonnegative number is a fourth power. The condition is
--   <i>not</i> checked. If a number passes the
--   <a>isPossibleFourthPower</a> test, its integer fourth root is
--   calculated.
isFourthPower' :: Integral a => a -> Bool

-- | Test whether a nonnegative number is a possible fourth power. The
--   condition is <i>not</i> checked. This eliminates about 99.958% of
--   numbers.
isPossibleFourthPower :: Integral a => a -> Bool


-- | Functions dealing with cubes. Moderately efficient calculation of
--   integer cube roots and testing for cubeness.
module Math.NumberTheory.Powers.Cubes

-- | Calculate the integer cube root of an integer <tt>n</tt>, that is the
--   largest integer <tt>r</tt> such that <tt>r^3 &lt;= n</tt>. Note that
--   this is not symmetric about <tt>0</tt>, for example
--   <tt>integerCubeRoot (-2) = (-2)</tt> while <tt>integerCubeRoot 2 =
--   1</tt>.
integerCubeRoot :: Integral a => a -> a

-- | Calculate the integer cube root of a nonnegative integer <tt>n</tt>,
--   that is, the largest integer <tt>r</tt> such that <tt>r^3 &lt;=
--   n</tt>. The precondition <tt>n &gt;= 0</tt> is not checked.
integerCubeRoot' :: Integral a => a -> a

-- | Returns <tt>Nothing</tt> if the argument is not a cube, <tt>Just
--   r</tt> if <tt>n == r^3</tt>.
exactCubeRoot :: Integral a => a -> Maybe a

-- | Test whether an integer is a cube.
isCube :: Integral a => a -> Bool

-- | Test whether a nonnegative integer is a cube. Before
--   <a>integerCubeRoot</a> is calculated, a few tests of remainders modulo
--   small primes weed out most non-cubes. For testing many numbers, most
--   of which aren't cubes, this is much faster than <tt>let r = cubeRoot n
--   in r*r*r == n</tt>. The condition <tt>n &gt;= 0</tt> is <i>not</i>
--   checked.
isCube' :: Integral a => a -> Bool

-- | Test whether a nonnegative number is possibly a cube. Only about 0.08%
--   of all numbers pass this test. The precondition <tt>n &gt;= 0</tt> is
--   <i>not</i> checked.
isPossibleCube :: Integral a => a -> Bool


-- | Generalised Möbius inversion for <a>Int</a> valued functions.
module Math.NumberTheory.MoebiusInversion.Int

-- | <tt>generalInversion g n</tt> evaluates the generalised Möbius
--   inversion of <tt>g</tt> at the argument <tt>n</tt>.
--   
--   The generalised Möbius inversion implemented here allows an efficient
--   calculation of isolated values of the function <tt>f : N -&gt; Z</tt>
--   if the function <tt>g</tt> defined by
--   
--   <pre>
--   g n = sum [f (n `quot` k) | k &lt;- [1 .. n]]
--   </pre>
--   
--   can be cheaply computed. By the generalised Möbius inversion formula,
--   then
--   
--   <pre>
--   f n = sum [moebius k * g (n `quot` k) | k &lt;- [1 .. n]]
--   </pre>
--   
--   which allows the computation in <i>O</i>(n) steps, if the values of
--   the Möbius function are known. A slightly different formula, used
--   here, does not need the values of the Möbius function and allows the
--   computation in <i>O</i>(n^0.75) steps, using <i>O</i>(n^0.5) memory.
--   
--   An example of a pair of such functions where the inversion allows a
--   more efficient computation than the direct approach is
--   
--   <pre>
--   f n = number of reduced proper fractions with denominator &lt;= n
--   g n = number of proper fractions with denominator &lt;= n
--   </pre>
--   
--   (a <i>proper fraction</i> is a fraction <tt>0 &lt; n/d &lt; 1</tt>).
--   Then <tt>f n</tt> is the cardinality of the Farey sequence of order
--   <tt>n</tt> (minus 1 or 2 if 0 and/or 1 are included in the Farey
--   sequence), or the sum of the totients of the numbers <tt>2 &lt;= k
--   &lt;= n</tt>. <tt>f n</tt> is not easily directly computable, but then
--   <tt>g n = n*(n-1)/2</tt> is very easy to compute, and hence the
--   inversion gives an efficient method of computing <tt>f n</tt>.
--   
--   For <a>Int</a> valued functions, unboxed arrays can be used for
--   greater efficiency. That bears the risk of overflow, however, so be
--   sure to use it only when it's safe.
--   
--   The value <tt>f n</tt> is then computed by <tt>generalInversion g
--   n</tt>. Note that when many values of <tt>f</tt> are needed, there are
--   far more efficient methods, this method is only appropriate to compute
--   isolated values of <tt>f</tt>.
generalInversion :: (Int -> Int) -> Int -> Int

-- | <tt>totientSum n</tt> is, for <tt>n &gt; 0</tt>, the sum of
--   <tt>[totient k | k &lt;- [1 .. n]]</tt>, computed via generalised
--   Möbius inversion. See
--   <a>http://mathworld.wolfram.com/TotientSummatoryFunction.html</a> for
--   the formula used for <tt>totientSum</tt>.
totientSum :: Int -> Int


-- | Generalised Möbius inversion
module Math.NumberTheory.MoebiusInversion

-- | <tt>generalInversion g n</tt> evaluates the generalised Möbius
--   inversion of <tt>g</tt> at the argument <tt>n</tt>.
--   
--   The generalised Möbius inversion implemented here allows an efficient
--   calculation of isolated values of the function <tt>f : N -&gt; Z</tt>
--   if the function <tt>g</tt> defined by
--   
--   <pre>
--   g n = sum [f (n `quot` k) | k &lt;- [1 .. n]]
--   </pre>
--   
--   can be cheaply computed. By the generalised Möbius inversion formula,
--   then
--   
--   <pre>
--   f n = sum [moebius k * g (n `quot` k) | k &lt;- [1 .. n]]
--   </pre>
--   
--   which allows the computation in <i>O</i>(n) steps, if the values of
--   the Möbius function are known. A slightly different formula, used
--   here, does not need the values of the Möbius function and allows the
--   computation in <i>O</i>(n^0.75) steps, using <i>O</i>(n^0.5) memory.
--   
--   An example of a pair of such functions where the inversion allows a
--   more efficient computation than the direct approach is
--   
--   <pre>
--   f n = number of reduced proper fractions with denominator &lt;= n
--   
--   g n = number of proper fractions with denominator &lt;= n
--   </pre>
--   
--   (a <i>proper fraction</i> is a fraction <tt>0 &lt; n/d &lt; 1</tt>).
--   Then <tt>f n</tt> is the cardinality of the Farey sequence of order
--   <tt>n</tt> (minus 1 or 2 if 0 and/or 1 are included in the Farey
--   sequence), or the sum of the totients of the numbers <tt>2 &lt;= k
--   &lt;= n</tt>. <tt>f n</tt> is not easily directly computable, but then
--   <tt>g n = n*(n-1)/2</tt> is very easy to compute, and hence the
--   inversion gives an efficient method of computing <tt>f n</tt>.
--   
--   Since the function arguments are used as array indices, the domain of
--   <tt>f</tt> is restricted to <a>Int</a>.
--   
--   The value <tt>f n</tt> is then computed by <tt>generalInversion g
--   n</tt>. Note that when many values of <tt>f</tt> are needed, there are
--   far more efficient methods, this method is only appropriate to compute
--   isolated values of <tt>f</tt>.
generalInversion :: (Int -> Integer) -> Int -> Integer

-- | <tt>totientSum n</tt> is, for <tt>n &gt; 0</tt>, the sum of
--   <tt>[totient k | k &lt;- [1 .. n]]</tt>, computed via generalised
--   Möbius inversion. See
--   <a>http://mathworld.wolfram.com/TotientSummatoryFunction.html</a> for
--   the formula used for <tt>totientSum</tt>.
totientSum :: Int -> Integer


-- | Calculating integer roots and determining perfect powers. The
--   algorithms are moderately efficient.
module Math.NumberTheory.Powers.General

-- | Calculate an integer root, <tt><a>integerRoot</a> k n</tt> computes
--   the (floor of) the <tt>k</tt>-th root of <tt>n</tt>, where <tt>k</tt>
--   must be positive. <tt>r = <a>integerRoot</a> k n</tt> means <tt>r^k
--   &lt;= n &lt; (r+1)^k</tt> if that is possible at all. It is impossible
--   if <tt>k</tt> is even and <tt>n &lt; 0</tt>, since then <tt>r^k &gt;=
--   0</tt> for all <tt>r</tt>, then, and if <tt>k &lt;= 0</tt>,
--   <tt><a>integerRoot</a></tt> raises an error. For <tt>k &lt; 5</tt>, a
--   specialised version is called which should be more efficient than the
--   general algorithm. However, it is not guaranteed that the rewrite
--   rules for those fire, so if <tt>k</tt> is known in advance, it is
--   safer to directly call the specialised versions.
integerRoot :: (Integral a, Integral b) => b -> a -> a

-- | <tt><a>exactRoot</a> k n</tt> returns <tt><a>Nothing</a></tt> if
--   <tt>n</tt> is not a <tt>k</tt>-th power, <tt><a>Just</a> r</tt> if
--   <tt>n == r^k</tt>. If <tt>k</tt> is divisible by <tt>4, 3</tt> or
--   <tt>2</tt>, a residue test is performed to avoid the expensive
--   calculation if it can thus be determined that <tt>n</tt> is not a
--   <tt>k</tt>-th power.
exactRoot :: (Integral a, Integral b) => b -> a -> Maybe a

-- | <tt><a>isKthPower</a> k n</tt> checks whether <tt>n</tt> is a
--   <tt>k</tt>-th power.
isKthPower :: (Integral a, Integral b) => b -> a -> Bool

-- | <tt><a>isPerfectPower</a> n</tt> checks whether <tt>n == r^k</tt> for
--   some <tt>k &gt; 1</tt>.
isPerfectPower :: Integral a => a -> Bool

-- | <tt><a>highestPower</a> n</tt> produces the pair <tt>(b,k)</tt> with
--   the largest exponent <tt>k</tt> such that <tt>n == b^k</tt>, except
--   for <tt><a>abs</a> n &lt;= 1</tt>, in which case arbitrarily large
--   exponents exist, and by an arbitrary decision <tt>(n,3)</tt> is
--   returned.
--   
--   First, by trial division with small primes, the range of possible
--   exponents is reduced (if <tt>p^e</tt> exactly divides <tt>n</tt>, then
--   <tt>k</tt> must be a divisor of <tt>e</tt>, if several small primes
--   divide <tt>n</tt>, <tt>k</tt> must divide the greatest common divisor
--   of their exponents, which mostly will be <tt>1</tt>, generally small;
--   if none of the small primes divides <tt>n</tt>, the range of possible
--   exponents is reduced since the base is necessarily large), if that has
--   not yet determined the result, the remaining factor is examined by
--   trying the divisors of the <tt>gcd</tt> of the prime exponents if some
--   have been found, otherwise by trying prime exponents recursively.
highestPower :: Integral a => a -> (a, Int)

-- | <tt><a>largePFPower</a> bd n</tt> produces the pair <tt>(b,k)</tt>
--   with the largest exponent <tt>k</tt> such that <tt>n == b^k</tt>,
--   where <tt>bd &gt; 1</tt> (it is expected that <tt>bd</tt> is much
--   larger, at least <tt>1000</tt> or so), <tt>n &gt; bd^2</tt> and
--   <tt>n</tt> has no prime factors <tt>p &lt;= bd</tt>, skipping the
--   trial division phase of <tt><a>highestPower</a></tt> when that is a
--   priori known to be superfluous. It is only present to avoid
--   duplication of work in factorisation and primality testing, it is not
--   expected to be generally useful. The assumptions are not checked, if
--   they are not satisfied, wrong results and wasted work may be the
--   consequence.
largePFPower :: Integer -> Integer -> (Integer, Int)


-- | <a>Jacobi symbol</a> is a generalization of the Legendre symbol,
--   useful for primality testing and integer factorization.
module Math.NumberTheory.Moduli.Jacobi

-- | Represents three possible values of <a>Jacobi symbol</a>.
data JacobiSymbol
MinusOne :: JacobiSymbol
Zero :: JacobiSymbol
One :: JacobiSymbol

-- | <a>Jacobi symbol</a> of two arguments. The lower argument
--   ("denominator") must be odd and positive, this condition is checked.
--   
--   If arguments have a common factor, the result is <a>Zero</a>,
--   otherwise it is <a>MinusOne</a> or <a>One</a>.
--   
--   <pre>
--   &gt;&gt;&gt; jacobi 1001 9911
--   Zero -- arguments have a common factor 11
--   
--   &gt;&gt;&gt; jacobi 1001 9907
--   MinusOne
--   </pre>
jacobi :: (Integral a, Bits a) => a -> a -> JacobiSymbol

-- | <i>Deprecated: Use <a>jacobi</a> instead</i>
jacobi' :: (Integral a, Bits a) => a -> a -> JacobiSymbol
instance GHC.Show.Show Math.NumberTheory.Moduli.Jacobi.JacobiSymbol
instance GHC.Classes.Ord Math.NumberTheory.Moduli.Jacobi.JacobiSymbol
instance GHC.Classes.Eq Math.NumberTheory.Moduli.Jacobi.JacobiSymbol
instance GHC.Base.Semigroup Math.NumberTheory.Moduli.Jacobi.JacobiSymbol
instance GHC.Base.Monoid Math.NumberTheory.Moduli.Jacobi.JacobiSymbol


-- | Chinese remainder theorem
module Math.NumberTheory.Moduli.Chinese

-- | Given a list <tt>[(r_1,m_1), ..., (r_n,m_n)]</tt> of
--   <tt>(residue,modulus)</tt> pairs, <tt>chineseRemainder</tt> calculates
--   the solution to the simultaneous congruences
--   
--   <pre>
--   r ≡ r_k (mod m_k)
--   </pre>
--   
--   if all moduli are positive and pairwise coprime. Otherwise the result
--   is <tt>Nothing</tt> regardless of whether a solution exists.
chineseRemainder :: [(Integer, Integer)] -> Maybe Integer

-- | <tt>chineseRemainder2 (r_1,m_1) (r_2,m_2)</tt> calculates the solution
--   of
--   
--   <pre>
--   r ≡ r_k (mod m_k)
--   </pre>
--   
--   if <tt>m_1</tt> and <tt>m_2</tt> are coprime.
chineseRemainder2 :: (Integer, Integer) -> (Integer, Integer) -> Integer


-- | Low level gcd and coprimality functions using the binary gcd
--   algorithm. Normally, accessing these via the higher level interface of
--   <a>Math.NumberTheory.GCD</a> should be sufficient.
module Math.NumberTheory.GCD.LowLevel

-- | Greatest common divisor of two <a>Int</a>s, calculated with the binary
--   gcd algorithm.

-- | <i>Deprecated: Use Math.NumberTheory.Euclidean.gcd</i>
gcdInt :: Int -> Int -> Int

-- | Greatest common divisor of two <a>Word</a>s, calculated with the
--   binary gcd algorithm.

-- | <i>Deprecated: Use Math.NumberTheory.Euclidean.gcd</i>
gcdWord :: Word -> Word -> Word

-- | Greatest common divisor of two <a>Int#</a>s, calculated with the
--   binary gcd algorithm.

-- | <i>Deprecated: Use Math.NumberTheory.Euclidean.gcd</i>
gcdInt# :: Int# -> Int# -> Int#

-- | Greatest common divisor of two <a>Word#</a>s, calculated with the
--   binary gcd algorithm.

-- | <i>Deprecated: Use Math.NumberTheory.Euclidean.gcd</i>
gcdWord# :: Word# -> Word# -> Word#

-- | Test whether two <a>Int</a>s are coprime, using an abbreviated binary
--   gcd algorithm.

-- | <i>Deprecated: Math.NumberTheory.Euclidean.</i>
coprimeInt :: Int -> Int -> Bool

-- | Test whether two <a>Word</a>s are coprime, using an abbreviated binary
--   gcd algorithm.

-- | <i>Deprecated: Math.NumberTheory.Euclidean.</i>
coprimeWord :: Word -> Word -> Bool

-- | Test whether two <a>Int#</a>s are coprime.

-- | <i>Deprecated: Math.NumberTheory.Euclidean.</i>
coprimeInt# :: Int# -> Int# -> Bool

-- | Test whether two <a>Word#</a>s are coprime.

-- | <i>Deprecated: Math.NumberTheory.Euclidean.</i>
coprimeWord# :: Word# -> Word# -> Bool


-- | This module exports GCD and coprimality test using the binary gcd
--   algorithm and GCD with the extended Euclidean algorithm.
--   
--   Efficiently counting the number of trailing zeros, the binary gcd
--   algorithm can perform considerably faster than the Euclidean algorithm
--   on average. For <a>Int</a>, GHC has a rewrite rule to use GMP's fast
--   gcd, depending on hardware and/or GMP version, that can be faster or
--   slower than the binary algorithm (on my 32-bit box, binary is faster,
--   on my 64-bit box, GMP). For <a>Word</a> and the sized
--   <tt>IntN/WordN</tt> types, there is no rewrite rule (yet) in GHC, and
--   the binary algorithm performs consistently (so far as my tests go)
--   much better (if this module's rewrite rules fire).
--   
--   When using this module, always compile with optimisations turned on to
--   benefit from GHC's primops and the rewrite rules.
module Math.NumberTheory.GCD

-- | Calculate the greatest common divisor using the binary gcd algorithm.
--   Depending on type and hardware, that can be considerably faster than
--   <tt><a>gcd</a></tt> but it may also be significantly slower.
--   
--   There are specialised functions for <tt><a>Int</a></tt> and
--   <tt><a>Word</a></tt> and rewrite rules for those and <tt>IntN</tt> and
--   <tt>WordN</tt>, <tt>N &lt;= 64</tt>, to use the specialised variants.
--   These types are worth benchmarking, others probably not.
--   
--   It is very slow for <a>Integer</a> (and probably every type except the
--   abovementioned), I recommend not using it for those.
--   
--   Relies on twos complement or sign and magnitude representaion for
--   signed types.

-- | <i>Deprecated: Use <a>gcd</a></i>
binaryGCD :: (Integral a, Bits a) => a -> a -> a

-- | Calculate the greatest common divisor of two numbers and coefficients
--   for the linear combination.
--   
--   For signed types satisfies:
--   
--   <pre>
--   case extendedGCD a b of
--     (d, u, v) -&gt; u*a + v*b == d
--                  &amp;&amp; d == gcd a b
--   </pre>
--   
--   For unsigned and bounded types the property above holds, but since
--   <tt>u</tt> and <tt>v</tt> must also be unsigned, the result may look
--   weird. E. g., on 64-bit architecture
--   
--   <pre>
--   extendedGCD (2 :: Word) (3 :: Word) == (1, 2^64-1, 1)
--   </pre>
--   
--   For unsigned and unbounded types (like <a>Natural</a>) the result is
--   undefined.
--   
--   For signed types we also have
--   
--   <pre>
--   abs u &lt; abs b || abs b &lt;= 1
--   
--   abs v &lt; abs a || abs a &lt;= 1
--   </pre>
--   
--   (except if one of <tt>a</tt> and <tt>b</tt> is <a>minBound</a> of a
--   signed type).

-- | <i>Deprecated: Use <a>extendedGCD</a></i>
extendedGCD :: Integral a => a -> a -> (a, a, a)

-- | Test whether two numbers are coprime using an abbreviated binary gcd
--   algorithm. A little bit faster than checking <tt>binaryGCD a b ==
--   1</tt> if one of the arguments is even, much faster if both are even.
--   
--   The remarks about performance at <a>binaryGCD</a> apply here too, use
--   this function only at the types with rewrite rules.
--   
--   Relies on twos complement or sign and magnitude representaion for
--   signed types.

-- | <i>Deprecated: Use <a>coprime</a></i>
coprime :: (Integral a, Bits a) => a -> a -> Bool


-- | Arithmetic on Montgomery elliptic curve.
module Math.NumberTheory.Curves.Montgomery

-- | We use the Montgomery form of elliptic curve: b Y² = X³ + a X² + X
--   (mod n). See Eq. (10.3.1.1) at p. 260 of <a>Speeding the Pollard and
--   Elliptic Curve Methods of Factorization</a> by P. L. Montgomery.
--   
--   Switching to projective space by substitutions Y = y / z, X = x / z,
--   we get b y² z = x³ + a x² z + x z² (mod n). The point on projective
--   elliptic curve is characterized by three coordinates, but it appears
--   that only x- and z-components matter for computations. By the same
--   reason there is no need to store coefficient b.
--   
--   That said, the chosen curve is represented by a24 = (a + 2) / 4 and
--   modulo n at type level, making points on different curves
--   incompatible.
data Point (a24 :: Nat) (n :: Nat)

-- | Extract x-coordinate.
pointX :: Point a24 n -> Integer

-- | Extract z-coordinate.
pointZ :: Point a24 n -> Integer
pointN :: forall a24 n. KnownNat n => Point a24 n -> Integer
pointA24 :: forall a24 n. KnownNat a24 => Point a24 n -> Integer

-- | Point on unknown curve.
data SomePoint
[SomePoint] :: (KnownNat a24, KnownNat n) => Point a24 n -> SomePoint

-- | <a>newPoint</a> <tt>s</tt> <tt>n</tt> creates a point on an elliptic
--   curve modulo <tt>n</tt>, uniquely determined by seed <tt>s</tt>. Some
--   choices of <tt>s</tt> and <tt>n</tt> produce ill-parametrized curves,
--   which is reflected by return value <a>Nothing</a>.
--   
--   We choose a curve by Suyama's parametrization. See Eq. (3)-(4) at p. 4
--   of <a>Implementing the Elliptic Curve Method of Factoring in
--   Reconfigurable Hardware</a> by K. Gaj, S. Kwon et al.
newPoint :: Integer -> Integer -> Maybe SomePoint

-- | If <tt>p0</tt> + <tt>p1</tt> = <tt>p2</tt>, then <a>add</a>
--   <tt>p0</tt> <tt>p1</tt> <tt>p2</tt> equals to <tt>p1</tt> +
--   <tt>p2</tt>. It is also required that z-coordinates of <tt>p0</tt>,
--   <tt>p1</tt> and <tt>p2</tt> are coprime with modulo of elliptic curve;
--   and x-coordinate of <tt>p0</tt> is non-zero. If preconditions do not
--   hold, return value is undefined.
--   
--   Remarkably such addition does not require <a>KnownNat</a> <tt>a24</tt>
--   constraint.
--   
--   Computations follow Algorithm 3 at p. 4 of <a>Implementing the
--   Elliptic Curve Method of Factoring in Reconfigurable Hardware</a> by
--   K. Gaj, S. Kwon et al.
add :: KnownNat n => Point a24 n -> Point a24 n -> Point a24 n -> Point a24 n

-- | Multiply by 2.
--   
--   Computations follow Algorithm 3 at p. 4 of <a>Implementing the
--   Elliptic Curve Method of Factoring in Reconfigurable Hardware</a> by
--   K. Gaj, S. Kwon et al.
double :: (KnownNat a24, KnownNat n) => Point a24 n -> Point a24 n

-- | Multiply by given number, using binary algorithm.
multiply :: (KnownNat a24, KnownNat n) => Word -> Point a24 n -> Point a24 n
instance GHC.Show.Show Math.NumberTheory.Curves.Montgomery.SomePoint
instance GHC.TypeNats.KnownNat n => GHC.Classes.Eq (Math.NumberTheory.Curves.Montgomery.Point a24 n)
instance (GHC.TypeNats.KnownNat a24, GHC.TypeNats.KnownNat n) => GHC.Show.Show (Math.NumberTheory.Curves.Montgomery.Point a24 n)


-- | Prime generation using a sieve. Currently, an enhanced sieve of
--   Eratosthenes is used, switching to an Atkin sieve is planned (if I get
--   around to implementing it and it's not slower).
--   
--   The sieve used is segmented, with a chunk size chosen to give good
--   (enough) cache locality while still getting something substantial done
--   per chunk. However, that means we must store data for primes up to the
--   square root of where sieving is done, thus sieving primes up to
--   <tt>n</tt> requires <tt><i>O</i>(sqrt n/log n)</tt> space.
module Math.NumberTheory.Primes.Sieve

-- | Ascending list of primes.
--   
--   <pre>
--   &gt;&gt;&gt; take 10 primes
--   [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
--   </pre>
--   
--   <a>primes</a> is a polymorphic list, so the results of computations
--   are not retained in memory. Make it monomorphic to take advantages of
--   memoization. Compare
--   
--   <pre>
--   &gt;&gt;&gt; :set +s
--   
--   &gt;&gt;&gt; primes !! 1000000 :: Int
--   15485867
--   (5.32 secs, 6,945,267,496 bytes)
--   
--   &gt;&gt;&gt; primes !! 1000000 :: Int
--   15485867
--   (5.19 secs, 6,945,267,496 bytes)
--   </pre>
--   
--   against
--   
--   <pre>
--   &gt;&gt;&gt; let primes' = primes :: [Int]
--   
--   &gt;&gt;&gt; primes' !! 1000000 :: Int
--   15485867
--   (5.29 secs, 6,945,269,856 bytes)
--   
--   &gt;&gt;&gt; primes' !! 1000000 :: Int
--   15485867
--   (0.02 secs, 336,232 bytes)
--   </pre>
primes :: (Ord a, Num a) => [a]

-- | <tt><a>sieveFrom</a> n</tt> creates the list of primes not less than
--   <tt>n</tt>.
sieveFrom :: Integer -> [Integer]

-- | Compact store of primality flags.
data PrimeSieve

-- | Sieve primes up to (and including) a bound (or 7, if bound is
--   smaller). For small enough bounds, this is more efficient than using
--   the segmented sieve.
--   
--   Since arrays are <a>Int</a>-indexed, overflow occurs when the sieve
--   size comes near <tt><a>maxBound</a> :: <a>Int</a></tt>, that
--   corresponds to an upper bound near <tt>15/8*<a>maxBound</a></tt>. On
--   <tt>32</tt>-bit systems, that is often within memory limits, so don't
--   give bounds larger than <tt>8*10^9</tt> there.
primeSieve :: Integer -> PrimeSieve

-- | List of primes in the form of a list of <a>PrimeSieve</a>s, more
--   compact than <a>primes</a>, thus it may be better to use
--   <tt><a>psieveList</a> &gt;&gt;= <a>primeList</a></tt> than keeping the
--   list of primes alive during the entire run.
psieveList :: [PrimeSieve]

-- | <tt><a>psieveFrom</a> n</tt> creates the list of <a>PrimeSieve</a>s
--   starting roughly at <tt>n</tt>. Due to the organisation of the sieve,
--   the list may contain a few primes less than <tt>n</tt>. This form uses
--   less memory than <tt>[<a>Integer</a>]</tt>, hence it may be preferable
--   to use this if it is to be reused.
psieveFrom :: Integer -> [PrimeSieve]

-- | Generate a list of primes for consumption from a <a>PrimeSieve</a>.
primeList :: forall a. Integral a => PrimeSieve -> [a]


-- | A <a>smooth number</a> is an integer, which can be represented as a
--   product of powers of elements from a given set (smooth basis). E. g.,
--   48 = 3 * 4 * 4 is smooth over a set {3, 4}, and 24 is not.
module Math.NumberTheory.SmoothNumbers

-- | An abstract representation of a smooth basis. It consists of a set of
--   coprime numbers ≥2.
data SmoothBasis a

-- | Build a <a>SmoothBasis</a> from a set of coprime numbers ≥2.
--   
--   <pre>
--   &gt;&gt;&gt; import qualified Data.Set as Set
--   
--   &gt;&gt;&gt; fromSet (Set.fromList [2, 3])
--   Just (SmoothBasis [2, 3])
--   
--   &gt;&gt;&gt; fromSet (Set.fromList [2, 4]) -- should be coprime
--   Nothing
--   
--   &gt;&gt;&gt; fromSet (Set.fromList [1, 3]) -- should be &gt;= 2
--   Nothing
--   </pre>
fromSet :: Euclidean a => Set a -> Maybe (SmoothBasis a)

-- | Build a <a>SmoothBasis</a> from a list of coprime numbers ≥2.
--   
--   <pre>
--   &gt;&gt;&gt; fromList [2, 3]
--   Just (SmoothBasis [2, 3])
--   
--   &gt;&gt;&gt; fromList [2, 2]
--   Just (SmoothBasis [2])
--   
--   &gt;&gt;&gt; fromList [2, 4] -- should be coprime
--   Nothing
--   
--   &gt;&gt;&gt; fromList [1, 3] -- should be &gt;= 2
--   Nothing
--   </pre>
fromList :: Euclidean a => [a] -> Maybe (SmoothBasis a)

-- | Build a <a>SmoothBasis</a> from a list of primes below given bound.
--   
--   <pre>
--   &gt;&gt;&gt; fromSmoothUpperBound 10
--   Just (SmoothBasis [2, 3, 5, 7])
--   
--   &gt;&gt;&gt; fromSmoothUpperBound 1
--   Nothing
--   </pre>
fromSmoothUpperBound :: Integral a => a -> Maybe (SmoothBasis a)

-- | Generate an infinite ascending list of <a>smooth numbers</a> over a
--   given smooth basis.
--   
--   <pre>
--   &gt;&gt;&gt; import Data.Maybe
--   
--   &gt;&gt;&gt; take 10 (smoothOver (fromJust (fromList [2, 5])))
--   [1, 2, 4, 5, 8, 10, 16, 20, 25, 32]
--   </pre>
smoothOver :: Integral a => SmoothBasis a -> [a]

-- | Generate an ascending list of <a>smooth numbers</a> over a given
--   smooth basis in a given range.
--   
--   It may appear inefficient for short, but distant ranges; consider
--   using <a>smoothOverInRangeBF</a> in such cases.
--   
--   <pre>
--   &gt;&gt;&gt; import Data.Maybe
--   
--   &gt;&gt;&gt; smoothOverInRange (fromJust (fromList [2, 5])) 100 200
--   [100, 125, 128, 160, 200]
--   </pre>
smoothOverInRange :: forall a. Integral a => SmoothBasis a -> a -> a -> [a]

-- | Generate an ascending list of <a>smooth numbers</a> over a given
--   smooth basis in a given range.
--   
--   It is inefficient for large or starting near 0 ranges; consider using
--   <a>smoothOverInRange</a> in such cases.
--   
--   Suffix BF stands for the brute force algorithm, involving a lot of
--   divisions.
--   
--   <pre>
--   &gt;&gt;&gt; import Data.Maybe
--   
--   &gt;&gt;&gt; smoothOverInRangeBF (fromJust (fromList [2, 5])) 100 200
--   [100, 125, 128, 160, 200]
--   </pre>
smoothOverInRangeBF :: forall a. Integral a => SmoothBasis a -> a -> a -> [a]
instance GHC.Show.Show a => GHC.Show.Show (Math.NumberTheory.SmoothNumbers.SmoothBasis a)
instance GHC.Classes.Eq a => GHC.Classes.Eq (Math.NumberTheory.SmoothNumbers.SmoothBasis a)


-- | Primality tests.
module Math.NumberTheory.Primes.Testing

-- | <tt>isPrime n</tt> tests whether <tt>n</tt> is a prime (negative or
--   positive). It is a combination of trial division and Baillie-PSW test.
--   
--   If <tt>isPrime n</tt> returns <tt>False</tt> then <tt>n</tt> is
--   definitely composite. There is a theoretical possibility that
--   <tt>isPrime n</tt> is <tt>True</tt>, but in fact <tt>n</tt> is not
--   prime. However, no such numbers are known and none exist below
--   <tt>2^64</tt>. If you have found one, please report it, because it is
--   a major discovery.
isPrime :: Integer -> Bool

-- | <tt><a>isCertifiedPrime</a> n</tt> tests primality of <tt>n</tt>,
--   first trial division by small primes is performed, then a Baillie PSW
--   test and finally a prime certificate is constructed and verified,
--   provided no step before found <tt>n</tt> to be composite. Constructing
--   prime certificates can take a <i>very</i> long time, so use this with
--   care.
isCertifiedPrime :: Integer -> Bool

-- | Primality test after Baillie, Pomerance, Selfridge and Wagstaff. The
--   Baillie-PSW test consists of a strong Fermat probable primality test
--   followed by a (strong) Lucas primality test. This implementation
--   assumes that the number <tt>n</tt> to test is odd and larger than
--   <tt>3</tt>. Even and small numbers have to be handled before. Also,
--   before applying this test, trial division by small primes should be
--   performed to identify many composites cheaply (although the
--   Baillie-PSW test is rather fast, about the same speed as a strong
--   Fermat test for four or five bases usually, it is, for large numbers,
--   much more costly than trial division by small primes, the primes less
--   than <tt>1000</tt>, say, so eliminating numbers with small prime
--   factors beforehand is more efficient).
--   
--   The Baillie-PSW test is very reliable, so far no composite numbers
--   passing it are known, and it is known (Gilchrist 2010) that no
--   Baillie-PSW pseudoprimes exist below <tt>2^64</tt>. However, a
--   heuristic argument by Pomerance indicates that there are likely
--   infinitely many Baillie-PSW pseudoprimes. On the other hand, according
--   to <a>http://mathworld.wolfram.com/Baillie-PSWPrimalityTest.html</a>
--   there is reason to believe that there are none with less than several
--   thousand digits, so that for most use cases the test can be considered
--   definitive.
bailliePSW :: Integer -> Bool

-- | Miller-Rabin probabilistic primality test. It consists of the trial
--   division test and several rounds of the strong Fermat test with
--   different bases. The choice of trial divisors and bases are
--   implementation details and may change in future silently.
--   
--   First argument stands for the number of rounds of strong Fermat test.
--   If it is 0, only trial division test is performed.
--   
--   If <tt>millerRabinV k n</tt> returns <tt>False</tt> then <tt>n</tt> is
--   definitely composite. Otherwise <tt>n</tt> may appear composite with
--   probability <tt>1/4^k</tt>.
millerRabinV :: Int -> Integer -> Bool

-- | <tt><a>isStrongFermatPP</a> n b</tt> tests whether non-negative
--   <tt>n</tt> is a strong Fermat probable prime for base <tt>b</tt>.
--   
--   Apart from primes, also some composite numbers have the tested
--   property, but those are rare. Very rare are composite numbers having
--   the property for many bases, so testing a large prime candidate with
--   several bases can identify composite numbers with high probability. An
--   odd number <tt>n &gt; 3</tt> is prime if and only if
--   <tt><a>isStrongFermatPP</a> n b</tt> holds for all <tt>b</tt> with
--   <tt>2 &lt;= b &lt;= (n-1)/2</tt>, but of course checking all those
--   bases would be less efficient than trial division, so one normally
--   checks only a relatively small number of bases, depending on the
--   desired degree of certainty. The probability that a randomly chosen
--   base doesn't identify a composite number <tt>n</tt> is less than
--   <tt>1/4</tt>, so five to ten tests give a reasonable level of
--   certainty in general.
--   
--   Please consult <a>Deterministic variants of the Miller-Rabin primality
--   test</a> for the best choice of bases.
isStrongFermatPP :: Integer -> Integer -> Bool

-- | <tt><a>isFermatPP</a> n b</tt> tests whether <tt>n</tt> is a Fermat
--   probable prime for the base <tt>b</tt>, that is, whether <tt>b^(n-1)
--   <a>mod</a> n == 1</tt>. This is a weaker but simpler condition.
--   However, more is lost in strength than is gained in simplicity, so for
--   primality testing, the strong check should be used. The remarks about
--   the choice of bases to test from <tt><a>isStrongFermatPP</a></tt>
--   apply with the modification that if <tt>a</tt> and <tt>b</tt> are
--   Fermat bases for <tt>n</tt>, then <tt>a*b</tt> <i>always</i> is a
--   Fermat base for <tt>n</tt> too. A <i>Charmichael number</i> is a
--   composite number <tt>n</tt> which is a Fermat probable prime for all
--   bases <tt>b</tt> coprime to <tt>n</tt>. By the above, only primes
--   <tt>p &lt;= n/2</tt> not dividing <tt>n</tt> need to be tested to
--   identify Carmichael numbers (however, testing all those primes would
--   be less efficient than determining Carmichaelness from the prime
--   factorisation; but testing an appropriate number of prime bases is
--   reasonable to find out whether it's worth the effort to undertake the
--   prime factorisation).
isFermatPP :: Integer -> Integer -> Bool

-- | <tt><a>trialDivisionPrimeTo</a> bound n</tt> tests whether <tt>n</tt>
--   is coprime to all primes <tt>&lt;= bound</tt>. If <tt>n &lt;=
--   bound^2</tt>, this is a full prime test, but very slow if <tt>n</tt>
--   has no small prime divisors.
trialDivisionPrimeTo :: Integer -> Integer -> Bool


-- | Certificates for primality or compositeness.
module Math.NumberTheory.Primes.Testing.Certificates

-- | A certificate of either compositeness or primality of an
--   <a>Integer</a>. Only numbers <tt>&gt; 1</tt> can be certified, trying
--   to create a certificate for other numbers raises an error.
data Certificate
Composite :: !CompositenessProof -> Certificate
Prime :: !PrimalityProof -> Certificate
argueCertificate :: Certificate -> Either CompositenessArgument PrimalityArgument

-- | A proof of compositeness of a positive number. The type is abstract to
--   ensure the validity of proofs.
data CompositenessProof

-- | The number whose compositeness is proved.
composite :: CompositenessProof -> Integer

-- | A proof of primality of a positive number. The type is abstract to
--   ensure the validity of proofs.
data PrimalityProof

-- | The number whose primality is proved.
cprime :: PrimalityProof -> Integer

-- | An argument for compositeness of a number (which must be <tt>&gt;
--   1</tt>). <a>CompositenessProof</a>s translate directly to
--   <tt>CompositenessArguments</tt>, correct arguments can be transformed
--   into proofs. This type allows the manipulation of proofs while
--   maintaining their correctness. The only way to access components of a
--   <a>CompositenessProof</a> except the composite is through this type.
data CompositenessArgument

-- | <tt>compo == firstDiv*secondDiv</tt>, where all are <tt>&gt; 1</tt>
Divisors :: Integer -> CompositenessArgument
[compo, firstDivisor, secondDivisor] :: CompositenessArgument -> Integer

-- | <tt>compo</tt> fails the strong Fermat test for <tt>fermatBase</tt>
Fermat :: Integer -> CompositenessArgument
[compo, fermatBase] :: CompositenessArgument -> Integer

-- | <tt>compo</tt> fails the Lucas-Selfridge test
Lucas :: Integer -> CompositenessArgument
[compo] :: CompositenessArgument -> Integer

-- | No particular reason given
Belief :: Integer -> CompositenessArgument
[compo] :: CompositenessArgument -> Integer

-- | An argument for primality of a number (which must be <tt>&gt; 1</tt>).
--   <a>PrimalityProof</a>s translate directly to
--   <tt>PrimalityArguments</tt>, correct arguments can be transformed into
--   proofs. This type allows the manipulation of proofs while maintaining
--   their correctness. The only way to access components of a
--   <a>PrimalityProof</a> except the prime is through this type.
data PrimalityArgument

-- | A suggested Pocklington certificate
Pock :: Integer -> Integer -> [(Integer, Int, Integer, PrimalityArgument)] -> PrimalityArgument
[aprime] :: PrimalityArgument -> Integer
[largeFactor, smallFactor] :: PrimalityArgument -> Integer
[factorList] :: PrimalityArgument -> [(Integer, Int, Integer, PrimalityArgument)]

-- | Primality should be provable by trial division to <tt>alimit</tt>
Division :: Integer -> PrimalityArgument
[aprime, alimit] :: PrimalityArgument -> Integer

-- | <tt>aprime</tt> is said to be obviously prime, that holds for primes
--   <tt>&lt; 30</tt>
Obvious :: Integer -> PrimalityArgument
[aprime] :: PrimalityArgument -> Integer

-- | Primality assumed
Assumption :: Integer -> PrimalityArgument
[aprime] :: PrimalityArgument -> Integer

-- | <tt><a>arguePrimality</a></tt> transforms a proof of primality into an
--   argument for primality.
arguePrimality :: PrimalityProof -> PrimalityArgument

-- | <tt><a>argueCompositeness</a></tt> transforms a proof of compositeness
--   into an argument for compositeness.
argueCompositeness :: CompositenessProof -> CompositenessArgument

-- | <tt><a>verifyPrimalityArgument</a></tt> checks the given argument and
--   constructs a proof from it, if it is valid. For the explicit
--   arguments, this is simple and resonably fast, for an
--   <a>Assumption</a>, the verification uses <a>certify</a> and hence may
--   take a long time.
verifyPrimalityArgument :: PrimalityArgument -> Maybe PrimalityProof

-- | <tt><a>verifyCompositenessArgument</a></tt> checks the given argument
--   and constructs a proof from it, if it is valid. For the explicit
--   arguments, this is simple and resonably fast, for a <a>Belief</a>, the
--   verification uses <a>certify</a> and hence may take a long time.
verifyCompositenessArgument :: CompositenessArgument -> Maybe CompositenessProof

-- | <tt><a>certify</a> n</tt> constructs, for <tt>n &gt; 1</tt>, a proof
--   of either primality or compositeness of <tt>n</tt>. This may take a
--   very long time if the number has no small(ish) prime divisors
certify :: Integer -> Certificate

-- | Check the validity of a <a>Certificate</a>. Since it should be
--   impossible to construct invalid certificates by the public interface,
--   this should never return <a>False</a>.
checkCertificate :: Certificate -> Bool

-- | Check the validity of a <a>CompositenessProof</a>. Since it should be
--   impossible to create invalid proofs by the public interface, this
--   should never return <a>False</a>.
checkCompositenessProof :: CompositenessProof -> Bool

-- | Check the validity of a <a>PrimalityProof</a>. Since it should be
--   impossible to create invalid proofs by the public interface, this
--   should never return <a>False</a>.
checkPrimalityProof :: PrimalityProof -> Bool


-- | Factorisation proving the primality of the found factors.
--   
--   For large numbers, this will be very slow in general. Use only if
--   you're paranoid or must be <i>really</i> sure.
module Math.NumberTheory.Primes.Factorisation.Certified

-- | <tt><a>certifiedFactorisation</a> n</tt> produces the prime
--   factorisation of <tt>n</tt>, proving the primality of the factors, but
--   doesn't report the proofs.
certifiedFactorisation :: Integer -> [(Integer, Int)]

-- | <tt><a>certificateFactorisation</a> n</tt> produces a
--   <a>provenFactorisation</a> with a default bound of <tt>100000</tt>.
certificateFactorisation :: Integer -> [((Integer, Int), PrimalityProof)]

-- | <tt><a>provenFactorisation</a> bound n</tt> constructs a the prime
--   factorisation of <tt>n</tt> (which must be positive) together with
--   proofs of primality of the factors, using trial division up to
--   <tt>bound</tt> (which is arbitrarily replaced by <tt>2000</tt> if the
--   supplied value is smaller) and elliptic curve factorisation for the
--   remaining factors if necessary.
--   
--   Construction of primality proofs can take a <i>very</i> long time, so
--   this will usually be slow (but should be faster than using
--   <a>factorise</a> and proving the primality of the factors from
--   scratch).
provenFactorisation :: Integer -> Integer -> [((Integer, Int), PrimalityProof)]


-- | Various functions related to prime factorisation. Many of these
--   functions use the prime factorisation of an <a>Integer</a>. If several
--   of them are used on the same <a>Integer</a>, it would be inefficient
--   to recalculate the factorisation, hence there are also functions
--   working on the canonical factorisation, these require that the number
--   be positive and in the case of the Carmichael function that the list
--   of prime factors with their multiplicities is ascending.
module Math.NumberTheory.Primes.Factorisation

-- | <tt><a>factorise</a> n</tt> produces the prime factorisation of
--   <tt>n</tt>. <tt><a>factorise</a> 0</tt> is an error and the
--   factorisation of <tt>1</tt> is empty. Uses a <a>StdGen</a> produced in
--   an arbitrary manner from the bit-pattern of <tt>n</tt>.
factorise :: Integer -> [(Integer, Int)]

-- | <tt><a>defaultStdGenFactorisation</a></tt> first strips off all small
--   prime factors and then, if the factorisation is not complete, proceeds
--   to curve factorisation. For negative numbers, a factor of <tt>-1</tt>
--   is included, the factorisation of <tt>1</tt> is empty. Since
--   <tt>0</tt> has no prime factorisation, a zero argument causes an
--   error.
defaultStdGenFactorisation :: StdGen -> Integer -> [(Integer, Int)]

-- | <tt><a>stepFactorisation</a></tt> is like <a>factorise'</a>, except
--   that it doesn't use a pseudo random generator but steps through the
--   curves in order. This strategy turns out to be surprisingly fast, on
--   average it doesn't seem to be slower than the <a>StdGen</a> based
--   variant.
stepFactorisation :: Integer -> [(Integer, Int)]

-- | Like <a>factorise</a>, but without input checking, hence <tt>n &gt;
--   1</tt> is required.
factorise' :: Integer -> [(Integer, Int)]

-- | Like <a>defaultStdGenFactorisation</a>, but without input checking, so
--   <tt>n</tt> must be larger than <tt>1</tt>.
defaultStdGenFactorisation' :: StdGen -> Integer -> [(Integer, Int)]

-- | <tt><a>trialDivisionTo</a> bound n</tt> produces a factorisation of
--   <tt>n</tt> using the primes <tt>&lt;= bound</tt>. If <tt>n</tt> has
--   prime divisors <tt>&gt; bound</tt>, the last entry in the list is the
--   product of all these. If <tt>n &lt;= bound^2</tt>, this is a full
--   factorisation, but very slow if <tt>n</tt> has large prime divisors.
trialDivisionTo :: Integer -> Integer -> [(Integer, Int)]

-- | <tt><a>smallFactors</a> bound n</tt> finds all prime divisors of <tt>n
--   &gt; 1</tt> up to <tt>bound</tt> by trial division and returns the
--   list of these together with their multiplicities, and a possible
--   remaining factor which may be composite.
smallFactors :: Integer -> Integer -> ([(Integer, Int)], Maybe Integer)

-- | A wrapper around <a>curveFactorisation</a> providing a few default
--   arguments. The primality test is <a>bailliePSW</a>, the <tt>prng</tt>
--   function - naturally - <a>randomR</a>. This function also requires
--   small prime factors to have been stripped before.
stdGenFactorisation :: Maybe Integer -> StdGen -> Maybe Int -> Integer -> [(Integer, Int)]

-- | <a>curveFactorisation</a> is the driver for the factorisation. Its
--   performance (and success) can be influenced by passing appropriate
--   arguments. If you know that <tt>n</tt> has no prime divisors below
--   <tt>b</tt>, any divisor found less than <tt>b*b</tt> must be prime,
--   thus giving <tt>Just (b*b)</tt> as the first argument allows skipping
--   the comparatively expensive primality test for those. If <tt>n</tt> is
--   such that all prime divisors must have a specific easy to test for
--   structure, a custom primality test can improve the performance
--   (normally, it will make very little difference, since <tt>n</tt> has
--   not many divisors, and many curves have to be tried to find one). More
--   influence has the pseudo random generator (a function <tt>prng</tt>
--   with <tt>6 &lt;= fst (prng k s) &lt;= k-2</tt> and an initial state
--   for the PRNG) used to generate the curves to try. A lucky choice here
--   can make a huge difference. So, if the default takes too long, try
--   another one; or you can improve your chances for a quick result by
--   running several instances in parallel.
--   
--   <a>curveFactorisation</a> <tt>n</tt> requires that small (&lt; 100000)
--   prime factors of <tt>n</tt> have been stripped before. Otherwise it is
--   likely to cycle forever. When in doubt, use
--   <a>defaultStdGenFactorisation</a>.
--   
--   <a>curveFactorisation</a> is unlikely to succeed if <tt>n</tt> has
--   more than one (really) large prime factor.
curveFactorisation :: forall g. Maybe Integer -> (Integer -> Bool) -> (Integer -> g -> (Integer, g)) -> g -> Maybe Int -> Integer -> [(Integer, Int)]

-- | <tt><a>montgomeryFactorisation</a> n b1 b2 s</tt> tries to find a
--   factor of <tt>n</tt> using the curve and point determined by the seed
--   <tt>s</tt> (<tt>6 &lt;= s &lt; n-1</tt>), multiplying the point by the
--   least common multiple of all numbers <tt>&lt;= b1</tt> and all primes
--   between <tt>b1</tt> and <tt>b2</tt>. The idea is that there's a good
--   chance that the order of the point in the curve over one prime factor
--   divides the multiplier, but the order over another factor doesn't, if
--   <tt>b1</tt> and <tt>b2</tt> are appropriately chosen. If they are too
--   small, none of the orders will probably divide the multiplier, if they
--   are too large, all probably will, so they should be chosen to fit the
--   expected size of the smallest factor.
--   
--   It is assumed that <tt>n</tt> has no small prime factors.
--   
--   The result is maybe a nontrivial divisor of <tt>n</tt>.
montgomeryFactorisation :: KnownNat n => Word -> Word -> Mod n -> Maybe Integer


-- | Number of primes not exceeding <tt>n</tt>, <tt>π(n)</tt>, and
--   <tt>n</tt>-th prime; also fast, but reasonably accurate approximations
--   to these.
module Math.NumberTheory.Primes.Counting

-- | <tt><a>primeCount</a> n == π(n)</tt> is the number of (positive)
--   primes not exceeding <tt>n</tt>.
--   
--   For efficiency, the calculations are done on 64-bit signed integers,
--   therefore <tt>n</tt> must not exceed <a>primeCountMaxArg</a>.
--   
--   Requires <tt><i>O</i>(n^0.5)</tt> space, the time complexity is
--   roughly <tt><i>O</i>(n^0.7)</tt>. For small bounds,
--   <tt><a>primeCount</a> n</tt> simply counts the primes not exceeding
--   <tt>n</tt>, for bounds from <tt>30000</tt> on, Meissel's algorithm is
--   used in the improved form due to D.H. Lehmer, cf.
--   <a>http://en.wikipedia.org/wiki/Prime_counting_function#Algorithms_for_evaluating_.CF.80.28x.29</a>.
primeCount :: Integer -> Integer

-- | Maximal allowed argument of <a>primeCount</a>. Currently 8e18.
primeCountMaxArg :: Integer

-- | <tt><a>nthPrime</a> n</tt> calculates the <tt>n</tt>-th prime.
--   Numbering of primes is <tt>1</tt>-based, so <tt><a>nthPrime</a> 1 ==
--   2</tt>.
--   
--   Requires <tt><i>O</i>((n*log n)^0.5)</tt> space, the time complexity
--   is roughly <tt><i>O</i>((n*log n)^0.7</tt>. The argument must be
--   strictly positive, and must not exceed <a>nthPrimeMaxArg</a>.
nthPrime :: Integer -> Integer

-- | Maximal allowed argument of <a>nthPrime</a>. Currently 1.5e17.
nthPrimeMaxArg :: Integer

-- | <tt><a>approxPrimeCount</a> n</tt> gives an approximation of the
--   number of primes not exceeding <tt>n</tt>. The approximation is fairly
--   good for <tt>n</tt> large enough.
approxPrimeCount :: Integral a => a -> a

-- | Following property holds:
--   
--   <pre>
--   approxPrimeCount n &gt;= primeCount n || n &gt;= approxPrimeCountOverestimateLimit
--   </pre>
approxPrimeCountOverestimateLimit :: Integral a => a

-- | <tt><a>nthPrimeApprox</a> n</tt> gives an approximation to the n-th
--   prime. The approximation is fairly good for <tt>n</tt> large enough.
nthPrimeApprox :: Integer -> Integer

-- | Following property holds:
--   
--   <pre>
--   nthPrimeApprox n &lt;= nthPrime n || n &gt;= nthPrimeApproxUnderestimateLimit
--   </pre>
nthPrimeApproxUnderestimateLimit :: Integer


module Math.NumberTheory.Primes


-- | Modular powers (a. k. a. modular exponentiation).
module Math.NumberTheory.Powers.Modular

-- | <tt>powMod</tt> <tt>b</tt> <tt>e</tt> <tt>m</tt> computes
--   (<tt>b^e</tt>) `mod` <tt>m</tt> in effective way. An exponent
--   <tt>e</tt> must be non-negative, a modulo <tt>m</tt> must be positive.
--   Otherwise the behaviour of <tt>powMod</tt> is undefined.
--   
--   <pre>
--   &gt;&gt;&gt; powMod 2 3 5
--   3
--   
--   &gt;&gt;&gt; powMod 3 12345678901234567890 1001
--   1
--   </pre>
--   
--   See also <a>powMod</a> and <a>powSomeMod</a>.
--   
--   For finite numeric types (<a>Int</a>, <a>Word</a>, etc.) modulo
--   <tt>m</tt> should be such that <tt>m^2</tt> 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
--   <a>powModInt</a> and <a>powModWord</a>.
--   
--   <pre>
--   &gt;&gt;&gt; powMod 3 101 (2^60-1 :: Integer)
--   1018105167100379328 -- correct
--   
--   &gt;&gt;&gt; powMod 3 101 (2^60-1 :: Int64)
--   1115647832265427613 -- incorrect due to overflow
--   
--   &gt;&gt;&gt; powModInt 3 101 (2^60-1 :: Int)
--   1018105167100379328 -- correct
--   </pre>
powMod :: (Integral a, Integral b) => a -> b -> a -> a

-- | Specialised version of <a>powMod</a>, able to handle large moduli
--   correctly.
--   
--   <pre>
--   &gt;&gt;&gt; powModWord 3 101 (2^60-1)
--   1018105167100379328
--   </pre>
powModWord :: Word -> Word -> Word -> Word

-- | Specialised version of <a>powMod</a>, able to handle large moduli
--   correctly.
--   
--   <pre>
--   &gt; powModInt 3 101 (2^60-1)
--   </pre>
--   
--   1018105167100379328
powModInt :: Int -> Int -> Int -> Int


-- | Calculating integer roots, modular powers and related things. This
--   module reexports the most needed functions from the implementation
--   modules. The implementation modules provide some additional functions,
--   in particular some unsafe functions which omit some tests for
--   performance reasons.
module Math.NumberTheory.Powers

-- | Calculate the integer square root of a nonnegative number <tt>n</tt>,
--   that is, the largest integer <tt>r</tt> with <tt>r*r &lt;= n</tt>.
--   Throws an error on negative input.
integerSquareRoot :: Integral a => a -> a

-- | Test whether the argument is a square. After a number is found to be
--   positive, first <a>isPossibleSquare</a> is checked, if it is, the
--   integer square root is calculated.
isSquare :: Integral a => a -> Bool

-- | Returns <a>Nothing</a> if the argument is not a square,
--   <tt><a>Just</a> r</tt> if <tt>r*r == n</tt> and <tt>r &gt;= 0</tt>.
--   Avoids the expensive calculation of the square root if <tt>n</tt> is
--   recognized as a non-square before, prevents repeated calculation of
--   the square root if only the roots of perfect squares are needed.
--   Checks for negativity and <a>isPossibleSquare</a>.
exactSquareRoot :: Integral a => a -> Maybe a

-- | Calculate the integer cube root of an integer <tt>n</tt>, that is the
--   largest integer <tt>r</tt> such that <tt>r^3 &lt;= n</tt>. Note that
--   this is not symmetric about <tt>0</tt>, for example
--   <tt>integerCubeRoot (-2) = (-2)</tt> while <tt>integerCubeRoot 2 =
--   1</tt>.
integerCubeRoot :: Integral a => a -> a

-- | Test whether an integer is a cube.
isCube :: Integral a => a -> Bool

-- | Returns <tt>Nothing</tt> if the argument is not a cube, <tt>Just
--   r</tt> if <tt>n == r^3</tt>.
exactCubeRoot :: Integral a => a -> Maybe a

-- | Calculate the integer fourth root of a nonnegative number, that is,
--   the largest integer <tt>r</tt> with <tt>r^4 &lt;= n</tt>. Throws an
--   error on negaitve input.
integerFourthRoot :: Integral a => a -> a

-- | Test whether an integer is a fourth power. First nonnegativity is
--   checked, then the unchecked test is called.
isFourthPower :: Integral a => a -> Bool

-- | Returns <tt>Nothing</tt> if <tt>n</tt> is not a fourth power, <tt>Just
--   r</tt> if <tt>n == r^4</tt> and <tt>r &gt;= 0</tt>.
exactFourthRoot :: Integral a => a -> Maybe a

-- | Calculate an integer root, <tt><a>integerRoot</a> k n</tt> computes
--   the (floor of) the <tt>k</tt>-th root of <tt>n</tt>, where <tt>k</tt>
--   must be positive. <tt>r = <a>integerRoot</a> k n</tt> means <tt>r^k
--   &lt;= n &lt; (r+1)^k</tt> if that is possible at all. It is impossible
--   if <tt>k</tt> is even and <tt>n &lt; 0</tt>, since then <tt>r^k &gt;=
--   0</tt> for all <tt>r</tt>, then, and if <tt>k &lt;= 0</tt>,
--   <tt><a>integerRoot</a></tt> raises an error. For <tt>k &lt; 5</tt>, a
--   specialised version is called which should be more efficient than the
--   general algorithm. However, it is not guaranteed that the rewrite
--   rules for those fire, so if <tt>k</tt> is known in advance, it is
--   safer to directly call the specialised versions.
integerRoot :: (Integral a, Integral b) => b -> a -> a

-- | <tt><a>isKthPower</a> k n</tt> checks whether <tt>n</tt> is a
--   <tt>k</tt>-th power.
isKthPower :: (Integral a, Integral b) => b -> a -> Bool

-- | <tt><a>exactRoot</a> k n</tt> returns <tt><a>Nothing</a></tt> if
--   <tt>n</tt> is not a <tt>k</tt>-th power, <tt><a>Just</a> r</tt> if
--   <tt>n == r^k</tt>. If <tt>k</tt> is divisible by <tt>4, 3</tt> or
--   <tt>2</tt>, a residue test is performed to avoid the expensive
--   calculation if it can thus be determined that <tt>n</tt> is not a
--   <tt>k</tt>-th power.
exactRoot :: (Integral a, Integral b) => b -> a -> Maybe a

-- | <tt><a>isPerfectPower</a> n</tt> checks whether <tt>n == r^k</tt> for
--   some <tt>k &gt; 1</tt>.
isPerfectPower :: Integral a => a -> Bool

-- | <tt><a>highestPower</a> n</tt> produces the pair <tt>(b,k)</tt> with
--   the largest exponent <tt>k</tt> such that <tt>n == b^k</tt>, except
--   for <tt><a>abs</a> n &lt;= 1</tt>, in which case arbitrarily large
--   exponents exist, and by an arbitrary decision <tt>(n,3)</tt> is
--   returned.
--   
--   First, by trial division with small primes, the range of possible
--   exponents is reduced (if <tt>p^e</tt> exactly divides <tt>n</tt>, then
--   <tt>k</tt> must be a divisor of <tt>e</tt>, if several small primes
--   divide <tt>n</tt>, <tt>k</tt> must divide the greatest common divisor
--   of their exponents, which mostly will be <tt>1</tt>, generally small;
--   if none of the small primes divides <tt>n</tt>, the range of possible
--   exponents is reduced since the base is necessarily large), if that has
--   not yet determined the result, the remaining factor is examined by
--   trying the divisors of the <tt>gcd</tt> of the prime exponents if some
--   have been found, otherwise by trying prime exponents recursively.
highestPower :: Integral a => a -> (a, Int)

-- | <tt>powMod</tt> <tt>b</tt> <tt>e</tt> <tt>m</tt> computes
--   (<tt>b^e</tt>) `mod` <tt>m</tt> in effective way. An exponent
--   <tt>e</tt> must be non-negative, a modulo <tt>m</tt> must be positive.
--   Otherwise the behaviour of <tt>powMod</tt> is undefined.
--   
--   <pre>
--   &gt;&gt;&gt; powMod 2 3 5
--   3
--   
--   &gt;&gt;&gt; powMod 3 12345678901234567890 1001
--   1
--   </pre>
--   
--   See also <a>powMod</a> and <a>powSomeMod</a>.
--   
--   For finite numeric types (<a>Int</a>, <a>Word</a>, etc.) modulo
--   <tt>m</tt> should be such that <tt>m^2</tt> 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
--   <a>powModInt</a> and <a>powModWord</a>.
--   
--   <pre>
--   &gt;&gt;&gt; powMod 3 101 (2^60-1 :: Integer)
--   1018105167100379328 -- correct
--   
--   &gt;&gt;&gt; powMod 3 101 (2^60-1 :: Int64)
--   1115647832265427613 -- incorrect due to overflow
--   
--   &gt;&gt;&gt; powModInt 3 101 (2^60-1 :: Int)
--   1018105167100379328 -- correct
--   </pre>
powMod :: (Integral a, Integral b) => a -> b -> a -> a


-- | Modular square roots.
module Math.NumberTheory.Moduli.Sqrt

-- | List all modular square roots.
--   
--   <pre>
--   &gt;&gt;&gt; :set -XDataKinds
--   
--   &gt;&gt;&gt; sqrtsMod (1 :: Mod 60)
--   &gt; [(1 `modulo` 60),(49 `modulo` 60),(41 `modulo` 60),(29 `modulo` 60),(31 `modulo` 60),(19 `modulo` 60),(11 `modulo` 60),(59 `modulo` 60)]
--   </pre>
sqrtsMod :: KnownNat m => Mod m -> [Mod m]

-- | List all square roots modulo a number, which factorisation is passed
--   as a second argument.
--   
--   <pre>
--   &gt;&gt;&gt; sqrtsModFactorisation 1 (factorise 60)
--   [1,49,41,29,31,19,11,59]
--   </pre>
sqrtsModFactorisation :: Integer -> [(Prime Integer, Word)] -> [Integer]

-- | List all square roots modulo power of a prime.
--   
--   <pre>
--   &gt;&gt;&gt; sqrtsModPrimePower 7 (fromJust (isPrime 3)) 2
--   [4,5]
--   
--   &gt;&gt;&gt; sqrtsModPrimePower 9 (fromJust (isPrime 3)) 3
--   [3,12,21,24,6,15]
--   </pre>
sqrtsModPrimePower :: Integer -> Prime Integer -> Word -> [Integer]

-- | List all square roots by prime modulo.
--   
--   <pre>
--   &gt;&gt;&gt; sqrtsModPrime 1 (fromJust (isPrime 5))
--   [1,4]
--   
--   &gt;&gt;&gt; sqrtsModPrime 0 (fromJust (isPrime 5))
--   [0]
--   
--   &gt;&gt;&gt; sqrtsModPrime 2 (fromJust (isPrime 5))
--   []
--   </pre>
sqrtsModPrime :: Integer -> Prime Integer -> [Integer]

-- | <tt>sqrtModP n prime</tt> calculates a modular square root of
--   <tt>n</tt> modulo <tt>prime</tt> if that exists. The second argument
--   <i>must</i> be a (positive) prime, otherwise the computation may not
--   terminate and if it does, may yield a wrong result. The precondition
--   is <i>not</i> checked.
--   
--   If <tt>prime</tt> is a prime and <tt>n</tt> a quadratic residue modulo
--   <tt>prime</tt>, the result is <tt>Just r</tt> where <tt>r^2 ≡ n (mod
--   prime)</tt>, if <tt>n</tt> is a quadratic nonresidue, the result is
--   <tt>Nothing</tt>.

-- | <i>Deprecated: Use <a>sqrtsModPrime</a> instead</i>
sqrtModP :: Integer -> Integer -> Maybe Integer

-- | <tt>sqrtModPList n prime</tt> computes the list of all square roots of
--   <tt>n</tt> modulo <tt>prime</tt>. <tt>prime</tt> <i>must</i> be a
--   (positive) prime. The precondition is <i>not</i> checked.

-- | <i>Deprecated: Use <a>sqrtsModPrime</a> instead</i>
sqrtModPList :: Integer -> Integer -> [Integer]

-- | <tt>sqrtModP' square prime</tt> finds a square root of <tt>square</tt>
--   modulo prime. <tt>prime</tt> <i>must</i> be a (positive) prime, and
--   <tt>square</tt> <i>must</i> be a positive quadratic residue modulo
--   <tt>prime</tt>, i.e. <tt>'jacobi square prime == 1</tt>. The
--   precondition is <i>not</i> checked.

-- | <i>Deprecated: Use <a>sqrtsModPrime</a> instead</i>
sqrtModP' :: Integer -> Integer -> Integer

-- | <tt>tonelliShanks square prime</tt> calculates a square root of
--   <tt>square</tt> modulo <tt>prime</tt>, where <tt>prime</tt> is a prime
--   of the form <tt>4*k + 1</tt> and <tt>square</tt> is a positive
--   quadratic residue modulo <tt>prime</tt>, using the Tonelli-Shanks
--   algorithm. No checks on the input are performed.

-- | <i>Deprecated: Use <a>sqrtsModPrime</a> instead</i>
tonelliShanks :: Integer -> Integer -> Integer

-- | <tt>sqrtModPP n (prime,expo)</tt> calculates a square root of
--   <tt>n</tt> modulo <tt>prime^expo</tt> if one exists. <tt>prime</tt>
--   <i>must</i> be a (positive) prime. <tt>expo</tt> must be positive,
--   <tt>n</tt> must be coprime to <tt>prime</tt>

-- | <i>Deprecated: Use <a>sqrtsModPrimePower</a> instead</i>
sqrtModPP :: Integer -> (Integer, Int) -> Maybe Integer

-- | <tt>sqrtModPPList n (prime,expo)</tt> calculates the list of all
--   square roots of <tt>n</tt> modulo <tt>prime^expo</tt>. The same
--   restriction as in <a>sqrtModPP</a> applies to the arguments.

-- | <i>Deprecated: Use <a>sqrtsModPrimePower</a> instead</i>
sqrtModPPList :: Integer -> (Integer, Int) -> [Integer]

-- | <tt>sqrtModF n primePowers</tt> calculates a square root of <tt>n</tt>
--   modulo <tt>product [p^k | (p,k) &lt;- primePowers]</tt> if one exists
--   and all primes are distinct. The list must be non-empty, <tt>n</tt>
--   must be coprime with all primes.

-- | <i>Deprecated: Use <a>sqrtsModFactorisation</a> or <a>sqrtsMod</a>
--   instead</i>
sqrtModF :: Integer -> [(Integer, Int)] -> Maybe Integer

-- | <tt>sqrtModFList n primePowers</tt> calculates all square roots of
--   <tt>n</tt> modulo <tt>product [p^k | (p,k) &lt;- primePowers]</tt> if
--   all primes are distinct. The list must be non-empty, <tt>n</tt> must
--   be coprime with all primes.

-- | <i>Deprecated: Use <a>sqrtsModFactorisation</a> or <a>sqrtsMod</a>
--   instead</i>
sqrtModFList :: Integer -> [(Integer, Int)] -> [Integer]


-- | This module exports functions for manipulating Gaussian integers,
--   including computing their prime factorisations.
module Math.NumberTheory.Quadratic.GaussianIntegers

-- | A Gaussian integer is a+bi, where a and b are both integers.
data GaussianInteger
(:+) :: !Integer -> !Integer -> GaussianInteger
[real] :: GaussianInteger -> !Integer
[imag] :: GaussianInteger -> !Integer
infix 6 :+

-- | The imaginary unit, where
--   
--   <pre>
--   ι .^ 2 == -1
--   </pre>
ι :: GaussianInteger

-- | Conjugate a Gaussian integer.
conjugate :: GaussianInteger -> GaussianInteger

-- | The square of the magnitude of a Gaussian integer.
norm :: GaussianInteger -> Integer

-- | Raise a Gaussian integer to a given power.

-- | <i>Deprecated: Use (^) instead.</i>
(.^) :: Integral a => GaussianInteger -> a -> GaussianInteger
infixr 8 .^

-- | Compute whether a given Gaussian integer is prime.
isPrime :: GaussianInteger -> Bool

-- | An infinite list of the Gaussian primes. Uses primes in Z to
--   exhaustively generate all Gaussian primes (up to associates), in order
--   of ascending magnitude.
primes :: [GaussianInteger]

-- | Compute the GCD of two Gaussian integers. Result is always in the
--   first quadrant.

-- | <i>Deprecated: Use <a>gcd</a> instead.</i>
gcdG :: GaussianInteger -> GaussianInteger -> GaussianInteger

-- | <i>Deprecated: Use <a>gcd</a> instead.</i>
gcdG' :: GaussianInteger -> GaussianInteger -> GaussianInteger

-- | Find a Gaussian integer whose norm is the given prime number of form
--   4k + 1 using <a>Hermite-Serret algorithm</a>.
findPrime :: Integer -> GaussianInteger

-- | <i>Deprecated: Use <a>findPrime</a> instead.</i>
findPrime' :: Integer -> GaussianInteger

-- | Compute the prime factorisation of a Gaussian integer. This is unique
--   up to units (+<i>- 1, +</i>- i). Unit factors are not included in the
--   result.
factorise :: GaussianInteger -> [(GaussianInteger, Int)]
instance GHC.Generics.Generic Math.NumberTheory.Quadratic.GaussianIntegers.GaussianInteger
instance GHC.Classes.Ord Math.NumberTheory.Quadratic.GaussianIntegers.GaussianInteger
instance GHC.Classes.Eq Math.NumberTheory.Quadratic.GaussianIntegers.GaussianInteger
instance Control.DeepSeq.NFData Math.NumberTheory.Quadratic.GaussianIntegers.GaussianInteger
instance GHC.Show.Show Math.NumberTheory.Quadratic.GaussianIntegers.GaussianInteger
instance GHC.Num.Num Math.NumberTheory.Quadratic.GaussianIntegers.GaussianInteger
instance Math.NumberTheory.Euclidean.Euclidean Math.NumberTheory.Quadratic.GaussianIntegers.GaussianInteger


-- | This module exports functions for manipulating Gaussian integers,
--   including computing their prime factorisations.

-- | <i>Deprecated: Use Math.NumberTheory.Quadratic.GaussianIntegers
--   instead</i>
module Math.NumberTheory.GaussianIntegers


-- | This module exports functions for manipulating Eisenstein integers,
--   including computing their prime factorisations.
module Math.NumberTheory.Quadratic.EisensteinIntegers

-- | An Eisenstein integer is a + bω, where a and b are both integers.
data EisensteinInteger
(:+) :: !Integer -> !Integer -> EisensteinInteger
[real] :: EisensteinInteger -> !Integer
[imag] :: EisensteinInteger -> !Integer
infix 6 :+

-- | The imaginary unit for Eisenstein integers, where
--   
--   <pre>
--   ω == (-1/2) + ((sqrt 3)/2)ι == exp(2*pi*ι/3)
--   </pre>
--   
--   and ι is the usual imaginary unit with ι² == -1.
ω :: EisensteinInteger

-- | Conjugate a Eisenstein integer.
conjugate :: EisensteinInteger -> EisensteinInteger

-- | The square of the magnitude of a Eisenstein integer.
norm :: EisensteinInteger -> Integer

-- | Produce a list of an <tt>EisensteinInteger</tt>'s associates.
associates :: EisensteinInteger -> [EisensteinInteger]

-- | List of all Eisenstein units, counterclockwise across all sextants,
--   starting with <tt>1</tt>.
ids :: [EisensteinInteger]

-- | Remove <tt>1 - ω</tt> factors from an <tt>EisensteinInteger</tt>, and
--   calculate that prime's multiplicity in the number's factorisation.
divideByThree :: EisensteinInteger -> (Int, EisensteinInteger)

-- | Compute the prime factorisation of a Eisenstein integer. This is
--   unique up to units (+<i>- 1, +</i>- ω, +/- ω²). * Unit factors are not
--   included in the result. * All prime factors are primary i.e. <tt>e ≡ 2
--   (modE 3)</tt>, for an Eisenstein prime factor <tt>e</tt>.
--   
--   <ul>
--   <li>This function works by factorising the norm of an Eisenstein
--   integer and then, for each prime factor, finding the Eisenstein prime
--   whose norm is said prime factor with <tt>findPrime</tt>.</li>
--   <li>This is only possible because the norm function of the Euclidean
--   Domain of Eisenstein integers is multiplicative: <tt>norm (e1 * e2) ==
--   norm e1 * norm e2</tt> for any two <tt>EisensteinInteger</tt>s <tt>e1,
--   e2</tt>.</li>
--   <li>In the previously mentioned work <a>Bandara, Sarada, "An
--   Exposition of the Eisenstein Integers" (2016)</a>, in Theorem 8.4 in
--   Chapter 8, a way is given to express any Eisenstein integer <tt>μ</tt>
--   as <tt>(-1)^a * ω^b * (1 - ω)^c * product [π_i^a_i | i &lt;-
--   [1..N]]</tt> where <tt>a, b, c, a_i</tt> are nonnegative integers,
--   <tt>N &gt; 1</tt> is an integer and <tt>π_i</tt> are primary primes
--   (for a primary Eisenstein prime <tt>p</tt>, <tt>p ≡ 2 (modE 3)</tt>,
--   see <tt>primary</tt> above).</li>
--   <li>Aplying <tt>norm</tt> to both sides of Theorem 8.4: <tt>norm μ =
--   norm ((-1)^a * ω^b * (1 - ω)^c * product [ π_i^a_i | i &lt;-
--   [1..N]])</tt> == <tt>norm μ = norm ((-1)^a) * norm (ω^b) * norm ((1 -
--   ω)^c) * norm (product [ π_i^a_i | i &lt;- [1..N]])</tt> == <tt>norm μ
--   = (norm (-1))^a * (norm ω)^b * (norm (1 - ω))^c * product [ norm
--   (π_i^a_i) | i &lt;- [1..N]]</tt> == <tt>norm μ = (norm (-1))^a * (norm
--   ω)^b * (norm (1 - ω))^c * product [ (norm π_i)^a_i) | i &lt;-
--   [1..N]]</tt> == <tt>norm μ = 1^a * 1^b * 3^c * product [ (norm
--   π_i)^a_i) | i &lt;- [1..N]]</tt> == <tt>norm μ = 3^c * product [ (norm
--   π_i)^a_i) | i &lt;- [1..N]]</tt> where <tt>a, b, c, a_i</tt> are
--   nonnegative integers, and <tt>N &gt; 1</tt> is an integer.</li>
--   <li>The remainder of the Eisenstein integer factorisation problem is
--   about finding appropriate <tt>[e_i | i &lt;- [1..M]</tt> such that
--   <tt>(nub . map norm) [e_i | i &lt;- [1..N]] == [π_i | i &lt;-
--   [1..N]]</tt> where <tt> 1 &lt; N &lt;= M</tt> are integers,
--   <tt>nub</tt> removes duplicates and <tt>==</tt> is equality on
--   sets.</li>
--   <li>The reason <tt>M &gt;= N</tt> is because the prime factors of an
--   Eisenstein integer may include a prime factor and its conjugate,
--   meaning the number may have more Eisenstein prime factors than its
--   norm has integer prime factors.</li>
--   </ul>
factorise :: EisensteinInteger -> [(EisensteinInteger, Int)]

-- | Find an Eisenstein integer whose norm is the given prime number in the
--   form <tt>3k + 1</tt> using a modification of the <a>Hermite-Serret
--   algorithm</a>.
--   
--   The maintainer <a>Andrew Lelechenko</a> derived the following: * Each
--   prime of form <tt>3n+1</tt> is actually of form <tt>6k+1</tt>. * One
--   has <tt>(z+3k)^2 ≡ z^2 + 6kz + 9k^2 ≡ z^2 + (6k+1)z - z + 9k^2 ≡ z^2 -
--   z + 9k^2 (mod 6k+1)</tt>.
--   
--   <ul>
--   <li>The goal is to solve <tt>z^2 - z + 1 ≡ 0 (mod 6k+1)</tt>. One has:
--   <tt>z^2 - z + 9k^2 ≡ 9k^2 - 1 (mod 6k+1)</tt> <tt>(z+3k)^2 ≡ 9k^2-1
--   (mod 6k+1)</tt> <tt>z+3k = sqrtMod(9k^2-1)</tt> <tt>z =
--   sqrtMod(9k^2-1) - 3k</tt></li>
--   <li>For example, let <tt>p = 7</tt>, then <tt>k = 1</tt>. Square root
--   of <tt>9*1^2-1 modulo 7</tt> is <tt>1</tt>.</li>
--   <li>And <tt>z = 1 - 3*1 = -2 ≡ 5 (mod 7)</tt>.</li>
--   <li>Truly, <tt>norm (5 :+ 1) = 25 - 5 + 1 = 21 ≡ 0 (mod 7)</tt>.</li>
--   </ul>
findPrime :: Integer -> EisensteinInteger

-- | Checks if a given <tt>EisensteinInteger</tt> is prime.
--   <tt>EisensteinInteger</tt>s whose norm is a prime congruent to
--   <tt>0</tt> or <tt>1</tt> modulo 3 are prime. See <a>Bandara, Sarada,
--   "An Exposition of the Eisenstein Integers" (2016)</a>, page 12.
isPrime :: EisensteinInteger -> Bool

-- | An infinite list of Eisenstein primes. Uses primes in Z to
--   exhaustively generate all Eisenstein primes in order of ascending
--   magnitude. * Every prime is in the first sextant, so the list contains
--   no associates. * Eisenstein primes from the whole complex plane can be
--   generated by applying <tt>associates</tt> to each prime in this list.
primes :: [EisensteinInteger]
instance GHC.Generics.Generic Math.NumberTheory.Quadratic.EisensteinIntegers.EisensteinInteger
instance GHC.Classes.Ord Math.NumberTheory.Quadratic.EisensteinIntegers.EisensteinInteger
instance GHC.Classes.Eq Math.NumberTheory.Quadratic.EisensteinIntegers.EisensteinInteger
instance GHC.Show.Show Math.NumberTheory.Quadratic.EisensteinIntegers.EisensteinInteger
instance GHC.Num.Num Math.NumberTheory.Quadratic.EisensteinIntegers.EisensteinInteger
instance Math.NumberTheory.Euclidean.Euclidean Math.NumberTheory.Quadratic.EisensteinIntegers.EisensteinInteger


-- | An abstract type class for unique factorisation domains.
module Math.NumberTheory.UniqueFactorisation

-- | Type of primes of a given unique factorisation domain.
--   
--   <tt>abs (unPrime n) == unPrime n</tt> must hold for all <tt>n</tt> of
--   type <tt>Prime t</tt>
type family Prime (f :: *) :: *

-- | The following invariant must hold for <tt>n /= 0</tt>:
--   
--   <pre>
--   abs n == abs (product (map (\(p, k) -&gt; unPrime p ^ k) (factorise n)))
--   </pre>
--   
--   The result of <a>factorise</a> should not contain zero powers and
--   should not change after multiplication of the argument by domain's
--   unit.
class UniqueFactorisation a
unPrime :: UniqueFactorisation a => Prime a -> a
factorise :: UniqueFactorisation a => a -> [(Prime a, Word)]
isPrime :: UniqueFactorisation a => a -> Maybe (Prime a)
isPrime :: (UniqueFactorisation a, Eq a, Num a) => a -> Maybe (Prime a)
instance GHC.Show.Show Math.NumberTheory.UniqueFactorisation.EisensteinPrime
instance GHC.Classes.Eq Math.NumberTheory.UniqueFactorisation.EisensteinPrime
instance GHC.Show.Show Math.NumberTheory.UniqueFactorisation.GaussianPrime
instance GHC.Classes.Eq Math.NumberTheory.UniqueFactorisation.GaussianPrime
instance Math.NumberTheory.UniqueFactorisation.UniqueFactorisation GHC.Types.Int
instance Math.NumberTheory.UniqueFactorisation.UniqueFactorisation GHC.Types.Word
instance Math.NumberTheory.UniqueFactorisation.UniqueFactorisation GHC.Integer.Type.Integer
instance Math.NumberTheory.UniqueFactorisation.UniqueFactorisation GHC.Natural.Natural
instance Math.NumberTheory.UniqueFactorisation.UniqueFactorisation Math.NumberTheory.Quadratic.GaussianIntegers.GaussianInteger
instance Math.NumberTheory.UniqueFactorisation.UniqueFactorisation Math.NumberTheory.Quadratic.EisensteinIntegers.EisensteinInteger


-- | Type for numbers, accompanied by their factorisation.
module Math.NumberTheory.Prefactored

-- | A container for a number and its pairwise coprime (but not
--   neccessarily prime) factorisation. It is designed to preserve
--   information about factors under multiplication. One can use this
--   representation to speed up prime factorisation and computation of
--   arithmetic functions.
--   
--   For instance, let <tt>p</tt> and <tt>q</tt> be big primes:
--   
--   <pre>
--   &gt;&gt;&gt; let p, q :: Integer
--   
--   &gt;&gt;&gt; p = 1000000000000000000000000000057
--   
--   &gt;&gt;&gt; q = 2000000000000000000000000000071
--   </pre>
--   
--   It would be difficult to compute prime factorisation of their product
--   as is: <a>factorise</a> would take ages. Things become different if we
--   simply change types of <tt>p</tt> and <tt>q</tt> to prefactored ones:
--   
--   <pre>
--   &gt;&gt;&gt; let p, q :: Prefactored Integer
--   
--   &gt;&gt;&gt; p = 1000000000000000000000000000057
--   
--   &gt;&gt;&gt; q = 2000000000000000000000000000071
--   </pre>
--   
--   Now prime factorisation is done instantly:
--   
--   <pre>
--   &gt;&gt;&gt; factorise (p * q)
--   [(PrimeNat 1000000000000000000000000000057, 1), (PrimeNat 2000000000000000000000000000071, 1)]
--   
--   &gt;&gt;&gt; factorise (p^2 * q^3)
--   [(PrimeNat 1000000000000000000000000000057, 2), (PrimeNat 2000000000000000000000000000071, 3)]
--   </pre>
--   
--   Moreover, we can instantly compute <tt>totient</tt> and its
--   iterations. It works fine, because output of <tt>totient</tt> is also
--   prefactored.
--   
--   <pre>
--   &gt;&gt;&gt; prefValue $ totient (p^2 * q^3)
--   8000000000000000000000000001752000000000000000000000000151322000000000000000000000006445392000000000000000000000135513014000000000000000000001126361040
--   
--   &gt;&gt;&gt; prefValue $ totient $ totient (p^2 * q^3)
--   2133305798262843681544648472180210822742702690942899511234131900112583590230336435053688694839034890779375223070157301188739881477320529552945446912000
--   </pre>
--   
--   Let us look under the hood:
--   
--   <pre>
--   &gt;&gt;&gt; prefFactors $ totient (p^2 * q^3)
--   Coprimes {unCoprimes = fromList [(2,4),(3,3),
--     (41666666666666666666666666669,1),(111111111111111111111111111115,1),
--     (1000000000000000000000000000057,1),(2000000000000000000000000000071,2)]}
--   
--   &gt;&gt;&gt; prefFactors $ totient $ totient (p^2 * q^3)
--   Coprimes {unCoprimes = fromList [(2,22),(3,8),(5,3),(39521,1),(199937,1),(6046667,1),
--     (227098769,1),(361696272343,1),(85331809838489,1),(22222222222222222222222222223,1),
--     (41666666666666666666666666669,1),(2000000000000000000000000000071,1)]}
--   </pre>
--   
--   Pairwise coprimality of factors is crucial, because it allows us to
--   process them independently, possibly even in parallel or concurrent
--   fashion.
--   
--   Following invariant is guaranteed to hold:
--   
--   <pre>
--   abs (prefValue x) = abs (product (map (uncurry (^)) (prefFactors x)))
--   </pre>
data Prefactored a

-- | Create <a>Prefactored</a> from a given number.
--   
--   <pre>
--   &gt;&gt;&gt; fromValue 123
--   Prefactored {prefValue = 123, prefFactors = Coprimes {unCoprimes = fromList [(123,1)]}}
--   </pre>
fromValue :: (Eq a, Num a) => a -> Prefactored a

-- | Create <a>Prefactored</a> from a given list of pairwise coprime (but
--   not neccesarily prime) factors with multiplicities.
--   
--   <pre>
--   &gt;&gt;&gt; fromFactors (splitIntoCoprimes [(140, 1), (165, 1)])
--   Prefactored {prefValue = 23100, prefFactors = Coprimes {unCoprimes = fromList [(5,2),(28,1),(33,1)]}}
--   
--   &gt;&gt;&gt; fromFactors (splitIntoCoprimes [(140, 2), (165, 3)])
--   Prefactored {prefValue = 88045650000, prefFactors = Coprimes {unCoprimes = fromList [(5,5),(28,2),(33,3)]}}
--   </pre>
fromFactors :: Num a => Coprimes a Word -> Prefactored a
instance GHC.Show.Show a => GHC.Show.Show (Math.NumberTheory.Prefactored.Prefactored a)
instance (Math.NumberTheory.Euclidean.Euclidean a, GHC.Classes.Ord a) => GHC.Num.Num (Math.NumberTheory.Prefactored.Prefactored a)
instance (GHC.Classes.Eq a, GHC.Num.Num a, Math.NumberTheory.UniqueFactorisation.UniqueFactorisation a) => Math.NumberTheory.UniqueFactorisation.UniqueFactorisation (Math.NumberTheory.Prefactored.Prefactored a)


-- | Polynomial modular equations.
module Math.NumberTheory.Moduli.Equations

-- | Find all solutions of ax + b ≡ 0 (mod m).
--   
--   <pre>
--   &gt;&gt;&gt; :set -XDataKinds
--   
--   &gt;&gt;&gt; solveLinear (6 :: Mod 10) 4 -- solving 6x + 4 ≡ 0 (mod 10)
--   [(1 `modulo` 10),(6 `modulo` 10)]
--   </pre>
solveLinear :: KnownNat m => Mod m -> Mod m -> [Mod m]

-- | Find all solutions of ax² + bx + c ≡ 0 (mod m).
--   
--   <pre>
--   &gt;&gt;&gt; :set -XDataKinds
--   
--   &gt;&gt;&gt; solveQuadratic (1 :: Mod 32) 0 (-17) -- solving x² - 17 ≡ 0 (mod 32)
--   [(9 `modulo` 32),(25 `modulo` 32),(7 `modulo` 32),(23 `modulo` 32)]
--   </pre>
solveQuadratic :: KnownNat m => Mod m -> Mod m -> Mod m -> [Mod m]


-- | Values of <a>Möbius function</a>.
module Math.NumberTheory.ArithmeticFunctions.Moebius

-- | Represents three possible values of <a>Möbius function</a>.
data Moebius

-- | −1
MoebiusN :: Moebius

-- | 0
MoebiusZ :: Moebius

-- | 1
MoebiusP :: Moebius

-- | Convert to any numeric type.
runMoebius :: Num a => Moebius -> a

-- | Evaluate the Möbius function over a block. Value of <tt>f</tt> at 0,
--   if zero falls into block, is undefined.
--   
--   Based on the sieving algorithm from p. 3 of <a>Computations of the
--   Mertens function and improved bounds on the Mertens conjecture</a> by
--   G. Hurst. It is approximately 5x faster than <a>sieveBlockUnboxed</a>.
--   
--   <pre>
--   &gt;&gt;&gt; sieveBlockMoebius 1 10
--   [MoebiusP, MoebiusN, MoebiusN, MoebiusZ, MoebiusN, MoebiusP, MoebiusN, MoebiusZ, MoebiusZ, MoebiusP]
--   </pre>
sieveBlockMoebius :: Word -> Word -> Vector Moebius
instance GHC.Show.Show Math.NumberTheory.ArithmeticFunctions.Moebius.Moebius
instance GHC.Classes.Ord Math.NumberTheory.ArithmeticFunctions.Moebius.Moebius
instance GHC.Classes.Eq Math.NumberTheory.ArithmeticFunctions.Moebius.Moebius
instance Data.Vector.Unboxed.Base.Unbox Math.NumberTheory.ArithmeticFunctions.Moebius.Moebius
instance Data.Vector.Generic.Mutable.Base.MVector Data.Vector.Unboxed.Base.MVector Math.NumberTheory.ArithmeticFunctions.Moebius.Moebius
instance Data.Vector.Generic.Base.Vector Data.Vector.Unboxed.Base.Vector Math.NumberTheory.ArithmeticFunctions.Moebius.Moebius
instance GHC.Base.Semigroup Math.NumberTheory.ArithmeticFunctions.Moebius.Moebius
instance GHC.Base.Monoid Math.NumberTheory.ArithmeticFunctions.Moebius.Moebius


-- | This module provides an interface for defining and manipulating
--   arithmetic functions. It also defines several most widespreaded
--   arithmetic functions.
module Math.NumberTheory.ArithmeticFunctions

-- | A typical arithmetic function operates on the canonical factorisation
--   of a number into prime's powers and consists of two rules. The first
--   one determines the values of the function on the powers of primes. The
--   second one determines how to combine these values into final result.
--   
--   In the following definition the first argument is the function on
--   prime's powers, the monoid instance determines a rule of combination
--   (typically <tt>Product</tt> or <tt>Sum</tt>), and the second argument
--   is convenient for unwrapping (typically, <tt>getProduct</tt> or
--   <tt>getSum</tt>).
data ArithmeticFunction n a
[ArithmeticFunction] :: Monoid m => (Prime n -> Word -> m) -> (m -> a) -> ArithmeticFunction n a

-- | Convert to function. The value on 0 is undefined.
runFunction :: UniqueFactorisation n => ArithmeticFunction n a -> n -> a

-- | Create a multiplicative function from the function on prime's powers.
--   See examples below.
multiplicative :: Num a => (Prime n -> Word -> a) -> ArithmeticFunction n a
divisors :: (UniqueFactorisation n, Num n, Ord n) => n -> Set n

-- | The set of all (positive) divisors of an argument.
divisorsA :: forall n. (UniqueFactorisation n, Num n, Ord n) => ArithmeticFunction n (Set n)
divisorsList :: (UniqueFactorisation n, Num n) => n -> [n]

-- | The unsorted list of all (positive) divisors of an argument, produced
--   in lazy fashion.
divisorsListA :: forall n. (UniqueFactorisation n, Num n) => ArithmeticFunction n [n]
divisorsSmall :: (UniqueFactorisation n, Prime n ~ Prime Int) => n -> IntSet

-- | Same as <a>divisors</a>, but with better performance on cost of type
--   restriction.
divisorsSmallA :: forall n. Prime n ~ Prime Int => ArithmeticFunction n IntSet
tau :: (UniqueFactorisation n, Num a) => n -> a

-- | The number of (positive) divisors of an argument.
--   
--   <pre>
--   tauA = multiplicative (\_ k -&gt; k + 1)
--   </pre>
tauA :: Num a => ArithmeticFunction n a
sigma :: (UniqueFactorisation n, Integral n) => Word -> n -> n

-- | The sum of the <tt>k</tt>-th powers of (positive) divisors of an
--   argument.
--   
--   <pre>
--   sigmaA = multiplicative (\p k -&gt; sum $ map (p ^) [0..k])
--   sigmaA 0 = tauA
--   </pre>
sigmaA :: forall n. (UniqueFactorisation n, Integral n) => Word -> ArithmeticFunction n n
totient :: (UniqueFactorisation n, Num n) => n -> n

-- | Calculates the totient of a positive number <tt>n</tt>, i.e. the
--   number of <tt>k</tt> with <tt>1 &lt;= k &lt;= n</tt> and
--   <tt><a>gcd</a> n k == 1</tt>, in other words, the order of the group
--   of units in <tt>ℤ/(n)</tt>.
totientA :: forall n. (UniqueFactorisation n, Num n) => ArithmeticFunction n n
jordan :: (UniqueFactorisation n, Num n) => Word -> n -> n

-- | Calculates the k-th Jordan function of an argument.
--   
--   <pre>
--   jordanA 1 = totientA
--   </pre>
jordanA :: forall n. (UniqueFactorisation n, Num n) => Word -> ArithmeticFunction n n
ramanujan :: Integer -> Integer

-- | Calculates the <a>Ramanujan tau function</a> of a positive number
--   <tt>n</tt>, using formulas given <a>here</a>
ramanujanA :: ArithmeticFunction Integer Integer
moebius :: UniqueFactorisation n => n -> Moebius

-- | Calculates the Möbius function of an argument.
moebiusA :: ArithmeticFunction n Moebius

-- | Represents three possible values of <a>Möbius function</a>.
data Moebius

-- | −1
MoebiusN :: Moebius

-- | 0
MoebiusZ :: Moebius

-- | 1
MoebiusP :: Moebius

-- | Convert to any numeric type.
runMoebius :: Num a => Moebius -> a
liouville :: (UniqueFactorisation n, Num a) => n -> a

-- | Calculates the Liouville function of an argument.
liouvilleA :: Num a => ArithmeticFunction n a

-- | Create an additive function from the function on prime's powers. See
--   examples below.
additive :: Num a => (Prime n -> Word -> a) -> ArithmeticFunction n a
smallOmega :: (UniqueFactorisation n, Num a) => n -> a

-- | Number of distinct prime factors.
--   
--   <pre>
--   smallOmegaA = additive (\_ _ -&gt; 1)
--   </pre>
smallOmegaA :: Num a => ArithmeticFunction n a
bigOmega :: UniqueFactorisation n => n -> Word

-- | Number of prime factors, counted with multiplicity.
--   
--   <pre>
--   bigOmegaA = additive (\_ k -&gt; k)
--   </pre>
bigOmegaA :: ArithmeticFunction n Word
carmichael :: (UniqueFactorisation n, Integral n) => n -> n

-- | Calculates the Carmichael function for a positive integer, that is,
--   the (smallest) exponent of the group of units in <tt>ℤ/(n)</tt>.
carmichaelA :: forall n. (UniqueFactorisation n, Integral n) => ArithmeticFunction n n
expMangoldt :: (UniqueFactorisation n, Num n) => n -> n

-- | The exponent of von Mangoldt function. Use <tt>log expMangoldtA</tt>
--   to recover von Mangoldt function itself.
expMangoldtA :: forall n. (UniqueFactorisation n, Num n) => ArithmeticFunction n n


-- | Primitive roots and cyclic groups.
module Math.NumberTheory.Moduli.PrimitiveRoot

-- | A multiplicative group of residues is called cyclic, if there is a
--   primitive root <tt>g</tt>, whose powers generates all elements. Any
--   cyclic multiplicative group of residues falls into one of the
--   following cases.
data CyclicGroup a

-- | Residues modulo 2.
CG2 :: CyclicGroup a

-- | Residues modulo 4.
CG4 :: CyclicGroup a

-- | Residues modulo <tt>p</tt>^<tt>k</tt> for <b>odd</b> prime <tt>p</tt>.
CGOddPrimePower :: Prime a -> Word -> CyclicGroup a

-- | Residues modulo 2<tt>p</tt>^<tt>k</tt> for <b>odd</b> prime
--   <tt>p</tt>.
CGDoubleOddPrimePower :: Prime a -> Word -> CyclicGroup a

-- | Check whether a multiplicative group of residues, characterized by its
--   modulo, is cyclic and, if yes, return its form.
--   
--   <pre>
--   &gt;&gt;&gt; cyclicGroupFromModulo 4
--   Just CG4
--   
--   &gt;&gt;&gt; cyclicGroupFromModulo (2 * 13 ^ 3)
--   Just (CGDoubleOddPrimePower (PrimeNat 13) 3)
--   
--   &gt;&gt;&gt; cyclicGroupFromModulo (4 * 13)
--   Nothing
--   </pre>
cyclicGroupFromModulo :: (Ord a, Integral a, UniqueFactorisation a) => a -> Maybe (CyclicGroup a)

-- | Extract modulo and its factorisation from a cyclic multiplicative
--   group of residues.
--   
--   <pre>
--   &gt;&gt;&gt; cyclicGroupToModulo CG4
--   Prefactored {prefValue = 4, prefFactors = Coprimes {unCoprimes = fromList [(2,2)]}}
--   </pre>
--   
--   <pre>
--   &gt;&gt;&gt; :set -XTypeFamilies
--   
--   &gt;&gt;&gt; cyclicGroupToModulo (CGDoubleOddPrimePower (PrimeNat 13) 3)
--   Prefactored {prefValue = 4394, prefFactors = Coprimes {unCoprimes = fromList [(2,1),(13,3)]}}
--   </pre>
cyclicGroupToModulo :: (Euclidean a, Ord a, UniqueFactorisation a) => CyclicGroup a -> Prefactored a

-- | Calculate the size of a given cyclic group.
groupSize :: (Euclidean a, Ord a, UniqueFactorisation a) => CyclicGroup a -> Prefactored a

-- | 'PrimitiveRoot m' is a type which is only inhabited by primitive roots
--   of n.
data PrimitiveRoot m

-- | Extract primitive root value.
unPrimitiveRoot :: PrimitiveRoot m -> MultMod m

-- | Get cyclic group structure.
getGroup :: PrimitiveRoot m -> CyclicGroup Natural

-- | Check whether a given modular residue is a <a>primitive root</a>.
--   
--   <pre>
--   &gt;&gt;&gt; :set -XDataKinds
--   
--   &gt;&gt;&gt; isPrimitiveRoot (1 :: Mod 13)
--   False
--   
--   &gt;&gt;&gt; isPrimitiveRoot (2 :: Mod 13)
--   True
--   </pre>
--   
--   Here is how to list all primitive roots:
--   
--   <pre>
--   &gt;&gt;&gt; mapMaybe isPrimitiveRoot [minBound .. maxBound] :: [Mod 13]
--   [(2 `modulo` 13), (6 `modulo` 13), (7 `modulo` 13), (11 `modulo` 13)]
--   </pre>
--   
--   This function is a convenient wrapper around <a>isPrimitiveRoot'</a>.
--   The latter provides better control and performance, if you need them.
isPrimitiveRoot :: KnownNat n => Mod n -> Maybe (PrimitiveRoot n)

-- | <a>isPrimitiveRoot'</a> <tt>cg</tt> <tt>a</tt> checks whether
--   <tt>a</tt> is a <a>primitive root</a> of a given cyclic multiplicative
--   group of residues <tt>cg</tt>.
--   
--   <pre>
--   &gt;&gt;&gt; let Just cg = cyclicGroupFromModulo 13
--   
--   &gt;&gt;&gt; isPrimitiveRoot' cg 1
--   False
--   
--   &gt;&gt;&gt; isPrimitiveRoot' cg 2
--   True
--   </pre>
isPrimitiveRoot' :: (Integral a, UniqueFactorisation a) => CyclicGroup a -> a -> Bool
instance GHC.TypeNats.KnownNat m => GHC.Show.Show (Math.NumberTheory.Moduli.PrimitiveRoot.PrimitiveRoot m)
instance GHC.Classes.Eq (Math.NumberTheory.Moduli.PrimitiveRoot.PrimitiveRoot m)
instance GHC.Generics.Generic (Math.NumberTheory.Moduli.PrimitiveRoot.CyclicGroup a)
instance GHC.Classes.Eq (Math.NumberTheory.Primes.Types.Prime a) => GHC.Classes.Eq (Math.NumberTheory.Moduli.PrimitiveRoot.CyclicGroup a)
instance GHC.Show.Show (Math.NumberTheory.Primes.Types.Prime a) => GHC.Show.Show (Math.NumberTheory.Moduli.PrimitiveRoot.CyclicGroup a)
instance Control.DeepSeq.NFData (Math.NumberTheory.Primes.Types.Prime a) => Control.DeepSeq.NFData (Math.NumberTheory.Moduli.PrimitiveRoot.CyclicGroup a)


module Math.NumberTheory.Moduli.DiscreteLogarithm

-- | Computes the discrete logarithm. Currently uses a combination of the
--   baby-step giant-step method and Pollard's rho algorithm, with Bach
--   reduction.
discreteLogarithm :: KnownNat m => PrimitiveRoot m -> MultMod m -> Natural


-- | Miscellaneous functions related to modular arithmetic.
module Math.NumberTheory.Moduli


-- | Bulk evaluation of arithmetic functions over continuous intervals
--   without factorisation.
module Math.NumberTheory.ArithmeticFunctions.SieveBlock

-- | <a>runFunctionOverBlock</a> <tt>f</tt> <tt>x</tt> <tt>l</tt> evaluates
--   an arithmetic function for integers between <tt>x</tt> and
--   <tt>x+l-1</tt> and returns a vector of length <tt>l</tt>. It
--   completely avoids factorisation, so it is asymptotically faster than
--   pointwise evaluation of <tt>f</tt>.
--   
--   Value of <tt>f</tt> at 0, if zero falls into block, is undefined.
--   
--   Beware that for underlying non-commutative monoids the results may
--   potentially differ from pointwise application via <a>runFunction</a>.
--   
--   This is a thin wrapper over <a>sieveBlock</a>, read more details
--   there.
--   
--   <pre>
--   &gt;&gt;&gt; runFunctionOverBlock carmichaelA 1 10
--   [1,1,2,2,4,2,6,2,6,4]
--   </pre>
runFunctionOverBlock :: ArithmeticFunction Word a -> Word -> Word -> Vector a

-- | A record, which specifies a function to evaluate over a block.
--   
--   For example, here is a configuration for the totient function:
--   
--   <pre>
--   SieveBlockConfig
--     { sbcEmpty                = 1
--     , sbcFunctionOnPrimePower = (\p a -&gt; (p - 1) * p ^ (a - 1)
--     , sbcAppend               = (*)
--     }
--   </pre>
data SieveBlockConfig a
SieveBlockConfig :: a -> (Word -> Word -> a) -> (a -> a -> a) -> SieveBlockConfig a

-- | value of a function on 1
[sbcEmpty] :: SieveBlockConfig a -> a

-- | how to evaluate a function on prime powers
[sbcFunctionOnPrimePower] :: SieveBlockConfig a -> Word -> Word -> a

-- | how to combine values of a function on coprime arguments
[sbcAppend] :: SieveBlockConfig a -> a -> a -> a

-- | Create a config for a multiplicative function from its definition on
--   prime powers.
multiplicativeSieveBlockConfig :: Num a => (Word -> Word -> a) -> SieveBlockConfig a

-- | Create a config for an additive function from its definition on prime
--   powers.
additiveSieveBlockConfig :: Num a => (Word -> Word -> a) -> SieveBlockConfig a

-- | Evaluate a function over a block in accordance to provided
--   configuration. Value of <tt>f</tt> at 0, if zero falls into block, is
--   undefined.
--   
--   Based on Algorithm M of <a>Parity of the number of primes in a given
--   interval and algorithms of the sublinear summation</a> by A. V.
--   Lelechenko. See Lemma 2 on p. 5 on its algorithmic complexity. For the
--   majority of use-cases its time complexity is O(x^(1+ε)).
--   
--   <a>sieveBlock</a> is similar to <a>sieveBlockUnboxed</a> up to flavour
--   of <a>Vector</a>, but is typically 7x-10x slower and consumes 3x
--   memory. Use unboxed version whenever possible.
--   
--   For example, following code lists smallest prime factors:
--   
--   <pre>
--   &gt;&gt;&gt; sieveBlock (SieveBlockConfig maxBound (\p _ -&gt; p) min) 2 10
--   [2,3,2,5,2,7,2,3,2,11]
--   </pre>
--   
--   And this is how to factorise all numbers in a block:
--   
--   <pre>
--   &gt;&gt;&gt; sieveBlock (SieveBlockConfig [] (\p k -&gt; [(p,k)]) (++)) 2 10
--   [[(2,1)],[(3,1)],[(2,2)],[(5,1)],[(2,1),(3,1)],[(7,1)],[(2,3)],[(3,2)],[(2,1),(5,1)],[(11,1)]]
--   </pre>
sieveBlock :: SieveBlockConfig a -> Word -> Word -> Vector a

-- | Evaluate a function over a block in accordance to provided
--   configuration. Value of <tt>f</tt> at 0, if zero falls into block, is
--   undefined.
--   
--   Based on Algorithm M of <a>Parity of the number of primes in a given
--   interval and algorithms of the sublinear summation</a> by A. V.
--   Lelechenko. See Lemma 2 on p. 5 on its algorithmic complexity. For the
--   majority of use-cases its time complexity is O(x^(1+ε)).
--   
--   For example, here is an analogue of divisor function <tt>tau</tt>:
--   
--   <pre>
--   &gt;&gt;&gt; sieveBlockUnboxed (SieveBlockConfig 1 (\_ a -&gt; a + 1) (*)) 1 10
--   [1,2,2,3,2,4,2,4,3,4]
--   </pre>
sieveBlockUnboxed :: Unbox a => SieveBlockConfig a -> Word -> Word -> Vector a

-- | Evaluate the Möbius function over a block. Value of <tt>f</tt> at 0,
--   if zero falls into block, is undefined.
--   
--   Based on the sieving algorithm from p. 3 of <a>Computations of the
--   Mertens function and improved bounds on the Mertens conjecture</a> by
--   G. Hurst. It is approximately 5x faster than <a>sieveBlockUnboxed</a>.
--   
--   <pre>
--   &gt;&gt;&gt; sieveBlockMoebius 1 10
--   [MoebiusP, MoebiusN, MoebiusN, MoebiusZ, MoebiusN, MoebiusP, MoebiusN, MoebiusZ, MoebiusZ, MoebiusP]
--   </pre>
sieveBlockMoebius :: Word -> Word -> Vector Moebius


-- | Values of <a>Mertens function</a>.
module Math.NumberTheory.ArithmeticFunctions.Mertens

-- | Compute individual values of Mertens function in O(n^(2/3)) time and
--   space.
--   
--   <pre>
--   &gt;&gt;&gt; map (mertens . (10 ^)) [0..9]
--   [1,-1,1,2,-23,-48,212,1037,1928,-222]
--   </pre>
--   
--   The implementation follows Theorem 3.1 from <a>Computations of the
--   Mertens function and improved bounds on the Mertens conjecture</a> by
--   G. Hurst, excluding segmentation of sieves.
mertens :: Word -> Int


-- | Riemann zeta-function.
module Math.NumberTheory.Zeta.Riemann

-- | Infinite sequence of approximate (up to given precision) values of
--   Riemann zeta-function at integer arguments, starting with
--   <tt>ζ(0)</tt>. Computations for odd arguments are performed in
--   accordance to <a>Computational strategies for the Riemann zeta
--   function</a> by J. M. Borwein, D. M. Bradley, R. E. Crandall, formula
--   (57).
--   
--   <pre>
--   &gt;&gt;&gt; take 5 (zetas 1e-14) :: [Double]
--   &gt; [-0.5,Infinity,1.6449340668482262,1.2020569031595942,1.0823232337111381]
--   </pre>
--   
--   Beware to force evaluation of <tt>zetas !! 1</tt>, if the type
--   <tt>a</tt> does not support infinite values (for instance,
--   <a>Fixed</a>).
zetas :: (Floating a, Ord a) => a -> [a]

-- | Infinite sequence of exact values of Riemann zeta-function at even
--   arguments, starting with <tt>ζ(0)</tt>. Note that due to numerical
--   errors conversion to <a>Double</a> may return values below 1:
--   
--   <pre>
--   &gt;&gt;&gt; approximateValue (zetasEven !! 25) :: Double
--   0.9999999999999996
--   </pre>
--   
--   Use your favorite type for long-precision arithmetic. For instance,
--   <a>Fixed</a> works fine:
--   
--   <pre>
--   &gt;&gt;&gt; import Data.Number.Fixed
--   
--   &gt;&gt;&gt; approximateValue (zetasEven !! 25) :: Fixed Prec50
--   1.00000000000000088817842111574532859293035196051773
--   </pre>
zetasEven :: [ExactPi]
zetasOdd :: forall a. (Floating a, Ord a) => a -> [a]


-- | Dirichlet beta-function.
module Math.NumberTheory.Zeta.Dirichlet

-- | Infinite sequence of approximate (up to given precision) values of
--   Dirichlet beta-function at integer arguments, starting with
--   <tt>β(0)</tt>. The algorithm used to compute <tt>β</tt> for even
--   arguments was derived from <a>An Euler-type formula for β(2n) and
--   closed-form expressions for a class of zeta series</a> by F. M. S.
--   Lima, formula (12).
--   
--   <pre>
--   &gt; take 5 (betas 1e-14) :: [Double]
--   [0.5,0.7853981633974483,0.9159655941772191,0.9689461462593693,0.988944551741105]
--   </pre>
betas :: (Floating a, Ord a) => a -> [a]

-- | Infinite sequence of approximate values of the Dirichlet <tt>β</tt>
--   function at positive even integer arguments, starting with
--   <tt>β(0)</tt>.
betasEven :: forall a. (Floating a, Ord a) => a -> [a]

-- | Infinite sequence of exact values of Dirichlet beta-function at odd
--   arguments, starting with <tt>β(1)</tt>.
--   
--   <pre>
--   &gt; approximateValue (betasOdd !! 25) :: Double
--   0.9999999999999987
--   </pre>
--   
--   Using <a>Fixed</a>:
--   
--   <pre>
--   &gt; approximateValue (betasOdd !! 25) :: Fixed Prec50
--   0.99999999999999999999999960726927497384196726751694z
--   </pre>
betasOdd :: [ExactPi]


-- | Interface to work with Riemann zeta-function and Dirichlet
--   beta-function.
module Math.NumberTheory.Zeta

-- | Joins two lists element-by-element together into one, starting with
--   the first one provided as argument.
--   
--   <pre>
--   &gt;&gt;&gt; take 10 $ intertwine [0, 2 ..] [1, 3 ..]
--   [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
--   </pre>
intertwine :: [a] -> [a] -> [a]

-- | Skips every even-indexed element from an infinite list. Do NOT use
--   with finite lists.
--   
--   <pre>
--   &gt;&gt;&gt; take 10 (skipEvens [0, 1 ..])
--   [1, 3, 5, 7, 9, 11, 13, 15, 17, 19]
--   </pre>
skipEvens :: [a] -> [a]

-- | Skips every odd-indexed element from an infinite list. Do NOT use with
--   finite lists.
--   
--   <pre>
--   &gt;&gt;&gt; take 10 (skipOdds [0, 1 ..])
--   [0, 2, 4, 6, 8, 10, 12, 14, 16, 18]
--   </pre>
skipOdds :: [a] -> [a]

-- | Sums every element of an infinite list up to a certain precision. I.e.
--   once an element falls below the given threshold it stops traversing
--   the list.
--   
--   <pre>
--   &gt;&gt;&gt; suminf 1e-14 (iterate (/ 10) 1)
--   1.1111111111111112
--   </pre>
suminf :: (Floating a, Ord a) => a -> [a] -> a
