Class: FasterPrime::MPQS

Inherits:
Object
  • Object
show all
Defined in:
lib/faster_prime/prime_factorization.rb

Overview

An implementation of MPQS factoring algorithm.

    1. Silverman,

“The Multiple Polynomial Quadratic Sieve.” (1987) www.ams.org/journals/mcom/1987-48-177/S0025-5718-1987-0866119-8/S0025-5718-1987-0866119-8.pdf

Defined Under Namespace

Classes: GaussianElimination

Constant Summary collapse

PARAMETERS =

MPQS parameters (N digits, factor base size, sieve interval, and large prime tolerance)

[
  [24,  100,    250, 1.5 ],
  [30,  200,  1_250, 1.5 ],
  [36,  400,  1_250, 1.75],
  [42,  900,  2_500, 2.0 ],
  [48, 1200,  7_000, 2.0 ],
  [54, 2000, 12_500, 2.2 ],
  [60, 3000, 17_500, 2.4 ],
  [66, 4500, 25_000, 2.6 ],
]
SMALL_PRIMES = 
[2, 3, 5, 7, 11]
S = 
65536

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(n) ⇒ MPQS

Returns a new instance of MPQS.



150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
 
# File 'lib/faster_prime/prime_factorization.rb', line 150

def initialize(n)
  @n = n

  if @n <= 1 || PrimalityTest.prime?(@n)
    @trivial_factor = 1
    return
  end
  SMALL_PRIMES.each do |p|
    if @n % p == 0
      @trivial_factor = p
      return
    end
  end

  m = Utils.integer_square_root(n)
  if m * m == n
    @trivial_factor = m
    return
  end

  # select parameters
  # @k: a multiplier
  @k = MPQS.k_table[@n % 4][SMALL_PRIMES.map {|q| Utils.kronecker(@n, q) }]
  @kn = @k * @n

  # @f: the size of the factor base
  # @m: the length of the sieve interval (2M+1)
  # @t: the large prime tolerance
  _digit, @f, @m, @t =
    PARAMETERS.find {|d,| Math.log(@n, 10) < d } || PARAMETERS.last

  # compute factor base
  # @ps: the factor base
  # @log_ps: log of each factor base
  # @sqrts: sqrt(kN) mod p
  @ps, @log_ps, @sqrts = [], [], []
  Sieve.each do |p|
    if @n % p == 0
      @trivial_factor = p
      return
    end
    if p >= 3 && Utils.kronecker(@kn, p) == 1
      @ps << p
      @log_ps << (Math.log(p) * S).floor
      @sqrts << Utils.mod_sqrt(@kn, p)
      break if @ps.size >= @f
    end
  end
  @pmax = @ps.last # max factor base

  # for enumerating coefficients
  @base_d = Math.sqrt(Math.sqrt(@kn) / (Math.sqrt(2) * @m)).floor | 3
  @diff_d = 0 # 0, 4, -4, 8, -8, 12, -12, ...

  # sieve table
  @ws = [0] * (2 * @m + 1)

  # a test value for sieve
  @test_value = (Math.log(@m * Math.sqrt(@kn) / @pmax ** @t) * S).round

  # relations found
  @incomplete_relations = {}

  # for gaussian elimination
  @elim = GaussianElimination.new(@ps.size + 2)
end

Class Method Details

.k_table ⇒ Object

Generates a table for finding k such that kN = 1 mod 4 and that kN is rich in small quadratic residues, i.e., (kN | p) = 1 for all p in SMALL_PRIMES.



130
131
132
133
134
135
136
137
138
139
140
141
142
 
# File 'lib/faster_prime/prime_factorization.rb', line 130

def self.k_table
  return @@k_table if defined?(@@k_table)

  k_table = { 1 => {}, 3 => {} }
  Sieve.each do |p|
    next if p <= SMALL_PRIMES.last
    state = SMALL_PRIMES.map {|q| Utils.kronecker(p, q) }
    k_table[p % 4][state] ||= p
    break if k_table.all? {|r, t| t.size == 2 ** SMALL_PRIMES.size }
  end

  @@k_table = k_table
end

.try_find_factor(n) ⇒ Object 



144
145
146
 
# File 'lib/faster_prime/prime_factorization.rb', line 144

def self.try_find_factor(n)
  new(n).try_find_factor
end

Instance Method Details

#check_factor_candidate(relations) ⇒ Object

Computes and checks a factor candidate



336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
 
# File 'lib/faster_prime/prime_factorization.rb', line 336

def check_factor_candidate(relations)
  # computes the factor candidate
  p1 = 1
  p2_vec = Hash.new(0) # XXX: use FactoredInteger
  relations.each do |h, qx_vec|
    p1 = p1 * h % @kn
    qx_vec.each {|p, e| p2_vec[p] += e }
  end
  p2 = 1
  p2_vec.each {|p, e| p2 = (p2 * Utils.mod_pow(p, e / 2, @kn)) % @kn }

  return if p1 == p2 || (p1 + p2) % @kn == 0

  factor = (p1 + p2).gcd(@n)

  # check the factor candidate
  if factor != 1 && factor != @n
    throw :factor_found, factor
  end
end

#check_relation(x) ⇒ Object

Tests whether Q(x) is actually a product of factor bases



286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
 
# File 'lib/faster_prime/prime_factorization.rb', line 286

