Coverage for /usr/lib/python3/dist-packages/sympy/polys/monomials.py: 41%

1"""Tools and arithmetics for monomials of distributed polynomials. """

4from itertools import combinations_with_replacement, product

5from textwrap import dedent

7from sympy.core import Mul, S, Tuple, sympify

8from sympy.polys.polyerrors import ExactQuotientFailed

9from sympy.polys.polyutils import PicklableWithSlots, dict_from_expr

10from sympy.utilities import public

11from sympy.utilities.iterables import is_sequence, iterable

13@public

14def itermonomials(variables, max_degrees, min_degrees=None):

15 r"""

16 ``max_degrees`` and ``min_degrees`` are either both integers or both lists.

17 Unless otherwise specified, ``min_degrees`` is either ``0`` or

18 ``[0, ..., 0]``.

20 A generator of all monomials ``monom`` is returned, such that

21 either

22 ``min_degree <= total_degree(monom) <= max_degree``,

23 or

24 ``min_degrees[i] <= degree_list(monom)[i] <= max_degrees[i]``,

25 for all ``i``.

27 Case I. ``max_degrees`` and ``min_degrees`` are both integers

28 =============================================================

30 Given a set of variables $V$ and a min_degree $N$ and a max_degree $M$

31 generate a set of monomials of degree less than or equal to $N$ and greater

32 than or equal to $M$. The total number of monomials in commutative

33 variables is huge and is given by the following formula if $M = 0$:

35 .. math::

36 \frac{(\#V + N)!}{\#V! N!}

38 For example if we would like to generate a dense polynomial of

39 a total degree $N = 50$ and $M = 0$, which is the worst case, in 5

40 variables, assuming that exponents and all of coefficients are 32-bit long

41 and stored in an array we would need almost 80 GiB of memory! Fortunately

42 most polynomials, that we will encounter, are sparse.

44 Consider monomials in commutative variables $x$ and $y$

45 and non-commutative variables $a$ and $b$::

47 >>> from sympy import symbols

48 >>> from sympy.polys.monomials import itermonomials

49 >>> from sympy.polys.orderings import monomial_key

50 >>> from sympy.abc import x, y

52 >>> sorted(itermonomials([x, y], 2), key=monomial_key('grlex', [y, x]))

53 [1, x, y, x**2, x*y, y**2]

55 >>> sorted(itermonomials([x, y], 3), key=monomial_key('grlex', [y, x]))

56 [1, x, y, x**2, x*y, y**2, x**3, x**2*y, x*y**2, y**3]

58 >>> a, b = symbols('a, b', commutative=False)

59 >>> set(itermonomials([a, b, x], 2))

60 {1, a, a**2, b, b**2, x, x**2, a*b, b*a, x*a, x*b}

62 >>> sorted(itermonomials([x, y], 2, 1), key=monomial_key('grlex', [y, x]))

63 [x, y, x**2, x*y, y**2]

65 Case II. ``max_degrees`` and ``min_degrees`` are both lists

66 ===========================================================

68 If ``max_degrees = [d_1, ..., d_n]`` and

69 ``min_degrees = [e_1, ..., e_n]``, the number of monomials generated

70 is:

72 .. math::

73 (d_1 - e_1 + 1) (d_2 - e_2 + 1) \cdots (d_n - e_n + 1)

75 Let us generate all monomials ``monom`` in variables $x$ and $y$

76 such that ``[1, 2][i] <= degree_list(monom)[i] <= [2, 4][i]``,

77 ``i = 0, 1`` ::

79 >>> from sympy import symbols

80 >>> from sympy.polys.monomials import itermonomials

81 >>> from sympy.polys.orderings import monomial_key

82 >>> from sympy.abc import x, y

84 >>> sorted(itermonomials([x, y], [2, 4], [1, 2]), reverse=True, key=monomial_key('lex', [x, y]))

85 [x**2*y**4, x**2*y**3, x**2*y**2, x*y**4, x*y**3, x*y**2]

86 """

87 n = len(variables)

88 if is_sequence(max_degrees):

89 if len(max_degrees) != n:

90 raise ValueError('Argument sizes do not match')

91 if min_degrees is None:

92 min_degrees = [0]*n

93 elif not is_sequence(min_degrees):

94 raise ValueError('min_degrees is not a list')

95 else:

96 if len(min_degrees) != n:

97 raise ValueError('Argument sizes do not match')

98 if any(i < 0 for i in min_degrees):

99 raise ValueError("min_degrees cannot contain negative numbers")

100 total_degree = False

101 else:

102 max_degree = max_degrees

103 if max_degree < 0:

104 raise ValueError("max_degrees cannot be negative")

105 if min_degrees is None:

106 min_degree = 0

107 else:

108 if min_degrees < 0:

109 raise ValueError("min_degrees cannot be negative")

110 min_degree = min_degrees

111 total_degree = True

112 if total_degree:

113 if min_degree > max_degree:

114 return

115 if not variables or max_degree == 0:

116 yield S.One

117 return

118 # Force to list in case of passed tuple or other incompatible collection

119 variables = list(variables) + [S.One]

120 if all(variable.is_commutative for variable in variables):

121 monomials_list_comm = []

122 for item in combinations_with_replacement(variables, max_degree):

123 powers = {variable: 0 for variable in variables}

124 for variable in item:

125 if variable != 1:

126 powers[variable] += 1

127 if sum(powers.values()) >= min_degree:

128 monomials_list_comm.append(Mul(*item))

129 yield from set(monomials_list_comm)

130 else:

131 monomials_list_non_comm = []

132 for item in product(variables, repeat=max_degree):

133 powers = {variable: 0 for variable in variables}

134 for variable in item:

135 if variable != 1:

136 powers[variable] += 1

137 if sum(powers.values()) >= min_degree:

138 monomials_list_non_comm.append(Mul(*item))

139 yield from set(monomials_list_non_comm)

140 else:

141 if any(min_degrees[i] > max_degrees[i] for i in range(n)):

142 raise ValueError('min_degrees[i] must be <= max_degrees[i] for all i')

143 power_lists = []

144 for var, min_d, max_d in zip(variables, min_degrees, max_degrees):

145 power_lists.append([var**i for i in range(min_d, max_d + 1)])

146 for powers in product(*power_lists):

147 yield Mul(*powers)

148

149def monomial_count(V, N):

150 r"""

151 Computes the number of monomials.

152

153 The number of monomials is given by the following formula:

154

155 .. math::

156

157 \frac{(\#V + N)!}{\#V! N!}

158

159 where `N` is a total degree and `V` is a set of variables.

160

161 Examples

162 ========

163

164 >>> from sympy.polys.monomials import itermonomials, monomial_count

165 >>> from sympy.polys.orderings import monomial_key

166 >>> from sympy.abc import x, y

167

168 >>> monomial_count(2, 2)

169 6

170

171 >>> M = list(itermonomials([x, y], 2))

172

173 >>> sorted(M, key=monomial_key('grlex', [y, x]))

174 [1, x, y, x**2, x*y, y**2]

175 >>> len(M)

176 6

177

178 """

179 from sympy.functions.combinatorial.factorials import factorial

180 return factorial(V + N) / factorial(V) / factorial(N)

181

182def monomial_mul(A, B):

183 """

184 Multiplication of tuples representing monomials.

185

186 Examples

187 ========

188

189 Lets multiply `x**3*y**4*z` with `x*y**2`::

190

191 >>> from sympy.polys.monomials import monomial_mul

192

193 >>> monomial_mul((3, 4, 1), (1, 2, 0))

194 (4, 6, 1)

195

196 which gives `x**4*y**5*z`.

197

198 """

199 return tuple([ a + b for a, b in zip(A, B) ])

200

201def monomial_div(A, B):

202 """

203 Division of tuples representing monomials.

204

205 Examples

206 ========

207

208 Lets divide `x**3*y**4*z` by `x*y**2`::

209

210 >>> from sympy.polys.monomials import monomial_div

211

212 >>> monomial_div((3, 4, 1), (1, 2, 0))

213 (2, 2, 1)

214

215 which gives `x**2*y**2*z`. However::

216

217 >>> monomial_div((3, 4, 1), (1, 2, 2)) is None

218 True

219

220 `x*y**2*z**2` does not divide `x**3*y**4*z`.

221

222 """

