-- 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.6.0.1


-- | 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)


-- | 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; (3 :: Mod 10) + (4 :: Mod 12)
--   error: Couldn't match type ‘12’ with ‘10’...
--   &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 :: KnownNat m => Mod m -> Integer

-- | The canonical representative of the residue class, always between 0
--   and <tt>m-1</tt> inclusively.
getNatVal :: KnownNat m => 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; invertMod (3 :: Mod 10)
--   Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10
--   &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; powMod (3 :: Mod 10) 4
--   (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 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; 2 `modulo` 10 + 4 `modulo` 15
--   (1 `modulo` 5)
--   &gt; 2 `modulo` 10 * 4 `modulo` 15
--   (3 `modulo` 5)
--   &gt; 2 `modulo` 10 + fromRational (3 % 7)
--   (1 `modulo` 10)
--   &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; invertSomeMod (3 `modulo` 10)
--   Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10
--   &gt; invertMod (4 `modulo` 10)
--   Nothing
--   &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; 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.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.Show.Show (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)


-- | Potentially faster power function for <a>Integer</a> base and
--   <a>Int</a> or <a>Word</a> exponent.

-- | <i>Deprecated: It is no faster than (^)</i>
module Math.NumberTheory.Powers.Integer

-- | Power of an <a>Integer</a> by the left-to-right repeated squaring
--   algorithm. This needs two multiplications in each step while the
--   right-to-left algorithm needs only one multiplication for 0-bits, but
--   here the two factors always have approximately the same size, which on
--   average gains a bit when the result is large.
--   
--   For small results, it is unlikely to be any faster than '(^)', quite
--   possibly slower (though the difference shouldn't be large), and for
--   exponents with few bits set, the same holds. But for exponents with
--   many bits set, the speedup can be significant.
--   
--   <i>Warning:</i> No check for the negativity of the exponent is
--   performed, a negative exponent is interpreted as a large positive
--   exponent.
integerPower :: Integer -> Int -> Integer

-- | Same as <a>integerPower</a>, but for exponents of type <a>Word</a>.
integerWordPower :: Integer -> Word -> Integer


-- | Internal functions dealing with square roots. End-users should not
--   import this module.
module Math.NumberTheory.Powers.Squares.Internal
karatsubaSqrt :: Integer -> (Integer, Integer)
isqrtA :: Integral a => a -> a


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

-- | Infinite zero-based table of factorials.
--   
--   <pre>
--   &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; :set +s
--   &gt; binomial !! 1000 !! 1000 :: Integer
--   1
--   (0.01 secs, 1,385,512 bytes)
--   &gt; binomial !! 1000 !! 1000 :: Integer
--   1
--   (0.01 secs, 1,381,616 bytes)
--   </pre>
--   
--   against
--   
--   <pre>
--   &gt; let binomial' = binomial :: [[Integer]]
--   &gt; binomial' !! 1000 !! 1000 :: Integer
--   1
--   (0.01 secs, 1,381,696 bytes)
--   &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; 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; 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; 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; 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; 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; 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; 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]


-- | Prime generation using a priority queue for the composites. The
--   algorithm is basically the one described in
--   <a>http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf</a>, but it uses
--   a more efficient heap for the priority queue and a larger wheel, thus
--   it is faster (in particular, the speed penalty for
--   <tt><a>Integer</a></tt> is much smaller) and uses less memory. It is
--   nevertheless very slow compared to a bit sieve. This module is mainly
--   intended for comparison and verification.
module Math.NumberTheory.Primes.Heap

-- | A list of primes. The sieve does not handle overflow, hence for
--   bounded types, garbage occurs near <tt><a>maxBound</a></tt>. If primes
--   that large are requested, use
--   
--   <pre>
--   <a>map</a> <a>fromInteger</a> $ <a>takeWhile</a> (&lt;= <a>fromIntegral</a> <a>maxBound</a>) <a>primes</a>
--   </pre>
--   
--   instead. Checking for overflow would be slower. The sieve is
--   specialised for <tt><a>Int</a></tt>, <tt><a>Word</a></tt> and
--   <tt><a>Integer</a></tt>, since these are the most commonly used. For
--   the fixed-width <tt>Int</tt> or <tt>Word</tt> types, sieving at one of
--   the specialised types and converting is faster. To ensure sharing of
--   the list of primes, bind it to a monomorphic variable, to make sure
--   that it is not shared, use <tt><a>sieveFrom</a></tt> with different
--   arguments.
primes :: Integral a => [a]

-- | <tt><a>sieveFrom</a> n</tt> generates the list of primes <tt>&gt;=
--   n</tt>. The remarks about overflow and performance in the
--   documentation of <tt><a>primes</a></tt> apply here too.
sieveFrom :: Integral a => 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 Moebius inversion for <a>Int</a> valued functions.
module Math.NumberTheory.MoebiusInversion.Int

-- | <tt>generalInversion g n</tt> evaluates the generalised Moebius
--   inversion of <tt>g</tt> at the argument <tt>n</tt>.
--   
--   The generalised Moebius 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 Moebius 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 Moebius function are known. A slightly different formula, used
--   here, does not need the values of the Moebius 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
--   Moebius inversion. See
--   <a>http://mathworld.wolfram.com/TotientSummatoryFunction.html</a> for
--   the formula used for <tt>totientSum</tt>.
totientSum :: Int -> Int


-- | Generalised Moebius inversion
module Math.NumberTheory.MoebiusInversion

-- | <tt>generalInversion g n</tt> evaluates the generalised Moebius
--   inversion of <tt>g</tt> at the argument <tt>n</tt>.
--   
--   The generalised Moebius 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 Moebius 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 Moebius function are known. A slightly different formula, used
--   here, does not need the values of the Moebius 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
--   Moebius inversion. See
--   <a>http://mathworld.wolfram.com/TotientSummatoryFunction.html</a> for
--   the formula used for <tt>totientSum</tt>.
totientSum :: Int -> Integer


-- | 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

-- | List of primes. Since the sieve uses unboxed arrays, overflow occurs
--   at some point. On 64-bit systems, that point is beyond the memory
--   limits, on 32-bit systems, it is at about <tt>1.7*10^18</tt>.
primes :: [Integer]

-- | <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 :: PrimeSieve -> [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


-- | 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)


-- | 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)


