Coverage for /usr/lib/python3/dist-packages/sympy/ntheory/qs.py: 9%

1from sympy.core.numbers import igcd, mod_inverse

2from sympy.core.power import integer_nthroot

3from sympy.ntheory.residue_ntheory import _sqrt_mod_prime_power

4from sympy.ntheory import isprime

5from math import log, sqrt

6import random

8rgen = random.Random()

10class SievePolynomial:

11 def __init__(self, modified_coeff=(), a=None, b=None):

12 """This class denotes the seive polynomial.

13 If ``g(x) = (a*x + b)**2 - N``. `g(x)` can be expanded

14 to ``a*x**2 + 2*a*b*x + b**2 - N``, so the coefficient

15 is stored in the form `[a**2, 2*a*b, b**2 - N]`. This

16 ensures faster `eval` method because we dont have to

17 perform `a**2, 2*a*b, b**2` every time we call the

18 `eval` method. As multiplication is more expensive

19 than addition, by using modified_coefficient we get

20 a faster seiving process.

22 Parameters

23 ==========

25 modified_coeff : modified_coefficient of sieve polynomial

26 a : parameter of the sieve polynomial

27 b : parameter of the sieve polynomial

28 """

29 self.modified_coeff = modified_coeff

30 self.a = a

31 self.b = b

33 def eval(self, x):

34 """

35 Compute the value of the sieve polynomial at point x.

37 Parameters

38 ==========

40 x : Integer parameter for sieve polynomial

41 """

42 ans = 0

43 for coeff in self.modified_coeff:

44 ans *= x

45 ans += coeff

46 return ans

49class FactorBaseElem:

50 """This class stores an element of the `factor_base`.

51 """

52 def __init__(self, prime, tmem_p, log_p):

53 """

54 Initialization of factor_base_elem.

56 Parameters

57 ==========

59 prime : prime number of the factor_base

60 tmem_p : Integer square root of x**2 = n mod prime

61 log_p : Compute Natural Logarithm of the prime

62 """

63 self.prime = prime

64 self.tmem_p = tmem_p

65 self.log_p = log_p

66 self.soln1 = None

67 self.soln2 = None

68 self.a_inv = None

69 self.b_ainv = None

72def _generate_factor_base(prime_bound, n):

73 """Generate `factor_base` for Quadratic Sieve. The `factor_base`

74 consists of all the points whose ``legendre_symbol(n, p) == 1``

75 and ``p < num_primes``. Along with the prime `factor_base` also stores

76 natural logarithm of prime and the residue n modulo p.

77 It also returns the of primes numbers in the `factor_base` which are

78 close to 1000 and 5000.

80 Parameters

81 ==========

83 prime_bound : upper prime bound of the factor_base

84 n : integer to be factored

85 """

86 from sympy.ntheory.generate import sieve

87 factor_base = []

88 idx_1000, idx_5000 = None, None

89 for prime in sieve.primerange(1, prime_bound):

