Coverage for /usr/lib/python3/dist-packages/sympy/polys/densetools.py: 12%

1"""Advanced tools for dense recursive polynomials in ``K[x]`` or ``K[X]``. """

4from sympy.polys.densearith import (

5 dup_add_term, dmp_add_term,

6 dup_lshift,

7 dup_add, dmp_add,

8 dup_sub, dmp_sub,

9 dup_mul, dmp_mul,

10 dup_sqr,

11 dup_div,

12 dup_rem, dmp_rem,

13 dmp_expand,

14 dup_mul_ground, dmp_mul_ground,

15 dup_quo_ground, dmp_quo_ground,

16 dup_exquo_ground, dmp_exquo_ground,

17)

18from sympy.polys.densebasic import (

19 dup_strip, dmp_strip,

20 dup_convert, dmp_convert,

21 dup_degree, dmp_degree,

22 dmp_to_dict,

23 dmp_from_dict,

24 dup_LC, dmp_LC, dmp_ground_LC,

25 dup_TC, dmp_TC,

26 dmp_zero, dmp_ground,

27 dmp_zero_p,

28 dup_to_raw_dict, dup_from_raw_dict,

29 dmp_zeros

30)

31from sympy.polys.polyerrors import (

32 MultivariatePolynomialError,

33 DomainError

34)

35from sympy.utilities import variations

37from math import ceil as _ceil, log as _log

39def dup_integrate(f, m, K):

40 """

41 Computes the indefinite integral of ``f`` in ``K[x]``.

43 Examples

44 ========

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

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

49 >>> R.dup_integrate(x**2 + 2*x, 1)

50 1/3*x**3 + x**2

51 >>> R.dup_integrate(x**2 + 2*x, 2)

52 1/12*x**4 + 1/3*x**3

54 """

55 if m <= 0 or not f:

56 return f

58 g = [K.zero]*m

60 for i, c in enumerate(reversed(f)):

61 n = i + 1

63 for j in range(1, m):

64 n *= i + j + 1

66 g.insert(0, K.exquo(c, K(n)))

68 return g

71def dmp_integrate(f, m, u, K):

72 """

73 Computes the indefinite integral of ``f`` in ``x_0`` in ``K[X]``.

75 Examples

76 ========

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

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

81 >>> R.dmp_integrate(x + 2*y, 1)

82 1/2*x**2 + 2*x*y

83 >>> R.dmp_integrate(x + 2*y, 2)

84 1/6*x**3 + x**2*y

86 """

87 if not u:

88 return dup_integrate(f, m, K)

90 if m <= 0 or dmp_zero_p(f, u):

91 return f

93 g, v = dmp_zeros(m, u - 1, K), u - 1

95 for i, c in enumerate(reversed(f)):

96 n = i + 1

98 for j in range(1, m):

99 n *= i + j + 1

100

101 g.insert(0, dmp_quo_ground(c, K(n), v, K))

102

103 return g

104

105

106def _rec_integrate_in(g, m, v, i, j, K):

107 """Recursive helper for :func:`dmp_integrate_in`."""

108 if i == j:

109 return dmp_integrate(g, m, v, K)

110

111 w, i = v - 1, i + 1

112

113 return dmp_strip([ _rec_integrate_in(c, m, w, i, j, K) for c in g ], v)

114

115

116def dmp_integrate_in(f, m, j, u, K):

117 """

118 Computes the indefinite integral of ``f`` in ``x_j`` in ``K[X]``.

119

120 Examples

121 ========

122

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

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

125

126 >>> R.dmp_integrate_in(x + 2*y, 1, 0)

127 1/2*x**2 + 2*x*y

128 >>> R.dmp_integrate_in(x + 2*y, 1, 1)

129 x*y + y**2

130

131 """

132 if j < 0 or j > u:

133 raise IndexError("0 <= j <= u expected, got u = %d, j = %d" % (u, j))

134

135 return _rec_integrate_in(f, m, u, 0, j, K)

136

137

138def dup_diff(f, m, K):

139 """

140 ``m``-th order derivative of a polynomial in ``K[x]``.

141

142 Examples

143 ========

144

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

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

147

148 >>> R.dup_diff(x**3 + 2*x**2 + 3*x + 4, 1)

149 3*x**2 + 4*x + 3

150 >>> R.dup_diff(x**3 + 2*x**2 + 3*x + 4, 2)

151 6*x + 4

152

153 """

