Coverage for /usr/lib/python3/dist-packages/sympy/polys/factortools.py: 7%

1"""Polynomial factorization routines in characteristic zero. """

3from sympy.core.random import _randint

5from sympy.polys.galoistools import (

6 gf_from_int_poly, gf_to_int_poly,

7 gf_lshift, gf_add_mul, gf_mul,

8 gf_div, gf_rem,

9 gf_gcdex,

10 gf_sqf_p,

11 gf_factor_sqf, gf_factor)

13from sympy.polys.densebasic import (

14 dup_LC, dmp_LC, dmp_ground_LC,

15 dup_TC,

16 dup_convert, dmp_convert,

17 dup_degree, dmp_degree,

18 dmp_degree_in, dmp_degree_list,

19 dmp_from_dict,

20 dmp_zero_p,

21 dmp_one,

22 dmp_nest, dmp_raise,

23 dup_strip,

24 dmp_ground,

25 dup_inflate,

26 dmp_exclude, dmp_include,

27 dmp_inject, dmp_eject,

28 dup_terms_gcd, dmp_terms_gcd)

30from sympy.polys.densearith import (

31 dup_neg, dmp_neg,

32 dup_add, dmp_add,

33 dup_sub, dmp_sub,

34 dup_mul, dmp_mul,

35 dup_sqr,

36 dmp_pow,

37 dup_div, dmp_div,

38 dup_quo, dmp_quo,

39 dmp_expand,

40 dmp_add_mul,

41 dup_sub_mul, dmp_sub_mul,

42 dup_lshift,

43 dup_max_norm, dmp_max_norm,

44 dup_l1_norm,

45 dup_mul_ground, dmp_mul_ground,

46 dup_quo_ground, dmp_quo_ground)

48from sympy.polys.densetools import (

49 dup_clear_denoms, dmp_clear_denoms,

50 dup_trunc, dmp_ground_trunc,

51 dup_content,

52 dup_monic, dmp_ground_monic,

53 dup_primitive, dmp_ground_primitive,

54 dmp_eval_tail,

55 dmp_eval_in, dmp_diff_eval_in,

56 dmp_compose,

57 dup_shift, dup_mirror)

59from sympy.polys.euclidtools import (

60 dmp_primitive,

61 dup_inner_gcd, dmp_inner_gcd)

63from sympy.polys.sqfreetools import (

64 dup_sqf_p,

65 dup_sqf_norm, dmp_sqf_norm,

66 dup_sqf_part, dmp_sqf_part)

68from sympy.polys.polyutils import _sort_factors

69from sympy.polys.polyconfig import query

71from sympy.polys.polyerrors import (

72 ExtraneousFactors, DomainError, CoercionFailed, EvaluationFailed)

74from sympy.utilities import subsets

76from math import ceil as _ceil, log as _log

79def dup_trial_division(f, factors, K):

80 """

81 Determine multiplicities of factors for a univariate polynomial

82 using trial division.

83 """

84 result = []

86 for factor in factors:

87 k = 0

89 while True:

90 q, r = dup_div(f, factor, K)

92 if not r:

93 f, k = q, k + 1

94 else:

95 break

97 result.append((factor, k))

99 return _sort_factors(result)

100

101

102def dmp_trial_division(f, factors, u, K):

103 """

104 Determine multiplicities of factors for a multivariate polynomial

105 using trial division.

106 """

107 result = []

108

109 for factor in factors:

110 k = 0

111

112 while True:

113 q, r = dmp_div(f, factor, u, K)

114

115 if dmp_zero_p(r, u):

116 f, k = q, k + 1

117 else:

118 break

119

120 result.append((factor, k))

121

122 return _sort_factors(result)

123

124

125def dup_zz_mignotte_bound(f, K):

126 """

127 The Knuth-Cohen variant of Mignotte bound for

128 univariate polynomials in `K[x]`.

129

130 Examples

131 ========

132

133 >>> from sympy.polys import ring, ZZ

134 >>> R, x = ring("x", ZZ)

135

136 >>> f = x**3 + 14*x**2 + 56*x + 64

137 >>> R.dup_zz_mignotte_bound(f)

138 152

139

140 By checking `factor(f)` we can see that max coeff is 8

141

142 Also consider a case that `f` is irreducible for example `f = 2*x**2 + 3*x + 4`

143 To avoid a bug for these cases, we return the bound plus the max coefficient of `f`

144

145 >>> f = 2*x**2 + 3*x + 4

146 >>> R.dup_zz_mignotte_bound(f)

147 6

148

149 Lastly,To see the difference between the new and the old Mignotte bound

150 consider the irreducible polynomial::

151

152 >>> f = 87*x**7 + 4*x**6 + 80*x**5 + 17*x**4 + 9*x**3 + 12*x**2 + 49*x + 26

153 >>> R.dup_zz_mignotte_bound(f)

154 744

155

156 The new Mignotte bound is 744 whereas the old one (SymPy 1.5.1) is 1937664.

157

158

159 References

160 ==========

161

162 ..[1] [Abbott2013]_

163

164 """

165 from sympy.functions.combinatorial.factorials import binomial

166 d = dup_degree(f)

167 delta = _ceil(d / 2)

168 delta2 = _ceil(delta / 2)

169

170 # euclidean-norm

171 eucl_norm = K.sqrt( sum( [cf**2 for cf in f] ) )

172

173 # biggest values of binomial coefficients (p. 538 of reference)

174 t1 = binomial(delta - 1, delta2)

175 t2 = binomial(delta - 1, delta2 - 1)

176

177 lc = K.abs(dup_LC(f, K)) # leading coefficient

178 bound = t1 * eucl_norm + t2 * lc # (p. 538 of reference)

179 bound += dup_max_norm(f, K) # add max coeff for irreducible polys

180 bound = _ceil(bound / 2) * 2 # round up to even integer

181

182 return bound

183

184def dmp_zz_mignotte_bound(f, u, K):

185 """Mignotte bound for multivariate polynomials in `K[X]`. """

186 a = dmp_max_norm(f, u, K)

187 b = abs(dmp_ground_LC(f, u, K))

188 n = sum(dmp_degree_list(f, u))

189

190 return K.sqrt(K(n + 1))*2**n*a*b

191

192

193def dup_zz_hensel_step(m, f, g, h, s, t, K):