-- | Jacobi symbol.
module Math.NumberTheory.Moduli.Jacobi

-- | Type for result of <a>jacobi</a>.
data JacobiSymbol
MinusOne :: JacobiSymbol
Zero :: JacobiSymbol
One :: JacobiSymbol

-- | Jacobi symbol of two numbers. The "denominator" must be odd and
--   positive, this condition is checked.
--   
--   If both numbers have a common prime factor, the result is <tt>0</tt>,
--   otherwise it is ±1.
jacobi :: (Integral a, Bits a) => a -> a -> JacobiSymbol

-- | Jacobi symbol of two numbers without validity check of the
--   "denominator".
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 Data.Semigroup.Semigroup Math.NumberTheory.Moduli.Jacobi.JacobiSymbol
instance GHC.Base.Monoid Math.NumberTheory.Moduli.Jacobi.JacobiSymbol


-- | 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


-- | 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

-- | A compact store of smallest prime factors.
data FactorSieve

-- | Test primality using a <a>FactorSieve</a>. If <tt>n</tt> is out of
--   bounds of the sieve, fall back to <a>isPrime</a>.
fsIsPrime :: FactorSieve -> 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


-- | 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>, including a factor of <tt>(-1)</tt> if <tt>n &lt; 0</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)]

-- | A compact store of smallest prime factors.
data FactorSieve

-- | <tt><a>factorSieve</a> n</tt> creates a store of smallest prime
--   factors of the numbers not exceeding <tt>n</tt>. If you need to
--   factorise many smallish numbers, this can give a big speedup since it
--   avoids many superfluous divisions. However, a too large sieve leads to
--   a slowdown due to cache misses. The prime factors are stored as
--   <a>Word16</a> for compactness, so <tt>n</tt> must be smaller than
--   <tt>2^32</tt>.
factorSieve :: Integer -> FactorSieve

