Coverage for /usr/lib/python3/dist-packages/sympy/polys/euclidtools.py: 9%

1"""Euclidean algorithms, GCDs, LCMs and polynomial remainder sequences. """

4from sympy.polys.densearith import (

5 dup_sub_mul,

6 dup_neg, dmp_neg,

7 dmp_add,

8 dmp_sub,

9 dup_mul, dmp_mul,

10 dmp_pow,

11 dup_div, dmp_div,

12 dup_rem,

13 dup_quo, dmp_quo,

14 dup_prem, dmp_prem,

15 dup_mul_ground, dmp_mul_ground,

16 dmp_mul_term,

17 dup_quo_ground, dmp_quo_ground,

18 dup_max_norm, dmp_max_norm)

19from sympy.polys.densebasic import (

20 dup_strip, dmp_raise,

21 dmp_zero, dmp_one, dmp_ground,

22 dmp_one_p, dmp_zero_p,

23 dmp_zeros,

24 dup_degree, dmp_degree, dmp_degree_in,

25 dup_LC, dmp_LC, dmp_ground_LC,

26 dmp_multi_deflate, dmp_inflate,

27 dup_convert, dmp_convert,

28 dmp_apply_pairs)

29from sympy.polys.densetools import (

30 dup_clear_denoms, dmp_clear_denoms,

31 dup_diff, dmp_diff,

32 dup_eval, dmp_eval, dmp_eval_in,

33 dup_trunc, dmp_ground_trunc,

34 dup_monic, dmp_ground_monic,

35 dup_primitive, dmp_ground_primitive,

36 dup_extract, dmp_ground_extract)

37from sympy.polys.galoistools import (

38 gf_int, gf_crt)

39from sympy.polys.polyconfig import query

40from sympy.polys.polyerrors import (

41 MultivariatePolynomialError,

42 HeuristicGCDFailed,

43 HomomorphismFailed,

44 NotInvertible,

45 DomainError)

50def dup_half_gcdex(f, g, K):

51 """

52 Half extended Euclidean algorithm in `F[x]`.

54 Returns ``(s, h)`` such that ``h = gcd(f, g)`` and ``s*f = h (mod g)``.

56 Examples

57 ========

59 >>> from sympy.polys import ring, QQ

60 >>> R, x = ring("x", QQ)

62 >>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15

63 >>> g = x**3 + x**2 - 4*x - 4

65 >>> R.dup_half_gcdex(f, g)

66 (-1/5*x + 3/5, x + 1)

68 """

69 if not K.is_Field:

70 raise DomainError("Cannot compute half extended GCD over %s" % K)

72 a, b = [K.one], []

74 while g:

75 q, r = dup_div(f, g, K)

76 f, g = g, r

77 a, b = b, dup_sub_mul(a, q, b, K)

79 a = dup_quo_ground(a, dup_LC(f, K), K)

80 f = dup_monic(f, K)

82 return a, f

85def dmp_half_gcdex(f, g, u, K):

86 """

87 Half extended Euclidean algorithm in `F[X]`.

89 Examples

90 ========

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

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

95 """

96 if not u:

97 return dup_half_gcdex(f, g, K)

98 else:

99 raise MultivariatePolynomialError(f, g)

100

101

102def dup_gcdex(f, g, K):

103 """

104 Extended Euclidean algorithm in `F[x]`.

105

106 Returns ``(s, t, h)`` such that ``h = gcd(f, g)`` and ``s*f + t*g = h``.

107

108 Examples

109 ========

110

111 >>> from sympy.polys import ring, QQ

112 >>> R, x = ring("x", QQ)

113

114 >>> f = x**4 - 2*x**3 - 6*x**2 + 12*x + 15

115 >>> g = x**3 + x**2 - 4*x - 4

116

117 >>> R.dup_gcdex(f, g)

118 (-1/5*x + 3/5, 1/5*x**2 - 6/5*x + 2, x + 1)

119

120 """

121 s, h = dup_half_gcdex(f, g, K)

122

123 F = dup_sub_mul(h, s, f, K)

124 t = dup_quo(F, g, K)

125

126 return s, t, h

127

128

129def dmp_gcdex(f, g, u, K):

130 """

131 Extended Euclidean algorithm in `F[X]`.

132

133 Examples

134 ========

135

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

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

138

139 """

140 if not u:

141 return dup_gcdex(f, g, K)

142 else:

143 raise MultivariatePolynomialError(f, g)

144

145

146def dup_invert(f, g, K):

147 """

148 Compute multiplicative inverse of `f` modulo `g` in `F[x]`.

149

150 Examples

151 ========

152

153 >>> from sympy.polys import ring, QQ

154 >>> R, x = ring("x", QQ)

155

156 >>> f = x**2 - 1

157 >>> g = 2*x - 1

158 >>> h = x - 1

159

160 >>> R.dup_invert(f, g)

161 -4/3

162

163 >>> R.dup_invert(f, h)

164 Traceback (most recent call last):

165 ...

166 NotInvertible: zero divisor

167

168 """

169 s, h = dup_half_gcdex(f, g, K)

170

171 if h == [K.one]:

172 return dup_rem(s, g, K)

173 else:

174 raise NotInvertible("zero divisor")

175

176

177def dmp_invert(f, g, u, K):

178 """

179 Compute multiplicative inverse of `f` modulo `g` in `F[X]`.

180

181 Examples

182 ========

183

184 >>> from sympy.polys import ring, QQ

185 >>> R, x = ring("x", QQ)

186

187 """

188 if not u:

189 return dup_invert(f, g, K)

190 else:

191 raise MultivariatePolynomialError(f, g)

192

193

194def dup_euclidean_prs(f, g, K):

195 """

196 Euclidean polynomial remainder sequence (PRS) in `K[x]`.

197

198 Examples

199 ========

200

201 >>> from sympy.polys import ring, QQ

202 >>> R, x = ring("x", QQ)

203

204 >>> f = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5

205 >>> g = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21

206

207 >>> prs = R.dup_euclidean_prs(f, g)

208

209 >>> prs[0]

210 x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5

211 >>> prs[1]

212 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21

213 >>> prs[2]

214 -5/9*x**4 + 1/9*x**2 - 1/3

215 >>> prs[3]

216 -117/25*x**2 - 9*x + 441/25

217 >>> prs[4]

218 233150/19773*x - 102500/6591

219 >>> prs[5]

220 -1288744821/543589225

221

222 """

223 prs = [f, g]

224 h = dup_rem(f, g, K)

225

226 while h:

227 prs.append(h)

228 f, g = g, h

229 h = dup_rem(f, g, K)

230

231 return prs

232

233

234def dmp_euclidean_prs(f, g, u, K):

235 """

236 Euclidean polynomial remainder sequence (PRS) in `K[X]`.

237

238 Examples

239 ========

240

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

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

243

244 """

245 if not u:

246 return dup_euclidean_prs(f, g, K)

247 else:

248 raise MultivariatePolynomialError(f, g)

249

250

251def dup_primitive_prs(f, g, K):