223 C = monomial_ldiv(A, B)

224

225 if all(c >= 0 for c in C):

226 return tuple(C)

227 else:

228 return None

229

230def monomial_ldiv(A, B):

231 """

232 Division of tuples representing monomials.

233

234 Examples

235 ========

236

237 Lets divide `x**3*y**4*z` by `x*y**2`::

238

239 >>> from sympy.polys.monomials import monomial_ldiv

240

241 >>> monomial_ldiv((3, 4, 1), (1, 2, 0))

242 (2, 2, 1)

243

244 which gives `x**2*y**2*z`.

245

246 >>> monomial_ldiv((3, 4, 1), (1, 2, 2))

247 (2, 2, -1)

248

249 which gives `x**2*y**2*z**-1`.

250

251 """

252 return tuple([ a - b for a, b in zip(A, B) ])

253

254def monomial_pow(A, n):

255 """Return the n-th pow of the monomial. """

256 return tuple([ a*n for a in A ])

257

258def monomial_gcd(A, B):

259 """

260 Greatest common divisor of tuples representing monomials.

261

262 Examples

263 ========

264

265 Lets compute GCD of `x*y**4*z` and `x**3*y**2`::

266

267 >>> from sympy.polys.monomials import monomial_gcd

268

269 >>> monomial_gcd((1, 4, 1), (3, 2, 0))

270 (1, 2, 0)

271

272 which gives `x*y**2`.

273

274 """

275 return tuple([ min(a, b) for a, b in zip(A, B) ])

276

277def monomial_lcm(A, B):

278 """

279 Least common multiple of tuples representing monomials.

280

281 Examples

282 ========

283

284 Lets compute LCM of `x*y**4*z` and `x**3*y**2`::

285

286 >>> from sympy.polys.monomials import monomial_lcm

287

288 >>> monomial_lcm((1, 4, 1), (3, 2, 0))

289 (3, 4, 1)

290

291 which gives `x**3*y**4*z`.

292

293 """

294 return tuple([ max(a, b) for a, b in zip(A, B) ])

295

296def monomial_divides(A, B):

297 """

298 Does there exist a monomial X such that XA == B?

299

300 Examples

301 ========

302

303 >>> from sympy.polys.monomials import monomial_divides

304 >>> monomial_divides((1, 2), (3, 4))

305 True

306 >>> monomial_divides((1, 2), (0, 2))

307 False

308 """

309 return all(a <= b for a, b in zip(A, B))

310

311def monomial_max(*monoms):

312 """

313 Returns maximal degree for each variable in a set of monomials.

314

315 Examples

316 ========

317

318 Consider monomials `x**3*y**4*z**5`, `y**5*z` and `x**6*y**3*z**9`.

319 We wish to find out what is the maximal degree for each of `x`, `y`

320 and `z` variables::

321

322 >>> from sympy.polys.monomials import monomial_max

323

324 >>> monomial_max((3,4,5), (0,5,1), (6,3,9))

325 (6, 5, 9)

326

327 """

328 M = list(monoms[0])

329

330 for N in monoms[1:]:

331 for i, n in enumerate(N):

332 M[i] = max(M[i], n)

333

334 return tuple(M)

335

336def monomial_min(*monoms):

337 """

338 Returns minimal degree for each variable in a set of monomials.

339

340 Examples

341 ========

342

343 Consider monomials `x**3*y**4*z**5`, `y**5*z` and `x**6*y**3*z**9`.

344 We wish to find out what is the minimal degree for each of `x`, `y`

345 and `z` variables::

346

347 >>> from sympy.polys.monomials import monomial_min

348

349 >>> monomial_min((3,4,5), (0,5,1), (6,3,9))

350 (0, 3, 1)

351

352 """

353 M = list(monoms[0])

354

355 for N in monoms[1:]:

356 for i, n in enumerate(N):

357 M[i] = min(M[i], n)

358

359 return tuple(M)

360

361def monomial_deg(M):