-- | <tt><a>sieveFactor</a> fs n</tt> finds the prime factorisation of
--   <tt>n</tt> using the <a>FactorSieve</a> <tt>fs</tt>. For negative
--   <tt>n</tt>, a factor of <tt>-1</tt> is included with multiplicity
--   <tt>1</tt>. After stripping any present factors <tt>2</tt>, the
--   remaining cofactor <tt>c</tt> (if larger than <tt>1</tt>) is
--   factorised with <tt>fs</tt>. This is most efficient of course if
--   <tt>c</tt> does not exceed the bound with which <tt>fs</tt> was
--   constructed. If it does, trial division is performed until either the
--   cofactor falls below the bound or the sieve is exhausted. In the
--   latter case, the elliptic curve method is used to finish the
--   factorisation.
sieveFactor :: FactorSieve -> 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)]

-- | <tt><a>curveFactorisation</a></tt> 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.
--   
--   <tt><a>curveFactorisation</a></tt> requires that small prime factors
--   have been stripped before. Also, it is unlikely to succeed if
--   <tt>n</tt> has more than one (really) large prime factor.
curveFactorisation :: 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

-- | A compact store of totients.
data TotientSieve

-- | <tt><a>totientSieve</a> n</tt> creates a store of the totients of the
--   numbers not exceeding <tt>n</tt>. A <a>TotientSieve</a> only stores
--   values for numbers coprime to <tt>30</tt> to reduce space usage. The
--   maximal admissible value for <tt>n</tt> is <tt><a>fromIntegral</a>
--   (<a>maxBound</a> :: <a>Word</a>)</tt>.
totientSieve :: Integer -> TotientSieve

-- | <tt><a>sieveTotient</a> ts n</tt> finds the totient <tt>π(n)</tt>,
--   i.e. the number of integers <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>, using the <a>TotientSieve</a>
--   <tt>ts</tt>. First, factors of <tt>2, 3</tt> and <tt>5</tt> are
--   handled individually, if the remaining cofactor of <tt>n</tt> is
--   within the sieve range, its totient is looked up, otherwise the
--   calculation falls back on factorisation, first by trial division and
--   if necessary, elliptic curves.
sieveTotient :: TotientSieve -> Integer -> Integer

-- | A compact store of values of the Carmichael function.
data CarmichaelSieve

-- | <tt><a>carmichaelSieve</a> n</tt> creates a store of values of the
--   Carmichael function for numbers not exceeding <tt>n</tt>. Like a
--   <a>TotientSieve</a>, a <a>CarmichaelSieve</a> only stores values for
--   numbers coprime to <tt>30</tt> to reduce space usage. The maximal
--   admissible value for <tt>n</tt> is <tt><a>fromIntegral</a>
--   (<a>maxBound</a> :: <a>Word</a>)</tt>.
carmichaelSieve :: Integer -> CarmichaelSieve

-- | <tt><a>sieveCarmichael</a> cs n</tt> finds the value of <tt>λ(n)</tt>
--   (or <tt>ψ(n)</tt>), the smallest positive integer <tt>e</tt> such that
--   for all <tt>a</tt> with <tt>gcd a n == 1</tt> the congruence <tt>a^e ≡
--   1 (mod n)</tt> holds, in other words, the (smallest) exponent of the
--   group of units in <tt>ℤ/(n)</tt>. The strategy is analogous to
--   <a>sieveTotient</a>.
sieveCarmichael :: CarmichaelSieve -> Integer -> Integer


module Math.NumberTheory.Primes


-- | 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)]


-- | 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.
gcdInt :: Int -> Int -> Int

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

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

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

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

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

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

-- | Test whether two <a>Word#</a>s are coprime.
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.
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).
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.
coprime :: (Integral a, Bits a) => a -> a -> Bool


-- | 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


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

-- | <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>.
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.
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.
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.
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>
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.
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.
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.
sqrtModFList :: Integer -> [(Integer, Int)] -> [Integer]


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


-- | This module exports functions for manipulating Gaussian integers,
--   including computing their prime factorisations.
module Math.NumberTheory.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

-- | 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

-- | Simultaneous <a>div</a> and <a>mod</a>.
divModG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)

-- | Gaussian integer division, truncating toward negative infinity.
divG :: GaussianInteger -> GaussianInteger -> GaussianInteger