252 """

253 Primitive polynomial remainder sequence (PRS) in `K[x]`.

254

255 Examples

256 ========

257

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

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

260

261 >>> f = x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5

262 >>> g = 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21

263

264 >>> prs = R.dup_primitive_prs(f, g)

265

266 >>> prs[0]

267 x**8 + x**6 - 3*x**4 - 3*x**3 + 8*x**2 + 2*x - 5

268 >>> prs[1]

269 3*x**6 + 5*x**4 - 4*x**2 - 9*x + 21

270 >>> prs[2]

271 -5*x**4 + x**2 - 3

272 >>> prs[3]

273 13*x**2 + 25*x - 49

274 >>> prs[4]

275 4663*x - 6150

276 >>> prs[5]

277 1

278

279 """

280 prs = [f, g]

281 _, h = dup_primitive(dup_prem(f, g, K), K)

282

283 while h:

284 prs.append(h)

285 f, g = g, h

286 _, h = dup_primitive(dup_prem(f, g, K), K)

287

288 return prs

289

290

291def dmp_primitive_prs(f, g, u, K):

292 """

293 Primitive polynomial remainder sequence (PRS) in `K[X]`.

294

295 Examples

296 ========

297

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

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

300

301 """

302 if not u:

303 return dup_primitive_prs(f, g, K)

304 else:

305 raise MultivariatePolynomialError(f, g)

306

307

308def dup_inner_subresultants(f, g, K):

309 """

310 Subresultant PRS algorithm in `K[x]`.

311

312 Computes the subresultant polynomial remainder sequence (PRS)

313 and the non-zero scalar subresultants of `f` and `g`.

314 By [1] Thm. 3, these are the constants '-c' (- to optimize

315 computation of sign).

316 The first subdeterminant is set to 1 by convention to match

317 the polynomial and the scalar subdeterminants.

318 If 'deg(f) < deg(g)', the subresultants of '(g,f)' are computed.

319

320 Examples

321 ========

322

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

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

325

326 >>> R.dup_inner_subresultants(x**2 + 1, x**2 - 1)

327 ([x**2 + 1, x**2 - 1, -2], [1, 1, 4])

328

329 References

330 ==========

331

332 .. [1] W.S. Brown, The Subresultant PRS Algorithm.

333 ACM Transaction of Mathematical Software 4 (1978) 237-249

334

335 """

336 n = dup_degree(f)

337 m = dup_degree(g)

338

339 if n < m:

340 f, g = g, f

341 n, m = m, n

342

343 if not f:

344 return [], []

345

346 if not g:

347 return [f], [K.one]

348

349 R = [f, g]

350 d = n - m

351

352 b = (-K.one)**(d + 1)

353

354 h = dup_prem(f, g, K)

355 h = dup_mul_ground(h, b, K)

356

357 lc = dup_LC(g, K)

358 c = lc**d

359

360 # Conventional first scalar subdeterminant is 1

361 S = [K.one, c]

362 c = -c

363

364 while h:

365 k = dup_degree(h)

366 R.append(h)

367

368 f, g, m, d = g, h, k, m - k

369

370 b = -lc * c**d

371

372 h = dup_prem(f, g, K)

373 h = dup_quo_ground(h, b, K)

374

375 lc = dup_LC(g, K)

376

377 if d > 1: # abnormal case

378 q = c**(d - 1)

379 c = K.quo((-lc)**d, q)

380 else:

381 c = -lc

383 S.append(-c)

385 return R, S

388def dup_subresultants(f, g, K):

389 """

390 Computes subresultant PRS of two polynomials in `K[x]`.

391

392 Examples

393 ========

394

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

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

397

398 >>> R.dup_subresultants(x**2 + 1, x**2 - 1)

399 [x**2 + 1, x**2 - 1, -2]

400

401 """

402 return dup_inner_subresultants(f, g, K)[0]

403

404

405def dup_prs_resultant(f, g, K):

406 """

407 Resultant algorithm in `K[x]` using subresultant PRS.

408

409 Examples

410 ========

411

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

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

414

415 >>> R.dup_prs_resultant(x**2 + 1, x**2 - 1)

416 (4, [x**2 + 1, x**2 - 1, -2])

417

418 """

419 if not f or not g:

420 return (K.zero, [])

421

422 R, S = dup_inner_subresultants(f, g, K)

423

424 if dup_degree(R[-1]) > 0:

425 return (K.zero, R)

426

427 return S[-1], R

428

429

430def dup_resultant(f, g, K, includePRS=False):

431 """

432 Computes resultant of two polynomials in `K[x]`.

433

434 Examples

435 ========

436

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

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

439

440 >>> R.dup_resultant(x**2 + 1, x**2 - 1)

441 4

442

443 """

444 if includePRS:

445 return dup_prs_resultant(f, g, K)

446 return dup_prs_resultant(f, g, K)[0]

447

448

449def dmp_inner_subresultants(f, g, u, K):

450 """

451 Subresultant PRS algorithm in `K[X]`.

452

453 Examples

454 ========

455

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

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

458

459 >>> f = 3*x**2*y - y**3 - 4

460 >>> g = x**2 + x*y**3 - 9

461

462 >>> a = 3*x*y**4 + y**3 - 27*y + 4

463 >>> b = -3*y**10 - 12*y**7 + y**6 - 54*y**4 + 8*y**3 + 729*y**2 - 216*y + 16

464

465 >>> prs = [f, g, a, b]

466 >>> sres = [[1], [1], [3, 0, 0, 0, 0], [-3, 0, 0, -12, 1, 0, -54, 8, 729, -216, 16]]

467

468 >>> R.dmp_inner_subresultants(f, g) == (prs, sres)

469 True

470

471 """

472 if not u:

473 return dup_inner_subresultants(f, g, K)

474

475 n = dmp_degree(f, u)

476 m = dmp_degree(g, u)

477

478 if n < m:

479 f, g = g, f

480 n, m = m, n

481

482 if dmp_zero_p(f, u):

483 return [], []

484

485 v = u - 1

486 if dmp_zero_p(g, u):

487 return [f], [dmp_ground(K.one, v)]

488

489 R = [f, g]

490 d = n - m

491

492 b = dmp_pow(dmp_ground(-K.one, v), d + 1, v, K)

493

494 h = dmp_prem(f, g, u, K)

495 h = dmp_mul_term(h, b, 0, u, K)

496

497 lc = dmp_LC(g, K)

498 c = dmp_pow(lc, d, v, K)

499

500 S = [dmp_ground(K.one, v), c]

501 c = dmp_neg(c, v, K)

502

503 while not dmp_zero_p(h, u):

504 k = dmp_degree(h, u)

505 R.append(h)

506

507 f, g, m, d = g, h, k, m - k

508

509 b = dmp_mul(dmp_neg(lc, v, K),

510 dmp_pow(c, d, v, K), v, K)

511

512 h = dmp_prem(f, g, u, K)

513 h = [ dmp_quo(ch, b, v, K) for ch in h ]

514

515 lc = dmp_LC(g, K)

516

517 if d > 1:

518 p = dmp_pow(dmp_neg(lc, v, K), d, v, K)

519 q = dmp_pow(c, d - 1, v, K)

520 c = dmp_quo(p, q, v, K)

521 else:

522 c = dmp_neg(lc, v, K)

523

524 S.append(dmp_neg(c, v, K))

525

526 return R, S

527

528

529def dmp_subresultants(f, g, u, K):