90 if pow(n, (prime - 1) // 2, prime) == 1:

91 if prime > 1000 and idx_1000 is None:

92 idx_1000 = len(factor_base) - 1

93 if prime > 5000 and idx_5000 is None:

94 idx_5000 = len(factor_base) - 1

95 residue = _sqrt_mod_prime_power(n, prime, 1)[0]

96 log_p = round(log(prime)*2**10)

97 factor_base.append(FactorBaseElem(prime, residue, log_p))

98 return idx_1000, idx_5000, factor_base

100

101def _initialize_first_polynomial(N, M, factor_base, idx_1000, idx_5000, seed=None):

102 """This step is the initialization of the 1st sieve polynomial.

103 Here `a` is selected as a product of several primes of the factor_base

104 such that `a` is about to ``sqrt(2*N) / M``. Other initial values of

105 factor_base elem are also initialized which includes a_inv, b_ainv, soln1,

106 soln2 which are used when the sieve polynomial is changed. The b_ainv

107 is required for fast polynomial change as we do not have to calculate

108 `2*b*mod_inverse(a, prime)` every time.

109 We also ensure that the `factor_base` primes which make `a` are between

110 1000 and 5000.

111

112 Parameters

113 ==========

114

115 N : Number to be factored

116 M : sieve interval

117 factor_base : factor_base primes

118 idx_1000 : index of prime number in the factor_base near 1000

119 idx_5000 : index of prime number in the factor_base near to 5000

120 seed : Generate pseudoprime numbers

121 """

122 if seed is not None:

123 rgen.seed(seed)

124 approx_val = sqrt(2*N) / M

125 # `a` is a parameter of the sieve polynomial and `q` is the prime factors of `a`

126 # randomly search for a combination of primes whose multiplication is close to approx_val

127 # This multiplication of primes will be `a` and the primes will be `q`

128 # `best_a` denotes that `a` is close to approx_val in the random search of combination

129 best_a, best_q, best_ratio = None, None, None

130 start = 0 if idx_1000 is None else idx_1000

131 end = len(factor_base) - 1 if idx_5000 is None else idx_5000

132 for _ in range(50):

133 a = 1

134 q = []

135 while(a < approx_val):

136 rand_p = 0

137 while(rand_p == 0 or rand_p in q):

138 rand_p = rgen.randint(start, end)

139 p = factor_base[rand_p].prime

140 a *= p

141 q.append(rand_p)

142 ratio = a / approx_val

143 if best_ratio is None or abs(ratio - 1) < abs(best_ratio - 1):

144 best_q = q

145 best_a = a

146 best_ratio = ratio

147

148 a = best_a

149 q = best_q

150

151 B = []

152 for idx, val in enumerate(q):

153 q_l = factor_base[val].prime

154 gamma = factor_base[val].tmem_p * mod_inverse(a // q_l, q_l) % q_l

155 if gamma > q_l / 2:

156 gamma = q_l - gamma

157 B.append(a//q_l*gamma)

158

159 b = sum(B)

160 g = SievePolynomial([a*a, 2*a*b, b*b - N], a, b)

161

162 for fb in factor_base:

163 if a % fb.prime == 0:

164 continue

165 fb.a_inv = mod_inverse(a, fb.prime)

166 fb.b_ainv = [2*b_elem*fb.a_inv % fb.prime for b_elem in B]

167 fb.soln1 = (fb.a_inv*(fb.tmem_p - b)) % fb.prime

168 fb.soln2 = (fb.a_inv*(-fb.tmem_p - b)) % fb.prime

169 return g, B

170

171

172def _initialize_ith_poly(N, factor_base, i, g, B):

173 """Initialization stage of ith poly. After we finish sieving 1`st polynomial

174 here we quickly change to the next polynomial from which we will again

175 start sieving. Suppose we generated ith sieve polynomial and now we

176 want to generate (i + 1)th polynomial, where ``1 <= i <= 2**(j - 1) - 1``

177 where `j` is the number of prime factors of the coefficient `a`

178 then this function can be used to go to the next polynomial. If

179 ``i = 2**(j - 1) - 1`` then go to _initialize_first_polynomial stage.

180

181 Parameters

182 ==========

183

184 N : number to be factored

185 factor_base : factor_base primes

186 i : integer denoting ith polynomial

187 g : (i - 1)th polynomial

188 B : array that stores a//q_l*gamma

189 """

190 from sympy.functions.elementary.integers import ceiling

191 v = 1

192 j = i

193 while(j % 2 == 0):

194 v += 1

195 j //= 2

196 if ceiling(i / (2**v)) % 2 == 1:

197 neg_pow = -1

198 else:

199 neg_pow = 1

200 b = g.b + 2*neg_pow*B[v - 1]

201 a = g.a

202 g = SievePolynomial([a*a, 2*a*b, b*b - N], a, b)

203 for fb in factor_base:

204 if a % fb.prime == 0:

205 continue

206 fb.soln1 = (fb.soln1 - neg_pow*fb.b_ainv[v - 1]) % fb.prime

207 fb.soln2 = (fb.soln2 - neg_pow*fb.b_ainv[v - 1]) % fb.prime

208

209 return g

210

211

212def _gen_sieve_array(M, factor_base):

213 """Sieve Stage of the Quadratic Sieve. For every prime in the factor_base

214 that does not divide the coefficient `a` we add log_p over the sieve_array

215 such that ``-M <= soln1 + i*p <= M`` and ``-M <= soln2 + i*p <= M`` where `i`

216 is an integer. When p = 2 then log_p is only added using

217 ``-M <= soln1 + i*p <= M``.

218

219 Parameters

220 ==========

221

222 M : sieve interval

223 factor_base : factor_base primes

224 """

225 sieve_array = [0]*(2*M + 1)

226 for factor in factor_base:

227 if factor.soln1 is None: #The prime does not divides a

228 continue

229 for idx in range((M + factor.soln1) % factor.prime, 2*M, factor.prime):

230 sieve_array[idx] += factor.log_p

231 if factor.prime == 2:

232 continue

233 #if prime is 2 then sieve only with soln_1_p

234 for idx in range((M + factor.soln2) % factor.prime, 2*M, factor.prime):

235 sieve_array[idx] += factor.log_p

236 return sieve_array

237

238

239def _check_smoothness(num, factor_base):

240 """Here we check that if `num` is a smooth number or not. If `a` is a smooth

241 number then it returns a vector of prime exponents modulo 2. For example

242 if a = 2 * 5**2 * 7**3 and the factor base contains {2, 3, 5, 7} then

243 `a` is a smooth number and this function returns ([1, 0, 0, 1], True). If

244 `a` is a partial relation which means that `a` a has one prime factor

245 greater than the `factor_base` then it returns `(a, False)` which denotes `a`

246 is a partial relation.

247

248 Parameters

249 ==========

250

251 a : integer whose smootheness is to be checked

252 factor_base : factor_base primes

253 """

254 vec = []

255 if num < 0:

256 vec.append(1)

257 num *= -1

258 else:

259 vec.append(0)

260 #-1 is not included in factor_base add -1 in vector

261 for factor in factor_base:

262 if num % factor.prime != 0:

263 vec.append(0)

264 continue

265 factor_exp = 0

266 while num % factor.prime == 0:

267 factor_exp += 1

268 num //= factor.prime

269 vec.append(factor_exp % 2)

270 if num == 1:

271 return vec, True

272 if isprime(num):

273 return num, False

274 return None, None

275

276

277def _trial_division_stage(N, M, factor_base, sieve_array, sieve_poly, partial_relations, ERROR_TERM):

278 """Trial division stage. Here we trial divide the values generetated

279 by sieve_poly in the sieve interval and if it is a smooth number then

280 it is stored in `smooth_relations`. Moreover, if we find two partial relations

281 with same large prime then they are combined to form a smooth relation.

282 First we iterate over sieve array and look for values which are greater

283 than accumulated_val, as these values have a high chance of being smooth

284 number. Then using these values we find smooth relations.

285 In general, let ``t**2 = u*p modN`` and ``r**2 = v*p modN`` be two partial relations

286 with the same large prime p. Then they can be combined ``(t*r/p)**2 = u*v modN``

287 to form a smooth relation.

288

289 Parameters

290 ==========

291

292 N : Number to be factored

293 M : sieve interval

294 factor_base : factor_base primes

295 sieve_array : stores log_p values

296 sieve_poly : polynomial from which we find smooth relations

297 partial_relations : stores partial relations with one large prime

298 ERROR_TERM : error term for accumulated_val

299 """

300 sqrt_n = sqrt(float(N))

301 accumulated_val = log(M * sqrt_n)*2**10 - ERROR_TERM

302 smooth_relations = []

303 proper_factor = set()

304 partial_relation_upper_bound = 128*factor_base[-1].prime

305 for idx, val in enumerate(sieve_array):

306 if val < accumulated_val:

307 continue

308 x = idx - M

309 v = sieve_poly.eval(x)

310 vec, is_smooth = _check_smoothness(v, factor_base)

311 if is_smooth is None:#Neither smooth nor partial

312 continue

313 u = sieve_poly.a*x + sieve_poly.b

314 # Update the partial relation

315 # If 2 partial relation with same large prime is found then generate smooth relation

316 if is_smooth is False:#partial relation found

317 large_prime = vec

318 #Consider the large_primes under 128*F

319 if large_prime > partial_relation_upper_bound:

320 continue

321 if large_prime not in partial_relations:

322 partial_relations[large_prime] = (u, v)

323 continue

324 else:

325 u_prev, v_prev = partial_relations[large_prime]

326 partial_relations.pop(large_prime)

327 try:

328 large_prime_inv = mod_inverse(large_prime, N)

329 except ValueError:#if large_prine divides N

330 proper_factor.add(large_prime)

331 continue

332 u = u*u_prev*large_prime_inv

333 v = v*v_prev // (large_prime*large_prime)

334 vec, is_smooth = _check_smoothness(v, factor_base)

335 #assert u*u % N == v % N

336 smooth_relations.append((u, v, vec))

337 return smooth_relations, proper_factor

338

339

340#LINEAR ALGEBRA STAGE

341def _build_matrix(smooth_relations):

342 """Build a 2D matrix from smooth relations.

343

344 Parameters

345 ==========

346

347 smooth_relations : Stores smooth relations

348 """

349 matrix = []

350 for s_relation in smooth_relations:

351 matrix.append(s_relation[2])

352 return matrix

353

354

355def _gauss_mod_2(A):

356 """Fast gaussian reduction for modulo 2 matrix.

357

358 Parameters

359 ==========

360

361 A : Matrix

362

363 Examples

364 ========

365

366 >>> from sympy.ntheory.qs import _gauss_mod_2

367 >>> _gauss_mod_2([[0, 1, 1], [1, 0, 1], [0, 1, 0], [1, 1, 1]])

368 ([[[1, 0, 1], 3]],

369 [True, True, True, False],

370 [[0, 1, 0], [1, 0, 0], [0, 0, 1], [1, 0, 1]])

371

372 Reference

373 ==========

374

375 .. [1] A fast algorithm for gaussian elimination over GF(2) and

376 its implementation on the GAPP. Cetin K.Koc, Sarath N.Arachchige"""

377 import copy

378 matrix = copy.deepcopy(A)

379 row = len(matrix)

380 col = len(matrix[0])

381 mark = [False]*row

382 for c in range(col):

383 for r in range(row):

384 if matrix[r][c] == 1:

385 break

386 mark[r] = True

387 for c1 in range(col):

388 if c1 == c:

389 continue

390 if matrix[r][c1] == 1:

391 for r2 in range(row):

392 matrix[r2][c1] = (matrix[r2][c1] + matrix[r2][c]) % 2

393 dependent_row = []

394 for idx, val in enumerate(mark):

395 if val == False:

396 dependent_row.append([matrix[idx], idx])

397 return dependent_row, mark, matrix

398

399

400def _find_factor(dependent_rows, mark, gauss_matrix, index, smooth_relations, N):

401 """Finds proper factor of N. Here, transform the dependent rows as a

402 combination of independent rows of the gauss_matrix to form the desired

403 relation of the form ``X**2 = Y**2 modN``. After obtaining the desired relation

404 we obtain a proper factor of N by `gcd(X - Y, N)`.

405

406 Parameters

407 ==========

408

409 dependent_rows : denoted dependent rows in the reduced matrix form

410 mark : boolean array to denoted dependent and independent rows

411 gauss_matrix : Reduced form of the smooth relations matrix

412 index : denoted the index of the dependent_rows

413 smooth_relations : Smooth relations vectors matrix

414 N : Number to be factored

415 """

416 idx_in_smooth = dependent_rows[index][1]

417 independent_u = [smooth_relations[idx_in_smooth][0]]

418 independent_v = [smooth_relations[idx_in_smooth][1]]

419 dept_row = dependent_rows[index][0]

420

421 for idx, val in enumerate(dept_row):

422 if val == 1:

423 for row in range(len(gauss_matrix)):

424 if gauss_matrix[row][idx] == 1 and mark[row] == True:

425 independent_u.append(smooth_relations[row][0])

426 independent_v.append(smooth_relations[row][1])

427 break

428

429 u = 1

430 v = 1

431 for i in independent_u:

432 u *= i

433 for i in independent_v:

434 v *= i

435 #assert u**2 % N == v % N

436 v = integer_nthroot(v, 2)[0]

437 return igcd(u - v, N)

438

439

440def qs(N, prime_bound, M, ERROR_TERM=25, seed=1234):

441 """Performs factorization using Self-Initializing Quadratic Sieve.

442 In SIQS, let N be a number to be factored, and this N should not be a

443 perfect power. If we find two integers such that ``X**2 = Y**2 modN`` and

444 ``X != +-Y modN``, then `gcd(X + Y, N)` will reveal a proper factor of N.

445 In order to find these integers X and Y we try to find relations of form

446 t**2 = u modN where u is a product of small primes. If we have enough of

447 these relations then we can form ``(t1*t2...ti)**2 = u1*u2...ui modN`` such that

448 the right hand side is a square, thus we found a relation of ``X**2 = Y**2 modN``.

449

450 Here, several optimizations are done like using multiple polynomials for

451 sieving, fast changing between polynomials and using partial relations.

452 The use of partial relations can speeds up the factoring by 2 times.

453

454 Parameters

455 ==========

456

457 N : Number to be Factored

458 prime_bound : upper bound for primes in the factor base

459 M : Sieve Interval

460 ERROR_TERM : Error term for checking smoothness

461 threshold : Extra smooth relations for factorization

462 seed : generate pseudo prime numbers

463

464 Examples

465 ========

466

467 >>> from sympy.ntheory import qs

468 >>> qs(25645121643901801, 2000, 10000)

469 {5394769, 4753701529}

470 >>> qs(9804659461513846513, 2000, 10000)

471 {4641991, 2112166839943}

472

473 References

474 ==========

475

476 .. [1] https://pdfs.semanticscholar.org/5c52/8a975c1405bd35c65993abf5a4edb667c1db.pdf

477 .. [2] https://www.rieselprime.de/ziki/Self-initializing_quadratic_sieve

478 """

479 ERROR_TERM*=2**10

480 rgen.seed(seed)

481 idx_1000, idx_5000, factor_base = _generate_factor_base(prime_bound, N)

482 smooth_relations = []

483 ith_poly = 0

484 partial_relations = {}

485 proper_factor = set()

486 threshold = 5*len(factor_base) // 100

487 while True:

488 if ith_poly == 0:

489 ith_sieve_poly, B_array = _initialize_first_polynomial(N, M, factor_base, idx_1000, idx_5000)

490 else:

491 ith_sieve_poly = _initialize_ith_poly(N, factor_base, ith_poly, ith_sieve_poly, B_array)

492 ith_poly += 1

493 if ith_poly >= 2**(len(B_array) - 1): # time to start with a new sieve polynomial

494 ith_poly = 0

495 sieve_array = _gen_sieve_array(M, factor_base)

496 s_rel, p_f = _trial_division_stage(N, M, factor_base, sieve_array, ith_sieve_poly, partial_relations, ERROR_TERM)

497 smooth_relations += s_rel

498 proper_factor |= p_f

499 if len(smooth_relations) >= len(factor_base) + threshold:

500 break

501 matrix = _build_matrix(smooth_relations)

502 dependent_row, mark, gauss_matrix = _gauss_mod_2(matrix)

503 N_copy = N

504 for index in range(len(dependent_row)):

505 factor = _find_factor(dependent_row, mark, gauss_matrix, index, smooth_relations, N)

506 if factor > 1 and factor < N:

507 proper_factor.add(factor)

508 while(N_copy % factor == 0):

509 N_copy //= factor

510 if isprime(N_copy):

511 proper_factor.add(N_copy)

512 break

513 if(N_copy == 1):

514 break

515 return proper_factor