194 """

195 One step in Hensel lifting in `Z[x]`.

196

197 Given positive integer `m` and `Z[x]` polynomials `f`, `g`, `h`, `s`

198 and `t` such that::

199

200 f = g*h (mod m)

201 s*g + t*h = 1 (mod m)

202

203 lc(f) is not a zero divisor (mod m)

204 lc(h) = 1

205

206 deg(f) = deg(g) + deg(h)

207 deg(s) < deg(h)

208 deg(t) < deg(g)

209

210 returns polynomials `G`, `H`, `S` and `T`, such that::

211

212 f = G*H (mod m**2)

213 S*G + T*H = 1 (mod m**2)

214

215 References

216 ==========

217

218 .. [1] [Gathen99]_

219

220 """

221 M = m**2

222

223 e = dup_sub_mul(f, g, h, K)

224 e = dup_trunc(e, M, K)

225

226 q, r = dup_div(dup_mul(s, e, K), h, K)

227

228 q = dup_trunc(q, M, K)

229 r = dup_trunc(r, M, K)

230

231 u = dup_add(dup_mul(t, e, K), dup_mul(q, g, K), K)

232 G = dup_trunc(dup_add(g, u, K), M, K)

233 H = dup_trunc(dup_add(h, r, K), M, K)

234

235 u = dup_add(dup_mul(s, G, K), dup_mul(t, H, K), K)

236 b = dup_trunc(dup_sub(u, [K.one], K), M, K)

237

238 c, d = dup_div(dup_mul(s, b, K), H, K)

239

240 c = dup_trunc(c, M, K)

241 d = dup_trunc(d, M, K)

242

243 u = dup_add(dup_mul(t, b, K), dup_mul(c, G, K), K)

244 S = dup_trunc(dup_sub(s, d, K), M, K)

245 T = dup_trunc(dup_sub(t, u, K), M, K)

246

247 return G, H, S, T

248

249

250def dup_zz_hensel_lift(p, f, f_list, l, K):

251 r"""

252 Multifactor Hensel lifting in `Z[x]`.

253

254 Given a prime `p`, polynomial `f` over `Z[x]` such that `lc(f)`

255 is a unit modulo `p`, monic pair-wise coprime polynomials `f_i`

256 over `Z[x]` satisfying::

257

258 f = lc(f) f_1 ... f_r (mod p)

259

260 and a positive integer `l`, returns a list of monic polynomials

261 `F_1,\ F_2,\ \dots,\ F_r` satisfying::

262

263 f = lc(f) F_1 ... F_r (mod p**l)

264

265 F_i = f_i (mod p), i = 1..r

266

267 References

268 ==========

269

270 .. [1] [Gathen99]_

271

272 """

273 r = len(f_list)

274 lc = dup_LC(f, K)

275

276 if r == 1:

277 F = dup_mul_ground(f, K.gcdex(lc, p**l)[0], K)

278 return [ dup_trunc(F, p**l, K) ]

279

280 m = p

281 k = r // 2

282 d = int(_ceil(_log(l, 2)))

283

284 g = gf_from_int_poly([lc], p)

285

286 for f_i in f_list[:k]:

287 g = gf_mul(g, gf_from_int_poly(f_i, p), p, K)

288

289 h = gf_from_int_poly(f_list[k], p)

290

291 for f_i in f_list[k + 1:]:

292 h = gf_mul(h, gf_from_int_poly(f_i, p), p, K)

293

294 s, t, _ = gf_gcdex(g, h, p, K)

295

296 g = gf_to_int_poly(g, p)

297 h = gf_to_int_poly(h, p)

298 s = gf_to_int_poly(s, p)

299 t = gf_to_int_poly(t, p)

300

301 for _ in range(1, d + 1):

302 (g, h, s, t), m = dup_zz_hensel_step(m, f, g, h, s, t, K), m**2

303

304 return dup_zz_hensel_lift(p, g, f_list[:k], l, K) \

305 + dup_zz_hensel_lift(p, h, f_list[k:], l, K)

306

307def _test_pl(fc, q, pl):

308 if q > pl // 2:

309 q = q - pl

310 if not q:

311 return True

312 return fc % q == 0

313

314def dup_zz_zassenhaus(f, K):

315 """Factor primitive square-free polynomials in `Z[x]`. """

316 n = dup_degree(f)

317

318 if n == 1:

319 return [f]

320

321 from sympy.ntheory import isprime

322

323 fc = f[-1]

324 A = dup_max_norm(f, K)

325 b = dup_LC(f, K)

326 B = int(abs(K.sqrt(K(n + 1))*2**n*A*b))

327 C = int((n + 1)**(2*n)*A**(2*n - 1))

328 gamma = int(_ceil(2*_log(C, 2)))

329 bound = int(2*gamma*_log(gamma))

330 a = []

331 # choose a prime number `p` such that `f` be square free in Z_p

332 # if there are many factors in Z_p, choose among a few different `p`

333 # the one with fewer factors

334 for px in range(3, bound + 1):

335 if not isprime(px) or b % px == 0:

336 continue

337

338 px = K.convert(px)

339

340 F = gf_from_int_poly(f, px)

341

342 if not gf_sqf_p(F, px, K):

343 continue

344 fsqfx = gf_factor_sqf(F, px, K)[1]

345 a.append((px, fsqfx))

346 if len(fsqfx) < 15 or len(a) > 4:

347 break

348 p, fsqf = min(a, key=lambda x: len(x[1]))

349

350 l = int(_ceil(_log(2*B + 1, p)))

351

352 modular = [gf_to_int_poly(ff, p) for ff in fsqf]

353

354 g = dup_zz_hensel_lift(p, f, modular, l, K)

355

356 sorted_T = range(len(g))

357 T = set(sorted_T)

358 factors, s = [], 1

359 pl = p**l

360

361 while 2*s <= len(T):

362 for S in subsets(sorted_T, s):

363 # lift the constant coefficient of the product `G` of the factors

364 # in the subset `S`; if it is does not divide `fc`, `G` does

365 # not divide the input polynomial

366

367 if b == 1:

368 q = 1

369 for i in S:

370 q = q*g[i][-1]

371 q = q % pl

372 if not _test_pl(fc, q, pl):

373 continue

374 else:

375 G = [b]

376 for i in S:

377 G = dup_mul(G, g[i], K)

378 G = dup_trunc(G, pl, K)

379 G = dup_primitive(G, K)[1]

380 q = G[-1]

381 if q and fc % q != 0:

382 continue

383

384 H = [b]

385 S = set(S)

386 T_S = T - S

387

388 if b == 1:

389 G = [b]

390 for i in S:

391 G = dup_mul(G, g[i], K)

392 G = dup_trunc(G, pl, K)

393

394 for i in T_S:

395 H = dup_mul(H, g[i], K)

396

397 H = dup_trunc(H, pl, K)

398

399 G_norm = dup_l1_norm(G, K)

400 H_norm = dup_l1_norm(H, K)

401

402 if G_norm*H_norm <= B:

403 T = T_S

404 sorted_T = [i for i in sorted_T if i not in S]

405

406 G = dup_primitive(G, K)[1]

407 f = dup_primitive(H, K)[1]

408

409 factors.append(G)

410 b = dup_LC(f, K)

411

412 break

413 else:

414 s += 1

415

416 return factors + [f]

417

418

419def dup_zz_irreducible_p(f, K):

420 """Test irreducibility using Eisenstein's criterion. """

421 lc = dup_LC(f, K)

422 tc = dup_TC(f, K)

423

424 e_fc = dup_content(f[1:], K)

425

426 if e_fc:

427 from sympy.ntheory import factorint

428 e_ff = factorint(int(e_fc))

429

430 for p in e_ff.keys():

431 if (lc % p) and (tc % p**2):

432 return True

433

434

435def dup_cyclotomic_p(f, K, irreducible=False):

436 """

437 Efficiently test if ``f`` is a cyclotomic polynomial.

438

439 Examples

440 ========

441

442 >>> from sympy.polys import ring, ZZ

443 >>> R, x = ring("x", ZZ)

444

445 >>> f = x**16 + x**14 - x**10 + x**8 - x**6 + x**2 + 1

446 >>> R.dup_cyclotomic_p(f)

447 False

448

449 >>> g = x**16 + x**14 - x**10 - x**8 - x**6 + x**2 + 1

450 >>> R.dup_cyclotomic_p(g)

451 True

452

453 References

454 ==========

455

456 Bradford, Russell J., and James H. Davenport. "Effective tests for

457 cyclotomic polynomials." In International Symposium on Symbolic and

458 Algebraic Computation, pp. 244-251. Springer, Berlin, Heidelberg, 1988.

459

460 """

461 if K.is_QQ:

462 try:

463 K0, K = K, K.get_ring()

464 f = dup_convert(f, K0, K)

465 except CoercionFailed:

466 return False

467 elif not K.is_ZZ:

468 return False

469

470 lc = dup_LC(f, K)

471 tc = dup_TC(f, K)

472

473 if lc != 1 or (tc != -1 and tc != 1):

474 return False

475

476 if not irreducible:

477 coeff, factors = dup_factor_list(f, K)

478

479 if coeff != K.one or factors != [(f, 1)]:

480 return False

481

482 n = dup_degree(f)

483 g, h = [], []

484

485 for i in range(n, -1, -2):

486 g.insert(0, f[i])

487

488 for i in range(n - 1, -1, -2):

489 h.insert(0, f[i])

490

491 g = dup_sqr(dup_strip(g), K)

492 h = dup_sqr(dup_strip(h), K)

493

494 F = dup_sub(g, dup_lshift(h, 1, K), K)

495

496 if K.is_negative(dup_LC(F, K)):

497 F = dup_neg(F, K)

498

499 if F == f:

500 return True

501

502 g = dup_mirror(f, K)

503

504 if K.is_negative(dup_LC(g, K)):

505 g = dup_neg(g, K)

506

507 if F == g and dup_cyclotomic_p(g, K):

508 return True

509

510 G = dup_sqf_part(F, K)

511

512 if dup_sqr(G, K) == F and dup_cyclotomic_p(G, K):

513 return True

514

515 return False

516

517

518def dup_zz_cyclotomic_poly(n, K):

519 """Efficiently generate n-th cyclotomic polynomial. """

520 from sympy.ntheory import factorint

521 h = [K.one, -K.one]

522

523 for p, k in factorint(n).items():

524 h = dup_quo(dup_inflate(h, p, K), h, K)

525 h = dup_inflate(h, p**(k - 1), K)

526

527 return h

528

529

530def _dup_cyclotomic_decompose(n, K):

531 from sympy.ntheory import factorint

532

533 H = [[K.one, -K.one]]

534

535 for p, k in factorint(n).items():

536 Q = [ dup_quo(dup_inflate(h, p, K), h, K) for h in H ]

537 H.extend(Q)

538

539 for i in range(1, k):

540 Q = [ dup_inflate(q, p, K) for q in Q ]

541 H.extend(Q)

542

543 return H

544

545

546def dup_zz_cyclotomic_factor(f, K):

547 """

548 Efficiently factor polynomials `x**n - 1` and `x**n + 1` in `Z[x]`.

549

550 Given a univariate polynomial `f` in `Z[x]` returns a list of factors

551 of `f`, provided that `f` is in the form `x**n - 1` or `x**n + 1` for

552 `n >= 1`. Otherwise returns None.

553

554 Factorization is performed using cyclotomic decomposition of `f`,

555 which makes this method much faster that any other direct factorization

556 approach (e.g. Zassenhaus's).

557

558 References

559 ==========

560

561 .. [1] [Weisstein09]_

562

563 """

564 lc_f, tc_f = dup_LC(f, K), dup_TC(f, K)

565

566 if dup_degree(f) <= 0:

567 return None

568

569 if lc_f != 1 or tc_f not in [-1, 1]:

570 return None

571

572 if any(bool(cf) for cf in f[1:-1]):

573 return None

574

575 n = dup_degree(f)

576 F = _dup_cyclotomic_decompose(n, K)

577

578 if not K.is_one(tc_f):

579 return F

580 else:

581 H = []

582

583 for h in _dup_cyclotomic_decompose(2*n, K):

584 if h not in F:

585 H.append(h)

586

587 return H

588

589

590def dup_zz_factor_sqf(f, K):

591 """Factor square-free (non-primitive) polynomials in `Z[x]`. """

592 cont, g = dup_primitive(f, K)

593

594 n = dup_degree(g)

595

596 if dup_LC(g, K) < 0:

597 cont, g = -cont, dup_neg(g, K)

598

599 if n <= 0:

600 return cont, []

601 elif n == 1:

602 return cont, [g]

603

604 if query('USE_IRREDUCIBLE_IN_FACTOR'):