530 """

531 Computes subresultant PRS of two polynomials in `K[X]`.

532

533 Examples

534 ========

535

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

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

538

539 >>> f = 3*x**2*y - y**3 - 4

540 >>> g = x**2 + x*y**3 - 9

541

542 >>> a = 3*x*y**4 + y**3 - 27*y + 4

543 >>> b = -3*y**10 - 12*y**7 + y**6 - 54*y**4 + 8*y**3 + 729*y**2 - 216*y + 16

544

545 >>> R.dmp_subresultants(f, g) == [f, g, a, b]

546 True

547

548 """

549 return dmp_inner_subresultants(f, g, u, K)[0]

550

551

552def dmp_prs_resultant(f, g, u, K):

553 """

554 Resultant algorithm in `K[X]` using subresultant PRS.

555

556 Examples

557 ========

558

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

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

561

562 >>> f = 3*x**2*y - y**3 - 4

563 >>> g = x**2 + x*y**3 - 9

564

565 >>> a = 3*x*y**4 + y**3 - 27*y + 4

566 >>> b = -3*y**10 - 12*y**7 + y**6 - 54*y**4 + 8*y**3 + 729*y**2 - 216*y + 16

567

568 >>> res, prs = R.dmp_prs_resultant(f, g)

569

570 >>> res == b # resultant has n-1 variables

571 False

572 >>> res == b.drop(x)

573 True

574 >>> prs == [f, g, a, b]

575 True

576

577 """

578 if not u:

579 return dup_prs_resultant(f, g, K)

580

581 if dmp_zero_p(f, u) or dmp_zero_p(g, u):

582 return (dmp_zero(u - 1), [])

583

584 R, S = dmp_inner_subresultants(f, g, u, K)

585

586 if dmp_degree(R[-1], u) > 0:

587 return (dmp_zero(u - 1), R)

588

589 return S[-1], R

590

591

592def dmp_zz_modular_resultant(f, g, p, u, K):

593 """

594 Compute resultant of `f` and `g` modulo a prime `p`.

595

596 Examples

597 ========

598

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

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

601

602 >>> f = x + y + 2

603 >>> g = 2*x*y + x + 3

604

605 >>> R.dmp_zz_modular_resultant(f, g, 5)

606 -2*y**2 + 1

607

608 """

609 if not u:

610 return gf_int(dup_prs_resultant(f, g, K)[0] % p, p)

611

612 v = u - 1

613

614 n = dmp_degree(f, u)

615 m = dmp_degree(g, u)

616

617 N = dmp_degree_in(f, 1, u)

618 M = dmp_degree_in(g, 1, u)

619

620 B = n*M + m*N

621

622 D, a = [K.one], -K.one

623 r = dmp_zero(v)

624

625 while dup_degree(D) <= B:

626 while True:

627 a += K.one

628

629 if a == p:

630 raise HomomorphismFailed('no luck')

631

632 F = dmp_eval_in(f, gf_int(a, p), 1, u, K)

633

634 if dmp_degree(F, v) == n:

635 G = dmp_eval_in(g, gf_int(a, p), 1, u, K)

636

637 if dmp_degree(G, v) == m:

638 break

639

640 R = dmp_zz_modular_resultant(F, G, p, v, K)

641 e = dmp_eval(r, a, v, K)

642

643 if not v:

644 R = dup_strip([R])

645 e = dup_strip([e])

646 else:

647 R = [R]

648 e = [e]

649

650 d = K.invert(dup_eval(D, a, K), p)

651 d = dup_mul_ground(D, d, K)

652 d = dmp_raise(d, v, 0, K)

653

654 c = dmp_mul(d, dmp_sub(R, e, v, K), v, K)

655 r = dmp_add(r, c, v, K)

656

657 r = dmp_ground_trunc(r, p, v, K)

658

659 D = dup_mul(D, [K.one, -a], K)

660 D = dup_trunc(D, p, K)

661

662 return r

663

664

665def _collins_crt(r, R, P, p, K):

666 """Wrapper of CRT for Collins's resultant algorithm. """

667 return gf_int(gf_crt([r, R], [P, p], K), P*p)

668

669

670def dmp_zz_collins_resultant(f, g, u, K):

671 """

672 Collins's modular resultant algorithm in `Z[X]`.

673

674 Examples

675 ========

676

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

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

679

680 >>> f = x + y + 2

681 >>> g = 2*x*y + x + 3

682

683 >>> R.dmp_zz_collins_resultant(f, g)

684 -2*y**2 - 5*y + 1

685

686 """

687

688 n = dmp_degree(f, u)

689 m = dmp_degree(g, u)

690

691 if n < 0 or m < 0:

692 return dmp_zero(u - 1)

693

694 A = dmp_max_norm(f, u, K)

695 B = dmp_max_norm(g, u, K)

696

697 a = dmp_ground_LC(f, u, K)

698 b = dmp_ground_LC(g, u, K)

699

700 v = u - 1

701

702 B = K(2)*K.factorial(K(n + m))*A**m*B**n

703 r, p, P = dmp_zero(v), K.one, K.one

704

705 from sympy.ntheory import nextprime

706

707 while P <= B:

708 p = K(nextprime(p))

709

710 while not (a % p) or not (b % p):

711 p = K(nextprime(p))

712

713 F = dmp_ground_trunc(f, p, u, K)

714 G = dmp_ground_trunc(g, p, u, K)

715

716 try:

717 R = dmp_zz_modular_resultant(F, G, p, u, K)

718 except HomomorphismFailed:

719 continue

720

721 if K.is_one(P):

722 r = R

723 else:

724 r = dmp_apply_pairs(r, R, _collins_crt, (P, p, K), v, K)

726 P *= p

728 return r

731def dmp_qq_collins_resultant(f, g, u, K0):

732 """

733 Collins's modular resultant algorithm in `Q[X]`.

734

735 Examples

736 ========

737

738 >>> from sympy.polys import ring, QQ

739 >>> R, x,y = ring("x,y", QQ)

740

741 >>> f = QQ(1,2)*x + y + QQ(2,3)

742 >>> g = 2*x*y + x + 3

743

744 >>> R.dmp_qq_collins_resultant(f, g)

745 -2*y**2 - 7/3*y + 5/6

746

747 """

748 n = dmp_degree(f, u)

749 m = dmp_degree(g, u)

750

751 if n < 0 or m < 0:

752 return dmp_zero(u - 1)

753

754 K1 = K0.get_ring()

755

756 cf, f = dmp_clear_denoms(f, u, K0, K1)

757 cg, g = dmp_clear_denoms(g, u, K0, K1)

758

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

760 g = dmp_convert(g, u, K0, K1)

761

762 r = dmp_zz_collins_resultant(f, g, u, K1)

763 r = dmp_convert(r, u - 1, K1, K0)

764

765 c = K0.convert(cf**m * cg**n, K1)

766

767 return dmp_quo_ground(r, c, u - 1, K0)

768

769

770def dmp_resultant(f, g, u, K, includePRS=False):