154 if m <= 0:

155 return f

156

157 n = dup_degree(f)

158

159 if n < m:

160 return []

161

162 deriv = []

163

164 if m == 1:

165 for coeff in f[:-m]:

166 deriv.append(K(n)*coeff)

167 n -= 1

168 else:

169 for coeff in f[:-m]:

170 k = n

171

172 for i in range(n - 1, n - m, -1):

173 k *= i

174

175 deriv.append(K(k)*coeff)

176 n -= 1

177

178 return dup_strip(deriv)

179

180

181def dmp_diff(f, m, u, K):

182 """

183 ``m``-th order derivative in ``x_0`` of a polynomial in ``K[X]``.

184

185 Examples

186 ========

187

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

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

190

191 >>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1

192

193 >>> R.dmp_diff(f, 1)

194 y**2 + 2*y + 3

195 >>> R.dmp_diff(f, 2)

196 0

197

198 """

199 if not u:

200 return dup_diff(f, m, K)

201 if m <= 0:

202 return f

203

204 n = dmp_degree(f, u)

205

206 if n < m:

207 return dmp_zero(u)

208

209 deriv, v = [], u - 1

210

211 if m == 1:

212 for coeff in f[:-m]:

213 deriv.append(dmp_mul_ground(coeff, K(n), v, K))

214 n -= 1

215 else:

216 for coeff in f[:-m]:

217 k = n

218

219 for i in range(n - 1, n - m, -1):

220 k *= i

221

222 deriv.append(dmp_mul_ground(coeff, K(k), v, K))

223 n -= 1

224

225 return dmp_strip(deriv, u)

226

227

228def _rec_diff_in(g, m, v, i, j, K):

229 """Recursive helper for :func:`dmp_diff_in`."""

230 if i == j:

231 return dmp_diff(g, m, v, K)

232

233 w, i = v - 1, i + 1

234

235 return dmp_strip([ _rec_diff_in(c, m, w, i, j, K) for c in g ], v)

236

237

238def dmp_diff_in(f, m, j, u, K):

239 """

240 ``m``-th order derivative in ``x_j`` of a polynomial in ``K[X]``.

241

242 Examples

243 ========

244

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

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

247

248 >>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1

249

250 >>> R.dmp_diff_in(f, 1, 0)

251 y**2 + 2*y + 3

252 >>> R.dmp_diff_in(f, 1, 1)

253 2*x*y + 2*x + 4*y + 3

254

255 """

256 if j < 0 or j > u:

257 raise IndexError("0 <= j <= %s expected, got %s" % (u, j))

258

259 return _rec_diff_in(f, m, u, 0, j, K)

260

261

262def dup_eval(f, a, K):

263 """

264 Evaluate a polynomial at ``x = a`` in ``K[x]`` using Horner scheme.

265

266 Examples

267 ========

268

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

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

271

272 >>> R.dup_eval(x**2 + 2*x + 3, 2)

273 11

274

275 """

276 if not a:

277 return K.convert(dup_TC(f, K))

278

279 result = K.zero

280

281 for c in f:

282 result *= a

283 result += c

284

285 return result

286

287

288def dmp_eval(f, a, u, K):

289 """

290 Evaluate a polynomial at ``x_0 = a`` in ``K[X]`` using the Horner scheme.

291

292 Examples

293 ========

294

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

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

297

298 >>> R.dmp_eval(2*x*y + 3*x + y + 2, 2)

299 5*y + 8

300

301 """

302 if not u:

303 return dup_eval(f, a, K)

304

305 if not a:

306 return dmp_TC(f, K)

307

308 result, v = dmp_LC(f, K), u - 1

309

310 for coeff in f[1:]:

311 result = dmp_mul_ground(result, a, v, K)

312 result = dmp_add(result, coeff, v, K)

313

314 return result

315

316

317def _rec_eval_in(g, a, v, i, j, K):

318 """Recursive helper for :func:`dmp_eval_in`."""

319 if i == j:

320 return dmp_eval(g, a, v, K)

321

322 v, i = v - 1, i + 1

323

324 return dmp_strip([ _rec_eval_in(c, a, v, i, j, K) for c in g ], v)

325

326

327def dmp_eval_in(f, a, j, u, K):

