"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/poly/uniutil.py" (19 Nov 2012, 58772 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "uniutil.py" see the Fossies "Dox" file reference documentation.

    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
   19 import nzmath.compatibility
   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 
  338 class PseudoDivisionProvider(object):
  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 
  616 class SubresultantGcdProvider(object):
  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 
  732 class PrimeCharacteristicFunctionsProvider(object):
  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 
  840     def squarefree_decomposition(self):
  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 
  884     def distinct_degree_factorization(self):
  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 
 1142 class RingElementProvider(ring.CommutativeRingElement):
 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 
 1183 class RingPolynomial(OrderProvider,
 1184                      univar.SortedPolynomial,
 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 
 1216     def getCoefficientRing(self):
 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 
 1279 class DomainPolynomial(PseudoDivisionProvider,
 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 
 1316     def to_field_polynomial(self):
 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 
 1340 class UniqueFactorizationDomainPolynomial(SubresultantGcdProvider,
 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 
 1362 class IntegerPolynomial(UniqueFactorizationDomainPolynomial):
 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 
 1384 class FieldPolynomial(DivisionProvider,
 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 
 1482 class FiniteFieldPolynomial(PrimeCharacteristicFunctionsProvider,
 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 
 1500 class FinitePrimeFieldPolynomial(FiniteFieldPolynomial):
 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             ## stripped to bare integer
 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             ## back to decorated datatype automatically
 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