771 """

772 Computes resultant of two polynomials in `K[X]`.

773

774 Examples

775 ========

776

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

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

779

780 >>> f = 3*x**2*y - y**3 - 4

781 >>> g = x**2 + x*y**3 - 9

782

783 >>> R.dmp_resultant(f, g)

784 -3*y**10 - 12*y**7 + y**6 - 54*y**4 + 8*y**3 + 729*y**2 - 216*y + 16

785

786 """

787 if not u:

788 return dup_resultant(f, g, K, includePRS=includePRS)

789

790 if includePRS:

791 return dmp_prs_resultant(f, g, u, K)

792

793 if K.is_Field:

794 if K.is_QQ and query('USE_COLLINS_RESULTANT'):

795 return dmp_qq_collins_resultant(f, g, u, K)

796 else:

797 if K.is_ZZ and query('USE_COLLINS_RESULTANT'):

798 return dmp_zz_collins_resultant(f, g, u, K)

799

800 return dmp_prs_resultant(f, g, u, K)[0]

801

802

803def dup_discriminant(f, K):

804 """

805 Computes discriminant of a polynomial in `K[x]`.

806

807 Examples

808 ========

809

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

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

812

813 >>> R.dup_discriminant(x**2 + 2*x + 3)

814 -8

815

816 """

817 d = dup_degree(f)

818

819 if d <= 0:

820 return K.zero

821 else:

822 s = (-1)**((d*(d - 1)) // 2)

823 c = dup_LC(f, K)

824

825 r = dup_resultant(f, dup_diff(f, 1, K), K)

826

827 return K.quo(r, c*K(s))

828

829

830def dmp_discriminant(f, u, K):

831 """

832 Computes discriminant of a polynomial in `K[X]`.

833

834 Examples

835 ========

836

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

838 >>> R, x,y,z,t = ring("x,y,z,t", ZZ)

839

840 >>> R.dmp_discriminant(x**2*y + x*z + t)

841 -4*y*t + z**2

842

843 """

844 if not u:

845 return dup_discriminant(f, K)

846

847 d, v = dmp_degree(f, u), u - 1

848

849 if d <= 0:

850 return dmp_zero(v)

851 else:

852 s = (-1)**((d*(d - 1)) // 2)

853 c = dmp_LC(f, K)

854

855 r = dmp_resultant(f, dmp_diff(f, 1, u, K), u, K)

856 c = dmp_mul_ground(c, K(s), v, K)

857

858 return dmp_quo(r, c, v, K)

859

860

861def _dup_rr_trivial_gcd(f, g, K):

862 """Handle trivial cases in GCD algorithm over a ring. """

863 if not (f or g):

864 return [], [], []

865 elif not f:

866 if K.is_nonnegative(dup_LC(g, K)):

867 return g, [], [K.one]

868 else:

869 return dup_neg(g, K), [], [-K.one]

870 elif not g:

871 if K.is_nonnegative(dup_LC(f, K)):

872 return f, [K.one], []

873 else:

874 return dup_neg(f, K), [-K.one], []

875

876 return None

877

878

879def _dup_ff_trivial_gcd(f, g, K):

880 """Handle trivial cases in GCD algorithm over a field. """

881 if not (f or g):

882 return [], [], []

883 elif not f:

884 return dup_monic(g, K), [], [dup_LC(g, K)]

885 elif not g:

886 return dup_monic(f, K), [dup_LC(f, K)], []

887 else:

888 return None

889

890

891def _dmp_rr_trivial_gcd(f, g, u, K):

892 """Handle trivial cases in GCD algorithm over a ring. """

893 zero_f = dmp_zero_p(f, u)

894 zero_g = dmp_zero_p(g, u)

895 if_contain_one = dmp_one_p(f, u, K) or dmp_one_p(g, u, K)

896

897 if zero_f and zero_g:

898 return tuple(dmp_zeros(3, u, K))

899 elif zero_f:

900 if K.is_nonnegative(dmp_ground_LC(g, u, K)):

901 return g, dmp_zero(u), dmp_one(u, K)

902 else:

903 return dmp_neg(g, u, K), dmp_zero(u), dmp_ground(-K.one, u)

904 elif zero_g:

905 if K.is_nonnegative(dmp_ground_LC(f, u, K)):

906 return f, dmp_one(u, K), dmp_zero(u)

907 else:

908 return dmp_neg(f, u, K), dmp_ground(-K.one, u), dmp_zero(u)

909 elif if_contain_one:

910 return dmp_one(u, K), f, g

911 elif query('USE_SIMPLIFY_GCD'):

912 return _dmp_simplify_gcd(f, g, u, K)

913 else:

914 return None

915

916

917def _dmp_ff_trivial_gcd(f, g, u, K):

918 """Handle trivial cases in GCD algorithm over a field. """

919 zero_f = dmp_zero_p(f, u)

920 zero_g = dmp_zero_p(g, u)

921

922 if zero_f and zero_g:

923 return tuple(dmp_zeros(3, u, K))

924 elif zero_f:

925 return (dmp_ground_monic(g, u, K),

926 dmp_zero(u),

927 dmp_ground(dmp_ground_LC(g, u, K), u))

928 elif zero_g:

929 return (dmp_ground_monic(f, u, K),

930 dmp_ground(dmp_ground_LC(f, u, K), u),

931 dmp_zero(u))

932 elif query('USE_SIMPLIFY_GCD'):

933 return _dmp_simplify_gcd(f, g, u, K)

934 else:

935 return None

936

937

938def _dmp_simplify_gcd(f, g, u, K):

939 """Try to eliminate `x_0` from GCD computation in `K[X]`. """

940 df = dmp_degree(f, u)

941 dg = dmp_degree(g, u)

942

943 if df > 0 and dg > 0:

944 return None

945

946 if not (df or dg):

947 F = dmp_LC(f, K)

948 G = dmp_LC(g, K)

949 else:

950 if not df:

951 F = dmp_LC(f, K)

952 G = dmp_content(g, u, K)

953 else:

954 F = dmp_content(f, u, K)

955 G = dmp_LC(g, K)

956

957 v = u - 1

958 h = dmp_gcd(F, G, v, K)

959

960 cff = [ dmp_quo(cf, h, v, K) for cf in f ]

961 cfg = [ dmp_quo(cg, h, v, K) for cg in g ]

962

963 return [h], cff, cfg

964

965

966def dup_rr_prs_gcd(f, g, K):

967 """

968 Computes polynomial GCD using subresultants over a ring.

969

970 Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``, ``cff = quo(f, h)``,

971 and ``cfg = quo(g, h)``.

972

973 Examples

974 ========

975

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

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

978

979 >>> R.dup_rr_prs_gcd(x**2 - 1, x**2 - 3*x + 2)

980 (x - 1, x + 1, x - 2)

981

982 """

983 result = _dup_rr_trivial_gcd(f, g, K)

984

985 if result is not None:

986 return result

987

988 fc, F = dup_primitive(f, K)

989 gc, G = dup_primitive(g, K)

990

991 c = K.gcd(fc, gc)

992

993 h = dup_subresultants(F, G, K)[-1]

994 _, h = dup_primitive(h, K)

995

996 c *= K.canonical_unit(dup_LC(h, K))

997

998 h = dup_mul_ground(h, c, K)

999

1000 cff = dup_quo(f, h, K)

1001 cfg = dup_quo(g, h, K)

1002

1003 return h, cff, cfg

1004

1005

1006def dup_ff_prs_gcd(f, g, K):

1007 """

1008 Computes polynomial GCD using subresultants over a field.

1009

1010 Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``, ``cff = quo(f, h)``,

1011 and ``cfg = quo(g, h)``.

1012

1013 Examples

1014 ========

1015

1016 >>> from sympy.polys import ring, QQ

1017 >>> R, x = ring("x", QQ)

1018

1019 >>> R.dup_ff_prs_gcd(x**2 - 1, x**2 - 3*x + 2)

1020 (x - 1, x + 1, x - 2)

1021

1022 """

1023 result = _dup_ff_trivial_gcd(f, g, K)

1024

1025 if result is not None:

1026 return result

1027

1028 h = dup_subresultants(f, g, K)[-1]

1029 h = dup_monic(h, K)

1030

1031 cff = dup_quo(f, h, K)

1032 cfg = dup_quo(g, h, K)

1033

1034 return h, cff, cfg

1035

1036

1037def dmp_rr_prs_gcd(f, g, u, K):

1038 """

1039 Computes polynomial GCD using subresultants over a ring.

1040

1041 Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``, ``cff = quo(f, h)``,

1042 and ``cfg = quo(g, h)``.

1043

1044 Examples

1045 ========

1046

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

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

1049

1050 >>> f = x**2 + 2*x*y + y**2

1051 >>> g = x**2 + x*y

1052

1053 >>> R.dmp_rr_prs_gcd(f, g)

1054 (x + y, x + y, x)

1055

1056 """

1057 if not u:

1058 return dup_rr_prs_gcd(f, g, K)

1059

1060 result = _dmp_rr_trivial_gcd(f, g, u, K)

1061

1062 if result is not None:

1063 return result

1064

1065 fc, F = dmp_primitive(f, u, K)

1066 gc, G = dmp_primitive(g, u, K)

1067

1068 h = dmp_subresultants(F, G, u, K)[-1]

1069 c, _, _ = dmp_rr_prs_gcd(fc, gc, u - 1, K)

1070

1071 if K.is_negative(dmp_ground_LC(h, u, K)):

1072 h = dmp_neg(h, u, K)

1073

1074 _, h = dmp_primitive(h, u, K)

1075 h = dmp_mul_term(h, c, 0, u, K)

1076

1077 cff = dmp_quo(f, h, u, K)

1078 cfg = dmp_quo(g, h, u, K)

1079

1080 return h, cff, cfg

1081

1082

1083def dmp_ff_prs_gcd(f, g, u, K):

1084 """

1085 Computes polynomial GCD using subresultants over a field.

1086

1087 Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``, ``cff = quo(f, h)``,

1088 and ``cfg = quo(g, h)``.

1089

1090 Examples

1091 ========

1092

1093 >>> from sympy.polys import ring, QQ

1094 >>> R, x,y, = ring("x,y", QQ)

1095

1096 >>> f = QQ(1,2)*x**2 + x*y + QQ(1,2)*y**2

1097 >>> g = x**2 + x*y

1098

1099 >>> R.dmp_ff_prs_gcd(f, g)

1100 (x + y, 1/2*x + 1/2*y, x)

1101

1102 """

1103 if not u:

1104 return dup_ff_prs_gcd(f, g, K)

1105

1106 result = _dmp_ff_trivial_gcd(f, g, u, K)

1107

1108 if result is not None:

1109 return result

1110

1111 fc, F = dmp_primitive(f, u, K)

1112 gc, G = dmp_primitive(g, u, K)

1113

1114 h = dmp_subresultants(F, G, u, K)[-1]

1115 c, _, _ = dmp_ff_prs_gcd(fc, gc, u - 1, K)

1116

1117 _, h = dmp_primitive(h, u, K)

1118 h = dmp_mul_term(h, c, 0, u, K)

1119 h = dmp_ground_monic(h, u, K)

1120

1121 cff = dmp_quo(f, h, u, K)

1122 cfg = dmp_quo(g, h, u, K)

1123

1124 return h, cff, cfg

1125

1126HEU_GCD_MAX = 6

1127

1128

1129def _dup_zz_gcd_interpolate(h, x, K):

1130 """Interpolate polynomial GCD from integer GCD. """

1131 f = []

1132

1133 while h:

1134 g = h % x

1135

1136 if g > x // 2:

1137 g -= x

1138

1139 f.insert(0, g)

1140 h = (h - g) // x

1141

1142 return f

1143

1144

1145def dup_zz_heu_gcd(f, g, K):

1146 """

1147 Heuristic polynomial GCD in `Z[x]`.

1148

1149 Given univariate polynomials `f` and `g` in `Z[x]`, returns

1150 their GCD and cofactors, i.e. polynomials ``h``, ``cff`` and ``cfg``

1151 such that::

1152

1153 h = gcd(f, g), cff = quo(f, h) and cfg = quo(g, h)

1154

1155 The algorithm is purely heuristic which means it may fail to compute

1156 the GCD. This will be signaled by raising an exception. In this case

1157 you will need to switch to another GCD method.

1158

1159 The algorithm computes the polynomial GCD by evaluating polynomials

1160 f and g at certain points and computing (fast) integer GCD of those

1161 evaluations. The polynomial GCD is recovered from the integer image

1162 by interpolation. The final step is to verify if the result is the

1163 correct GCD. This gives cofactors as a side effect.

1164

1165 Examples

1166 ========

1167

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

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

1170

1171 >>> R.dup_zz_heu_gcd(x**2 - 1, x**2 - 3*x + 2)

1172 (x - 1, x + 1, x - 2)

1173

1174 References

1175 ==========

1176

1177 .. [1] [Liao95]_

1178

1179 """

1180 result = _dup_rr_trivial_gcd(f, g, K)

1181

1182 if result is not None:

1183 return result

1184

1185 df = dup_degree(f)

1186 dg = dup_degree(g)

1187

1188 gcd, f, g = dup_extract(f, g, K)

1189

1190 if df == 0 or dg == 0:

1191 return [gcd], f, g

1192

1193 f_norm = dup_max_norm(f, K)

1194 g_norm = dup_max_norm(g, K)

1195

1196 B = K(2*min(f_norm, g_norm) + 29)

1197

1198 x = max(min(B, 99*K.sqrt(B)),

1199 2*min(f_norm // abs(dup_LC(f, K)),

1200 g_norm // abs(dup_LC(g, K))) + 2)

1201

1202 for i in range(0, HEU_GCD_MAX):

1203 ff = dup_eval(f, x, K)

1204 gg = dup_eval(g, x, K)

1205

1206 if ff and gg:

1207 h = K.gcd(ff, gg)

1208

1209 cff = ff // h

1210 cfg = gg // h

1211

1212 h = _dup_zz_gcd_interpolate(h, x, K)

1213 h = dup_primitive(h, K)[1]

1214

1215 cff_, r = dup_div(f, h, K)

1216

1217 if not r:

1218 cfg_, r = dup_div(g, h, K)

1219

1220 if not r:

1221 h = dup_mul_ground(h, gcd, K)

1222 return h, cff_, cfg_

1223

1224 cff = _dup_zz_gcd_interpolate(cff, x, K)

1225

1226 h, r = dup_div(f, cff, K)

1227

1228 if not r:

1229 cfg_, r = dup_div(g, h, K)

1230

1231 if not r:

1232 h = dup_mul_ground(h, gcd, K)

1233 return h, cff, cfg_

1234

1235 cfg = _dup_zz_gcd_interpolate(cfg, x, K)

1236

1237 h, r = dup_div(g, cfg, K)

1238

1239 if not r:

1240 cff_, r = dup_div(f, h, K)

1241

1242 if not r:

1243 h = dup_mul_ground(h, gcd, K)

1244 return h, cff_, cfg

1245

1246 x = 73794*x * K.sqrt(K.sqrt(x)) // 27011

1247

1248 raise HeuristicGCDFailed('no luck')

1249

1250

1251def _dmp_zz_gcd_interpolate(h, x, v, K):

1252 """Interpolate polynomial GCD from integer GCD. """

1253 f = []

1254

1255 while not dmp_zero_p(h, v):

1256 g = dmp_ground_trunc(h, x, v, K)

1257 f.insert(0, g)

1258

1259 h = dmp_sub(h, g, v, K)

1260 h = dmp_quo_ground(h, x, v, K)

1261

1262 if K.is_negative(dmp_ground_LC(f, v + 1, K)):

1263 return dmp_neg(f, v + 1, K)

1264 else:

1265 return f

1266

1267

1268def dmp_zz_heu_gcd(f, g, u, K):

1269 """

1270 Heuristic polynomial GCD in `Z[X]`.

1271

1272 Given univariate polynomials `f` and `g` in `Z[X]`, returns

1273 their GCD and cofactors, i.e. polynomials ``h``, ``cff`` and ``cfg``

1274 such that::

1275

1276 h = gcd(f, g), cff = quo(f, h) and cfg = quo(g, h)

1277

1278 The algorithm is purely heuristic which means it may fail to compute

1279 the GCD. This will be signaled by raising an exception. In this case

1280 you will need to switch to another GCD method.

1281

1282 The algorithm computes the polynomial GCD by evaluating polynomials

1283 f and g at certain points and computing (fast) integer GCD of those

1284 evaluations. The polynomial GCD is recovered from the integer image

1285 by interpolation. The evaluation process reduces f and g variable by

1286 variable into a large integer. The final step is to verify if the

1287 interpolated polynomial is the correct GCD. This gives cofactors of

1288 the input polynomials as a side effect.

1289

1290 Examples

1291 ========

1292

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

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

1295

1296 >>> f = x**2 + 2*x*y + y**2

1297 >>> g = x**2 + x*y

1298

1299 >>> R.dmp_zz_heu_gcd(f, g)

1300 (x + y, x + y, x)

1301

1302 References

1303 ==========

1304

1305 .. [1] [Liao95]_

1306

1307 """

1308 if not u:

1309 return dup_zz_heu_gcd(f, g, K)

1310

1311 result = _dmp_rr_trivial_gcd(f, g, u, K)

1312

1313 if result is not None:

1314 return result

1315

1316 gcd, f, g = dmp_ground_extract(f, g, u, K)

1317

1318 f_norm = dmp_max_norm(f, u, K)

1319 g_norm = dmp_max_norm(g, u, K)

1320

1321 B = K(2*min(f_norm, g_norm) + 29)

1322

1323 x = max(min(B, 99*K.sqrt(B)),

1324 2*min(f_norm // abs(dmp_ground_LC(f, u, K)),

1325 g_norm // abs(dmp_ground_LC(g, u, K))) + 2)

1326

1327 for i in range(0, HEU_GCD_MAX):

1328 ff = dmp_eval(f, x, u, K)

1329 gg = dmp_eval(g, x, u, K)

1330

1331 v = u - 1

1332

1333 if not (dmp_zero_p(ff, v) or dmp_zero_p(gg, v)):

1334 h, cff, cfg = dmp_zz_heu_gcd(ff, gg, v, K)

1335

1336 h = _dmp_zz_gcd_interpolate(h, x, v, K)

1337 h = dmp_ground_primitive(h, u, K)[1]

1338

1339 cff_, r = dmp_div(f, h, u, K)

1340

1341 if dmp_zero_p(r, u):

1342 cfg_, r = dmp_div(g, h, u, K)

1343

1344 if dmp_zero_p(r, u):

1345 h = dmp_mul_ground(h, gcd, u, K)

1346 return h, cff_, cfg_

1347

1348 cff = _dmp_zz_gcd_interpolate(cff, x, v, K)

1349

1350 h, r = dmp_div(f, cff, u, K)

1351

1352 if dmp_zero_p(r, u):

1353 cfg_, r = dmp_div(g, h, u, K)

1354

1355 if dmp_zero_p(r, u):

1356 h = dmp_mul_ground(h, gcd, u, K)

1357 return h, cff, cfg_

1358

1359 cfg = _dmp_zz_gcd_interpolate(cfg, x, v, K)

1360

1361 h, r = dmp_div(g, cfg, u, K)

1362

1363 if dmp_zero_p(r, u):

1364 cff_, r = dmp_div(f, h, u, K)

1365

1366 if dmp_zero_p(r, u):

1367 h = dmp_mul_ground(h, gcd, u, K)

1368 return h, cff_, cfg

1369

1370 x = 73794*x * K.sqrt(K.sqrt(x)) // 27011

1371

1372 raise HeuristicGCDFailed('no luck')

1373

1374

1375def dup_qq_heu_gcd(f, g, K0):

1376 """

1377 Heuristic polynomial GCD in `Q[x]`.

1378

1379 Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``,

1380 ``cff = quo(f, h)``, and ``cfg = quo(g, h)``.

1381

1382 Examples

1383 ========

1384

1385 >>> from sympy.polys import ring, QQ

1386 >>> R, x = ring("x", QQ)

1387

1388 >>> f = QQ(1,2)*x**2 + QQ(7,4)*x + QQ(3,2)

1389 >>> g = QQ(1,2)*x**2 + x

1390

1391 >>> R.dup_qq_heu_gcd(f, g)

1392 (x + 2, 1/2*x + 3/4, 1/2*x)

1393

1394 """

1395 result = _dup_ff_trivial_gcd(f, g, K0)

1396

1397 if result is not None:

1398 return result

1399

1400 K1 = K0.get_ring()

1401

1402 cf, f = dup_clear_denoms(f, K0, K1)

1403 cg, g = dup_clear_denoms(g, K0, K1)

1404

1405 f = dup_convert(f, K0, K1)

1406 g = dup_convert(g, K0, K1)

1407

1408 h, cff, cfg = dup_zz_heu_gcd(f, g, K1)

1409

1410 h = dup_convert(h, K1, K0)

1411

1412 c = dup_LC(h, K0)

1413 h = dup_monic(h, K0)

1414

1415 cff = dup_convert(cff, K1, K0)

1416 cfg = dup_convert(cfg, K1, K0)

1417

1418 cff = dup_mul_ground(cff, K0.quo(c, cf), K0)

1419 cfg = dup_mul_ground(cfg, K0.quo(c, cg), K0)

1420

1421 return h, cff, cfg

1422

1423

1424def dmp_qq_heu_gcd(f, g, u, K0):

1425 """

1426 Heuristic polynomial GCD in `Q[X]`.

1427

1428 Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``,

1429 ``cff = quo(f, h)``, and ``cfg = quo(g, h)``.

1430

1431 Examples

1432 ========

1433

1434 >>> from sympy.polys import ring, QQ

1435 >>> R, x,y, = ring("x,y", QQ)

1436

1437 >>> f = QQ(1,4)*x**2 + x*y + y**2

1438 >>> g = QQ(1,2)*x**2 + x*y

1439

1440 >>> R.dmp_qq_heu_gcd(f, g)

1441 (x + 2*y, 1/4*x + 1/2*y, 1/2*x)

1442

1443 """

1444 result = _dmp_ff_trivial_gcd(f, g, u, K0)

1445

1446 if result is not None:

1447 return result

1448

1449 K1 = K0.get_ring()

1450

1451 cf, f = dmp_clear_denoms(f, u, K0, K1)

1452 cg, g = dmp_clear_denoms(g, u, K0, K1)

1453

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

1455 g = dmp_convert(g, u, K0, K1)

1456

1457 h, cff, cfg = dmp_zz_heu_gcd(f, g, u, K1)

1458

1459 h = dmp_convert(h, u, K1, K0)

1460

1461 c = dmp_ground_LC(h, u, K0)

1462 h = dmp_ground_monic(h, u, K0)

1463

1464 cff = dmp_convert(cff, u, K1, K0)

1465 cfg = dmp_convert(cfg, u, K1, K0)

1466

1467 cff = dmp_mul_ground(cff, K0.quo(c, cf), u, K0)

1468 cfg = dmp_mul_ground(cfg, K0.quo(c, cg), u, K0)

1469

1470 return h, cff, cfg

1471

1472

1473def dup_inner_gcd(f, g, K):

1474 """

1475 Computes polynomial GCD and cofactors of `f` and `g` in `K[x]`.

1476

1477 Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``,

1478 ``cff = quo(f, h)``, and ``cfg = quo(g, h)``.

1479

1480 Examples

1481 ========

1482

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

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

1485

1486 >>> R.dup_inner_gcd(x**2 - 1, x**2 - 3*x + 2)

1487 (x - 1, x + 1, x - 2)

1488

1489 """

1490 if not K.is_Exact:

1491 try:

1492 exact = K.get_exact()

1493 except DomainError:

1494 return [K.one], f, g

1495

1496 f = dup_convert(f, K, exact)

1497 g = dup_convert(g, K, exact)

1498

1499 h, cff, cfg = dup_inner_gcd(f, g, exact)

1500

1501 h = dup_convert(h, exact, K)

1502 cff = dup_convert(cff, exact, K)

1503 cfg = dup_convert(cfg, exact, K)

1504

1505 return h, cff, cfg

1506 elif K.is_Field:

1507 if K.is_QQ and query('USE_HEU_GCD'):

1508 try:

1509 return dup_qq_heu_gcd(f, g, K)

1510 except HeuristicGCDFailed:

1511 pass

1512

1513 return dup_ff_prs_gcd(f, g, K)

1514 else:

1515 if K.is_ZZ and query('USE_HEU_GCD'):

1516 try:

1517 return dup_zz_heu_gcd(f, g, K)

1518 except HeuristicGCDFailed:

1519 pass

1520

1521 return dup_rr_prs_gcd(f, g, K)

1522

1523

1524def _dmp_inner_gcd(f, g, u, K):

1525 """Helper function for `dmp_inner_gcd()`. """

1526 if not K.is_Exact:

1527 try:

1528 exact = K.get_exact()

1529 except DomainError:

1530 return dmp_one(u, K), f, g

1531

1532 f = dmp_convert(f, u, K, exact)

1533 g = dmp_convert(g, u, K, exact)

1534

1535 h, cff, cfg = _dmp_inner_gcd(f, g, u, exact)

1536

1537 h = dmp_convert(h, u, exact, K)

1538 cff = dmp_convert(cff, u, exact, K)

1539 cfg = dmp_convert(cfg, u, exact, K)

1540

1541 return h, cff, cfg

1542 elif K.is_Field:

1543 if K.is_QQ and query('USE_HEU_GCD'):

1544 try:

1545 return dmp_qq_heu_gcd(f, g, u, K)

1546 except HeuristicGCDFailed:

1547 pass

1548

1549 return dmp_ff_prs_gcd(f, g, u, K)

1550 else:

1551 if K.is_ZZ and query('USE_HEU_GCD'):

1552 try:

1553 return dmp_zz_heu_gcd(f, g, u, K)

1554 except HeuristicGCDFailed:

1555 pass

1556

1557 return dmp_rr_prs_gcd(f, g, u, K)

1558

1559

1560def dmp_inner_gcd(f, g, u, K):

1561 """

1562 Computes polynomial GCD and cofactors of `f` and `g` in `K[X]`.

1563

1564 Returns ``(h, cff, cfg)`` such that ``a = gcd(f, g)``,

1565 ``cff = quo(f, h)``, and ``cfg = quo(g, h)``.

1566

1567 Examples

1568 ========

1569

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

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

1572

1573 >>> f = x**2 + 2*x*y + y**2

1574 >>> g = x**2 + x*y

1575

1576 >>> R.dmp_inner_gcd(f, g)

1577 (x + y, x + y, x)

1578

1579 """

1580 if not u:

1581 return dup_inner_gcd(f, g, K)

1582

1583 J, (f, g) = dmp_multi_deflate((f, g), u, K)

1584 h, cff, cfg = _dmp_inner_gcd(f, g, u, K)

1585

1586 return (dmp_inflate(h, J, u, K),

1587 dmp_inflate(cff, J, u, K),

1588 dmp_inflate(cfg, J, u, K))

1589

1590

1591def dup_gcd(f, g, K):

1592 """

1593 Computes polynomial GCD of `f` and `g` in `K[x]`.

1594

1595 Examples

1596 ========

1597

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

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

1600

1601 >>> R.dup_gcd(x**2 - 1, x**2 - 3*x + 2)

1602 x - 1

1603

1604 """

1605 return dup_inner_gcd(f, g, K)[0]

1606

1607

1608def dmp_gcd(f, g, u, K):

1609 """

1610 Computes polynomial GCD of `f` and `g` in `K[X]`.

1611

1612 Examples

1613 ========

1614

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

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

1617

1618 >>> f = x**2 + 2*x*y + y**2

1619 >>> g = x**2 + x*y

1620

1621 >>> R.dmp_gcd(f, g)

1622 x + y

1623

1624 """

1625 return dmp_inner_gcd(f, g, u, K)[0]

1626

1627

1628def dup_rr_lcm(f, g, K):

1629 """

1630 Computes polynomial LCM over a ring in `K[x]`.

1631

1632 Examples

1633 ========

1634

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

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

1637

1638 >>> R.dup_rr_lcm(x**2 - 1, x**2 - 3*x + 2)

1639 x**3 - 2*x**2 - x + 2

1640

1641 """

1642 fc, f = dup_primitive(f, K)

1643 gc, g = dup_primitive(g, K)

1644

1645 c = K.lcm(fc, gc)

1646

1647 h = dup_quo(dup_mul(f, g, K),

1648 dup_gcd(f, g, K), K)

1649

1650 return dup_mul_ground(h, c, K)

1651

1652

1653def dup_ff_lcm(f, g, K):

1654 """

1655 Computes polynomial LCM over a field in `K[x]`.

1656

1657 Examples

1658 ========

1659

1660 >>> from sympy.polys import ring, QQ

1661 >>> R, x = ring("x", QQ)

1662

1663 >>> f = QQ(1,2)*x**2 + QQ(7,4)*x + QQ(3,2)

1664 >>> g = QQ(1,2)*x**2 + x

1665

1666 >>> R.dup_ff_lcm(f, g)

1667 x**3 + 7/2*x**2 + 3*x

1668

1669 """

1670 h = dup_quo(dup_mul(f, g, K),

1671 dup_gcd(f, g, K), K)

1672

1673 return dup_monic(h, K)

1674

1675

1676def dup_lcm(f, g, K):

1677 """

1678 Computes polynomial LCM of `f` and `g` in `K[x]`.

1679

1680 Examples

1681 ========

1682

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

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

1685

1686 >>> R.dup_lcm(x**2 - 1, x**2 - 3*x + 2)

1687 x**3 - 2*x**2 - x + 2

1688

1689 """

1690 if K.is_Field:

1691 return dup_ff_lcm(f, g, K)

1692 else:

1693 return dup_rr_lcm(f, g, K)

1694

1695

1696def dmp_rr_lcm(f, g, u, K):

1697 """

1698 Computes polynomial LCM over a ring in `K[X]`.

1699

1700 Examples

1701 ========

1702

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

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

1705

1706 >>> f = x**2 + 2*x*y + y**2

1707 >>> g = x**2 + x*y

1708

1709 >>> R.dmp_rr_lcm(f, g)

1710 x**3 + 2*x**2*y + x*y**2

1711

1712 """

1713 fc, f = dmp_ground_primitive(f, u, K)

1714 gc, g = dmp_ground_primitive(g, u, K)

1715

1716 c = K.lcm(fc, gc)

1717

1718 h = dmp_quo(dmp_mul(f, g, u, K),

1719 dmp_gcd(f, g, u, K), u, K)

1720

1721 return dmp_mul_ground(h, c, u, K)

1722

1723

1724def dmp_ff_lcm(f, g, u, K):

1725 """

1726 Computes polynomial LCM over a field in `K[X]`.

1727

1728 Examples

1729 ========

1730

1731 >>> from sympy.polys import ring, QQ

1732 >>> R, x,y, = ring("x,y", QQ)

1733

1734 >>> f = QQ(1,4)*x**2 + x*y + y**2

1735 >>> g = QQ(1,2)*x**2 + x*y

1736

1737 >>> R.dmp_ff_lcm(f, g)

1738 x**3 + 4*x**2*y + 4*x*y**2

1739

1740 """

1741 h = dmp_quo(dmp_mul(f, g, u, K),

1742 dmp_gcd(f, g, u, K), u, K)

1743

1744 return dmp_ground_monic(h, u, K)

1745

1746

1747def dmp_lcm(f, g, u, K):

1748 """

1749 Computes polynomial LCM of `f` and `g` in `K[X]`.

1750

1751 Examples

1752 ========

1753

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

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

1756

1757 >>> f = x**2 + 2*x*y + y**2

1758 >>> g = x**2 + x*y

1759

1760 >>> R.dmp_lcm(f, g)

1761 x**3 + 2*x**2*y + x*y**2

1762

1763 """

1764 if not u:

1765 return dup_lcm(f, g, K)

1766

1767 if K.is_Field:

1768 return dmp_ff_lcm(f, g, u, K)

1769 else:

1770 return dmp_rr_lcm(f, g, u, K)

1771

1772

1773def dmp_content(f, u, K):

1774 """

1775 Returns GCD of multivariate coefficients.

1776

1777 Examples

1778 ========

1779

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

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

1782

1783 >>> R.dmp_content(2*x*y + 6*x + 4*y + 12)

1784 2*y + 6

1785

1786 """

1787 cont, v = dmp_LC(f, K), u - 1

1788

1789 if dmp_zero_p(f, u):

1790 return cont

1791

1792 for c in f[1:]:

1793 cont = dmp_gcd(cont, c, v, K)

1794

1795 if dmp_one_p(cont, v, K):

1796 break

1797

1798 if K.is_negative(dmp_ground_LC(cont, v, K)):

1799 return dmp_neg(cont, v, K)

1800 else:

1801 return cont

1802

1803

1804def dmp_primitive(f, u, K):

1805 """

1806 Returns multivariate content and a primitive polynomial.

1807

1808 Examples

1809 ========

1810

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

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

1813

1814 >>> R.dmp_primitive(2*x*y + 6*x + 4*y + 12)

1815 (2*y + 6, x + 2)

1816

1817 """

1818 cont, v = dmp_content(f, u, K), u - 1

1819

1820 if dmp_zero_p(f, u) or dmp_one_p(cont, v, K):

1821 return cont, f

1822 else:

1823 return cont, [ dmp_quo(c, cont, v, K) for c in f ]

1824

1825

1826def dup_cancel(f, g, K, include=True):

1827 """

1828 Cancel common factors in a rational function `f/g`.

1829

1830 Examples

1831 ========

1832

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

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

1835

1836 >>> R.dup_cancel(2*x**2 - 2, x**2 - 2*x + 1)

1837 (2*x + 2, x - 1)

1838

1839 """

1840 return dmp_cancel(f, g, 0, K, include=include)

1841

1842

1843def dmp_cancel(f, g, u, K, include=True):

1844 """

1845 Cancel common factors in a rational function `f/g`.

1846

1847 Examples

1848 ========

1849

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

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

1852

1853 >>> R.dmp_cancel(2*x**2 - 2, x**2 - 2*x + 1)

1854 (2*x + 2, x - 1)

1855

1856 """

1857 K0 = None

1858

1859 if K.is_Field and K.has_assoc_Ring:

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

1861

1862 cq, f = dmp_clear_denoms(f, u, K0, K, convert=True)

1863 cp, g = dmp_clear_denoms(g, u, K0, K, convert=True)

1864 else:

1865 cp, cq = K.one, K.one

1866

1867 _, p, q = dmp_inner_gcd(f, g, u, K)

1868

1869 if K0 is not None:

1870 _, cp, cq = K.cofactors(cp, cq)

1871

1872 p = dmp_convert(p, u, K, K0)

1873 q = dmp_convert(q, u, K, K0)

1874

1875 K = K0

1876

1877 p_neg = K.is_negative(dmp_ground_LC(p, u, K))

1878 q_neg = K.is_negative(dmp_ground_LC(q, u, K))

1879

1880 if p_neg and q_neg:

1881 p, q = dmp_neg(p, u, K), dmp_neg(q, u, K)

1882 elif p_neg:

1883 cp, p = -cp, dmp_neg(p, u, K)

1884 elif q_neg:

1885 cp, q = -cp, dmp_neg(q, u, K)

1886

1887 if not include:

1888 return cp, cq, p, q

1889

1890 p = dmp_mul_ground(p, cp, u, K)

1891 q = dmp_mul_ground(q, cq, u, K)

1892

1893 return p, q