Coverage for /usr/lib/python3/dist-packages/sympy/polys/polyutils.py: 37%

1"""Useful utilities for higher level polynomial classes. """

4from sympy.core import (S, Add, Mul, Pow, Eq, Expr,

5 expand_mul, expand_multinomial)

6from sympy.core.exprtools import decompose_power, decompose_power_rat

7from sympy.core.numbers import _illegal

8from sympy.polys.polyerrors import PolynomialError, GeneratorsError

9from sympy.polys.polyoptions import build_options

12import re

14_gens_order = {

15 'a': 301, 'b': 302, 'c': 303, 'd': 304,

16 'e': 305, 'f': 306, 'g': 307, 'h': 308,

17 'i': 309, 'j': 310, 'k': 311, 'l': 312,

18 'm': 313, 'n': 314, 'o': 315, 'p': 216,

19 'q': 217, 'r': 218, 's': 219, 't': 220,

20 'u': 221, 'v': 222, 'w': 223, 'x': 124,

21 'y': 125, 'z': 126,

22}

24_max_order = 1000

25_re_gen = re.compile(r"^(.*?)(\d*)$", re.MULTILINE)

28def _nsort(roots, separated=False):

29 """Sort the numerical roots putting the real roots first, then sorting

30 according to real and imaginary parts. If ``separated`` is True, then

31 the real and imaginary roots will be returned in two lists, respectively.

33 This routine tries to avoid issue 6137 by separating the roots into real

34 and imaginary parts before evaluation. In addition, the sorting will raise

35 an error if any computation cannot be done with precision.

36 """

37 if not all(r.is_number for r in roots):

38 raise NotImplementedError

39 # see issue 6137:

40 # get the real part of the evaluated real and imaginary parts of each root

41 key = [[i.n(2).as_real_imag()[0] for i in r.as_real_imag()] for r in roots]

42 # make sure the parts were computed with precision

43 if len(roots) > 1 and any(i._prec == 1 for k in key for i in k):

44 raise NotImplementedError("could not compute root with precision")

45 # insert a key to indicate if the root has an imaginary part

46 key = [(1 if i else 0, r, i) for r, i in key]

47 key = sorted(zip(key, roots))

48 # return the real and imaginary roots separately if desired

49 if separated:

50 r = []

51 i = []

52 for (im, _, _), v in key:

53 if im:

54 i.append(v)

55 else:

56 r.append(v)

57 return r, i

58 _, roots = zip(*key)

59 return list(roots)

62def _sort_gens(gens, **args):

63 """Sort generators in a reasonably intelligent way. """

64 opt = build_options(args)

66 gens_order, wrt = {}, None

68 if opt is not None:

69 gens_order, wrt = {}, opt.wrt

71 for i, gen in enumerate(opt.sort):

72 gens_order[gen] = i + 1

74 def order_key(gen):

75 gen = str(gen)

77 if wrt is not None:

78 try:

79 return (-len(wrt) + wrt.index(gen), gen, 0)

80 except ValueError:

81 pass

83 name, index = _re_gen.match(gen).groups()

85 if index:

86 index = int(index)

87 else:

88 index = 0

90 try:

91 return ( gens_order[name], name, index)

92 except KeyError:

93 pass

95 try:

96 return (_gens_order[name], name, index)

97 except KeyError:

98 pass

100 return (_max_order, name, index)

101

102 try:

103 gens = sorted(gens, key=order_key)

104 except TypeError: # pragma: no cover

105 pass

106

107 return tuple(gens)

108

109

110def _unify_gens(f_gens, g_gens):

111 """Unify generators in a reasonably intelligent way. """

112 f_gens = list(f_gens)

113 g_gens = list(g_gens)

114

115 if f_gens == g_gens:

116 return tuple(f_gens)

117

118 gens, common, k = [], [], 0

119

120 for gen in f_gens:

121 if gen in g_gens:

122 common.append(gen)

123

124 for i, gen in enumerate(g_gens):

125 if gen in common:

126 g_gens[i], k = common[k], k + 1

127

128 for gen in common:

129 i = f_gens.index(gen)

130

131 gens.extend(f_gens[:i])

132 f_gens = f_gens[i + 1:]

133

134 i = g_gens.index(gen)

135

