1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
-- |
-- 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 prs n
| n < 0 = trialDivisionWith prs (-n)
| n == 0 = error "trialDivision of 0"
| n == 1 = []
| otherwise = go n (integerSquareRoot n) prs
where
go !m !r (p:ps)
| r < p = [(m,1)]
| otherwise =
case splitOff p m of
(0,_) -> go m r ps
(k,q) -> (p,k) : if q == 1
then []
else go q (integerSquareRoot q) ps
go m _ _ = [(m,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 bd
| bd < 100 = trialDivisionTo 100
| bd < 10000000 = trialDivisionWith (map unPrime $ primeList $ primeSieve bd)
| otherwise = trialDivisionWith (takeWhile (<= bd) $ map unPrime $ psieveList >>= 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 prs n
| n < 0 = trialDivisionPrimeWith prs (-n)
| n < 2 = False
| otherwise = go n (integerSquareRoot n) prs
where
go !m !r (p:ps) = r < p || m `rem` p /= 0 && go m r ps
go _ _ _ = 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 bd
| bd < 100 = trialDivisionPrimeTo 100
| bd < 10000000 = trialDivisionPrimeWith (map unPrime $ primeList $ primeSieve bd)
| otherwise = trialDivisionPrimeWith (takeWhile (<= bd) $ map unPrime $ psieveList >>= primeList)
|