328 """

329 Evaluate a polynomial at ``x_j = a`` in ``K[X]`` using the Horner scheme.

330

331 Examples

332 ========

333

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

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

336

337 >>> f = 2*x*y + 3*x + y + 2

338

339 >>> R.dmp_eval_in(f, 2, 0)

340 5*y + 8

341 >>> R.dmp_eval_in(f, 2, 1)

342 7*x + 4

343

344 """

345 if j < 0 or j > u:

346 raise IndexError("0 <= j <= %s expected, got %s" % (u, j))

347

348 return _rec_eval_in(f, a, u, 0, j, K)

349

350

351def _rec_eval_tail(g, i, A, u, K):

352 """Recursive helper for :func:`dmp_eval_tail`."""

353 if i == u:

354 return dup_eval(g, A[-1], K)

355 else:

356 h = [ _rec_eval_tail(c, i + 1, A, u, K) for c in g ]

357

358 if i < u - len(A) + 1:

359 return h

360 else:

361 return dup_eval(h, A[-u + i - 1], K)

362

363

364def dmp_eval_tail(f, A, u, K):

365 """

366 Evaluate a polynomial at ``x_j = a_j, ...`` in ``K[X]``.

367

368 Examples

369 ========

370

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

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

373

374 >>> f = 2*x*y + 3*x + y + 2

375

376 >>> R.dmp_eval_tail(f, [2])

377 7*x + 4

378 >>> R.dmp_eval_tail(f, [2, 2])

379 18

380

381 """

382 if not A:

383 return f

384

385 if dmp_zero_p(f, u):

386 return dmp_zero(u - len(A))

387

388 e = _rec_eval_tail(f, 0, A, u, K)

389

390 if u == len(A) - 1:

391 return e

392 else:

393 return dmp_strip(e, u - len(A))

394

395

396def _rec_diff_eval(g, m, a, v, i, j, K):

397 """Recursive helper for :func:`dmp_diff_eval`."""

398 if i == j:

399 return dmp_eval(dmp_diff(g, m, v, K), a, v, K)

400

401 v, i = v - 1, i + 1

402

403 return dmp_strip([ _rec_diff_eval(c, m, a, v, i, j, K) for c in g ], v)

404

405

406def dmp_diff_eval_in(f, m, a, j, u, K):

407 """

408 Differentiate and evaluate a polynomial in ``x_j`` at ``a`` in ``K[X]``.

409

410 Examples

411 ========

412

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

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

415

416 >>> f = x*y**2 + 2*x*y + 3*x + 2*y**2 + 3*y + 1

417

418 >>> R.dmp_diff_eval_in(f, 1, 2, 0)

419 y**2 + 2*y + 3

420 >>> R.dmp_diff_eval_in(f, 1, 2, 1)

421 6*x + 11

422

423 """

424 if j > u:

425 raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))

426 if not j:

427 return dmp_eval(dmp_diff(f, m, u, K), a, u, K)

428

429 return _rec_diff_eval(f, m, a, u, 0, j, K)

430

431

432def dup_trunc(f, p, K):

433 """

434 Reduce a ``K[x]`` polynomial modulo a constant ``p`` in ``K``.

435

436 Examples

437 ========

438

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

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

441

442 >>> R.dup_trunc(2*x**3 + 3*x**2 + 5*x + 7, ZZ(3))

443 -x**3 - x + 1

444

445 """

446 if K.is_ZZ:

447 g = []

448

449 for c in f:

450 c = c % p

451

452 if c > p // 2:

453 g.append(c - p)

454 else:

455 g.append(c)

456 else:

457 g = [ c % p for c in f ]

458

459 return dup_strip(g)

460

461

462def dmp_trunc(f, p, u, K):

463 """

464 Reduce a ``K[X]`` polynomial modulo a polynomial ``p`` in ``K[Y]``.

465

466 Examples

467 ========

468

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

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

471

472 >>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3

473 >>> g = (y - 1).drop(x)

474

475 >>> R.dmp_trunc(f, g)

476 11*x**2 + 11*x + 5

477

478 """

479 return dmp_strip([ dmp_rem(c, p, u - 1, K) for c in f ], u)

480

481

482def dmp_ground_trunc(f, p, u, K):