136 gens.extend(g_gens[:i])

137 g_gens = g_gens[i + 1:]

138

139 gens.append(gen)

140

141 gens.extend(f_gens)

142 gens.extend(g_gens)

143

144 return tuple(gens)

145

146

147def _analyze_gens(gens):

148 """Support for passing generators as `*gens` and `[gens]`. """

149 if len(gens) == 1 and hasattr(gens[0], '__iter__'):

150 return tuple(gens[0])

151 else:

152 return tuple(gens)

153

154

155def _sort_factors(factors, **args):

156 """Sort low-level factors in increasing 'complexity' order. """

157 def order_if_multiple_key(factor):

158 (f, n) = factor

159 return (len(f), n, f)

160

161 def order_no_multiple_key(f):

162 return (len(f), f)

163

164 if args.get('multiple', True):

165 return sorted(factors, key=order_if_multiple_key)

166 else:

167 return sorted(factors, key=order_no_multiple_key)

168

169illegal_types = [type(obj) for obj in _illegal]

170finf = [float(i) for i in _illegal[1:3]]

171def _not_a_coeff(expr):

172 """Do not treat NaN and infinities as valid polynomial coefficients. """

173 if type(expr) in illegal_types or expr in finf:

174 return True

175 if isinstance(expr, float) and float(expr) != expr:

176 return True # nan

177 return # could be

178

179

180def _parallel_dict_from_expr_if_gens(exprs, opt):

181 """Transform expressions into a multinomial form given generators. """

182 k, indices = len(opt.gens), {}

183

184 for i, g in enumerate(opt.gens):

185 indices[g] = i

186

187 polys = []

188

189 for expr in exprs:

190 poly = {}

191

192 if expr.is_Equality:

193 expr = expr.lhs - expr.rhs

194

195 for term in Add.make_args(expr):

196 coeff, monom = [], [0]*k

197

198 for factor in Mul.make_args(term):

199 if not _not_a_coeff(factor) and factor.is_Number:

200 coeff.append(factor)

201 else:

202 try:

203 if opt.series is False:

204 base, exp = decompose_power(factor)

205

206 if exp < 0:

207 exp, base = -exp, Pow(base, -S.One)

208 else:

209 base, exp = decompose_power_rat(factor)

210

211 monom[indices[base]] = exp

212 except KeyError:

213 if not factor.has_free(*opt.gens):

214 coeff.append(factor)

215 else:

216 raise PolynomialError("%s contains an element of "

217 "the set of generators." % factor)

218

219 monom = tuple(monom)

220

221 if monom in poly:

222 poly[monom] += Mul(*coeff)

223 else:

224 poly[monom] = Mul(*coeff)

225

226 polys.append(poly)

227

228 return polys, opt.gens

229

230

231def _parallel_dict_from_expr_no_gens(exprs, opt):

232 """Transform expressions into a multinomial form and figure out generators. """

233 if opt.domain is not None:

234 def _is_coeff(factor):

235 return factor in opt.domain

236 elif opt.extension is True:

237 def _is_coeff(factor):

238 return factor.is_algebraic

239 elif opt.greedy is not False:

240 def _is_coeff(factor):

241 return factor is S.ImaginaryUnit

242 else:

243 def _is_coeff(factor):

244 return factor.is_number

245

246 gens, reprs = set(), []

247

248 for expr in exprs:

249 terms = []

250

251 if expr.is_Equality:

252 expr = expr.lhs - expr.rhs

253

254 for term in Add.make_args(expr):

255 coeff, elements = [], {}

256

257 for factor in Mul.make_args(term):

258 if not _not_a_coeff(factor) and (factor.is_Number or _is_coeff(factor)):

259 coeff.append(factor)

260 else:

261 if opt.series is False:

262 base, exp = decompose_power(factor)

263

264 if exp < 0:

265 exp, base = -exp, Pow(base, -S.One)

266 else:

267 base, exp = decompose_power_rat(factor)

268

269 elements[base] = elements.setdefault(base, 0) + exp

270 gens.add(base)

271

272 terms.append((coeff, elements))

273

274 reprs.append(terms)

275

276 gens = _sort_gens(gens, opt=opt)

277 k, indices = len(gens), {}

278

279 for i, g in enumerate(gens):

280 indices[g] = i