605 if dup_zz_irreducible_p(g, K):

606 return cont, [g]

607

608 factors = None

609

610 if query('USE_CYCLOTOMIC_FACTOR'):

611 factors = dup_zz_cyclotomic_factor(g, K)

612

613 if factors is None:

614 factors = dup_zz_zassenhaus(g, K)

615

616 return cont, _sort_factors(factors, multiple=False)

617

618

619def dup_zz_factor(f, K):

620 """

621 Factor (non square-free) polynomials in `Z[x]`.

622

623 Given a univariate polynomial `f` in `Z[x]` computes its complete

624 factorization `f_1, ..., f_n` into irreducibles over integers::

625

626 f = content(f) f_1**k_1 ... f_n**k_n

627

628 The factorization is computed by reducing the input polynomial

629 into a primitive square-free polynomial and factoring it using

630 Zassenhaus algorithm. Trial division is used to recover the

631 multiplicities of factors.

632

633 The result is returned as a tuple consisting of::

634

635 (content(f), [(f_1, k_1), ..., (f_n, k_n))

636

637 Examples

638 ========

639

640 Consider the polynomial `f = 2*x**4 - 2`::

641

642 >>> from sympy.polys import ring, ZZ

643 >>> R, x = ring("x", ZZ)

644

645 >>> R.dup_zz_factor(2*x**4 - 2)

646 (2, [(x - 1, 1), (x + 1, 1), (x**2 + 1, 1)])

647

648 In result we got the following factorization::

649

650 f = 2 (x - 1) (x + 1) (x**2 + 1)

651

652 Note that this is a complete factorization over integers,

653 however over Gaussian integers we can factor the last term.

654

655 By default, polynomials `x**n - 1` and `x**n + 1` are factored

656 using cyclotomic decomposition to speedup computations. To

657 disable this behaviour set cyclotomic=False.

658

659 References

660 ==========

661

662 .. [1] [Gathen99]_

663

664 """

665 cont, g = dup_primitive(f, K)

666

667 n = dup_degree(g)

668

669 if dup_LC(g, K) < 0:

670 cont, g = -cont, dup_neg(g, K)

671

672 if n <= 0:

673 return cont, []

674 elif n == 1:

675 return cont, [(g, 1)]

676

677 if query('USE_IRREDUCIBLE_IN_FACTOR'):

678 if dup_zz_irreducible_p(g, K):

679 return cont, [(g, 1)]

680

681 g = dup_sqf_part(g, K)

682 H = None

683

684 if query('USE_CYCLOTOMIC_FACTOR'):

685 H = dup_zz_cyclotomic_factor(g, K)

686

687 if H is None:

688 H = dup_zz_zassenhaus(g, K)

689

690 factors = dup_trial_division(f, H, K)

691 return cont, factors

692

693

694def dmp_zz_wang_non_divisors(E, cs, ct, K):

695 """Wang/EEZ: Compute a set of valid divisors. """

696 result = [ cs*ct ]

697

698 for q in E:

699 q = abs(q)

700

701 for r in reversed(result):

702 while r != 1:

703 r = K.gcd(r, q)

704 q = q // r

705

706 if K.is_one(q):

707 return None

708

709 result.append(q)

710

711 return result[1:]

712

713

714def dmp_zz_wang_test_points(f, T, ct, A, u, K):

715 """Wang/EEZ: Test evaluation points for suitability. """

716 if not dmp_eval_tail(dmp_LC(f, K), A, u - 1, K):

717 raise EvaluationFailed('no luck')

718

719 g = dmp_eval_tail(f, A, u, K)

720

721 if not dup_sqf_p(g, K):

722 raise EvaluationFailed('no luck')

723

724 c, h = dup_primitive(g, K)

725

726 if K.is_negative(dup_LC(h, K)):

727 c, h = -c, dup_neg(h, K)

728

729 v = u - 1

730

731 E = [ dmp_eval_tail(t, A, v, K) for t, _ in T ]

732 D = dmp_zz_wang_non_divisors(E, c, ct, K)

733

734 if D is not None:

735 return c, h, E

736 else:

737 raise EvaluationFailed('no luck')

738

739

740def dmp_zz_wang_lead_coeffs(f, T, cs, E, H, A, u, K):

741 """Wang/EEZ: Compute correct leading coefficients. """

742 C, J, v = [], [0]*len(E), u - 1

743

744 for h in H:

745 c = dmp_one(v, K)

746 d = dup_LC(h, K)*cs

747

748 for i in reversed(range(len(E))):

749 k, e, (t, _) = 0, E[i], T[i]

750

751 while not (d % e):

752 d, k = d//e, k + 1

753

754 if k != 0:

755 c, J[i] = dmp_mul(c, dmp_pow(t, k, v, K), v, K), 1

756

757 C.append(c)

758

759 if not all(J):

760 raise ExtraneousFactors # pragma: no cover

761

762 CC, HH = [], []

763

764 for c, h in zip(C, H):

765 d = dmp_eval_tail(c, A, v, K)

766 lc = dup_LC(h, K)

767

768 if K.is_one(cs):

769 cc = lc//d

770 else:

771 g = K.gcd(lc, d)

772 d, cc = d//g, lc//g

773 h, cs = dup_mul_ground(h, d, K), cs//d

774

775 c = dmp_mul_ground(c, cc, v, K)

776

777 CC.append(c)

778 HH.append(h)

779

780 if K.is_one(cs):

781 return f, HH, CC

782

783 CCC, HHH = [], []

784

785 for c, h in zip(CC, HH):

786 CCC.append(dmp_mul_ground(c, cs, v, K))

787 HHH.append(dmp_mul_ground(h, cs, 0, K))

788

789 f = dmp_mul_ground(f, cs**(len(H) - 1), u, K)

790

791 return f, HHH, CCC

792

793

794def dup_zz_diophantine(F, m, p, K):

795 """Wang/EEZ: Solve univariate Diophantine equations. """

796 if len(F) == 2:

797 a, b = F

798

799 f = gf_from_int_poly(a, p)

800 g = gf_from_int_poly(b, p)

801

802 s, t, G = gf_gcdex(g, f, p, K)

803

804 s = gf_lshift(s, m, K)

805 t = gf_lshift(t, m, K)

806

807 q, s = gf_div(s, f, p, K)

808

809 t = gf_add_mul(t, q, g, p, K)

810

811 s = gf_to_int_poly(s, p)