483 """

484 Reduce a ``K[X]`` polynomial modulo a constant ``p`` in ``K``.

485

486 Examples

487 ========

488

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

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

491

492 >>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3

493

494 >>> R.dmp_ground_trunc(f, ZZ(3))

495 -x**2 - x*y - y

496

497 """

498 if not u:

499 return dup_trunc(f, p, K)

500

501 v = u - 1

502

503 return dmp_strip([ dmp_ground_trunc(c, p, v, K) for c in f ], u)

504

505

506def dup_monic(f, K):

507 """

508 Divide all coefficients by ``LC(f)`` in ``K[x]``.

509

510 Examples

511 ========

512

513 >>> from sympy.polys import ring, ZZ, QQ

514

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

516 >>> R.dup_monic(3*x**2 + 6*x + 9)

517 x**2 + 2*x + 3

518

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

520 >>> R.dup_monic(3*x**2 + 4*x + 2)

521 x**2 + 4/3*x + 2/3

522

523 """

524 if not f:

525 return f

526

527 lc = dup_LC(f, K)

528

529 if K.is_one(lc):

530 return f

531 else:

532 return dup_exquo_ground(f, lc, K)

533

534

535def dmp_ground_monic(f, u, K):

536 """

537 Divide all coefficients by ``LC(f)`` in ``K[X]``.

538

539 Examples

540 ========

541

542 >>> from sympy.polys import ring, ZZ, QQ

543

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

545 >>> f = 3*x**2*y + 6*x**2 + 3*x*y + 9*y + 3

546

547 >>> R.dmp_ground_monic(f)

548 x**2*y + 2*x**2 + x*y + 3*y + 1

549

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

551 >>> f = 3*x**2*y + 8*x**2 + 5*x*y + 6*x + 2*y + 3

552

553 >>> R.dmp_ground_monic(f)

554 x**2*y + 8/3*x**2 + 5/3*x*y + 2*x + 2/3*y + 1

555

556 """

557 if not u:

558 return dup_monic(f, K)

559

560 if dmp_zero_p(f, u):

561 return f

562

563 lc = dmp_ground_LC(f, u, K)

564

565 if K.is_one(lc):

566 return f

567 else:

568 return dmp_exquo_ground(f, lc, u, K)

569

570

571def dup_content(f, K):

572 """

573 Compute the GCD of coefficients of ``f`` in ``K[x]``.

574

575 Examples

576 ========

577

578 >>> from sympy.polys import ring, ZZ, QQ

579

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

581 >>> f = 6*x**2 + 8*x + 12

582

583 >>> R.dup_content(f)

584 2

585

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

587 >>> f = 6*x**2 + 8*x + 12

588

589 >>> R.dup_content(f)

590 2

591

592 """

593 from sympy.polys.domains import QQ

594

595 if not f:

596 return K.zero

597

598 cont = K.zero

599

600 if K == QQ:

601 for c in f:

602 cont = K.gcd(cont, c)

603 else:

604 for c in f:

605 cont = K.gcd(cont, c)

606

607 if K.is_one(cont):

608 break

609

610 return cont

611

612

613def dmp_ground_content(f, u, K):

614 """

615 Compute the GCD of coefficients of ``f`` in ``K[X]``.

616

617 Examples

618 ========

619

620 >>> from sympy.polys import ring, ZZ, QQ

621

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

623 >>> f = 2*x*y + 6*x + 4*y + 12

624

625 >>> R.dmp_ground_content(f)

626 2

627

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

629 >>> f = 2*x*y + 6*x + 4*y + 12

630

631 >>> R.dmp_ground_content(f)

632 2

633

634 """

635 from sympy.polys.domains import QQ

636

637 if not u:

638 return dup_content(f, K)

639

640 if dmp_zero_p(f, u):

641 return K.zero

642

643 cont, v = K.zero, u - 1

644

645 if K == QQ:

646 for c in f:

647 cont = K.gcd(cont, dmp_ground_content(c, v, K))

648 else:

649 for c in f:

650 cont = K.gcd(cont, dmp_ground_content(c, v, K))

651

652 if K.is_one(cont):

653 break

654

655 return cont

656

657

658def dup_primitive(f, K):

