    -- |
    -- Module:      Math.NumberTheory.Primes.Factorisation.TrialDivision
    -- Copyright:   (c) 2011 Daniel Fischer
    -- Licence:     MIT
    -- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>
    --
    -- Factorisation and primality testing using trial division.
    --
    -- Used to create and check certificates.
    -- Currently not exposed because it's not that useful, is it?
    -- But the trial...To functions are exported from other modules.
    {-# LANGUAGE BangPatterns #-}
    module Math.NumberTheory.Primes.Factorisation.TrialDivision
        ( trialDivisionWith
        , trialDivisionTo
        , trialDivisionPrimeTo
        ) where

    import Math.NumberTheory.Primes.Sieve.Eratosthenes (primeList, primeSieve, psieveList)
    import Math.NumberTheory.Roots
    import Math.NumberTheory.Primes.Types
    import Math.NumberTheory.Utils

    -- | Factorise an 'Integer' using a given list of numbers considered prime.
    --   If the list is not a list of primes containing all relevant primes, the
    --   result could be surprising.
    trialDivisionWith :: [Integer] -> Integer -> [(Integer, Word)]
    trialDivisionWith :: [Integer] -> Integer -> [(Integer, Word)]
    trialDivisionWith [Integer]
    prs Integer
    n
        | Integer
    n Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    < Integer
    0     = [Integer] -> Integer -> [(Integer, Word)]
    trialDivisionWith [Integer]
    prs (-Integer
    n)
        | Integer
    n Integer -> Integer -> Bool
    forall a. Eq a => a -> a -> Bool
    == Integer
    0    = [Char] -> [(Integer, Word)]
    forall a. HasCallStack => [Char] -> a
    error [Char]
    "trialDivision of 0"
        | Integer
    n Integer -> Integer -> Bool
    forall a. Eq a => a -> a -> Bool
    == Integer
    1    = []
        | Bool
    otherwise = Integer -> Integer -> [Integer] -> [(Integer, Word)]
    forall {a}.
    (GcdDomain a, Integral a) =>
    a -> a -> [a] -> [(a, Word)]
    go Integer
    n (Integer -> Integer
    forall a. Integral a => a -> a
    integerSquareRoot Integer
    n) [Integer]
    prs
          where
            go :: a -> a -> [a] -> [(a, Word)]
    go !a
    m !a
    r (a
    p:[a]
    ps)
                | a
    r a -> a -> Bool
    forall a. Ord a => a -> a -> Bool
    < a
    p     = [(a
    m,Word
    1)]
                | Bool
    otherwise =
                    case a -> a -> (Word, a)
    forall a. (Eq a, GcdDomain a) => a -> a -> (Word, a)
    splitOff a
    p a
    m of
                      (Word
    0,a
    _) -> a -> a -> [a] -> [(a, Word)]
    go a
    m a
    r [a]
    ps
                      (Word
    k,a
    q) -> (a
    p,Word
    k) (a, Word) -> [(a, Word)] -> [(a, Word)]
    forall a. a -> [a] -> [a]
    : if a
    q a -> a -> Bool
    forall a. Eq a => a -> a -> Bool
    == a
    1
                                         then []
                                         else a -> a -> [a] -> [(a, Word)]
    go a
    q (a -> a
    forall a. Integral a => a -> a
    integerSquareRoot a
    q) [a]
    ps
            go a
    m a
    _ [a]
    _    = [(a
    m,Word
    1)]

    -- | @'trialDivisionTo' bound n@ produces a factorisation of @n@ using the
    --   primes @<= bound@. If @n@ has prime divisors @> bound@, the last entry
    --   in the list is the product of all these. If @n <= bound^2@, this is a
    --   full factorisation, but very slow if @n@ has large prime divisors.
    trialDivisionTo :: Integer -> Integer -> [(Integer, Word)]
    trialDivisionTo :: Integer -> Integer -> [(Integer, Word)]
    trialDivisionTo Integer
    bd
        | Integer
    bd Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    < Integer
    100      = Integer -> Integer -> [(Integer, Word)]
    trialDivisionTo Integer
    100
        | Integer
    bd Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    < Integer
    10000000 = [Integer] -> Integer -> [(Integer, Word)]
    trialDivisionWith ((Prime Integer -> Integer) -> [Prime Integer] -> [Integer]
    forall a b. (a -> b) -> [a] -> [b]
    map Prime Integer -> Integer
    forall a. Prime a -> a
    unPrime ([Prime Integer] -> [Integer]) -> [Prime Integer] -> [Integer]
    forall a b. (a -> b) -> a -> b
    $ PrimeSieve -> [Prime Integer]
    forall a. Integral a => PrimeSieve -> [Prime a]
    primeList (PrimeSieve -> [Prime Integer]) -> PrimeSieve -> [Prime Integer]
    forall a b. (a -> b) -> a -> b
    $ Integer -> PrimeSieve
    primeSieve Integer
    bd)
        | Bool
    otherwise     = [Integer] -> Integer -> [(Integer, Word)]
    trialDivisionWith ((Integer -> Bool) -> [Integer] -> [Integer]
    forall a. (a -> Bool) -> [a] -> [a]
    takeWhile (Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Integer
    bd) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
    forall a b. (a -> b) -> a -> b
    $ (Prime Integer -> Integer) -> [Prime Integer] -> [Integer]
    forall a b. (a -> b) -> [a] -> [b]
    map Prime Integer -> Integer
    forall a. Prime a -> a
    unPrime ([Prime Integer] -> [Integer]) -> [Prime Integer] -> [Integer]
    forall a b. (a -> b) -> a -> b
    $ [PrimeSieve]
    psieveList [PrimeSieve] -> (PrimeSieve -> [Prime Integer]) -> [Prime Integer]
    forall a b. [a] -> (a -> [b]) -> [b]
    forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
    >>= PrimeSieve -> [Prime Integer]
    forall a. Integral a => PrimeSieve -> [Prime a]
    primeList)

    -- | Check whether a number is coprime to all of the numbers in the list
    --   (assuming that list contains only numbers > 1 and is ascending).
    trialDivisionPrimeWith :: [Integer] -> Integer -> Bool
    trialDivisionPrimeWith :: [Integer] -> Integer -> Bool
    trialDivisionPrimeWith [Integer]
    prs Integer
    n
        | Integer
    n Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    < Integer
    0     = [Integer] -> Integer -> Bool
    trialDivisionPrimeWith [Integer]
    prs (-Integer
    n)
        | Integer
    n Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    < Integer
    2     = Bool
    False
        | Bool
    otherwise = Integer -> Integer -> [Integer] -> Bool
    forall {t}. Integral t => t -> t -> [t] -> Bool
    go Integer
    n (Integer -> Integer
    forall a. Integral a => a -> a
    integerSquareRoot Integer
    n) [Integer]
    prs
          where
            go :: t -> t -> [t] -> Bool
    go !t
    m !t
    r (t
    p:[t]
    ps) = t
    r t -> t -> Bool
    forall a. Ord a => a -> a -> Bool
    < t
    p Bool -> Bool -> Bool
    || t
    m t -> t -> t
    forall a. Integral a => a -> a -> a
    `rem` t
    p t -> t -> Bool
    forall a. Eq a => a -> a -> Bool
    /= t
    0 Bool -> Bool -> Bool
    && t -> t -> [t] -> Bool
    go t
    m t
    r [t]
    ps
            go t
    _ t
    _ [t]
    _ = Bool
    True

    -- | @'trialDivisionPrimeTo' bound n@ tests whether @n@ is coprime to all primes @<= bound@.
    --   If @n <= bound^2@, this is a full prime test, but very slow if @n@ has no small prime divisors.
    trialDivisionPrimeTo :: Integer -> Integer -> Bool
    trialDivisionPrimeTo :: Integer -> Integer -> Bool
    trialDivisionPrimeTo Integer
    bd
        | Integer
    bd Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    < Integer
    100      = Integer -> Integer -> Bool
    trialDivisionPrimeTo Integer
    100
        | Integer
    bd Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    < Integer
    10000000 = [Integer] -> Integer -> Bool
    trialDivisionPrimeWith ((Prime Integer -> Integer) -> [Prime Integer] -> [Integer]
    forall a b. (a -> b) -> [a] -> [b]
    map Prime Integer -> Integer
    forall a. Prime a -> a
    unPrime ([Prime Integer] -> [Integer]) -> [Prime Integer] -> [Integer]
    forall a b. (a -> b) -> a -> b
    $ PrimeSieve -> [Prime Integer]
    forall a. Integral a => PrimeSieve -> [Prime a]
    primeList (PrimeSieve -> [Prime Integer]) -> PrimeSieve -> [Prime Integer]
    forall a b. (a -> b) -> a -> b
    $ Integer -> PrimeSieve
    primeSieve Integer
    bd)
        | Bool
    otherwise     = [Integer] -> Integer -> Bool
    trialDivisionPrimeWith ((Integer -> Bool) -> [Integer] -> [Integer]
    forall a. (a -> Bool) -> [a] -> [a]
    takeWhile (Integer -> Integer -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Integer
    bd) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
    forall a b. (a -> b) -> a -> b
    $ (Prime Integer -> Integer) -> [Prime Integer] -> [Integer]
    forall a b. (a -> b) -> [a] -> [b]
    map Prime Integer -> Integer
    forall a. Prime a -> a
    unPrime ([Prime Integer] -> [Integer]) -> [Prime Integer] -> [Integer]
    forall a b. (a -> b) -> a -> b
    $ [PrimeSieve]
    psieveList [PrimeSieve] -> (PrimeSieve -> [Prime Integer]) -> [Prime Integer]
    forall a b. [a] -> (a -> [b]) -> [b]
    forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
    >>= PrimeSieve -> [Prime Integer]
    forall a. Integral a => PrimeSieve -> [Prime a]
    primeList)