812 t = gf_to_int_poly(t, p)

813

814 result = [s, t]

815 else:

816 G = [F[-1]]

817

818 for f in reversed(F[1:-1]):

819 G.insert(0, dup_mul(f, G[0], K))

820

821 S, T = [], [[1]]

822

823 for f, g in zip(F, G):

824 t, s = dmp_zz_diophantine([g, f], T[-1], [], 0, p, 1, K)

825 T.append(t)

826 S.append(s)

827

828 result, S = [], S + [T[-1]]

829

830 for s, f in zip(S, F):

831 s = gf_from_int_poly(s, p)

832 f = gf_from_int_poly(f, p)

833

834 r = gf_rem(gf_lshift(s, m, K), f, p, K)

835 s = gf_to_int_poly(r, p)

837 result.append(s)

839 return result

842def dmp_zz_diophantine(F, c, A, d, p, u, K):

843 """Wang/EEZ: Solve multivariate Diophantine equations. """

844 if not A:

845 S = [ [] for _ in F ]

846 n = dup_degree(c)

847

848 for i, coeff in enumerate(c):

849 if not coeff:

850 continue

851

852 T = dup_zz_diophantine(F, n - i, p, K)

853

854 for j, (s, t) in enumerate(zip(S, T)):

855 t = dup_mul_ground(t, coeff, K)

856 S[j] = dup_trunc(dup_add(s, t, K), p, K)

857 else:

858 n = len(A)

859 e = dmp_expand(F, u, K)

860

861 a, A = A[-1], A[:-1]

862 B, G = [], []

863

864 for f in F:

865 B.append(dmp_quo(e, f, u, K))

866 G.append(dmp_eval_in(f, a, n, u, K))

867

868 C = dmp_eval_in(c, a, n, u, K)

869

870 v = u - 1

871

872 S = dmp_zz_diophantine(G, C, A, d, p, v, K)

873 S = [ dmp_raise(s, 1, v, K) for s in S ]

874

875 for s, b in zip(S, B):

876 c = dmp_sub_mul(c, s, b, u, K)

877

878 c = dmp_ground_trunc(c, p, u, K)

879

880 m = dmp_nest([K.one, -a], n, K)

881 M = dmp_one(n, K)

882

883 for k in K.map(range(0, d)):

884 if dmp_zero_p(c, u):

885 break

886

887 M = dmp_mul(M, m, u, K)

888 C = dmp_diff_eval_in(c, k + 1, a, n, u, K)

889

890 if not dmp_zero_p(C, v):

891 C = dmp_quo_ground(C, K.factorial(k + 1), v, K)

892 T = dmp_zz_diophantine(G, C, A, d, p, v, K)

893

894 for i, t in enumerate(T):

895 T[i] = dmp_mul(dmp_raise(t, 1, v, K), M, u, K)

896

897 for i, (s, t) in enumerate(zip(S, T)):

898 S[i] = dmp_add(s, t, u, K)

899

900 for t, b in zip(T, B):

901 c = dmp_sub_mul(c, t, b, u, K)

902

903 c = dmp_ground_trunc(c, p, u, K)

904

905 S = [ dmp_ground_trunc(s, p, u, K) for s in S ]

906

907 return S

908

909

910def dmp_zz_wang_hensel_lifting(f, H, LC, A, p, u, K):

911 """Wang/EEZ: Parallel Hensel lifting algorithm. """

912 S, n, v = [f], len(A), u - 1

913

914 H = list(H)

915

916 for i, a in enumerate(reversed(A[1:])):

917 s = dmp_eval_in(S[0], a, n - i, u - i, K)

918 S.insert(0, dmp_ground_trunc(s, p, v - i, K))

919

920 d = max(dmp_degree_list(f, u)[1:])

921

922 for j, s, a in zip(range(2, n + 2), S, A):

923 G, w = list(H), j - 1

924

925 I, J = A[:j - 2], A[j - 1:]

926

927 for i, (h, lc) in enumerate(zip(H, LC)):

928 lc = dmp_ground_trunc(dmp_eval_tail(lc, J, v, K), p, w - 1, K)

929 H[i] = [lc] + dmp_raise(h[1:], 1, w - 1, K)

930

931 m = dmp_nest([K.one, -a], w, K)

932 M = dmp_one(w, K)

933

934 c = dmp_sub(s, dmp_expand(H, w, K), w, K)

935

936 dj = dmp_degree_in(s, w, w)

937

938 for k in K.map(range(0, dj)):

939 if dmp_zero_p(c, w):

940 break

941

942 M = dmp_mul(M, m, w, K)

943 C = dmp_diff_eval_in(c, k + 1, a, w, w, K)

944

945 if not dmp_zero_p(C, w - 1):

946 C = dmp_quo_ground(C, K.factorial(k + 1), w - 1, K)

947 T = dmp_zz_diophantine(G, C, I, d, p, w - 1, K)

948

949 for i, (h, t) in enumerate(zip(H, T)):

950 h = dmp_add_mul(h, dmp_raise(t, 1, w - 1, K), M, w, K)

951 H[i] = dmp_ground_trunc(h, p, w, K)

952

953 h = dmp_sub(s, dmp_expand(H, w, K), w, K)

954 c = dmp_ground_trunc(h, p, w, K)

955

956 if dmp_expand(H, u, K) != f:

957 raise ExtraneousFactors # pragma: no cover

958 else:

959 return H

960

961

962def dmp_zz_wang(f, u, K, mod=None, seed=None):

963 r"""

964 Factor primitive square-free polynomials in `Z[X]`.

965

966 Given a multivariate polynomial `f` in `Z[x_1,...,x_n]`, which is

967 primitive and square-free in `x_1`, computes factorization of `f` into

968 irreducibles over integers.

969

970 The procedure is based on Wang's Enhanced Extended Zassenhaus

971 algorithm. The algorithm works by viewing `f` as a univariate polynomial

972 in `Z[x_2,...,x_n][x_1]`, for which an evaluation mapping is computed::

973

974 x_2 -> a_2, ..., x_n -> a_n

975

976 where `a_i`, for `i = 2, \dots, n`, are carefully chosen integers. The

977 mapping is used to transform `f` into a univariate polynomial in `Z[x_1]`,

978 which can be factored efficiently using Zassenhaus algorithm. The last

979 step is to lift univariate factors to obtain true multivariate

980 factors. For this purpose a parallel Hensel lifting procedure is used.

981

982 The parameter ``seed`` is passed to _randint and can be used to seed randint

983 (when an integer) or (for testing purposes) can be a sequence of numbers.

984

985 References

986 ==========

987

988 .. [1] [Wang78]_

989 .. [2] [Geddes92]_

990

991 """