659 """

660 Compute content and the primitive form of ``f`` in ``K[x]``.

661

662 Examples

663 ========

664

665 >>> from sympy.polys import ring, ZZ, QQ

666

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

668 >>> f = 6*x**2 + 8*x + 12

669

670 >>> R.dup_primitive(f)

671 (2, 3*x**2 + 4*x + 6)

672

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

674 >>> f = 6*x**2 + 8*x + 12

675

676 >>> R.dup_primitive(f)

677 (2, 3*x**2 + 4*x + 6)

678

679 """

680 if not f:

681 return K.zero, f

682

683 cont = dup_content(f, K)

684

685 if K.is_one(cont):

686 return cont, f

687 else:

688 return cont, dup_quo_ground(f, cont, K)

689

690

691def dmp_ground_primitive(f, u, K):

692 """

693 Compute content and the primitive form of ``f`` in ``K[X]``.

694

695 Examples

696 ========

697

698 >>> from sympy.polys import ring, ZZ, QQ

699

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

701 >>> f = 2*x*y + 6*x + 4*y + 12

702

703 >>> R.dmp_ground_primitive(f)

704 (2, x*y + 3*x + 2*y + 6)

705

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

707 >>> f = 2*x*y + 6*x + 4*y + 12

708

709 >>> R.dmp_ground_primitive(f)

710 (2, x*y + 3*x + 2*y + 6)

711

712 """

713 if not u:

714 return dup_primitive(f, K)

715

716 if dmp_zero_p(f, u):

717 return K.zero, f

718

719 cont = dmp_ground_content(f, u, K)

720

721 if K.is_one(cont):

722 return cont, f

723 else:

724 return cont, dmp_quo_ground(f, cont, u, K)

725

726

727def dup_extract(f, g, K):

728 """

729 Extract common content from a pair of polynomials in ``K[x]``.

730

731 Examples

732 ========

733

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

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

736

737 >>> R.dup_extract(6*x**2 + 12*x + 18, 4*x**2 + 8*x + 12)

738 (2, 3*x**2 + 6*x + 9, 2*x**2 + 4*x + 6)

739

740 """

741 fc = dup_content(f, K)

742 gc = dup_content(g, K)

743

744 gcd = K.gcd(fc, gc)

745

746 if not K.is_one(gcd):

747 f = dup_quo_ground(f, gcd, K)

748 g = dup_quo_ground(g, gcd, K)

749

750 return gcd, f, g

751

752

753def dmp_ground_extract(f, g, u, K):

754 """

755 Extract common content from a pair of polynomials in ``K[X]``.

756

757 Examples

758 ========

759

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

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

762

763 >>> R.dmp_ground_extract(6*x*y + 12*x + 18, 4*x*y + 8*x + 12)

764 (2, 3*x*y + 6*x + 9, 2*x*y + 4*x + 6)

765

766 """

767 fc = dmp_ground_content(f, u, K)

768 gc = dmp_ground_content(g, u, K)

769

770 gcd = K.gcd(fc, gc)

771

772 if not K.is_one(gcd):

773 f = dmp_quo_ground(f, gcd, u, K)

774 g = dmp_quo_ground(g, gcd, u, K)

775

776 return gcd, f, g

777

778

779def dup_real_imag(f, K):

780 """

781 Return bivariate polynomials ``f1`` and ``f2``, such that ``f = f1 + f2*I``.

782

783 Examples

784 ========

785

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

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

788

789 >>> R.dup_real_imag(x**3 + x**2 + x + 1)

790 (x**3 + x**2 - 3*x*y**2 + x - y**2 + 1, 3*x**2*y + 2*x*y - y**3 + y)

791

792 """

793 if not K.is_ZZ and not K.is_QQ:

794 raise DomainError("computing real and imaginary parts is not supported over %s" % K)

795

796 f1 = dmp_zero(1)

797 f2 = dmp_zero(1)

798

799 if not f:

800 return f1, f2

801

802 g = [[[K.one, K.zero]], [[K.one], []]]

803 h = dmp_ground(f[0], 2)

804

805 for c in f[1:]:

806 h = dmp_mul(h, g, 2, K)

807 h = dmp_add_term(h, dmp_ground(c, 1), 0, 2, K)

808

809 H = dup_to_raw_dict(h)

810

811 for k, h in H.items():

812 m = k % 4

813

814 if not m:

815 f1 = dmp_add(f1, h, 1, K)

816 elif m == 1:

817 f2 = dmp_add(f2, h, 1, K)

818 elif m == 2:

819 f1 = dmp_sub(f1, h, 1, K)

820 else:

821 f2 = dmp_sub(f2, h, 1, K)

822

823 return f1, f2

824

825

826def dup_mirror(f, K):

827 """

828 Evaluate efficiently the composition ``f(-x)`` in ``K[x]``.

829

830 Examples

831 ========

832

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

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

835

836 >>> R.dup_mirror(x**3 + 2*x**2 - 4*x + 2)

837 -x**3 + 2*x**2 + 4*x + 2

838

839 """

840 f = list(f)

841

842 for i in range(len(f) - 2, -1, -2):

843 f[i] = -f[i]

844

845 return f

846

847

848def dup_scale(f, a, K):

849 """

850 Evaluate efficiently composition ``f(a*x)`` in ``K[x]``.

851

852 Examples

853 ========

854

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

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

857

858 >>> R.dup_scale(x**2 - 2*x + 1, ZZ(2))

859 4*x**2 - 4*x + 1

860

861 """

862 f, n, b = list(f), len(f) - 1, a

863

864 for i in range(n - 1, -1, -1):

865 f[i], b = b*f[i], b*a

866

867 return f

868

869

870def dup_shift(f, a, K):

871 """

872 Evaluate efficiently Taylor shift ``f(x + a)`` in ``K[x]``.

873

874 Examples

875 ========

876

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

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

879

880 >>> R.dup_shift(x**2 - 2*x + 1, ZZ(2))

881 x**2 + 2*x + 1

882

883 """

884 f, n = list(f), len(f) - 1

885

886 for i in range(n, 0, -1):

887 for j in range(0, i):

888 f[j + 1] += a*f[j]

889

890 return f

891

892

893def dup_transform(f, p, q, K):

894 """

895 Evaluate functional transformation ``q**n * f(p/q)`` in ``K[x]``.

896

897 Examples

898 ========

899

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

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

902

903 >>> R.dup_transform(x**2 - 2*x + 1, x**2 + 1, x - 1)

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

905

906 """

907 if not f:

908 return []

909

910 n = len(f) - 1

911 h, Q = [f[0]], [[K.one]]

912

913 for i in range(0, n):

914 Q.append(dup_mul(Q[-1], q, K))

915

916 for c, q in zip(f[1:], Q[1:]):

917 h = dup_mul(h, p, K)

918 q = dup_mul_ground(q, c, K)

919 h = dup_add(h, q, K)

920

921 return h

922

923

924def dup_compose(f, g, K):

925 """

926 Evaluate functional composition ``f(g)`` in ``K[x]``.

927

928 Examples

929 ========

930

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

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

933

934 >>> R.dup_compose(x**2 + x, x - 1)

935 x**2 - x

936

937 """

938 if len(g) <= 1:

939 return dup_strip([dup_eval(f, dup_LC(g, K), K)])

940

941 if not f:

942 return []

943

944 h = [f[0]]

945

946 for c in f[1:]:

947 h = dup_mul(h, g, K)

948 h = dup_add_term(h, c, 0, K)

949

950 return h

951

952

953def dmp_compose(f, g, u, K):

954 """

955 Evaluate functional composition ``f(g)`` in ``K[X]``.

956

957 Examples

958 ========

959

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

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

962

963 >>> R.dmp_compose(x*y + 2*x + y, y)

964 y**2 + 3*y

965

966 """

967 if not u:

968 return dup_compose(f, g, K)

969

970 if dmp_zero_p(f, u):

971 return f

972

973 h = [f[0]]

974

975 for c in f[1:]:

976 h = dmp_mul(h, g, u, K)

977 h = dmp_add_term(h, c, 0, u, K)

978

979 return h

980

981

982def _dup_right_decompose(f, s, K):

983 """Helper function for :func:`_dup_decompose`."""

984 n = len(f) - 1

985 lc = dup_LC(f, K)

986

987 f = dup_to_raw_dict(f)

988 g = { s: K.one }

989

990 r = n // s

991

992 for i in range(1, s):

993 coeff = K.zero

994

995 for j in range(0, i):

996 if not n + j - i in f:

997 continue

998

999 if not s - j in g:

1000 continue

1001

1002 fc, gc = f[n + j - i], g[s - j]

1003 coeff += (i - r*j)*fc*gc

1004

1005 g[s - i] = K.quo(coeff, i*r*lc)