362 """

363 Returns the total degree of a monomial.

364

365 Examples

366 ========

367

368 The total degree of `xy^2` is 3:

369

370 >>> from sympy.polys.monomials import monomial_deg

371 >>> monomial_deg((1, 2))

372 3

373 """

374 return sum(M)

375

376def term_div(a, b, domain):

377 """Division of two terms in over a ring/field. """

378 a_lm, a_lc = a

379 b_lm, b_lc = b

380

381 monom = monomial_div(a_lm, b_lm)

382

383 if domain.is_Field:

384 if monom is not None:

385 return monom, domain.quo(a_lc, b_lc)

386 else:

387 return None

388 else:

389 if not (monom is None or a_lc % b_lc):

390 return monom, domain.quo(a_lc, b_lc)

391 else:

392 return None

393

394class MonomialOps:

395 """Code generator of fast monomial arithmetic functions. """

396

397 def __init__(self, ngens):

398 self.ngens = ngens

399

400 def _build(self, code, name):

401 ns = {}

402 exec(code, ns)

403 return ns[name]

404

405 def _vars(self, name):

406 return [ "%s%s" % (name, i) for i in range(self.ngens) ]

407

408 def mul(self):

409 name = "monomial_mul"

410 template = dedent("""\

411 def %(name)s(A, B):

412 (%(A)s,) = A

413 (%(B)s,) = B

414 return (%(AB)s,)

415 """)

416 A = self._vars("a")

417 B = self._vars("b")

418 AB = [ "%s + %s" % (a, b) for a, b in zip(A, B) ]

419 code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}

420 return self._build(code, name)

421

422 def pow(self):

423 name = "monomial_pow"

424 template = dedent("""\

425 def %(name)s(A, k):

426 (%(A)s,) = A

427 return (%(Ak)s,)

428 """)

429 A = self._vars("a")

430 Ak = [ "%s*k" % a for a in A ]

431 code = template % {"name": name, "A": ", ".join(A), "Ak": ", ".join(Ak)}

432 return self._build(code, name)

433

434 def mulpow(self):

435 name = "monomial_mulpow"

436 template = dedent("""\

437 def %(name)s(A, B, k):

438 (%(A)s,) = A

439 (%(B)s,) = B

440 return (%(ABk)s,)

441 """)

442 A = self._vars("a")

443 B = self._vars("b")

444 ABk = [ "%s + %s*k" % (a, b) for a, b in zip(A, B) ]

445 code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "ABk": ", ".join(ABk)}

446 return self._build(code, name)

447

448 def ldiv(self):

449 name = "monomial_ldiv"

450 template = dedent("""\

451 def %(name)s(A, B):

452 (%(A)s,) = A

453 (%(B)s,) = B

454 return (%(AB)s,)

455 """)

456 A = self._vars("a")

457 B = self._vars("b")

458 AB = [ "%s - %s" % (a, b) for a, b in zip(A, B) ]

459 code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}

460 return self._build(code, name)

461

462 def div(self):

463 name = "monomial_div"

464 template = dedent("""\

465 def %(name)s(A, B):

466 (%(A)s,) = A

467 (%(B)s,) = B

468 %(RAB)s

469 return (%(R)s,)

470 """)

471 A = self._vars("a")

472 B = self._vars("b")

473 RAB = [ "r%(i)s = a%(i)s - b%(i)s\n if r%(i)s < 0: return None" % {"i": i} for i in range(self.ngens) ]

474 R = self._vars("r")

475 code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "RAB": "\n ".join(RAB), "R": ", ".join(R)}

476 return self._build(code, name)

477

478 def lcm(self):

479 name = "monomial_lcm"

480 template = dedent("""\

481 def %(name)s(A, B):

482 (%(A)s,) = A

483 (%(B)s,) = B

484 return (%(AB)s,)

485 """)

486 A = self._vars("a")

487 B = self._vars("b")

488 AB = [ "%s if %s >= %s else %s" % (a, a, b, b) for a, b in zip(A, B) ]

489 code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}

490 return self._build(code, name)

491

492 def gcd(self):

493 name = "monomial_gcd"

494 template = dedent("""\

495 def %(name)s(A, B):

496 (%(A)s,) = A

497 (%(B)s,) = B

498 return (%(AB)s,)

499 """)