992 from sympy.ntheory import nextprime

993

994 randint = _randint(seed)

995

996 ct, T = dmp_zz_factor(dmp_LC(f, K), u - 1, K)

997

998 b = dmp_zz_mignotte_bound(f, u, K)

999 p = K(nextprime(b))

1000

1001 if mod is None:

1002 if u == 1:

1003 mod = 2

1004 else:

1005 mod = 1

1006

1007 history, configs, A, r = set(), [], [K.zero]*u, None

1008

1009 try:

1010 cs, s, E = dmp_zz_wang_test_points(f, T, ct, A, u, K)

1011

1012 _, H = dup_zz_factor_sqf(s, K)

1013

1014 r = len(H)

1015

1016 if r == 1:

1017 return [f]

1018

1019 configs = [(s, cs, E, H, A)]

1020 except EvaluationFailed:

1021 pass

1022

1023 eez_num_configs = query('EEZ_NUMBER_OF_CONFIGS')

1024 eez_num_tries = query('EEZ_NUMBER_OF_TRIES')

1025 eez_mod_step = query('EEZ_MODULUS_STEP')

1026

1027 while len(configs) < eez_num_configs:

1028 for _ in range(eez_num_tries):

1029 A = [ K(randint(-mod, mod)) for _ in range(u) ]

1030

1031 if tuple(A) not in history:

1032 history.add(tuple(A))

1033 else:

1034 continue

1035

1036 try:

1037 cs, s, E = dmp_zz_wang_test_points(f, T, ct, A, u, K)

1038 except EvaluationFailed:

1039 continue

1040

1041 _, H = dup_zz_factor_sqf(s, K)

1042

1043 rr = len(H)

1044

1045 if r is not None:

1046 if rr != r: # pragma: no cover

1047 if rr < r:

1048 configs, r = [], rr

1049 else:

1050 continue

1051 else:

1052 r = rr

1053

1054 if r == 1:

1055 return [f]

1056

1057 configs.append((s, cs, E, H, A))

1058

1059 if len(configs) == eez_num_configs:

1060 break

1061 else:

1062 mod += eez_mod_step

1063

1064 s_norm, s_arg, i = None, 0, 0

1065

1066 for s, _, _, _, _ in configs:

1067 _s_norm = dup_max_norm(s, K)

1068

1069 if s_norm is not None:

1070 if _s_norm < s_norm:

1071 s_norm = _s_norm

1072 s_arg = i

1073 else:

1074 s_norm = _s_norm

1075

1076 i += 1

1077

1078 _, cs, E, H, A = configs[s_arg]

1079 orig_f = f

1080

1081 try:

1082 f, H, LC = dmp_zz_wang_lead_coeffs(f, T, cs, E, H, A, u, K)

1083 factors = dmp_zz_wang_hensel_lifting(f, H, LC, A, p, u, K)

1084 except ExtraneousFactors: # pragma: no cover

1085 if query('EEZ_RESTART_IF_NEEDED'):

1086 return dmp_zz_wang(orig_f, u, K, mod + 1)

1087 else:

1088 raise ExtraneousFactors(

1089 "we need to restart algorithm with better parameters")

1090

1091 result = []

1092

1093 for f in factors:

1094 _, f = dmp_ground_primitive(f, u, K)

1095

1096 if K.is_negative(dmp_ground_LC(f, u, K)):

1097 f = dmp_neg(f, u, K)

1098

1099 result.append(f)

1100

1101 return result

1102

1103

1104def dmp_zz_factor(f, u, K):

1105 r"""

1106 Factor (non square-free) polynomials in `Z[X]`.

1107

1108 Given a multivariate polynomial `f` in `Z[x]` computes its complete

1109 factorization `f_1, \dots, f_n` into irreducibles over integers::

1110

1111 f = content(f) f_1**k_1 ... f_n**k_n

1112

1113 The factorization is computed by reducing the input polynomial

1114 into a primitive square-free polynomial and factoring it using

1115 Enhanced Extended Zassenhaus (EEZ) algorithm. Trial division

1116 is used to recover the multiplicities of factors.

1117

1118 The result is returned as a tuple consisting of::

1119

1120 (content(f), [(f_1, k_1), ..., (f_n, k_n))

1121

1122 Consider polynomial `f = 2*(x**2 - y**2)`::

1123

1124 >>> from sympy.polys import ring, ZZ

1125 >>> R, x,y = ring("x,y", ZZ)

1126

1127 >>> R.dmp_zz_factor(2*x**2 - 2*y**2)

1128 (2, [(x - y, 1), (x + y, 1)])

1129

1130 In result we got the following factorization::

1131

1132 f = 2 (x - y) (x + y)

1133

1134 References

1135 ==========

1136

1137 .. [1] [Gathen99]_

1138

1139 """

1140 if not u:

1141 return dup_zz_factor(f, K)

1142

1143 if dmp_zero_p(f, u):

1144 return K.zero, []

1145

1146 cont, g = dmp_ground_primitive(f, u, K)

1147

1148 if dmp_ground_LC(g, u, K) < 0:

1149 cont, g = -cont, dmp_neg(g, u, K)

1150

1151 if all(d <= 0 for d in dmp_degree_list(g, u)):

1152 return cont, []

1153

1154 G, g = dmp_primitive(g, u, K)

1155

1156 factors = []

1157

1158 if dmp_degree(g, u) > 0:

1159 g = dmp_sqf_part(g, u, K)

1160 H = dmp_zz_wang(g, u, K)

1161 factors = dmp_trial_division(f, H, u, K)

1162

1163 for g, k in dmp_zz_factor(G, u - 1, K)[1]:

1164 factors.insert(0, ([g], k))

1165

1166 return cont, _sort_factors(factors)

1167

1168

1169def dup_qq_i_factor(f, K0):

1170 """Factor univariate polynomials into irreducibles in `QQ_I[x]`. """

1171 # Factor in QQ<I>

1172 K1 = K0.as_AlgebraicField()

1173 f = dup_convert(f, K0, K1)

1174 coeff, factors = dup_factor_list(f, K1)