1006

1007 return dup_from_raw_dict(g, K)

1008

1009

1010def _dup_left_decompose(f, h, K):

1011 """Helper function for :func:`_dup_decompose`."""

1012 g, i = {}, 0

1013

1014 while f:

1015 q, r = dup_div(f, h, K)

1016

1017 if dup_degree(r) > 0:

1018 return None

1019 else:

1020 g[i] = dup_LC(r, K)

1021 f, i = q, i + 1

1022

1023 return dup_from_raw_dict(g, K)

1024

1025

1026def _dup_decompose(f, K):

1027 """Helper function for :func:`dup_decompose`."""

1028 df = len(f) - 1

1029

1030 for s in range(2, df):

1031 if df % s != 0:

1032 continue

1033

1034 h = _dup_right_decompose(f, s, K)

1035

1036 if h is not None:

1037 g = _dup_left_decompose(f, h, K)

1038

1039 if g is not None:

1040 return g, h

1041

1042 return None

1043

1044

1045def dup_decompose(f, K):

1046 """

1047 Computes functional decomposition of ``f`` in ``K[x]``.

1048

1049 Given a univariate polynomial ``f`` with coefficients in a field of

1050 characteristic zero, returns list ``[f_1, f_2, ..., f_n]``, where::

1051

1052 f = f_1 o f_2 o ... f_n = f_1(f_2(... f_n))

1053

1054 and ``f_2, ..., f_n`` are monic and homogeneous polynomials of at

1055 least second degree.

1056

1057 Unlike factorization, complete functional decompositions of

1058 polynomials are not unique, consider examples:

1059

1060 1. ``f o g = f(x + b) o (g - b)``

1061 2. ``x**n o x**m = x**m o x**n``

1062 3. ``T_n o T_m = T_m o T_n``

1063

1064 where ``T_n`` and ``T_m`` are Chebyshev polynomials.

1065

1066 Examples

1067 ========

1068

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

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

1071

1072 >>> R.dup_decompose(x**4 - 2*x**3 + x**2)

1073 [x**2, x**2 - x]

1074

1075 References

1076 ==========

1077

1078 .. [1] [Kozen89]_

1079

1080 """

1081 F = []

1082

1083 while True:

1084 result = _dup_decompose(f, K)

1085

1086 if result is not None:

1087 f, h = result

1088 F = [h] + F

1089 else:

1090 break

1091

1092 return [f] + F

1093

1094

1095def dmp_lift(f, u, K):

1096 """

1097 Convert algebraic coefficients to integers in ``K[X]``.

1098

1099 Examples

1100 ========

1101

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

1103 >>> from sympy import I

1104

1105 >>> K = QQ.algebraic_field(I)

1106 >>> R, x = ring("x", K)

1107

1108 >>> f = x**2 + K([QQ(1), QQ(0)])*x + K([QQ(2), QQ(0)])

1109

1110 >>> R.dmp_lift(f)

1111 x**8 + 2*x**6 + 9*x**4 - 8*x**2 + 16

1112

1113 """

1114 if K.is_GaussianField:

1115 K1 = K.as_AlgebraicField()

1116 f = dmp_convert(f, u, K, K1)

1117 K = K1

1118

1119 if not K.is_Algebraic:

1120 raise DomainError(

1121 'computation can be done only in an algebraic domain')

1122

1123 F, monoms, polys = dmp_to_dict(f, u), [], []

1124

1125 for monom, coeff in F.items():

1126 if not coeff.is_ground:

1127 monoms.append(monom)

1128

1129 perms = variations([-1, 1], len(monoms), repetition=True)

1130

1131 for perm in perms:

1132 G = dict(F)

1133

1134 for sign, monom in zip(perm, monoms):

1135 if sign == -1:

1136 G[monom] = -G[monom]

1137

1138 polys.append(dmp_from_dict(G, u, K))

1139

1140 return dmp_convert(dmp_expand(polys, u, K), u, K, K.dom)

1141

1142

1143def dup_sign_variations(f, K):

1144 """

1145 Compute the number of sign variations of ``f`` in ``K[x]``.

1146

1147 Examples

1148 ========

1149

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

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

1152

1153 >>> R.dup_sign_variations(x**4 - x**2 - x + 1)

1154 2

1155

1156 """

1157 prev, k = K.zero, 0

