NZMATH  1.2.0
About: NZMATH is a Python based number theory oriented calculation system.
  Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

uniutil.py
Go to the documentation of this file.
1 """
2 Utilities for univar.
3 
4 The module provides higher level interfaces to univar classes and
5 functions.
6 """
7 
8 from __future__ import division
9 import logging
10 import nzmath.arith1 as arith1
11 import nzmath.bigrandom as bigrandom
12 import nzmath.gcd as gcd
13 import nzmath.rational as rational
14 import nzmath.ring as ring
15 import nzmath.poly.univar as univar
16 import nzmath.poly.termorder as termorder
17 import nzmath.poly.ring as poly_ring
18 import nzmath.poly.ratfunc as ratfunc
20 
21 
22 _log = logging.getLogger("nzmath.poly.uniutil")
23 
24 
25 _MIXIN_MSG = "%s is mix-in"
26 
27 
28 class OrderProvider(object):
29  """
30  OrderProvider provides order and related operations.
31  """
32  def __init__(self, order=termorder.ascending_order):
33  """
34  Do not instantiate OrderProvider.
35  This initializer should be called from descendant:
36  OrderProvider.__init__(self, order)
37  where order is default to termorder.ascending_order.
38  """
39  if type(self) is OrderProvider:
40  raise NotImplementedError(_MIXIN_MSG % self.__class__.__name__)
41  self.order = order
42 
43  def shift_degree_to(self, degree):
44  """
45  Return polynomial whose degree is the given degree.
46 
47  More precisely, let f(X) = a_0 + ... + a_n * X**n, then
48  order.shift_degree_to(f, m) returns:
49  - zero polynomial, if f is zero polynomial
50  - a_(n-m) + ... + a_n * X**m, if 0 <= m < n
51  - a_0 * X**(m-n) + ... + a_n * X**m, if m >= n
52  """
53  original_degree = self.order.degree(self)
54  difference = degree - original_degree
55  if difference == 0:
56  return self
57  elif difference < 0:
58  return self.construct_with_default([(d + difference, c) for (d, c) in self if d + difference >= 0])
59  else:
60  return self.upshift_degree(difference)
61 
62  def split_at(self, degree):
63  """
64  Return tuple of two polynomials, which are splitted at the
65  given degree. The term of the given degree, if exists,
66  belongs to the lower degree polynomial.
67  """
68  lower, upper = [], []
69  for d, c in self:
70  if self.order.cmp(d, degree) <= 0:
71  lower.append((d, c))
72  else:
73  upper.append((d, c))
74  return (self.construct_with_default(lower),
75  self.construct_with_default(upper))
76 
77 
78 class DivisionProvider(object):
79  """
80  DivisionProvider provides all kind of divisions for univariate
81  polynomials. It is assumed that the coefficient ring of the
82  polynomials is a field.
83 
84  The class should be used as a mix-in.
85 
86  REQUIRE: OrderProvider
87  """
88  def __init__(self):
89  """
90  Do not instantiate DivisionProvider.
91  This initializer should be called from descendant.
92  """
93  if type(self) is DivisionProvider:
94  raise NotImplementedError(_MIXIN_MSG % self.__class__.__name__)
95  self._reduced = {}
96 
97  def __divmod__(self, other):
98  """
99  divmod(self, other)
100 
101  The method does, as the built-in divmod, return a tuple of
102  (self // other, self % other).
103  """
104  degree, lc = self.order.leading_term(other)
105  quotient, remainder = [], self
106  while self.order.degree(remainder) >= degree:
107  rdegree, rlc = self.order.leading_term(remainder)
108  q = rdegree - degree, rlc / lc
109  remainder = remainder - other.term_mul(q)
110  quotient.append(q)
111  quotient = self.construct_with_default(quotient)
112  return quotient, remainder
113 
114  def __floordiv__(self, other):
115  """
116  self // other
117  """
118  degree, lc = self.order.leading_term(other)
119  quotient, remainder = [], self
120  while self.order.degree(remainder) >= degree:
121  rdegree, rlc = self.order.leading_term(remainder)
122  q = rdegree - degree, rlc / lc
123  remainder = remainder - other.term_mul(q)
124  quotient.append(q)
125  return self.construct_with_default(quotient)
126 
127  def __mod__(self, other):
128  """
129  self % other
130  """
131  degree, lc = self.order.leading_term(other)
132  remainder = self
133  rdegree, rlc = self.order.leading_term(remainder)
134  if rdegree > degree + 5 and degree >= 1:
135  return other.mod(remainder)
136  ilc = ring.inverse(lc)
137  while rdegree >= degree:
138  q = rdegree - degree, rlc * ilc
139  remainder = remainder - other.term_mul(q)
140  rdegree, rlc = self.order.leading_term(remainder)
141  return remainder
142 
143  def mod(self, dividend):
144  """
145  Return dividend % self.
146 
147  self should have attribute _reduced to cache reduced monomials.
148  """
149  degree, lc = self.order.leading_term(self)
150  div_deg = self.order.degree(dividend)
151  if div_deg < degree:
152  return dividend
153  upperbound = min(degree * 2, div_deg) + 1
154  populate = False
155  if not self._reduced:
156  populate = True
157  else:
158  i = min(self._reduced.keys())
159  while i in self._reduced.keys():
160  i += 1
161  if i < upperbound:
162  populate = True
163  if populate:
164  self._populate_reduced(degree, lc, upperbound)
165  if div_deg > degree * 2:
166  dividend_degrees = sorted(dividend.iterbases(), reverse=True)
167  dividend_degrees = [deg for deg in dividend_degrees if deg > degree]
168  self._populate_reduced_more(dividend_degrees)
169  accum = self.construct_with_default(())
170  lowers = []
171  for d, c in dividend:
172  if c:
173  if d < degree:
174  lowers.append((d, c))
175  else:
176  accum += self._reduced[d].scalar_mul(c)
177  return accum + self.construct_with_default(lowers)
178 
179  def mod_pow(self, polynom, index):
180  """
181  Return polynom ** index % self.
182  """
183  if not self:
184  raise ZeroDivisionError("polynomial division or modulo by zero.")
185  polynom %= self
186  if index == 1 or not polynom:
187  return polynom
188  elif polynom.degree() == 0:
189  return self.getRing().createElement([(0, polynom[0]**index)])
190  acoefficient = polynom.itercoefficients().next()
191  one = ring.getRing(acoefficient).one
192  power_product = self.construct_with_default([(0, one)])
193  if index:
194  power_of_2 = polynom
195  while index:
196  if index & 1:
197  power_product = self.mod(power_product * power_of_2)
198  power_of_2 = self.mod(power_of_2.square())
199  index //= 2
200  return power_product
201 
202  def _populate_reduced(self, degree, lc, upperbound):
203  """
204  Populate self._reduced.
205 
206  degree, lc is of self, and self._reduced is populated up to
207  the given upperbound.
208  """
209  one = ring.getRing(self.itercoefficients().next()).one
210  if not self._reduced:
211  minimum = degree
212  redux = self.construct_with_default([(degree - 1, one)])
213  else:
214  minimum = max(self._reduced.keys()) + 1
215  redux = self._reduced[minimum - 1]
216  moniced = self.scalar_mul(one / lc)
217  for i in range(minimum, upperbound):
218  redux = redux.term_mul((1, 1))
219  if self.order.degree(redux) == degree:
220  redux -= moniced.scalar_mul(self.order.leading_coefficient(redux))
221  self._reduced[i] = redux
222 
223  def _populate_reduced_more(self, degrees):
224  """
225  Populate self._reduced more for much higher degree dividend.
226  This method has to be called after _populate_reduced.
227 
228  self._reduced is populated so that it will include all of
229  degrees > degree * 2. The degrees are recommended to be
230  sorted in descending order.
231  """
232  degree = self.order.degree(self)
233  if not [deg for deg in degrees if deg not in self._reduced]:
234  # no need to populate more
235  return
236  # self._reduced has every key from degree up to maxreduced.
237  # The default (prepared by _populate_reduced) is degree*2.
238  maxreduced = degree * 2
239  while (maxreduced + 1) in self._reduced:
240  maxreduced += 1
241  # use binary chain multiplied by a stride
242  stride = maxreduced
243  if len(degrees) > 1:
244  # try to use wider stride
245  common_degree = gcd.gcd_of_list(degrees)[0]
246  if common_degree > maxreduced:
247  # common_degree is better than maxreduced.
248  if common_degree not in self._reduced:
249  self._populate_reduced_more([common_degree])
250  stride = common_degree
251  redux = self._reduced[stride]
252  maximum = max(degrees)
253  binary = {stride:redux}
254  while stride * 2 <= maximum:
255  stride += stride
256  if stride not in self._reduced:
257  redux = self.mod(redux.square())
258  else:
259  redux = self._reduced[stride]
260  binary[stride] = redux
261  binarykeys = sorted(binary.keys(), reverse=True)
262  for deg in (d for d in degrees if d > maxreduced):
263  pickup = []
264  rest = deg
265  for key in binarykeys:
266  if rest < key:
267  continue
268  rest -= key
269  pickup.append(key)
270  if rest < maxreduced or rest in self._reduced:
271  break
272  assert pickup
273  total = pickup.pop() # start with the smallest degree picked up
274  prod = binary[total]
275  for picked in reversed(pickup):
276  total += picked
277  prod = self.mod(prod * binary[picked])
278  self._reduced[total] = prod
279  if rest in self._reduced:
280  final = self.mod(prod * self._reduced[rest])
281  else: # rest < degree
282  final = self.mod(prod.term_mul((rest, 1)))
283  self._reduced[deg] = final
284 
285  def __truediv__(self, other):
286  """
287  self / other
288 
289  The result is a rational function.
290  """
291  return ratfunc.RationalFunction(self, other)
292 
293  def scalar_exact_division(self, scale):
294  """
295  Return quotient by a scalar which can divide each coefficient
296  exactly.
297  """
298  return self.coefficients_map(lambda c: c.exact_division(scale))
299 
300  def gcd(self, other):
301  """
302  Return a greatest common divisor of self and other.
303  Returned polynomial is always monic.
304  """
305  if self.order.degree(self) < self.order.degree(other):
306  divident, divisor = other, self
307  else:
308  divident, divisor = self, other
309  while divisor:
310  divident, divisor = divisor, divident % divisor
311  lc = self.order.leading_coefficient(divident)
312  if lc and lc != ring.getRing(lc).one:
313  divident = divident.scalar_exact_division(lc)
314  return divident
315 
316  def extgcd(self, other):
317  """
318  Return a tuple (u, v, d); they are the greatest common divisor
319  d of two polynomials self and other and u, v such that
320  d = self * u + other * v.
321 
322  see nzmath.gcd.extgcd
323  """
324  order = termorder.ascending_order
325  polynomial_ring = self.getRing()
326  zero, one = polynomial_ring.zero, polynomial_ring.one
327  a, b, g, u, v, w = one, zero, self, zero, one, other
328  while w:
329  q = g // w
330  a, b, g, u, v, w = u, v, w, a - q*u, b - q*v, g - q*w
331  lc = self.order.leading_coefficient(g)
332  if lc and lc != 1:
333  linv = lc.inverse()
334  a, b, g = linv * a, linv * b, linv * g
335  return (a, b, g)
336 
337 
339  """
340  PseudoDivisionProvider provides pseudo divisions for univariate
341  polynomials. It is assumed that the coefficient ring of the
342  polynomials is a domain.
343 
344  The class should be used as a mix-in.
345  """
346  def pseudo_divmod(self, other):
347  """
348  self.pseudo_divmod(other) -> (Q, R)
349 
350  Q, R are polynomials such that
351  d**(deg(self) - deg(other) + 1) * self == other * Q + R,
352  where d is the leading coefficient of other.
353  """
354  order = termorder.ascending_order
355  if hasattr(self, 'order'):
356  assert self.order is order
357  degree, lc = order.leading_term(other)
358  # step 1
359  quotient, remainder = self.construct_with_default(()), self
360  rdegree, rlc = order.leading_term(remainder)
361  e = rdegree - degree + 1
362  if e <= 0:
363  return quotient, remainder
364  while rdegree >= degree:
365  # step 3
366  canceller = self.construct_with_default([(rdegree - degree, rlc)])
367  quotient = quotient.scalar_mul(lc) + canceller
368  remainder = remainder.scalar_mul(lc) - canceller * other
369  e -= 1
370  rdegree, rlc = order.leading_term(remainder)
371  # step 2
372  q = lc ** e
373  quotient = quotient.scalar_mul(q)
374  remainder = remainder.scalar_mul(q)
375  return quotient, remainder
376 
377  def pseudo_floordiv(self, other):
378  """
379  self.pseudo_floordiv(other) -> Q
380 
381  Q is a polynomial such that
382  d**(deg(self) - deg(other) + 1) * self == other * Q + R,
383  where d is the leading coefficient of other and R is a
384  polynomial.
385  """
386  order = termorder.ascending_order
387  if hasattr(self, 'order'):
388  assert self.order is order
389  degree, lc = order.leading_term(other)
390  # step 1
391  quotient, remainder = self.construct_with_default(()), self
392  rdegree, rlc = order.leading_term(remainder)
393  e = order.degree(remainder) - degree + 1
394  if e <= 0:
395  return quotient
396  while rdegree >= degree:
397  # step 3
398  canceller = self.construct_with_default([(rdegree - degree, rlc)])
399  quotient = quotient.scalar_mul(lc) + canceller
400  remainder = remainder.scalar_mul(lc) - canceller * other
401  e -= 1
402  rdegree, rlc = order.leading_term(remainder)
403  # step 2
404  return quotient.scalar_mul(lc ** e)
405 
406  def pseudo_mod(self, other):
407  """
408  self.pseudo_mod(other) -> R
409 
410  R is a polynomial such that
411  d**(deg(self) - deg(other) + 1) * self == other * Q + R,
412  where d is the leading coefficient of other and Q a
413  polynomial.
414  """
415  order = termorder.ascending_order
416  if hasattr(self, 'order'):
417  assert self.order is order
418  degree, lc = order.leading_term(other)
419  # step 1
420  remainder = self
421  rdegree, rlc = order.leading_term(remainder)
422  e = order.degree(remainder) - degree + 1
423  if e <= 0:
424  return remainder
425  while rdegree >= degree:
426  # step 3
427  canceller = self.construct_with_default([(rdegree - degree, rlc)])
428  remainder = remainder.scalar_mul(lc) - canceller * other
429  e -= 1
430  rdegree, rlc = order.leading_term(remainder)
431  # step 2
432  return remainder.scalar_mul(lc ** e)
433 
434  def __truediv__(self, other):
435  """
436  self / other
437 
438  Return the result as a rational function.
439  """
440  order = termorder.ascending_order
441  return ratfunc.RationalFunction(self, other)
442 
443  def exact_division(self, other):
444  """
445  Return quotient of exact division.
446  """
447  order = termorder.ascending_order
448  if hasattr(self, 'order'):
449  assert self.order is order
450  quotient, remainder = self.pseudo_divmod(other)
451  if not remainder:
452  deg_other, lc = order.leading_term(other)
453  deg_self = order.degree(self)
454  extra_factor = lc ** (deg_self - deg_other + 1)
455  return quotient.scalar_exact_division(extra_factor)
456  raise ValueError("division is not exact")
457 
458  def scalar_exact_division(self, scale):
459  """
460  Return quotient by a scalar which can divide each coefficient
461  exactly.
462  """
463  return self.coefficients_map(lambda c: c.exact_division(scale))
464 
465  def monic_divmod(self, other):
466  """
467  self.monic_divmod(other) -> (Q, R)
468 
469  Q, R are polynomials such that
470  self == other * Q + R.
471 
472  The leading coefficient of other MUST be one.
473  """
474  order = termorder.ascending_order
475  if hasattr(self, 'order'):
476  assert self.order is order
477 
478  degree = order.degree(other)
479  quotient, remainder = self.construct_with_default(()), self
480  rdegree, rlc = order.leading_term(remainder)
481  if rdegree < degree:
482  return quotient, remainder
483  while rdegree >= degree:
484  # step 3
485  canceller = self.construct_with_default([(rdegree - degree, rlc)])
486  quotient += canceller
487  remainder -= canceller * other
488  rdegree, rlc = order.leading_term(remainder)
489 
490  return quotient, remainder
491 
492  def monic_floordiv(self, other):
493  """
494  self.monic_floordiv(other) -> Q
495 
496  Q is a polynomial such that
497  self == other * Q + R,
498  where R is a polynomial whose degree is smaller than other's.
499 
500  The leading coefficient of other MUST be one.
501  """
502  order = termorder.ascending_order
503  if hasattr(self, 'order'):
504  assert self.order is order
505  degree = order.degree(other)
506  # step 1
507  quotient, remainder = self.construct_with_default(()), self
508  rdegree, rlc = order.leading_term(remainder)
509  if rdegree < degree:
510  return quotient
511  while rdegree >= degree:
512  # step 3
513  canceller = self.construct_with_default([(rdegree - degree, rlc)])
514  quotient += canceller
515  remainder -= canceller * other
516  rdegree, rlc = order.leading_term(remainder)
517 
518  return quotient
519 
520  def monic_mod(self, other):
521  """
522  self.monic_mod(other) -> R
523 
524  R is a polynomial such that
525  self == other * Q + R,
526  where Q is another polynomial.
527 
528  The leading coefficient of other MUST be one.
529  """
530  order = termorder.ascending_order
531  if hasattr(self, 'order'):
532  assert self.order is order
533  degree = order.degree(other)
534  remainder = self
535  rdegree, rlc = order.leading_term(remainder)
536  if rdegree < degree:
537  return remainder
538  while rdegree >= degree:
539  canceller = self.construct_with_default([(rdegree - degree, rlc)])
540  remainder -= canceller * other
541  rdegree, rlc = order.leading_term(remainder)
542 
543  return remainder
544 
545  def monic_pow(self, index, mod):
546  """
547  Return self**index % mod.
548 
549  The leading coefficient of mod MUST be one.
550  """
551  if index < 0:
552  raise ValueError("negative index is not allowed.")
553  # special indices
554  elif index == 0:
555  return self.getRing().one
556  elif index == 1:
557  return self.monic_mod(mod)
558  elif index == 2:
559  return self.square().monic_mod(mod)
560  # special polynomials
561  if not self:
562  return self
563  # general
564  power_product = self.getRing().one
565  power_of_2 = self.monic_mod(mod)
566  while index:
567  if index & 1:
568  power_product = (power_product * power_of_2).monic_mod(mod)
569  index //= 2
570  if index:
571  power_of_2 = power_of_2.square().monic_mod(mod)
572  return power_product
573 
574 
575 class ContentProvider(object):
576  """
577  ContentProvider provides content and primitive part.
578 
579  The coefficient ring must be a unique factorization domain.
580  """
581  def content(self):
582  """
583  Return content of the polynomial.
584  """
585  coeffring = None
586  isquotfield = None
587  cont = 0
588  num, den = 0, 1
589  for c in self.itercoefficients():
590  if isquotfield is None:
591  coeffring = ring.getRing(c)
592  if coeffring.isfield() and isinstance(coeffring, ring.QuotientField):
593  isquotfield = True
594  basedomain = coeffring.basedomain
595  else:
596  isquotfield = False
597  if isquotfield:
598  num = basedomain.gcd(num, c.numerator)
599  den = basedomain.lcm(den, c.denominator)
600  else:
601  if not cont:
602  cont = c
603  else:
604  cont = coeffring.gcd(cont, c)
605  if isquotfield:
606  cont = coeffring.createElement(num, den)
607  return cont
608 
609  def primitive_part(self):
610  """
611  Return the primitive part of the polynomial.
612  """
613  return self.scalar_exact_division(self.content())
614 
615 
617  """
618  SubresultantGcdProvider provides gcd method using subresultant
619  algorithm.
620 
621  REQUIRE: PseudoDivisionProvider, ContentProvider
622  """
623  def resultant(self, other):
624  """
625  Return the resultant of self and other.
626  """
627  order = termorder.ascending_order
628  f, g = self, other
629  negative = False
630  if order.degree(f) < order.degree(g):
631  f, g = g, f
632  if (order.degree(f) & 1) and (order.degree(g) & 1):
633  negative = not negative
634  coeff_ring = self.getCoefficientRing()
635  a = b = coeff_ring.one
636  while order.degree(g) > 0:
637  delta = order.degree(f) - order.degree(g)
638  if (order.degree(f) & 1) and (order.degree(g) & 1):
639  negative = not negative
640  h = f.pseudo_mod(g)
641  f, g = g, h.scalar_exact_division(a * (b ** delta))
642  a = order.leading_coefficient(f)
643  b = ((a ** delta) * b).exact_division(b ** delta)
644  if not g:
645  return coeff_ring.zero
646  scalar = g[0]
647  degree_f = order.degree(f)
648  if negative:
649  return (-b * scalar ** degree_f).exact_division(b ** degree_f)
650  else:
651  return (b * scalar ** degree_f).exact_division(b ** degree_f)
652 
653  def subresultant_gcd(self, other):
654  """
655  Return the greatest common divisor of given polynomials. They
656  must be in the polynomial ring and its coefficient ring must
657  be a UFD.
658 
659  Reference: Algorithm 3.3.1 of Cohen's book
660  """
661  order = termorder.ascending_order
662  divident = self
663  divisor = other
664  polynomial_ring = self.getRing()
665  coeff_ring = polynomial_ring.getCoefficientRing()
666  one = coeff_ring.one
667  # step 1
668  if order.degree(divisor) > order.degree(divident):
669  divident, divisor = divisor, divident
670  if not divisor:
671  return divident
672  content_gcd = coeff_ring.gcd(divident.content(), divisor.content())
673  divident = divident.primitive_part()
674  divisor = divisor.primitive_part()
675  g = h = one
676 
677  while 1:
678  # step 2
679  delta = order.degree(divident) - order.degree(divisor)
680  quotient, remainder = divident.pseudo_divmod(divisor)
681  if not remainder:
682  return divisor.primitive_part().scalar_mul(content_gcd)
683  if order.degree(remainder) == 0:
684  return polynomial_ring.createElement(content_gcd)
685 
686  # step 3
687  divident, divisor = divisor, remainder
688  if delta == 0 and g != one:
689  divisor = divisor.scalar_exact_division(g)
690  elif delta == 1 and (g != one or h != one):
691  divisor = divisor.scalar_exact_division(g * h)
692  elif delta > 1 and (g != one or h != one):
693  divisor = divisor.scalar_exact_division(g * h**delta)
694  g = divident.leading_coefficient()
695  if delta == 1 and h != g:
696  h = g
697  elif delta > 1 and (g != one or h != one):
698  h = g * (h * g) ** (delta - 1)
699 
700  def subresultant_extgcd(self, other):
701  """
702  Return (A, B, P) s.t. A*self+B*other=P,
703  where P is the greatest common divisor of given polynomials.
704  They must be in the polynomial ring and its coefficient ring must
705  be a UFD.
706 
707  Reference: Kida's paper p.18
708  """
709  order = termorder.ascending_order
710  coeff_ring = self.getCoefficientRing()
711  P_1, P_2 = self, other
712  A_1, A_2 = coeff_ring.one, coeff_ring.zero
713  if order.degree(P_1) < order.degree(P_2):
714  P_1, P_2 = P_2, P_1
715  A_1, A_2 = A_2, A_1
716  if not P_2:
717  return (A_1, A_2, P_1)
718  a = b = coeff_ring.one
719  while P_2:
720  delta = order.degree(P_1) - order.degree(P_2)
721  b = a * (b ** delta)
722  Q, R = P_1.pseudo_divmod(P_2)
723  a = order.leading_coefficient(P_2)
724  ad = a ** delta
725  P_1, P_2 = P_2, R.scalar_exact_division(b)
726  A_1, A_2 = A_2, (ad * a * A_1 - Q * A_2).scalar_exact_division(b)
727  b = (ad * b) // (b ** delta)
728  B_1 = (P_1 - A_1 * self).exact_division(other)
729  return (A_1, B_1, P_1)
730 
731 
733  """
734  PrimeCharacteristicFunctionsProvider provides efficient powering
735  and factorization for polynomials whose coefficient ring has prime
736  characteristic.
737 
738  - A client of this mix-in class should use DivisionProvider also.
739  - A client of this mix-in class must have attribute ch, which
740  stores the prime characteristic of the coefficient ring.
741  """
742  def __init__(self, ch):
743  """
744  Do not instantiate PrimeCharacteristicFunctionsProvider.
745  This initializer should be called from descendant.
746  """
747  if type(self) is PrimeCharacteristicFunctionsProvider:
748  raise NotImplementedError(_MIXIN_MSG % self.__class__.__name__)
749  self.ch = ch
750 
751  def __pow__(self, index, mod=None):
752  """
753  self ** index
754 
755  A client of the mix-in class should write the name in front of
756  the main class so that the __pow__ defined here would be
757  preceded in method resolution order. For example:
758 
759  class Client (PrimeCharacteristicFunctionsProvider
760  DivisionProvider,
761  Polynomial):
762  """
763  if not isinstance(index, (int, long)):
764  raise TypeError("index must be an integer.")
765  if index < 0:
766  raise ValueError("index must be a non-negative integer.")
767  if mod is not None:
768  return mod.mod_pow(self, index)
769  if index > 0:
770  if not self:
771  return self
772  q = 1
773  while index % self.ch == 0:
774  q *= self.ch
775  index //= self.ch
776  if q > 1:
777  powered = self.construct_with_default([(d * q, c ** q) for (d, c) in self])
778  else:
779  powered = self
780  if index == 1:
781  return powered
782  acoefficient = self.itercoefficients().next()
783  one = ring.getRing(acoefficient).one
784  power_product = self.construct_with_default([(0, one)])
785  if index:
786  power_of_2 = powered
787  while index:
788  if index & 1:
789  power_product *= power_of_2
790  power_of_2 = power_of_2.square()
791  index //= 2
792  return power_product
793 
794  def mod_pow(self, polynom, index):
795  """
796  Return polynom ** index % self.
797  """
798  if not self:
799  raise ZeroDivisionError("polynomial division or modulo by zero.")
800  polynom %= self
801  if index == 1 or not polynom:
802  return polynom
803  elif polynom.degree() == 0:
804  return self.getRing().createElement([(0, polynom[0]**index)])
805  acoefficient = polynom.itercoefficients().next()
806  cardfq = card(ring.getRing(acoefficient))
807  final_product = self.getRing().one
808  qpow = polynom
809  # q-th power expansion
810  while index:
811  index, small = divmod(index, cardfq)
812  if small == 1:
813  final_product *= qpow
814  elif small:
815  final_product *= self._small_index_mod_pow(qpow, small)
816  # c ** q = c for any Fq element c, thus
817  qpow = self.mod(qpow.bases_map(lambda d: d * cardfq))
818 
819  return final_product
820 
821  def _small_index_mod_pow(self, polynom, index):
822  """
823  Return polynom ** index % self for small index (0 < index < q).
824  """
825  # binary powering
826  power_product = self.getRing().one
827  if index > 0:
828  power_of_2 = polynom
829  while index:
830  if index & 1:
831  power_product = self.mod(power_product * power_of_2)
832  if index == 1:
833  break
834  power_of_2 = self.mod(power_of_2.square())
835  index //= 2
836  else:
837  raise ValueError("this private method requires index > 0")
838  return power_product
839 
841  """
842  Return the square free decomposition of the polynomial. The
843  return value is a dict whose keys are integers and values are
844  corresponding powered factors. For example, if
845  A = A1 * A2**2,
846  the result is {1: A1, 2: A2}.
847 
848  gcd method is required.
849  """
850  result = {}
851  if self.order.degree(self) == 1:
852  return {1: self}
853  f = self
854  df = f.differentiate()
855  if df:
856  b = f.gcd(df)
857  a = f.exact_division(b)
858  i = 1
859  while self.order.degree(a) > 0:
860  c = a.gcd(b)
861  b = b.exact_division(c)
862  if a != c:
863  r = a.exact_division(c)
864  if self.order.degree(r) > 0:
865  result[i] = r
866  a = c
867  i += 1
868  f = b
869  if self.order.degree(f) > 0:
870  f = f.pthroot()
871  subresult = f.squarefree_decomposition()
872  for i, g in subresult.iteritems():
873  result[i*self.ch] = g
874  return result
875 
876  def pthroot(self):
877  """
878  Return a polynomial obtained by sending X**p to X, where p is
879  the characteristic. If the polynomial does not consist of
880  p-th powered terms only, result is nonsense.
881  """
882  return self.construct_with_default([(d // self.ch, c) for (d, c) in self if c])
883 
885  """
886  Return the distinct degree factorization of the polynomial.
887  The return value is a dict whose keys are integers and values
888  are corresponding product of factors of the degree. For
889  example, if A = A1 * A2, and all irreducible factors of A1
890  having degree 1 and all irreducible factors of A2 having
891  degree 2, then the result is:
892  {1: A1, 2: A2}.
893 
894  The given polynomial must be square free, and its coefficient
895  ring must be a finite field.
896  """
897  Fq = ring.getRing(self.itercoefficients().next())
898  q = card(Fq)
899  f = self
900  x = f.construct_with_default([(1, Fq.one)])
901  w = x
902  i = 1
903  result = {}
904  while i*2 <= self.order.degree(f):
905  w = pow(w, q, f)
906  result[i] = f.gcd(w - x)
907  if self.order.degree(result[i]) > 0:
908  f = f.exact_division(result[i])
909  w = w % f
910  else:
911  del result[i]
912  i += 1
913  if self.order.degree(f) != 0:
914  result[self.order.degree(f)] = f
915  return result
916 
917  def split_same_degrees(self, degree):
918  """
919  Return the irreducible factors of the polynomial. The
920  polynomial must be a product of irreducible factors of the
921  given degree.
922  """
923  r = self.order.degree(self) // degree
924  Fq = ring.getRing(self.itercoefficients().next())
925  q = card(Fq)
926  p = Fq.getCharacteristic()
927  if degree == 1:
928  result = []
929  X = self.construct_with_default([(1, Fq.one)])
930  f = self
931  while not f[0]:
932  f = f // X
933  result.append(X)
934  if self.order.degree(f) >= 1:
935  result.append(f)
936  else:
937  result = [self]
938  while len(result) < r:
939  # g is a random polynomial
940  rand_coeff = {}
941  for i in range(2 * degree):
942  rand_coeff[i] = Fq.createElement(bigrandom.randrange(q))
943  if not rand_coeff[2 * degree - 1]:
944  rand_coeff[2 * degree - 1] = Fq.one
945  randpoly = self.construct_with_default(rand_coeff)
946  if p == 2:
947  G = self.construct_with_default(())
948  for i in range(degree):
949  G = G + randpoly
950  randpoly = self.mod(randpoly.square())
951  else:
952  one = self.construct_with_default([(0, Fq.one)])
953  G = pow(randpoly, (q**degree - 1) >> 1, self) - one
954  subresult = []
955  while result:
956  h = result.pop()
957  if self.order.degree(h) == degree:
958  subresult.append(h)
959  continue
960  z = h.gcd(G)
961  if 0 < self.order.degree(z) < self.order.degree(h):
962  subresult.append(z)
963  subresult.append(h.exact_division(z))
964  else:
965  subresult.append(h)
966  result = subresult
967  return result
968 
969  def factor(self):
970  """
971  Factor the polynomial.
972 
973  The returned value is a list of tuples whose first component
974  is a factor and second component is its multiplicity.
975  """
976  result = []
977  lc = self.order.leading_coefficient(self)
978  if lc != ring.getRing(lc).one:
979  self = self.scalar_exact_division(lc)
980  result.append((self.getRing().one*lc, 1))
981  squarefreefactors = self.squarefree_decomposition()
982  for m, f in squarefreefactors.iteritems():
983  distinct_degree_factors = f.distinct_degree_factorization()
984  for d, g in distinct_degree_factors.iteritems():
985  if d == self.order.degree(g):
986  result.append((g, m))
987  else:
988  for irred in g.split_same_degrees(d):
989  result.append((irred, m))
990  return result
991 
992  def isirreducible(self):
993  """
994  f.isirreducible() -> bool
995 
996  If f is irreducible return True, otherwise False.
997  """
998  if 0 <= self.order.degree(self) <= 1:
999  # degree 1 polynomials and constants are irreducible
1000  return True
1001  elif not self[0]:
1002  # degree >= 2 polynomials are reducible if constant term is zero
1003  return False
1004  squareFreeFactors = self.squarefree_decomposition()
1005  if len(squareFreeFactors) != 1:
1006  return False
1007  m, f = squareFreeFactors.popitem()
1008  if m != 1:
1009  return False
1010  distinctDegreeFactors = f.distinct_degree_factorization()
1011  if len(distinctDegreeFactors) != 1:
1012  return False
1013  d, g = distinctDegreeFactors.popitem()
1014  if d != self.order.degree(g):
1015  return False
1016  return True
1017 
1018 
1019 class KaratsubaProvider(object):
1020  """
1021  define Karatsuba method multiplication and squaring.
1022 
1023  REQUIRE: OrderProvider
1024  """
1025  def ring_mul_karatsuba(self, other):
1026  """
1027  Multiplication of two polynomials in the same ring.
1028 
1029  Computation is carried out by Karatsuba method.
1030  """
1031  # zero
1032  if not self or not other:
1033  return self.construct_with_default(())
1034  # monomial
1035  if len(self) == 1:
1036  return other.term_mul(self)
1037  if len(other) == 1:
1038  return self.term_mul(other)
1039  # binomial
1040  if len(self) == 2:
1041  p, q = [other.term_mul(term) for term in self]
1042  return p + q
1043  if len(other) == 2:
1044  p, q = [self.term_mul(term) for term in other]
1045  return p + q
1046  # suppose self is black and other is red.
1047  black_least_degree = self.order.tail_degree(self)
1048  black_most_degree = self.order.degree(self)
1049  red_least_degree = self.order.tail_degree(other)
1050  red_most_degree = self.order.degree(other)
1051  least_degree = min(black_least_degree, red_least_degree)
1052  most_degree = max(black_most_degree, red_most_degree)
1053  assert least_degree < most_degree
1054  half_degree = (least_degree + most_degree) >> 1
1055 
1056  if black_least_degree > half_degree:
1057  return self.downshift_degree(black_least_degree).ring_mul_karatsuba(other).upshift_degree(black_least_degree)
1058  if red_least_degree > half_degree:
1059  return self.ring_mul_karatsuba(other.downshift_degree(red_least_degree)).upshift_degree(red_least_degree)
1060 
1061  black = self.downshift_degree(least_degree)
1062  red = other.downshift_degree(least_degree)
1063  club, spade = black.split_at(half_degree - least_degree)
1064  dia, heart = red.split_at(half_degree - least_degree)
1065  weaker = club.ring_mul_karatsuba(dia)
1066  stronger = spade.ring_mul_karatsuba(heart)
1067  karatsuba = (club + spade).ring_mul_karatsuba(dia + heart) - weaker - stronger
1068  return (weaker.upshift_degree(least_degree * 2) +
1069  karatsuba.upshift_degree(least_degree + half_degree) +
1070  stronger.upshift_degree(half_degree * 2))
1071 
1072  def square_karatsuba(self):
1073  """
1074  Return the square of self.
1075  """
1076  # zero
1077  if not self:
1078  return self
1079 
1080  polynomial = self.construct_with_default
1081  data_length = len(self)
1082  # monomial
1083  if data_length == 1:
1084  d, c = iter(self).next()
1085  if d:
1086  return polynomial([(d*2, c**2)])
1087  else:
1088  return polynomial([(0, c**2)])
1089  # binomial
1090  if data_length == 2:
1091  (d1, c1), (d2, c2) = [(d, c) for (d, c) in self]
1092  if "_sorted" in self._init_kwds:
1093  del self._init_kwds["_sorted"]
1094  return polynomial({d1*2:c1**2, d1+d2:c1*c2*2, d2*2:c2**2}, _sorted=False)
1095  # general (Karatsuba)
1096  least_degree = self.order.tail_degree(self)
1097  if least_degree:
1098  chopped = self.downshift_degree(least_degree)
1099  else:
1100  chopped = self
1101  half_degree = self.order.degree(self) >> 1
1102  fst, snd = chopped.split_at(half_degree)
1103  fst_squared = fst.square()
1104  snd_squared = snd.square()
1105  karatsuba = (fst + snd).square() - fst_squared - snd_squared
1106  result = (fst_squared +
1107  karatsuba.upshift_degree(half_degree) +
1108  snd_squared.upshift_degree(half_degree * 2))
1109  if least_degree:
1110  return result.upshift_degree(least_degree)
1111  else:
1112  return result
1113 
1114 
1115 class VariableProvider(object):
1116  """
1117  VariableProvider provides the variable name and other cariable
1118  related methods.
1119  """
1120  def __init__(self, varname):
1121  """
1122  Do not instantiate VariableProvider.
1123  This initializer should be called from descendant.
1124  """
1125  if type(self) is VariableProvider:
1126  raise NotImplementedError(_MIXIN_MSG % self.__class__.__name__)
1127  self.variable = varname
1128 
1129  def getVariable(self):
1130  """
1131  Get variable name
1132  """
1133  return self.variable
1134 
1135  def getVariableList(self):
1136  """
1137  Get variable name as list.
1138  """
1139  return [self.variable]
1140 
1141 
1143  """
1144  Provides interfaces for ring.CommutativeRingElement.
1145  """
1146  def __init__(self):
1147  """
1148  Do not instantiate RingElementProvider.
1149  This initializer should be called from descendant:
1150  RingElementProvider.__init__(self)
1151  """
1152  if type(self) is RingElementProvider:
1153  raise NotImplementedError(_MIXIN_MSG % self.__class__.__name__)
1154  ring.CommutativeRingElement.__init__(self)
1155  self._coefficient_ring = None
1156  self._ring = None
1157 
1158  def getRing(self):
1159  """
1160  Return an object of a subclass of Ring, to which the element
1161  belongs.
1162  """
1163  if self._coefficient_ring is None or self._ring is None:
1164  myring = None
1165  for c in self.itercoefficients():
1166  cring = ring.getRing(c)
1167  if not myring or myring != cring and myring.issubring(cring):
1168  myring = cring
1169  elif not cring.issubring(myring):
1170  myring = myring.getCommonSuperring(cring)
1171  if not myring:
1172  myring = rational.theIntegerRing
1173  self.set_coefficient_ring(myring)
1174  return self._ring
1175 
1176  def set_coefficient_ring(self, coeffring):
1177  if self._coefficient_ring is None:
1178  self._coefficient_ring = coeffring
1179  # variable names are ignored now
1180  self._ring = poly_ring.PolynomialRing.getInstance(self._coefficient_ring)
1181 
1182 
1185  RingElementProvider):
1186  """
1187  General polynomial with commutative ring coefficients.
1188  """
1189  def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds):
1190  """
1191  Initialize the polynomial.
1192 
1193  - coefficients: initializer for polynomial coefficients
1194  - coeffring: commutative ring
1195  """
1196  if coeffring is None:
1197  raise TypeError("argument `coeffring' is required")
1198  coefficients = dict(coefficients)
1199  if coefficients and coefficients.itervalues().next() not in coeffring:
1200  coefficients = [(d, coeffring.createElement(c)) for (d, c) in coefficients.iteritems()]
1201  _sorted = False
1202  kwds["coeffring"] = coeffring
1203  univar.SortedPolynomial.__init__(self, coefficients, _sorted, **kwds)
1204  OrderProvider.__init__(self, termorder.ascending_order)
1205  RingElementProvider.__init__(self)
1206  self.set_coefficient_ring(coeffring)
1207 
1208  def getRing(self):
1209  """
1210  Return an object of a subclass of Ring, to which the element
1211  belongs.
1212  """
1213  # short-cut self._ring is None case
1214  return self._ring
1215 
1217  """
1218  Return an object of a subclass of Ring, to which the all
1219  coefficients belong.
1220  """
1221  # short-cut self._coefficient_ring is None case
1222  return self._coefficient_ring
1223 
1224  def __repr__(self):
1225  return "%s(%s, %s)" % (self.__class__.__name__, repr(self.terms()), repr(self._coefficient_ring))
1226 
1227  def __add__(self, other):
1228  try:
1229  return univar.SortedPolynomial.__add__(self, other)
1230  except AttributeError:
1231  one = self.getRing().one
1232  try:
1233  return univar.SortedPolynomial.__add__(self, other * one)
1234  except Exception:
1235  return NotImplemented
1236 
1237  def __radd__(self, other):
1238  one = self.getRing().one
1239  try:
1240  return other * one + self
1241  except Exception:
1242  return NotImplemented
1243 
1244  def __sub__(self, other):
1245  try:
1246  return univar.SortedPolynomial.__sub__(self, other)
1247  except AttributeError:
1248  one = self.getRing().one
1249  try:
1250  return univar.SortedPolynomial.__sub__(self, other * one)
1251  except Exception:
1252  return NotImplemented
1253 
1254  def __rsub__(self, other):
1255  one = self.getRing().one
1256  try:
1257  return other * one - self
1258  except Exception:
1259  return NotImplemented
1260 
1261  def ismonic(self):
1262  """
1263  Return True if the polynomial is monic, i.e. its leading
1264  coefficient is one. False otherwise.
1265  """
1266  return self.leading_coefficient() == self._coefficient_ring.one
1267 
1268  def __getitem__(self, degree):
1269  """
1270  Return the coefficient of specified degree.
1271  If there is no term of degree, return 0.
1272  """
1273  result = univar.SortedPolynomial.__getitem__(self, degree)
1274  if result is 0:
1275  result = self._coefficient_ring.zero
1276  return result
1277 
1278 
1280  RingPolynomial):
1281  """
1282  Polynomial with domain coefficients.
1283  """
1284  def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds):
1285  """
1286  Initialize the polynomial.
1287 
1288  - coefficients: initializer for polynomial coefficients
1289  - coeffring: domain
1290  """
1291  if coeffring is None:
1292  raise TypeError("argument `coeffring' is required")
1293  elif not coeffring.isdomain():
1294  raise TypeError("coefficient ring has to be a domain")
1295  RingPolynomial.__init__(self, coefficients, coeffring, _sorted, **kwds)
1296  PseudoDivisionProvider.__init__(self)
1297 
1298  def discriminant(self):
1299  """
1300  Return discriminant of the polynomial.
1301  """
1302  df = self.differentiate()
1303  lc = self.leading_coefficient()
1304  if self.degree() & 3 in (2, 3):
1305  sign = -1
1306  else:
1307  sign = 1
1308  if self.getCoefficientRing().getCharacteristic() == 0:
1309  if lc != 1:
1310  return sign * self.resultant(df) / lc
1311  else:
1312  return sign * self.resultant(df)
1313  else:
1314  return sign * self.resultant(df) * lc**(self.degree() - df.degree() - 2)
1315 
1317  """
1318  Return a FieldPolynomial object obtained by embedding the
1319  polynomial ring over the domain D to over the quatient field
1320  of D.
1321  """
1322  field = self.getCoefficientRing().getQuotientField()
1323  return polynomial([(d, field.createElement(c)) for d, c in self.iterterms()], field)
1324 
1325  def __pow__(self, index, mod=None):
1326  """
1327  self ** index (% mod)
1328 
1329  It overrides the method from SortedPolynomial. The mod MUST
1330  be monic, otherwise the method raises ValueError.
1331  """
1332  if mod is None:
1333  return RingPolynomial.__pow__(self, index)
1334  elif mod.ismonic():
1335  return self.monic_pow(index, mod)
1336  else:
1337  raise ValueError("non-monic modulus")
1338 
1339 
1341  ContentProvider,
1342  DomainPolynomial):
1343  """
1344  Polynomial with unique factorization domain coefficients.
1345  """
1346  def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds):
1347  """
1348  Initialize the polynomial.
1349 
1350  - coefficients: initializer for polynomial coefficients
1351  - coeffring: unique factorization domain
1352  """
1353  if coeffring is None:
1354  raise TypeError("argument `coeffring' is required")
1355  elif not coeffring.isufd():
1356  raise TypeError("coefficient ring has to be a UFD")
1357  DomainPolynomial.__init__(self, coefficients, coeffring, _sorted, **kwds)
1358  ContentProvider.__init__(self)
1359  SubresultantGcdProvider.__init__(self)
1360 
1361 
1363  """
1364  Polynomial with integer coefficients.
1365 
1366  This class is required because special initialization must be done
1367  for built-in int/long.
1368  """
1369  def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds):
1370  dc = dict(coefficients)
1371  coefficients = [(d, rational.IntegerIfIntOrLong(c)) for (d, c) in dc.iteritems()]
1372  UniqueFactorizationDomainPolynomial.__init__(self, coefficients, coeffring, _sorted, **kwds)
1373 
1374  def normalize(self):
1375  """
1376  returns the unique normalized polynomial g
1377  which is associate to self (so g=u*self for some unit in coeffring).
1378  For IntegerPolynomial, g is positive.
1379  """
1380  if self.leading_coefficient() < 0:
1381  return -self
1382  return self
1383 
1385  ContentProvider,
1386  RingPolynomial):
1387  """
1388  Polynomial with field coefficients.
1389  """
1390  def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds):
1391  """
1392  Initialize the polynomial.
1393 
1394  - coefficients: initializer for polynomial coefficients
1395  - coeffring: field
1396  """
1397  if coeffring is None:
1398  raise TypeError("argument `coeffring' is required")
1399  elif not coeffring.isfield():
1400  raise TypeError("coefficient ring has to be a field")
1401  RingPolynomial.__init__(self, coefficients, coeffring, _sorted, **kwds)
1402  DivisionProvider.__init__(self)
1403  ContentProvider.__init__(self)
1404 
1405  def __pow__(self, index, mod=None):
1406  """
1407  self ** index (% mod).
1408  """
1409  if mod is None:
1410  return RingPolynomial.__pow__(self, index)
1411 
1412  if index < 0:
1413  raise ValueError("negative index is not allowed.")
1414  # special indices
1415  elif index == 0:
1416  return self.getRing().one
1417  elif index == 1:
1418  return self % mod
1419  elif index == 2:
1420  return self.square() % mod
1421  # special polynomials
1422  if not self:
1423  return self
1424  # general
1425  power_product = self.getRing().one
1426  power_of_2 = self % mod
1427  while index:
1428  if index & 1:
1429  power_product = mod.mod(power_product * power_of_2)
1430  index //= 2
1431  if index:
1432  power_of_2 = mod.mod(power_of_2.square())
1433  return power_product
1434 
1435  def resultant(self, other):
1436  """
1437  Return the resultant of self and other.
1438  """
1439  order = termorder.ascending_order
1440  f, g = self, other
1441  negative = False
1442  if order.degree(f) < order.degree(g):
1443  f, g = g, f
1444  if (order.degree(f) & 1) and (order.degree(g) & 1):
1445  negative = not negative
1446  coeff_ring = self.getCoefficientRing()
1447  a = b = coeff_ring.one
1448  while order.degree(g) > 0:
1449  delta = order.degree(f) - order.degree(g)
1450  if (order.degree(f) & 1) and (order.degree(g) & 1):
1451  negative = not negative
1452  h = f % g
1453  h *= order.leading_coefficient(g)**(order.degree(f) - order.degree(g) + 1)
1454  f, g = g, h.scalar_exact_division(a * (b ** delta))
1455  a = order.leading_coefficient(f)
1456  b = ((a ** delta) * b) // (b ** delta)
1457  if not g:
1458  return coeff_ring.zero
1459  scalar = g[0]
1460  degree_f = order.degree(f)
1461  if negative:
1462  return -(b * scalar ** degree_f) // (b ** degree_f)
1463  else:
1464  return (b * scalar ** degree_f) // (b ** degree_f)
1465 
1466  def discriminant(self):
1467  """
1468  Return discriminant of the polynomial.
1469  """
1470  df = self.differentiate()
1471  lc = self.leading_coefficient()
1472  if self.degree() & 3 in (2, 3):
1473  sign = -1
1474  else:
1475  sign = 1
1476  if self.getCoefficientRing().getCharacteristic() == 0:
1477  return sign * self.resultant(df) / lc
1478  else:
1479  return sign * self.resultant(df) * lc**(self.degree() - df.degree() - 2)
1480 
1481 
1483  FieldPolynomial):
1484  """
1485  Fq polynomial
1486  """
1487  def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds):
1488  """
1489  Initialize the polynomial.
1490 
1491  - coefficients: initializer for polynomial coefficients
1492  - coeffring: finite prime field
1493  """
1494  if coeffring is None:
1495  raise TypeError("argument `coeffring' is required")
1496  FieldPolynomial.__init__(self, coefficients, coeffring, _sorted, **kwds)
1497  PrimeCharacteristicFunctionsProvider.__init__(self, coeffring.getCharacteristic())
1498 
1499 
1501  """
1502  Fp polynomial
1503  """
1504  def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds):
1505  """
1506  Initialize the polynomial.
1507 
1508  - coefficients: initializer for polynomial coefficients
1509  - coeffring: finite prime field
1510  """
1511  FiniteFieldPolynomial.__init__(self, coefficients, coeffring, _sorted, **kwds)
1512 
1513 
1514  # IMPLEMENTATION REMARK:
1515  # The most time consuming part of computation is bunch of
1516  # object creations. Thus, here, the creation of finite field
1517  # elements is avoided during summing up coefficients.
1518 
1519  def ring_mul(self, other):
1520  """
1521  Multiplication of two polynomials in the same ring.
1522  """
1523  try:
1524  mul_coeff = {}
1525 
1526  stripped_self = [(ds, cs.n) for (ds, cs) in self]
1527  stripped_other = [(do, co.n) for (do, co) in other]
1528  for ds, cs in stripped_self:
1529  for do, co in stripped_other:
1530  term_degree = ds + do
1531  if term_degree in mul_coeff:
1532  mul_coeff[term_degree] = (mul_coeff[term_degree] + cs*co ) % self.ch
1533  else:
1534  mul_coeff[term_degree] = cs*co
1535 
1536  return self.construct_with_default(mul_coeff)
1537  except AttributeError:
1538  # maybe fail due to lacking attribute .n sometimes
1539  _log.debug("fall back to good old ring_mul")
1540  return univar.PolynomialInterface.ring_mul(self, other)
1541 
1542  def square(self):
1543  """
1544  Return the square of self.
1545  """
1546  # zero
1547  if not self:
1548  return self
1549  # monomial, binomial
1550  elif len(self.sorted) <= 2:
1551  return FieldPolynomial.square(self)
1552 
1553  # general
1554  try:
1555  result = {}
1556  # stripped to bare integer
1557  stripped_self = [(ds, cs.n) for (ds, cs) in self]
1558  for i, term in enumerate(stripped_self):
1559  di, ci = term[0], term[1]
1560  result[di*2] = result.get(di*2, 0) + pow(ci, 2, self.ch)
1561  for j in range(i + 1, len(stripped_self)):
1562  dj, cj = self.sorted[j][0], self.sorted[j][1]
1563  result[di + dj] = result.get(di + dj, 0) + 2 * ci * cj
1564  # back to decorated datatype automatically
1565  return self.construct_with_default(result)
1566  except AttributeError:
1567  _log.debug("fall back to old square")
1568  return FieldPolynomial.square(self)
1569 
1570 
1571 def inject_variable(polynom, variable):
1572  """
1573  Inject variable into polynom temporarily. The variable name will
1574  be lost after any arithmetic operations on polynom, though the
1575  class name of polynom will remain prefixed with 'Var'. If one need
1576  variable name permanently, he/she should define a class inheriting
1577  VariableProvider.
1578  """
1579  baseclasses = polynom.__class__.__bases__
1580  if VariableProvider not in baseclasses:
1581  polynom.__class__ = type("Var" + polynom.__class__.__name__,
1582  (polynom.__class__, VariableProvider,),
1583  {})
1584  polynom.variable = variable
1585 
1586 
1587 special_ring_table = {rational.IntegerRing: IntegerPolynomial}
1588 
1589 
1590 def polynomial(coefficients, coeffring):
1591  """
1592  Return a polynomial.
1593  - coefficients has to be a initializer for dict, whose keys are
1594  degrees and values are coefficients at degrees.
1595  - coeffring has to be an object inheriting ring.Ring.
1596 
1597  One can override the way to choose a polynomial type from a
1598  coefficient ring, by setting:
1599  special_ring_table[coeffring_type] = polynomial_type
1600  before the function call.
1601  """
1602  if type(coeffring) in special_ring_table:
1603  poly_type = special_ring_table[type(coeffring)]
1604  elif coeffring.isfield():
1605  poly_type = FieldPolynomial
1606  elif coeffring.isufd():
1607  poly_type = UniqueFactorizationDomainPolynomial
1608  elif coeffring.isdomain():
1609  poly_type = DomainPolynomial
1610  else:
1611  poly_type = RingPolynomial
1612  return poly_type(coefficients, coeffring)
1613 
1614 
1615 def OneVariableDensePolynomial(coefficient, variable, coeffring=None):
1616  """
1617  OneVariableDensePolynomial(coefficient, variable [,coeffring])
1618 
1619  - coefficient has to be a sequence of coefficients in ascending order
1620  of degree.
1621  - variable has to be a character string.
1622  - coeffring has to be, if specified, an object inheriting ring.Ring.
1623 
1624  This function is provided for backward compatible way of defining
1625  univariate polynomial. The argument variable is ignored.
1626  """
1627  _coefficients = dict(enumerate(coefficient))
1628  if coeffring is None:
1629  coeffring = init_coefficient_ring(_coefficients)
1630  return polynomial(_coefficients, coeffring)
1631 
1632 def OneVariableSparsePolynomial(coefficient, variable, coeffring=None):
1633  """
1634  OneVariableSparsePolynomial(coefficient, variable [,coeffring])
1635 
1636  - coefficient has to be a dictionary of degree-coefficient pairs.
1637  - variable has to be a character string.
1638  - coeffring has to be, if specified, an object inheriting ring.Ring.
1639 
1640  This function is provided for backward compatible way of defining
1641  univariate polynomial. The argument variable is ignored.
1642  """
1643  _coefficients = dict(coefficient)
1644  if coeffring is None:
1645  coeffring = init_coefficient_ring(_coefficients)
1646  return polynomial(_coefficients, coeffring)
1647 
1648 def init_coefficient_ring(coefficients):
1649  """
1650  Return a ring to which all coefficients belong. The argument
1651  coefficients is a dictionary whose values are the coefficients.
1652  """
1653  myRing = None
1654  for c in coefficients.itervalues():
1655  cring = ring.getRing(c)
1656  if not myRing or myRing != cring and myRing.issubring(cring):
1657  myRing = cring
1658  elif not cring.issubring(myRing):
1659  myRing = myRing.getCommonSuperring(cring)
1660  return myRing
nzmath.poly.uniutil.FiniteFieldPolynomial.__init__
def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds)
Definition: uniutil.py:1487
nzmath.poly.uniutil.RingElementProvider.getRing
def getRing(self)
Definition: uniutil.py:1158
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider
Definition: uniutil.py:732
nzmath.poly.univar.SortedPolynomial.degree
def degree(self)
Definition: univar.py:710
nzmath.bigrandom
Definition: bigrandom.py:1
nzmath.poly.uniutil.RingElementProvider.__init__
def __init__(self)
Definition: uniutil.py:1146
nzmath.poly.uniutil.KaratsubaProvider.square_karatsuba
def square_karatsuba(self)
Definition: uniutil.py:1072
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.ring
Definition: ring.py:1
nzmath.poly.uniutil.PseudoDivisionProvider
Definition: uniutil.py:338
nzmath.poly.uniutil.DivisionProvider
Definition: uniutil.py:78
nzmath.poly.uniutil.DomainPolynomial.to_field_polynomial
def to_field_polynomial(self)
Definition: uniutil.py:1316
nzmath.poly.uniutil.RingPolynomial.ismonic
def ismonic(self)
Definition: uniutil.py:1261
nzmath.poly.uniutil.DivisionProvider.scalar_exact_division
def scalar_exact_division(self, scale)
Definition: uniutil.py:293
nzmath.poly.uniutil.DivisionProvider.__truediv__
def __truediv__(self, other)
Definition: uniutil.py:285
nzmath.poly.uniutil.VariableProvider
Definition: uniutil.py:1115
nzmath.poly.uniutil.FinitePrimeFieldPolynomial.__init__
def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds)
Definition: uniutil.py:1504
nzmath.poly.uniutil.SubresultantGcdProvider.subresultant_extgcd
def subresultant_extgcd(self, other)
Definition: uniutil.py:700
nzmath.poly.uniutil.UniqueFactorizationDomainPolynomial.__init__
def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds)
Definition: uniutil.py:1346
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider._small_index_mod_pow
def _small_index_mod_pow(self, polynom, index)
Definition: uniutil.py:821
nzmath.poly.uniutil.RingElementProvider.set_coefficient_ring
def set_coefficient_ring(self, coeffring)
Definition: uniutil.py:1176
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.squarefree_decomposition
def squarefree_decomposition(self)
Definition: uniutil.py:840
nzmath.poly.univar.PolynomialInterface.differentiate
def differentiate(self)
Definition: univar.py:90
nzmath.poly.uniutil.OrderProvider.shift_degree_to
def shift_degree_to(self, degree)
Definition: uniutil.py:43
nzmath.poly.uniutil.FieldPolynomial.resultant
def resultant(self, other)
Definition: uniutil.py:1435
nzmath.poly.uniutil.SubresultantGcdProvider
Definition: uniutil.py:616
nzmath.poly.termorder
Definition: termorder.py:1
nzmath.poly.ring
Definition: ring.py:1
nzmath.poly.uniutil.DivisionProvider.gcd
def gcd(self, other)
Definition: uniutil.py:300
nzmath.poly.univar
Definition: univar.py:1
nzmath.poly.uniutil.VariableProvider.getVariableList
def getVariableList(self)
Definition: uniutil.py:1135
nzmath.gcd
Definition: gcd.py:1
nzmath.poly.uniutil.RingPolynomial.__getitem__
def __getitem__(self, degree)
Definition: uniutil.py:1268
nzmath.poly.uniutil.DomainPolynomial.discriminant
def discriminant(self)
Definition: uniutil.py:1298
nzmath.poly.uniutil.OrderProvider
Definition: uniutil.py:28
nzmath.poly.uniutil.KaratsubaProvider.ring_mul_karatsuba
def ring_mul_karatsuba(self, other)
Definition: uniutil.py:1025
nzmath.poly.uniutil.DivisionProvider._populate_reduced_more
def _populate_reduced_more(self, degrees)
Definition: uniutil.py:223
nzmath.poly.uniutil.DomainPolynomial.__pow__
def __pow__(self, index, mod=None)
Definition: uniutil.py:1325
nzmath.poly.uniutil.DivisionProvider.mod
def mod(self, dividend)
Definition: uniutil.py:143
nzmath.poly.formalsum.FormalSumContainerInterface.terms
def terms(self)
Definition: formalsum.py:192
nzmath.poly.univar.PolynomialInterface.square
def square(self)
Definition: univar.py:84
nzmath.rational
Definition: rational.py:1
nzmath.poly.uniutil.FiniteFieldPolynomial
Definition: uniutil.py:1483
nzmath.poly.uniutil.VariableProvider.getVariable
def getVariable(self)
Definition: uniutil.py:1129
nzmath.poly.uniutil.DivisionProvider.__floordiv__
def __floordiv__(self, other)
Definition: uniutil.py:114
nzmath.poly.uniutil.RingPolynomial.__rsub__
def __rsub__(self, other)
Definition: uniutil.py:1254
nzmath.poly.uniutil.RingPolynomial.__init__
def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds)
Definition: uniutil.py:1189
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.distinct_degree_factorization
def distinct_degree_factorization(self)
Definition: uniutil.py:884
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.pthroot
def pthroot(self)
Definition: uniutil.py:876
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.split_same_degrees
def split_same_degrees(self, degree)
Definition: uniutil.py:917
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.isirreducible
def isirreducible(self)
Definition: uniutil.py:992
nzmath.poly.uniutil.FinitePrimeFieldPolynomial
Definition: uniutil.py:1500
nzmath.poly.uniutil.PseudoDivisionProvider.scalar_exact_division
def scalar_exact_division(self, scale)
Definition: uniutil.py:458
nzmath.poly.uniutil.PseudoDivisionProvider.pseudo_mod
def pseudo_mod(self, other)
Definition: uniutil.py:406
nzmath.poly.uniutil.ContentProvider.content
def content(self)
Definition: uniutil.py:581
nzmath.poly.uniutil.PseudoDivisionProvider.pseudo_floordiv
def pseudo_floordiv(self, other)
Definition: uniutil.py:377
nzmath.poly.univar.SortedPolynomial.leading_coefficient
def leading_coefficient(self)
Definition: univar.py:720
nzmath.poly.uniutil.FieldPolynomial.discriminant
def discriminant(self)
Definition: uniutil.py:1466
nzmath.poly.uniutil.KaratsubaProvider
Definition: uniutil.py:1019
nzmath.poly.uniutil.RingPolynomial.__repr__
def __repr__(self)
Definition: uniutil.py:1224
nzmath.poly.uniutil.ContentProvider
Definition: uniutil.py:575
nzmath.ring.exact_division
def exact_division(self, other)
Definition: ring.py:852
nzmath.rational.IntegerRing
Definition: rational.py:789
nzmath.poly.uniutil.VariableProvider.variable
variable
Definition: uniutil.py:1127
nzmath.poly.uniutil.RingPolynomial.__sub__
def __sub__(self, other)
Definition: uniutil.py:1244
nzmath.poly.uniutil.PseudoDivisionProvider.monic_mod
def monic_mod(self, other)
Definition: uniutil.py:520
nzmath.poly.uniutil.DivisionProvider.__divmod__
def __divmod__(self, other)
Definition: uniutil.py:97
nzmath.poly.uniutil.IntegerPolynomial.normalize
def normalize(self)
Definition: uniutil.py:1374
nzmath.ring.QuotientField
Definition: ring.py:234
nzmath.poly.ratfunc
Definition: ratfunc.py:1
nzmath.poly.uniutil.FieldPolynomial.__init__
def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds)
Definition: uniutil.py:1390
nzmath.poly.uniutil.RingPolynomial.__add__
def __add__(self, other)
Definition: uniutil.py:1227
nzmath.poly.uniutil.OneVariableSparsePolynomial
def OneVariableSparsePolynomial(coefficient, variable, coeffring=None)
Definition: uniutil.py:1632
nzmath.compatibility.cmp
cmp
Definition: compatibility.py:20
nzmath.arith1
Definition: arith1.py:1
nzmath.poly.uniutil.PseudoDivisionProvider.__truediv__
def __truediv__(self, other)
Definition: uniutil.py:434
nzmath.ring.CommutativeRingElement
Definition: ring.py:291
nzmath.poly.uniutil.DivisionProvider.__init__
def __init__(self)
Definition: uniutil.py:88
nzmath.poly.uniutil.polynomial
def polynomial(coefficients, coeffring)
Definition: uniutil.py:1590
nzmath.poly.uniutil.DivisionProvider._populate_reduced
def _populate_reduced(self, degree, lc, upperbound)
Definition: uniutil.py:202
nzmath.poly.uniutil.PseudoDivisionProvider.pseudo_divmod
def pseudo_divmod(self, other)
Definition: uniutil.py:346
nzmath.poly.uniutil.FinitePrimeFieldPolynomial.square
def square(self)
Definition: uniutil.py:1542
nzmath.poly.formalsum.FormalSumContainerInterface.iterterms
def iterterms(self)
Definition: formalsum.py:174
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.__pow__
def __pow__(self, index, mod=None)
Definition: uniutil.py:751
nzmath.poly.uniutil.DivisionProvider._reduced
_reduced
Definition: uniutil.py:95
nzmath.compatibility
Definition: compatibility.py:1
nzmath.poly.uniutil.PseudoDivisionProvider.monic_divmod
def monic_divmod(self, other)
Definition: uniutil.py:465
nzmath.poly.uniutil.PseudoDivisionProvider.monic_pow
def monic_pow(self, index, mod)
Definition: uniutil.py:545
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.ch
ch
Definition: uniutil.py:749
nzmath.poly.uniutil.FieldPolynomial.__pow__
def __pow__(self, index, mod=None)
Definition: uniutil.py:1405
nzmath.poly.formalsum.FormalSumContainerInterface.construct_with_default
def construct_with_default(self, maindata)
Definition: formalsum.py:236
nzmath.poly.uniutil.IntegerPolynomial
Definition: uniutil.py:1362
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.mod_pow
def mod_pow(self, polynom, index)
Definition: uniutil.py:794
nzmath.poly.uniutil.DomainPolynomial
Definition: uniutil.py:1280
nzmath.poly.univar.SortedPolynomial
Definition: univar.py:360
nzmath.poly.uniutil.SubresultantGcdProvider.subresultant_gcd
def subresultant_gcd(self, other)
Definition: uniutil.py:653
nzmath.poly.univar.SortedPolynomial.sorted
sorted
Definition: univar.py:374
nzmath.poly.uniutil.IntegerPolynomial.__init__
def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds)
Definition: uniutil.py:1369
nzmath.poly.uniutil.RingPolynomial.__radd__
def __radd__(self, other)
Definition: uniutil.py:1237
nzmath.poly.uniutil.OneVariableDensePolynomial
def OneVariableDensePolynomial(coefficient, variable, coeffring=None)
Definition: uniutil.py:1615
nzmath.poly.uniutil.FinitePrimeFieldPolynomial.ring_mul
def ring_mul(self, other)
Definition: uniutil.py:1519
nzmath.poly.uniutil.DivisionProvider.mod_pow
def mod_pow(self, polynom, index)
Definition: uniutil.py:179
nzmath.poly.uniutil.DivisionProvider.__mod__
def __mod__(self, other)
Definition: uniutil.py:127
nzmath.poly.uniutil.RingElementProvider._ring
_ring
Definition: uniutil.py:1156
nzmath.poly.uniutil.UniqueFactorizationDomainPolynomial
Definition: uniutil.py:1342
nzmath.poly.ratfunc.RationalFunction
Definition: ratfunc.py:10
nzmath.poly.uniutil.FieldPolynomial
Definition: uniutil.py:1386
nzmath.poly.uniutil.init_coefficient_ring
def init_coefficient_ring(coefficients)
Definition: uniutil.py:1648
nzmath.poly.uniutil.ContentProvider.primitive_part
def primitive_part(self)
Definition: uniutil.py:609
nzmath.poly.uniutil.OrderProvider.order
order
Definition: uniutil.py:41
nzmath.poly.uniutil.PseudoDivisionProvider.exact_division
def exact_division(self, other)
Definition: uniutil.py:443
nzmath.poly.uniutil.RingPolynomial.getRing
def getRing(self)
Definition: uniutil.py:1208
nzmath.poly.uniutil.OrderProvider.split_at
def split_at(self, degree)
Definition: uniutil.py:62
nzmath.poly.uniutil.RingPolynomial.getCoefficientRing
def getCoefficientRing(self)
Definition: uniutil.py:1216
nzmath.poly.uniutil.RingElementProvider
Definition: uniutil.py:1142
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.__init__
def __init__(self, ch)
Definition: uniutil.py:742
nzmath.poly.uniutil.RingPolynomial
Definition: uniutil.py:1185
nzmath.poly.uniutil.DomainPolynomial.__init__
def __init__(self, coefficients, coeffring=None, _sorted=False, **kwds)
Definition: uniutil.py:1284
nzmath.poly.uniutil.PrimeCharacteristicFunctionsProvider.factor
def factor(self)
Definition: uniutil.py:969
nzmath.poly.uniutil.RingElementProvider._coefficient_ring
_coefficient_ring
Definition: uniutil.py:1155
nzmath.poly.uniutil.OrderProvider.__init__
def __init__(self, order=termorder.ascending_order)
Definition: uniutil.py:32
nzmath.poly.uniutil.SubresultantGcdProvider.resultant
def resultant(self, other)
Definition: uniutil.py:623
nzmath.poly.uniutil.DivisionProvider.extgcd
def extgcd(self, other)
Definition: uniutil.py:316
nzmath.compatibility.card
def card(virtualset)
Definition: compatibility.py:28
nzmath.poly.uniutil.VariableProvider.__init__
def __init__(self, varname)
Definition: uniutil.py:1120
nzmath.poly.uniutil.PseudoDivisionProvider.monic_floordiv
def monic_floordiv(self, other)
Definition: uniutil.py:492
nzmath.poly.uniutil.inject_variable
def inject_variable(polynom, variable)
Definition: uniutil.py:1571