-- |
-- 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 (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 (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= PrimeSieve -> [Prime Integer]
forall a. Integral a => PrimeSieve -> [Prime a]
primeList)