281

282 polys = []

283

284 for terms in reprs:

285 poly = {}

286

287 for coeff, term in terms:

288 monom = [0]*k

289

290 for base, exp in term.items():

291 monom[indices[base]] = exp

292

293 monom = tuple(monom)

294

295 if monom in poly:

296 poly[monom] += Mul(*coeff)

297 else:

298 poly[monom] = Mul(*coeff)

299

300 polys.append(poly)

301

302 return polys, tuple(gens)

303

304

305def _dict_from_expr_if_gens(expr, opt):

306 """Transform an expression into a multinomial form given generators. """

307 (poly,), gens = _parallel_dict_from_expr_if_gens((expr,), opt)

308 return poly, gens

309

310

311def _dict_from_expr_no_gens(expr, opt):

312 """Transform an expression into a multinomial form and figure out generators. """

313 (poly,), gens = _parallel_dict_from_expr_no_gens((expr,), opt)

314 return poly, gens

315

316

317def parallel_dict_from_expr(exprs, **args):

318 """Transform expressions into a multinomial form. """

319 reps, opt = _parallel_dict_from_expr(exprs, build_options(args))

320 return reps, opt.gens

321

322

323def _parallel_dict_from_expr(exprs, opt):

324 """Transform expressions into a multinomial form. """

325 if opt.expand is not False:

326 exprs = [ expr.expand() for expr in exprs ]

327

328 if any(expr.is_commutative is False for expr in exprs):

329 raise PolynomialError('non-commutative expressions are not supported')

330

331 if opt.gens:

332 reps, gens = _parallel_dict_from_expr_if_gens(exprs, opt)

333 else:

334 reps, gens = _parallel_dict_from_expr_no_gens(exprs, opt)

335

336 return reps, opt.clone({'gens': gens})

337

338

339def dict_from_expr(expr, **args):

340 """Transform an expression into a multinomial form. """

341 rep, opt = _dict_from_expr(expr, build_options(args))

342 return rep, opt.gens

343

344

345def _dict_from_expr(expr, opt):

346 """Transform an expression into a multinomial form. """

347 if expr.is_commutative is False:

348 raise PolynomialError('non-commutative expressions are not supported')

349

350 def _is_expandable_pow(expr):

351 return (expr.is_Pow and expr.exp.is_positive and expr.exp.is_Integer

352 and expr.base.is_Add)

353

354 if opt.expand is not False:

355 if not isinstance(expr, (Expr, Eq)):

356 raise PolynomialError('expression must be of type Expr')

357 expr = expr.expand()

358 # TODO: Integrate this into expand() itself

359 while any(_is_expandable_pow(i) or i.is_Mul and

360 any(_is_expandable_pow(j) for j in i.args) for i in

361 Add.make_args(expr)):

362

363 expr = expand_multinomial(expr)

364 while any(i.is_Mul and any(j.is_Add for j in i.args) for i in Add.make_args(expr)):

365 expr = expand_mul(expr)

366

367 if opt.gens:

368 rep, gens = _dict_from_expr_if_gens(expr, opt)

369 else:

370 rep, gens = _dict_from_expr_no_gens(expr, opt)

371

372 return rep, opt.clone({'gens': gens})

373

374

375def expr_from_dict(rep, *gens):

376 """Convert a multinomial form into an expression. """

377 result = []

378

379 for monom, coeff in rep.items():

380 term = [coeff]

381 for g, m in zip(gens, monom):

382 if m:

383 term.append(Pow(g, m))

384

385 result.append(Mul(*term))

386

387 return Add(*result)

388

389parallel_dict_from_basic = parallel_dict_from_expr

390dict_from_basic = dict_from_expr

391basic_from_dict = expr_from_dict

392

393

394def _dict_reorder(rep, gens, new_gens):

395 """Reorder levels using dict representation. """

396 gens = list(gens)

397

398 monoms = rep.keys()

399 coeffs = rep.values()

400

401 new_monoms = [ [] for _ in range(len(rep)) ]

402 used_indices = set()

403

404 for gen in new_gens:

405 try:

406 j = gens.index(gen)

407 used_indices.add(j)

408

409 for M, new_M in zip(monoms, new_monoms):

410 new_M.append(M[j])

411 except ValueError:

412 for new_M in new_monoms:

413 new_M.append(0)

414

415 for i, _ in enumerate(gens):

416 if i not in used_indices:

417 for monom in monoms:

418 if monom[i]:

419 raise GeneratorsError("unable to drop generators")

420

421 return map(tuple, new_monoms), coeffs

422

423

424class PicklableWithSlots:

425 """

426 Mixin class that allows to pickle objects with ``__slots__``.

427

428 Examples

429 ========

430

431 First define a class that mixes :class:`PicklableWithSlots` in::

432

433 >>> from sympy.polys.polyutils import PicklableWithSlots

434 >>> class Some(PicklableWithSlots):

435 ... __slots__ = ('foo', 'bar')

436 ...

437 ... def __init__(self, foo, bar):

438 ... self.foo = foo

439 ... self.bar = bar

440

441 To make :mod:`pickle` happy in doctest we have to use these hacks::

442

443 >>> import builtins

444 >>> builtins.Some = Some

445 >>> from sympy.polys import polyutils

446 >>> polyutils.Some = Some

447

448 Next lets see if we can create an instance, pickle it and unpickle::

449

450 >>> some = Some('abc', 10)

451 >>> some.foo, some.bar

452 ('abc', 10)

453

454 >>> from pickle import dumps, loads

455 >>> some2 = loads(dumps(some))

456

457 >>> some2.foo, some2.bar

458 ('abc', 10)

459

460 """

461

462 __slots__ = ()

463

464 def __getstate__(self, cls=None):

465 if cls is None:

466 # This is the case for the instance that gets pickled

467 cls = self.__class__

468

469 d = {}

470

471 # Get all data that should be stored from super classes

472 for c in cls.__bases__:

473 # XXX: Python 3.11 defines object.__getstate__ and it does not

474 # accept any arguments so we need to make sure not to call it with

475 # an argument here. To be compatible with Python < 3.11 we need to

476 # be careful not to assume that c or object has a __getstate__

477 # method though.

478 getstate = getattr(c, "__getstate__", None)

479 objstate = getattr(object, "__getstate__", None)

480 if getstate is not None and getstate is not objstate:

481 d.update(getstate(self, c))

482

483 # Get all information that should be stored from cls and return the dict

484 for name in cls.__slots__:

485 if hasattr(self, name):

486 d[name] = getattr(self, name)

487

488 return d

489

490 def __setstate__(self, d):

491 # All values that were pickled are now assigned to a fresh instance

492 for name, value in d.items():

493 try:

494 setattr(self, name, value)

495 except AttributeError: # This is needed in cases like Rational :> Half

496 pass

497

498

499class IntegerPowerable:

500 r"""

501 Mixin class for classes that define a `__mul__` method, and want to be

502 raised to integer powers in the natural way that follows. Implements

503 powering via binary expansion, for efficiency.

504

505 By default, only integer powers $\geq 2$ are supported. To support the

506 first, zeroth, or negative powers, override the corresponding methods,

507 `_first_power`, `_zeroth_power`, `_negative_power`, below.

508 """

509

510 def __pow__(self, e, modulo=None):

511 if e < 2:

512 try:

513 if e == 1:

514 return self._first_power()

515 elif e == 0:

516 return self._zeroth_power()

517 else:

518 return self._negative_power(e, modulo=modulo)

519 except NotImplementedError:

520 return NotImplemented

521 else:

522 bits = [int(d) for d in reversed(bin(e)[2:])]

523 n = len(bits)

524 p = self

525 first = True

526 for i in range(n):

527 if bits[i]:

528 if first:

529 r = p

530 first = False

531 else:

532 r *= p

533 if modulo is not None:

534 r %= modulo

535 if i < n - 1:

536 p *= p

537 if modulo is not None:

538 p %= modulo

539 return r

540

541 def _negative_power(self, e, modulo=None):

542 """

543 Compute inverse of self, then raise that to the abs(e) power.

544 For example, if the class has an `inv()` method,

545 return self.inv() ** abs(e) % modulo

546 """

547 raise NotImplementedError

548

549 def _zeroth_power(self):

550 """Return unity element of algebraic struct to which self belongs."""

551 raise NotImplementedError

552

553 def _first_power(self):

554 """Return a copy of self."""

555 raise NotImplementedError