1175 factors = [(dup_convert(fac, K1, K0), i) for fac, i in factors]

1176 coeff = K0.convert(coeff, K1)

1177 return coeff, factors

1178

1179

1180def dup_zz_i_factor(f, K0):

1181 """Factor univariate polynomials into irreducibles in `ZZ_I[x]`. """

1182 # First factor in QQ_I

1183 K1 = K0.get_field()

1184 f = dup_convert(f, K0, K1)

1185 coeff, factors = dup_qq_i_factor(f, K1)

1186

1187 new_factors = []

1188 for fac, i in factors:

1189 # Extract content

1190 fac_denom, fac_num = dup_clear_denoms(fac, K1)

1191 fac_num_ZZ_I = dup_convert(fac_num, K1, K0)

1192 content, fac_prim = dmp_ground_primitive(fac_num_ZZ_I, 0, K1)

1193

1194 coeff = (coeff * content ** i) // fac_denom ** i

1195 new_factors.append((fac_prim, i))

1196

1197 factors = new_factors

1198 coeff = K0.convert(coeff, K1)

1199 return coeff, factors

1200

1201

1202def dmp_qq_i_factor(f, u, K0):

1203 """Factor multivariate polynomials into irreducibles in `QQ_I[X]`. """

1204 # Factor in QQ<I>

1205 K1 = K0.as_AlgebraicField()

1206 f = dmp_convert(f, u, K0, K1)

1207 coeff, factors = dmp_factor_list(f, u, K1)

1208 factors = [(dmp_convert(fac, u, K1, K0), i) for fac, i in factors]

1209 coeff = K0.convert(coeff, K1)

1210 return coeff, factors

1211

1212

1213def dmp_zz_i_factor(f, u, K0):

1214 """Factor multivariate polynomials into irreducibles in `ZZ_I[X]`. """

1215 # First factor in QQ_I

1216 K1 = K0.get_field()

1217 f = dmp_convert(f, u, K0, K1)

1218 coeff, factors = dmp_qq_i_factor(f, u, K1)

1219

1220 new_factors = []

1221 for fac, i in factors:

1222 # Extract content

1223 fac_denom, fac_num = dmp_clear_denoms(fac, u, K1)

1224 fac_num_ZZ_I = dmp_convert(fac_num, u, K1, K0)

1225 content, fac_prim = dmp_ground_primitive(fac_num_ZZ_I, u, K1)

1226

1227 coeff = (coeff * content ** i) // fac_denom ** i

1228 new_factors.append((fac_prim, i))

1229

1230 factors = new_factors

1231 coeff = K0.convert(coeff, K1)

1232 return coeff, factors

1233

1234

1235def dup_ext_factor(f, K):

1236 """Factor univariate polynomials over algebraic number fields. """

1237 n, lc = dup_degree(f), dup_LC(f, K)

1238

1239 f = dup_monic(f, K)

1240

1241 if n <= 0:

1242 return lc, []

1243 if n == 1:

1244 return lc, [(f, 1)]

1245

1246 f, F = dup_sqf_part(f, K), f

1247 s, g, r = dup_sqf_norm(f, K)

1248

1249 factors = dup_factor_list_include(r, K.dom)

1250

1251 if len(factors) == 1:

