-
Notifications
You must be signed in to change notification settings - Fork 0
/
qs.hs
104 lines (92 loc) · 3.4 KB
/
qs.hs
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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
module QuadraticSieve where
import Control.Monad
import Data.Maybe
import qualified Data.IntSet as IntSet
import Data.List
import Data.Function
import Math.NumberTheory.Primes
import Math.NumberTheory.Moduli.Jacobi
import Math.NumberTheory.SmoothNumbers
-- Smoothness bound for @n@, computed as @exp(0.5 * sqrt(ln(n) * ln(ln(n))))@.
smoothnessBound :: Integer -> Integer
smoothnessBound =
ceiling . exp . (* 0.5) . sqrt . liftM2 (*) log (log . log) . fromInteger
-- Primes smaller than the smoothness bound for @n@ and for which @n@ is a quadratic residue.
factorBase :: Integer -> [Integer]
factorBase n =
takeWhile (<= smoothnessBound n)
$ unPrime (nextPrime 2)
: [ prime | prime <- unPrime <$> [nextPrime 3 ..], jacobi n prime == One ]
type Relation = (Integer, Integer)
-- Find smooth relations.
--
-- In a smooth relation @(x, x^2 - n)@, @x^2 - n@ is smooth over a given base.
smoothRelations :: [Integer] -> Integer -> [Relation]
smoothRelations base n =
let start = ceiling . sqrt $ fromInteger n
smoothBasis = fromJust $ fromList base
size = 1 + length base
in take
size
[ (x, y) | x <- [start ..], let y = x ^ 2 - n, isSmooth smoothBasis y ]
-- Find maximal @p@ such that @prime ^ p@ divides @n@
factorPower :: Integer -> Integer -> Integer
factorPower prime n =
if n `mod` prime == 0 then 1 + factorPower prime (n `div` prime) else 0
-- Give the exponent vector of @n@.
--
-- Exponent vector of @n@, over a given base @base@,
-- contains the powers of the primes in @base@ that factor @n@.
--
-- >>> exponentVector [2,3,5,7] 40
-- [3,0,1,0]
exponentVector :: [Integer] -> Integer -> [Integer]
exponentVector base n = snd $ foldr
(\prime (rem, vector) ->
let power = factorPower prime rem
in (rem `div` (prime ^ power), power : vector)
)
(n, [])
base
-- Find coefficients of a linear combination of @rows@.
combination :: [[Integer]] -> IntSet.IntSet
combination rows =
let coeffs = zip (map IntSet.singleton [0 ..]) rows
combineVecs
:: (IntSet.IntSet, [Integer])
-> (IntSet.IntSet, [Integer])
-> (IntSet.IntSet, [Integer])
combineVecs (s1, x) (s2, y) =
(IntSet.union s1 s2, map (`mod` 2) $ zipWith (+) x y)
sortCombinations = sortBy (compare `on` snd)
in fst . head . sortCombinations $ foldr
(\x cs -> x : [ combineVecs x y | y <- cs ] ++ cs)
[]
coeffs
-- Take only the elements of @xs@ specified by the indexes in @indexes@.
filterByIndexes :: IntSet.IntSet -> [a] -> [a]
filterByIndexes indexes xs =
reverse
. fst
. foldl'
(\(result, rest) i ->
let rest' = drop i rest in (head rest' : result, rest')
)
([], xs)
. diff
. IntSet.toAscList
$ indexes
where
diff [] = []
diff is@(i : _) = i : zipWith (-) (drop 1 is) is
-- Factor a given composite number @n@ using the quadratic sieve.
factor n =
let base = factorBase n
smoothRels = smoothRelations base n
exponentVecs = map (exponentVector base . snd) smoothRels
dependentSubset = (filterByIndexes =<< combination) exponentVecs
primePowers = foldr (zipWith (+)) [0 ..] dependentSubset
y = foldr (\(prime, exp) acc -> (acc * prime ^ exp) `mod` n) 1
$ zip base primePowers
x = product . map fst $ smoothRels
in gcd (x - y) n