-- | Gaussian integer remainder, satisfying
--   
--   <pre>
--   (x `divG` y)*y + (x `modG` y) == x
--   </pre>
modG :: GaussianInteger -> GaussianInteger -> GaussianInteger

-- | Simultaneous <a>quot</a> and <a>rem</a>.
quotRemG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)

-- | Gaussian integer division, truncating toward zero.
quotG :: GaussianInteger -> GaussianInteger -> GaussianInteger

-- | Gaussian integer remainder, satisfying
--   
--   <pre>
--   (x `quotG` y)*y + (x `remG` y) == x
--   </pre>
remG :: GaussianInteger -> GaussianInteger -> GaussianInteger

-- | Raise a Gaussian integer to a given power.
(.^) :: (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, but not quite in order of
--   ascending magnitude.
primes :: [GaussianInteger]

-- | Compute the GCD of two Gaussian integers. Enforces the precondition
--   that each integer must be in the first quadrant (or zero).
gcdG :: GaussianInteger -> GaussianInteger -> GaussianInteger

-- | Compute the GCD of two Gauss integers. Does not check the
--   precondition.
gcdG' :: GaussianInteger -> GaussianInteger -> GaussianInteger

-- | Find a Gaussian integer whose norm is the given prime number. Checks
--   the precondition that p is prime and that p <a>mod</a> 4 /= 3.
findPrime :: Integer -> GaussianInteger

-- | Find a Gaussian integer whose norm is the given prime number. Does not
--   check the precondition.
findPrime' :: Integer -> GaussianInteger

-- | Compute the prime factorisation of a Gaussian integer. This is unique
--   up to units (+<i>- 1, +</i>- i).
factorise :: GaussianInteger -> [(GaussianInteger, Int)]
instance GHC.Classes.Eq Math.NumberTheory.GaussianIntegers.GaussianInteger
instance GHC.Show.Show Math.NumberTheory.GaussianIntegers.GaussianInteger
instance GHC.Num.Num Math.NumberTheory.GaussianIntegers.GaussianInteger


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

-- | Type of primes of a given unique factorisation domain. When the domain
--   has exactly one unit, <tt>Prime t = t</tt>, but when units are
--   multiple more restricted types (or at least newtypes) should be
--   specified.
--   
--   <tt>abs (unPrime n) == unPrime n</tt> must hold for all <tt>n</tt> of
--   type <tt>Prime t</tt>

-- | 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)]
instance GHC.Show.Show Math.NumberTheory.UniqueFactorisation.GaussianPrime
instance GHC.Classes.Eq Math.NumberTheory.UniqueFactorisation.GaussianPrime
instance GHC.Show.Show Math.NumberTheory.UniqueFactorisation.BigPrime
instance GHC.Classes.Ord Math.NumberTheory.UniqueFactorisation.BigPrime
instance GHC.Classes.Eq Math.NumberTheory.UniqueFactorisation.BigPrime
instance GHC.Show.Show Math.NumberTheory.UniqueFactorisation.SmallPrime
instance GHC.Classes.Ord Math.NumberTheory.UniqueFactorisation.SmallPrime
instance GHC.Classes.Eq Math.NumberTheory.UniqueFactorisation.SmallPrime
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.GaussianIntegers.GaussianInteger


-- | 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 <a>Product</a> or <a>Sum</a>), and the second argument is
--   convenient for unwrapping (typically, <a>getProduct</a> or
--   <a>getSum</a>).
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, Integral 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, Integral n) => ArithmeticFunction n n
jordan :: (UniqueFactorisation n, Integral n) => Word -> n -> n

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

-- | Calculates the Moebius function of an argument.
moebiusA :: Num a => ArithmeticFunction n 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


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

-- | 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; take 5 (zetas 1e-14) :: [Double]
--   [-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 convertation to <a>Double</a> may return values below 1:
--   
--   <pre>
--   &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; approximateValue (zetasEven !! 25) :: Fixed Prec50
--   1.00000000000000088817842111574532859293035196051773
--   </pre>
zetasEven :: [ExactPi]

-- | Approximates an exact or approximate value, converting it to a
--   <a>Floating</a> type. This uses the value of <a>pi</a> supplied by the
--   destination type, to provide the appropriate precision.
approximateValue :: Floating a => ExactPi -> a