1252 return lc, [(f, n//dup_degree(f))]

1253

1254 H = s*K.unit

1255

1256 for i, (factor, _) in enumerate(factors):

1257 h = dup_convert(factor, K.dom, K)

1258 h, _, g = dup_inner_gcd(h, g, K)

1259 h = dup_shift(h, H, K)

1260 factors[i] = h

1261

1262 factors = dup_trial_division(F, factors, K)

1263 return lc, factors

1264

1265

1266def dmp_ext_factor(f, u, K):

1267 """Factor multivariate polynomials over algebraic number fields. """

1268 if not u:

1269 return dup_ext_factor(f, K)

1270

1271 lc = dmp_ground_LC(f, u, K)

1272 f = dmp_ground_monic(f, u, K)

1273

1274 if all(d <= 0 for d in dmp_degree_list(f, u)):

1275 return lc, []

1276

1277 f, F = dmp_sqf_part(f, u, K), f

1278 s, g, r = dmp_sqf_norm(f, u, K)

1279

1280 factors = dmp_factor_list_include(r, u, K.dom)

1281

1282 if len(factors) == 1:

1283 factors = [f]

1284 else:

1285 H = dmp_raise([K.one, s*K.unit], u, 0, K)

1286

1287 for i, (factor, _) in enumerate(factors):

1288 h = dmp_convert(factor, u, K.dom, K)

1289 h, _, g = dmp_inner_gcd(h, g, u, K)

1290 h = dmp_compose(h, H, u, K)

1291 factors[i] = h

1292

1293 return lc, dmp_trial_division(F, factors, u, K)

1294

1295

1296def dup_gf_factor(f, K):

1297 """Factor univariate polynomials over finite fields. """

1298 f = dup_convert(f, K, K.dom)

1299

1300 coeff, factors = gf_factor(f, K.mod, K.dom)

1301

1302 for i, (f, k) in enumerate(factors):

1303 factors[i] = (dup_convert(f, K.dom, K), k)

1304

1305 return K.convert(coeff, K.dom), factors

1306

1307

1308def dmp_gf_factor(f, u, K):

1309 """Factor multivariate polynomials over finite fields. """

1310 raise NotImplementedError('multivariate polynomials over finite fields')

1311

1312

1313def dup_factor_list(f, K0):

1314 """Factor univariate polynomials into irreducibles in `K[x]`. """

1315 j, f = dup_terms_gcd(f, K0)

1316 cont, f = dup_primitive(f, K0)

1317

1318 if K0.is_FiniteField:

1319 coeff, factors = dup_gf_factor(f, K0)

1320 elif K0.is_Algebraic:

1321 coeff, factors = dup_ext_factor(f, K0)

1322 elif K0.is_GaussianRing:

1323 coeff, factors = dup_zz_i_factor(f, K0)

1324 elif K0.is_GaussianField:

1325 coeff, factors = dup_qq_i_factor(f, K0)

1326 else:

1327 if not K0.is_Exact:

1328 K0_inexact, K0 = K0, K0.get_exact()

1329 f = dup_convert(f, K0_inexact, K0)

1330 else:

1331 K0_inexact = None

1332

1333 if K0.is_Field:

1334 K = K0.get_ring()

1335

1336 denom, f = dup_clear_denoms(f, K0, K)

1337 f = dup_convert(f, K0, K)

1338 else:

1339 K = K0

1340

1341 if K.is_ZZ:

1342 coeff, factors = dup_zz_factor(f, K)

1343 elif K.is_Poly:

1344 f, u = dmp_inject(f, 0, K)

1345

1346 coeff, factors = dmp_factor_list(f, u, K.dom)

1347

1348 for i, (f, k) in enumerate(factors):

1349 factors[i] = (dmp_eject(f, u, K), k)

1350

1351 coeff = K.convert(coeff, K.dom)

1352 else: # pragma: no cover

1353 raise DomainError('factorization not supported over %s' % K0)

1354

1355 if K0.is_Field:

1356 for i, (f, k) in enumerate(factors):

1357 factors[i] = (dup_convert(f, K, K0), k)

1358

1359 coeff = K0.convert(coeff, K)

1360 coeff = K0.quo(coeff, denom)

1361

1362 if K0_inexact:

1363 for i, (f, k) in enumerate(factors):

1364 max_norm = dup_max_norm(f, K0)

1365 f = dup_quo_ground(f, max_norm, K0)

1366 f = dup_convert(f, K0, K0_inexact)

1367 factors[i] = (f, k)

1368 coeff = K0.mul(coeff, K0.pow(max_norm, k))

1369

1370 coeff = K0_inexact.convert(coeff, K0)

1371 K0 = K0_inexact

1372

1373 if j:

1374 factors.insert(0, ([K0.one, K0.zero], j))

1375

1376 return coeff*cont, _sort_factors(factors)

1377

1378

1379def dup_factor_list_include(f, K):

1380 """Factor univariate polynomials into irreducibles in `K[x]`. """

1381 coeff, factors = dup_factor_list(f, K)

1382

1383 if not factors:

1384 return [(dup_strip([coeff]), 1)]

1385 else:

1386 g = dup_mul_ground(factors[0][0], coeff, K)

1387 return [(g, factors[0][1])] + factors[1:]

1388

1389

1390def dmp_factor_list(f, u, K0):

1391 """Factor multivariate polynomials into irreducibles in `K[X]`. """

1392 if not u:

1393 return dup_factor_list(f, K0)

1394

1395 J, f = dmp_terms_gcd(f, u, K0)

1396 cont, f = dmp_ground_primitive(f, u, K0)

1397

1398 if K0.is_FiniteField: # pragma: no cover

1399 coeff, factors = dmp_gf_factor(f, u, K0)

1400 elif K0.is_Algebraic:

1401 coeff, factors = dmp_ext_factor(f, u, K0)

1402 elif K0.is_GaussianRing:

1403 coeff, factors = dmp_zz_i_factor(f, u, K0)

1404 elif K0.is_GaussianField:

1405 coeff, factors = dmp_qq_i_factor(f, u, K0)

1406 else:

1407 if not K0.is_Exact:

1408 K0_inexact, K0 = K0, K0.get_exact()

1409 f = dmp_convert(f, u, K0_inexact, K0)

1410 else:

1411 K0_inexact = None

1412

1413 if K0.is_Field:

1414 K = K0.get_ring()

1415

1416 denom, f = dmp_clear_denoms(f, u, K0, K)

1417 f = dmp_convert(f, u, K0, K)

1418 else:

1419 K = K0

1420

1421 if K.is_ZZ:

1422 levels, f, v = dmp_exclude(f, u, K)

1423 coeff, factors = dmp_zz_factor(f, v, K)

1424

1425 for i, (f, k) in enumerate(factors):

1426 factors[i] = (dmp_include(f, levels, v, K), k)

1427 elif K.is_Poly:

1428 f, v = dmp_inject(f, u, K)

1429

1430 coeff, factors = dmp_factor_list(f, v, K.dom)

1431

1432 for i, (f, k) in enumerate(factors):

1433 factors[i] = (dmp_eject(f, v, K), k)

1434

1435 coeff = K.convert(coeff, K.dom)

1436 else: # pragma: no cover

1437 raise DomainError('factorization not supported over %s' % K0)

1438

1439 if K0.is_Field:

1440 for i, (f, k) in enumerate(factors):

1441 factors[i] = (dmp_convert(f, u, K, K0), k)

1442

1443 coeff = K0.convert(coeff, K)

1444 coeff = K0.quo(coeff, denom)

1445

1446 if K0_inexact:

1447 for i, (f, k) in enumerate(factors):

1448 max_norm = dmp_max_norm(f, u, K0)

1449 f = dmp_quo_ground(f, max_norm, u, K0)

1450 f = dmp_convert(f, u, K0, K0_inexact)

1451 factors[i] = (f, k)

1452 coeff = K0.mul(coeff, K0.pow(max_norm, k))

1453

1454 coeff = K0_inexact.convert(coeff, K0)

1455 K0 = K0_inexact

1456

1457 for i, j in enumerate(reversed(J)):

1458 if not j:

1459 continue

1460

1461 term = {(0,)*(u - i) + (1,) + (0,)*i: K0.one}

1462 factors.insert(0, (dmp_from_dict(term, u, K0), j))

1463

1464 return coeff*cont, _sort_factors(factors)

1465

1466

1467def dmp_factor_list_include(f, u, K):

1468 """Factor multivariate polynomials into irreducibles in `K[X]`. """

1469 if not u:

1470 return dup_factor_list_include(f, K)

1471

1472 coeff, factors = dmp_factor_list(f, u, K)

1473

1474 if not factors:

1475 return [(dmp_ground(coeff, u), 1)]

1476 else:

1477 g = dmp_mul_ground(factors[0][0], coeff, u, K)

1478 return [(g, factors[0][1])] + factors[1:]

1479

1480

1481def dup_irreducible_p(f, K):

1482 """

1483 Returns ``True`` if a univariate polynomial ``f`` has no factors

1484 over its domain.

1485 """

1486 return dmp_irreducible_p(f, 0, K)

1487

1488

1489def dmp_irreducible_p(f, u, K):

1490 """

1491 Returns ``True`` if a multivariate polynomial ``f`` has no factors

1492 over its domain.

1493 """

1494 _, factors = dmp_factor_list(f, u, K)

1495

1496 if not factors:

1497 return True

1498 elif len(factors) > 1:

1499 return False

1500 else:

1501 _, k = factors[0]

1502 return k == 1