def check_relation(x)
  h, qx = compute_q(x)
  return if qx == 0
  es, l = exponent_bitvector(qx)

  # discard this x if the residue L is too big
  return if l > @pmax ** @t

  if l == 1
    # complete relation found:
    # Q(x) = p0^e0 * p1^e1 * ... * pk^ek (pi in the factor base)
    qx_vec = Hash.new(0)
    PrimeFactorization.prime_factorization(qx) {|p, e| qx_vec[p] += e }
    collect_relation(es, h, qx_vec)

  elsif @incomplete_relations[l]
    # large prime procedure:
    # make a complete relation by multiplying two incomplete relations
    es2, h2, qx2 = @incomplete_relations[l]

    # XXX: use FactoredInteger
    qx_vec = Hash.new(0)
    PrimeFactorization.prime_factorization(qx) {|p, e| qx_vec[p] += e }
    PrimeFactorization.prime_factorization(qx2) {|p, e| qx_vec[p] += e }

    collect_relation(es ^ es2, h * h2 % @kn, qx_vec)

  else
    @incomplete_relations[l] = [es, h, qx]
  end
end

#collect_relation(es, h, q) ⇒ Object

Adds a complete relation found to gaussian elimination list



319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
 
# File 'lib/faster_prime/prime_factorization.rb', line 319

def collect_relation(es, h, q)
  #p "%0#{ @ps.size + 1 }b" % es
  @elim[es] = [h, q]

  if @elim.size > 0.9 * @f && @elim.size % [1, @f / 100].max == 0
    @elim.eliminate do |relations|
      # factor candidate found
      check_factor_candidate(relations)
    end
  end

  if @elim.size > @f * 1.5
    raise Failed, "failed to find a factor: #{ @n }"
  end
end

#compute_q(x) ⇒ Object

Computes Q(x)



358
359
360
361
362
363
364
365
366
 
# File 'lib/faster_prime/prime_factorization.rb', line 358

def compute_q(x)
  h = (2 * @a * x + @b) * Utils.mod_inv(2 * @d, @kn) % @kn

  q = h * h % @kn
  q -= @kn if @real_root_1 <= x && x <= @real_root_2
  # assert: q == @a*x*x + @b*x + @c

  return h, q
end

#exponent_bitvector(n) ⇒ Object

Factorizes n in factor base and returns an exponent bitvector



369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
 
# File 'lib/faster_prime/prime_factorization.rb', line 369

def exponent_bitvector(n)
  vec = 0
  if n < 0
    vec = 1
    n = -n
  end
  bit = 2
  ([2] + @ps).each do |p|
    e = false
    while n % p == 0
      n /= p
      e = !e
    end
    vec |= bit if e
    bit <<= 1
  end
  return vec, n
end

#select_q ⇒ Object

Selects Q(x) = Ax^2 + Bx + C for quadratic sieve



235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
 
# File 'lib/faster_prime/prime_factorization.rb', line 235

def select_q
  # find a prime D such that (kN | D) = 1 and D = 3 mod 4
  begin
    @d = @base_d + @diff_d
    @diff_d = @diff_d > 0 ? -@diff_d : 4 - @diff_d
  end until Utils.kronecker(@kn, @d) == 1 && PrimalityTest.probable_prime?(@d) != false && !@ps.include?(@d)

  # select coefficients (see the paper)
  @a = @d * @d

  h0 = Utils.mod_pow(@kn, @d / 4, @d)
  h1 = (@kn * h0) % @d
  # assert: h1*h1 = kN mod D
  # assert: h0*h1 = 1 mod D
  h2 = (Utils.mod_inv(2, @d) * h0 * ((@kn - h1 * h1) / @d)) % @d
  # assert: kN - h1^2 = 0 mod D

  @b = (h1 + h2 * @d) % @a
  @b -= @a if @b[0] == 0
  # assert: B^2 = kN mod 4A

  #@c = (@b * @b - @kn) / (4 * @a) # not needed

  # compute real roots of Q(x) = 0
  s = Utils.integer_square_root(@kn)
  @real_root_1 = (-@b - s) / (2 * @a)
  @real_root_2 = (-@b + s) / (2 * @a)
  # x is between real roots if @real_root_1 <= x && x <= @real_root_2 
end

#sieve ⇒ Object

Performs the quadratic sieve



277
278
279
280
281
282
283
 
# File 'lib/faster_prime/prime_factorization.rb', line 277

def sieve
  @ws.fill(0)
  @ps.zip(@log_ps, @roots) do |p, lp, (s1, s2)|
    ((@m + s1) % p).step(2 * @m, p) {|i| @ws[i] += lp }
    ((@m + s2) % p).step(2 * @m, p) {|j| @ws[j] += lp }
  end
end

#solve_roots ⇒ Object

Computes the roots of Q(x) = 0 mod p (p in the factor base)



266
267
268
269
270
271
272
273
274
 
# File 'lib/faster_prime/prime_factorization.rb', line 266

def solve_roots
  @roots = []
  @ps.zip(@sqrts) do |p, s|
    a_inv = Utils.mod_inv(2 * @a, p)
    x1 = (-@b + s) * a_inv % p
    x2 = (-@b - s) * a_inv % p
    @roots << [x1, x2]
  end
end

#try_find_factor ⇒ Object

try to find any factor (not necessarily a prime)



218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
 
# File 'lib/faster_prime/prime_factorization.rb', line 218

def try_find_factor
  return @trivial_factor if defined?(@trivial_factor)

  catch(:factor_found) do
    while true
      p :foo
      select_q
      solve_roots
      sieve
      @ws.each_with_index do |w, i|
        check_relation(i - @m) if w > @test_value
      end
    end
  end
end