500 A = self._vars("a")

501 B = self._vars("b")

502 AB = [ "%s if %s <= %s else %s" % (a, a, b, b) for a, b in zip(A, B) ]

503 code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}

504 return self._build(code, name)

505

506@public

507class Monomial(PicklableWithSlots):

508 """Class representing a monomial, i.e. a product of powers. """

509

510 __slots__ = ('exponents', 'gens')

511

512 def __init__(self, monom, gens=None):

513 if not iterable(monom):

514 rep, gens = dict_from_expr(sympify(monom), gens=gens)

515 if len(rep) == 1 and list(rep.values())[0] == 1:

516 monom = list(rep.keys())[0]

517 else:

518 raise ValueError("Expected a monomial got {}".format(monom))

519

520 self.exponents = tuple(map(int, monom))

521 self.gens = gens

522

523 def rebuild(self, exponents, gens=None):

524 return self.__class__(exponents, gens or self.gens)

525

526 def __len__(self):

527 return len(self.exponents)

528

529 def __iter__(self):

530 return iter(self.exponents)

531

532 def __getitem__(self, item):

533 return self.exponents[item]

534

535 def __hash__(self):

536 return hash((self.__class__.__name__, self.exponents, self.gens))

537

538 def __str__(self):

539 if self.gens:

540 return "*".join([ "%s**%s" % (gen, exp) for gen, exp in zip(self.gens, self.exponents) ])

541 else:

542 return "%s(%s)" % (self.__class__.__name__, self.exponents)

543

544 def as_expr(self, *gens):

545 """Convert a monomial instance to a SymPy expression. """

546 gens = gens or self.gens

547

548 if not gens:

549 raise ValueError(

550 "Cannot convert %s to an expression without generators" % self)

551

552 return Mul(*[ gen**exp for gen, exp in zip(gens, self.exponents) ])

553

554 def __eq__(self, other):

555 if isinstance(other, Monomial):

556 exponents = other.exponents

557 elif isinstance(other, (tuple, Tuple)):

558 exponents = other

559 else:

560 return False

561

562 return self.exponents == exponents

563

564 def __ne__(self, other):

565 return not self == other

566

567 def __mul__(self, other):

568 if isinstance(other, Monomial):

569 exponents = other.exponents

570 elif isinstance(other, (tuple, Tuple)):

571 exponents = other

572 else:

573 raise NotImplementedError

574

575 return self.rebuild(monomial_mul(self.exponents, exponents))

576

577 def __truediv__(self, other):

578 if isinstance(other, Monomial):

579 exponents = other.exponents

580 elif isinstance(other, (tuple, Tuple)):

581 exponents = other

582 else:

583 raise NotImplementedError

584

585 result = monomial_div(self.exponents, exponents)

586

587 if result is not None:

588 return self.rebuild(result)

589 else:

590 raise ExactQuotientFailed(self, Monomial(other))

591

592 __floordiv__ = __truediv__

593

594 def __pow__(self, other):

595 n = int(other)

596

597 if not n:

598 return self.rebuild([0]*len(self))

599 elif n > 0:

600 exponents = self.exponents

601

602 for i in range(1, n):

603 exponents = monomial_mul(exponents, self.exponents)

604

605 return self.rebuild(exponents)

606 else:

607 raise ValueError("a non-negative integer expected, got %s" % other)

608

609 def gcd(self, other):

610 """Greatest common divisor of monomials. """

611 if isinstance(other, Monomial):

612 exponents = other.exponents

613 elif isinstance(other, (tuple, Tuple)):

614 exponents = other

615 else:

616 raise TypeError(

617 "an instance of Monomial class expected, got %s" % other)

618

619 return self.rebuild(monomial_gcd(self.exponents, exponents))

620

621 def lcm(self, other):

622 """Least common multiple of monomials. """

623 if isinstance(other, Monomial):

624 exponents = other.exponents

625 elif isinstance(other, (tuple, Tuple)):

626 exponents = other

627 else:

628 raise TypeError(

629 "an instance of Monomial class expected, got %s" % other)

630

631 return self.rebuild(monomial_lcm(self.exponents, exponents))