1158

1159 for coeff in f:

1160 if K.is_negative(coeff*prev):

1161 k += 1

1162

1163 if coeff:

1164 prev = coeff

1165

1166 return k

1167

1168

1169def dup_clear_denoms(f, K0, K1=None, convert=False):

1170 """

1171 Clear denominators, i.e. transform ``K_0`` to ``K_1``.

1172

1173 Examples

1174 ========

1175

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

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

1178

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

1180

1181 >>> R.dup_clear_denoms(f, convert=False)

1182 (6, 3*x + 2)

1183 >>> R.dup_clear_denoms(f, convert=True)

1184 (6, 3*x + 2)

1185

1186 """

1187 if K1 is None:

1188 if K0.has_assoc_Ring:

1189 K1 = K0.get_ring()

1190 else:

1191 K1 = K0

1192

1193 common = K1.one

1194

1195 for c in f:

1196 common = K1.lcm(common, K0.denom(c))

1197

1198 if not K1.is_one(common):

1199 f = dup_mul_ground(f, common, K0)

1200

1201 if not convert:

1202 return common, f

1203 else:

1204 return common, dup_convert(f, K0, K1)

1205

1206

1207def _rec_clear_denoms(g, v, K0, K1):

1208 """Recursive helper for :func:`dmp_clear_denoms`."""

1209 common = K1.one

1210

1211 if not v:

1212 for c in g:

1213 common = K1.lcm(common, K0.denom(c))

1214 else:

1215 w = v - 1

1216

1217 for c in g:

1218 common = K1.lcm(common, _rec_clear_denoms(c, w, K0, K1))

1219

1220 return common

1221

1222

1223def dmp_clear_denoms(f, u, K0, K1=None, convert=False):

1224 """

1225 Clear denominators, i.e. transform ``K_0`` to ``K_1``.

1226

1227 Examples

1228 ========

1229

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

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

1232

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

1234

1235 >>> R.dmp_clear_denoms(f, convert=False)

1236 (6, 3*x + 2*y + 6)

1237 >>> R.dmp_clear_denoms(f, convert=True)

1238 (6, 3*x + 2*y + 6)

1239

1240 """

1241 if not u:

1242 return dup_clear_denoms(f, K0, K1, convert=convert)

1243

1244 if K1 is None:

1245 if K0.has_assoc_Ring:

1246 K1 = K0.get_ring()

1247 else:

1248 K1 = K0

1249

1250 common = _rec_clear_denoms(f, u, K0, K1)

1251

1252 if not K1.is_one(common):

1253 f = dmp_mul_ground(f, common, u, K0)

1254

1255 if not convert:

1256 return common, f

1257 else:

1258 return common, dmp_convert(f, u, K0, K1)

1259

1260

1261def dup_revert(f, n, K):

1262 """

1263 Compute ``f**(-1)`` mod ``x**n`` using Newton iteration.

1264

1265 This function computes first ``2**n`` terms of a polynomial that

1266 is a result of inversion of a polynomial modulo ``x**n``. This is

1267 useful to efficiently compute series expansion of ``1/f``.

1268

1269 Examples

1270 ========

1271

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

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

1274

1275 >>> f = -QQ(1,720)*x**6 + QQ(1,24)*x**4 - QQ(1,2)*x**2 + 1

1276

1277 >>> R.dup_revert(f, 8)

1278 61/720*x**6 + 5/24*x**4 + 1/2*x**2 + 1

1279

1280 """

1281 g = [K.revert(dup_TC(f, K))]

1282 h = [K.one, K.zero, K.zero]

1283

1284 N = int(_ceil(_log(n, 2)))

1285

1286 for i in range(1, N + 1):

1287 a = dup_mul_ground(g, K(2), K)

1288 b = dup_mul(f, dup_sqr(g, K), K)

1289 g = dup_rem(dup_sub(a, b, K), h, K)

1290 h = dup_lshift(h, dup_degree(h), K)

1291

1292 return g

1293

1294

1295def dmp_revert(f, g, u, K):

1296 """

1297 Compute ``f**(-1)`` mod ``x**n`` using Newton iteration.

1298

1299 Examples

1300 ========

1301

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

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

1304

1305 """

1306 if not u:

1307 return dup_revert(f, g, K)

1308 else:

1309 raise MultivariatePolynomialError(f, g)