"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/algfield.py" (19 Nov 2012, 27554 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 "algfield.py" see the Fossies "Dox" file reference documentation.

    1 from __future__ import division
    2 
    3 import math
    4 import nzmath.arith1 as arith1
    5 import nzmath.equation as equation
    6 import nzmath.gcd as gcd
    7 import nzmath.lattice as lattice
    8 import nzmath.vector as vector
    9 import nzmath.matrix as matrix
   10 import nzmath.factor.misc as misc
   11 import nzmath.finitefield as finitefield
   12 import nzmath.poly.uniutil as uniutil
   13 import nzmath.poly.multiutil as multiutil
   14 import nzmath.rational as rational
   15 import nzmath.ring as ring
   16 import nzmath.squarefree as squarefree
   17 import nzmath.round2 as round2
   18 
   19 
   20 class NumberField (ring.Field):
   21     """
   22     A class of number field.
   23     """
   24     def __init__(self, polynomial, precompute=False):
   25         """
   26         Initialize a number field with given polynomial coefficients
   27         (in ascending order).
   28         """
   29         ring.Field.__init__(self)
   30         self.polynomial = polynomial
   31         self.degree = len(polynomial) - 1
   32         if precompute:
   33             getConj()
   34             disc()
   35             signature()
   36             integer_ring()
   37 
   38     def __repr__(self):
   39         return_str = '%s(%s)' % (self.__class__.__name__, self.polynomial)
   40         return return_str
   41 
   42     def __mul__(self, other):
   43         """
   44         Output composite field of self and other.
   45         """
   46         common_options = {"coeffring": rational.theIntegerRing,
   47                           "number_of_variables": 2}
   48         flist = [((d, 0), c) for (d, c) in enumerate(self.polynomial)]
   49         f = multiutil.polynomial(flist, **common_options)
   50         g = zpoly(other.polynomial)
   51         diff = multiutil.polynomial({(1, 0): 1, (0, 1):-1}, **common_options)
   52         compos = f.resultant(g(diff), 0)
   53         return NumberField([compos[i] for i in range(compos.degree() + 1)])
   54 
   55     def getConj(self):
   56         """
   57         Return (approximate) solutions of self.polynomial.
   58         We can discriminate the conjugate field of self by these values.
   59         """
   60         if not hasattr(self, "conj"):
   61             conj = equation.SimMethod(self.polynomial)
   62             self.conj = conj
   63         return self.conj
   64 
   65     def disc(self):
   66         """
   67         Compute the discriminant of self.polynomial.
   68         (The output is not field disc of self but disc of self.polynomial.)
   69         """
   70         """
   71         Obasis, disc = round2.round2(self)
   72         A = [algfield.MatAlgNumber(obasis[0][i]) for i in range(degree)]
   73         return disc(A) <--- real disc of self field.
   74         """
   75 
   76         if not hasattr(self, "poldisc"):
   77             degree = self.degree
   78             traces = []
   79             for i in range(degree):
   80                 for j in range(degree):
   81                     s = self.basis(i)*self.basis(j)
   82                     traces.append(s.trace())
   83             M = matrix.RingSquareMatrix(degree, degree, traces)
   84             self.poldisc = M.determinant()
   85         return self.poldisc
   86 
   87     def integer_ring(self):
   88         """
   89         Return the integer ring of self (using round2 module)
   90         """
   91         if not hasattr(self, "integer_basis"):
   92             int_basis, self.discriminant = round2.round2(self.polynomial)
   93             self.integer_basis = matrix.createMatrix(
   94                 self.degree, [vector.Vector(int_basis[j]) for j in range(
   95                 len(int_basis))])
   96         return self.integer_basis
   97 
   98     def field_discriminant(self):
   99         """
  100         Return the (field) discriminant of self (using round2 module)
  101         """
  102         if not hasattr(self, "discriminant"):
  103             int_basis, self.discriminant = round2.round2(self.polynomial)
  104             self.integer_basis = matrix.createMatrix(
  105                 self.degree, [vector.Vector(int_basis[j]) for j in range(
  106                 len(int_basis))])
  107         return self.discriminant
  108 
  109     def basis(self, j):
  110         """
  111         Return the j-th basis of self.
  112         """
  113         base = [0] * self.degree
  114         base[j] = 1
  115         return BasicAlgNumber([base, 1], self.polynomial)
  116 
  117     def signature(self):
  118         """
  119         Using Strum's algorithm, compute the signature of self.
  120         Algorithm 4.1.11 in Cohen's Book
  121         """
  122         if not hasattr(self, "sign"):
  123             degree = self.degree
  124             #Step 1.
  125             if degree == 0:
  126                 return (0, 0)
  127             # no check for degree 1?
  128 
  129             minpoly = zpoly(self.polynomial)
  130             d_minpoly = minpoly.differentiate()
  131             A = minpoly.primitive_part()
  132             B = d_minpoly.primitive_part()
  133             g = 1
  134             h = 1
  135             pos_at_inf = A.leading_coefficient() > 0
  136             pos_at_neg = pos_at_inf == (degree & 1)
  137             r_1 = 1
  138 
  139             #Step 2.
  140             while True:
  141                 deg = A.degree() - B.degree()
  142                 residue = A.pseudo_mod(B)
  143                 if not residue:
  144                     raise ValueError("not squarefree")
  145                 if B.leading_coefficient() > 0 or deg & 1:
  146                     residue = - residue
  147                 #Step 3.
  148                 degree_res = residue.degree()
  149                 pos_at_inf_of_res = residue.leading_coefficient() > 0
  150 
  151                 if pos_at_inf_of_res != pos_at_inf:
  152                     pos_at_inf = not pos_at_inf
  153                     r_1 -= 1
  154 
  155                 if pos_at_inf_of_res != (pos_at_neg == (degree_res & 1 == 0)):
  156                     pos_at_neg = not pos_at_neg
  157                     r_1 += 1
  158 
  159                 #Step 4.
  160                 if degree_res == 0:
  161                     return (r_1, (degree - r_1) >> 1)
  162 
  163                 A, B = B, residue.scalar_exact_division(g*(h**deg))
  164                 g = abs(A.leading_coefficient())
  165                 if deg == 1:
  166                     h = g
  167                 elif deg > 1:
  168                     h = g**deg // h**(deg - 1)
  169         return self.sign
  170 
  171     def POLRED(self):
  172         """
  173         Given a polynomial f i.e. a field self, output some polynomials
  174         defining subfield of self, where self is a field defined by f.
  175         Algorithm 4.4.11 in Cohen's book.
  176         """
  177         n = self.degree
  178         appr = self.getConj()
  179 
  180         #Step 1.
  181         Basis = self.integer_ring()
  182         BaseList = []
  183         for i in range(n):
  184             AlgInt = MatAlgNumber(Basis[i].compo, self.polynomial)
  185             BaseList.append(AlgInt)
  186 
  187         #Step 2.
  188         traces = []
  189         if self.signature()[1] == 0:
  190             for i in range(n):
  191                 for j in range(n):
  192                     s = BaseList[i]*BaseList[j]
  193                     traces.append(s.trace())
  194         else:
  195             sigma = self.getConj()
  196             f = []
  197             for i in range(n):
  198                 f.append(zpoly(Basis[i]))
  199             for i in range(n):
  200                 for j in range(n):
  201                     m = 0 + 0j
  202                     for k in range(n):
  203                         m += f[i](sigma[k])*f[j](sigma[k].conjugate())
  204                     traces.append(m.real)
  205 
  206         #Step 3.
  207         M = matrix.createMatrix(n, n, traces)
  208         S = matrix.unitMatrix(n)
  209         L = lattice.LLL(S, M)[0]
  210 
  211 
  212         #Step 4.
  213         Ch_Basis = []
  214         for i in range(n):
  215             base_cor = changetype(0, self.polynomial).ch_matrix()
  216             for v in range(n):
  217                 base_cor += BaseList[v]*L.compo[v][i]
  218             Ch_Basis.append(base_cor)
  219 
  220         C = []
  221         a = Ch_Basis[0]
  222         for i in range(n):
  223             coeff = Ch_Basis[i].ch_basic().getCharPoly()
  224             C.append(zpoly(coeff))
  225 
  226         #Step 5.
  227         P = []
  228         for i in range(n):
  229             diff_C = C[i].differentiate()
  230             gcd_C = C[i].subresultant_gcd(diff_C)
  231             P.append(C[i].exact_division(gcd_C))
  232 
  233         return P
  234 
  235     def isIntBasis(self):
  236         """
  237         Determine whether self.basis is integral basis of self field.
  238         """
  239         D = self.disc()
  240         if squarefree.trial_division(abs(D)):
  241             return True
  242         else:
  243             if D & 3 == 0:
  244                 if squarefree.trial_division(abs(D) >> 2) and (D >> 2) & 3 != 1:
  245                     return True
  246         # compare with real integer ring
  247         return abs(self.integer_ring().determinant()) == 1
  248 
  249     def isGaloisField(self, other = None):
  250         """
  251         Determine whether self/other is Galois field.
  252         """
  253         if self.signature[0] == 0 or self.signature[1] == 0:
  254             return "Can not determined"
  255         else:
  256             return False
  257 
  258     def isFieldElement(self, A):
  259         """
  260         Determine whether A is field element of self field or not.
  261         """
  262         poly = A.polynomial
  263         if poly == self.polynomial:
  264             return True
  265         else:
  266             if poly == self.POLRED():
  267                 return True
  268             else:
  269                 return False
  270 
  271     def getCharacteristic(self):
  272         """
  273         Return characteristic of the field (it is always zero).
  274         """
  275         return 0
  276 
  277     def createElement(self, *seed):
  278         """
  279         createElement returns an element of the field with seed.
  280         """
  281         if len(seed) == 1:
  282             if isinstance(seed[0], list):
  283                 if isinstance(seed[0][0], list):
  284                     return BasicAlgNumber(seed[0], self.polynomial)
  285                 else:
  286                     return MatAlgNumber(seed[0], self.polynomial)
  287             else:
  288                 raise NotImplementedError
  289         else:
  290             raise NotImplementedError
  291 
  292     def issubring(self, other):
  293         """
  294         Report whether another ring contains the field as a subring.
  295         """
  296         if self is other or self.isSubField(other):
  297             return True
  298         raise NotImplementedError("don't know how to tell")
  299 
  300     def issuperring(self, other):
  301         """
  302         Report whether the field is a superring of another ring.
  303         """
  304         if self is other:
  305             return True
  306         elif other.issubring(rational.theRationalField):
  307             return True
  308         raise NotImplementedError("don't know how to tell")
  309 
  310     def __eq__(self, other):
  311         """
  312         Equality test.
  313         """
  314         if self.issubring(other) and self.issuperring(other):
  315             return True
  316         return False
  317 
  318     def __hash__(self):
  319         raise NotImplementedError
  320 
  321     # properties
  322     def _getOne(self):
  323         "getter for one"
  324         if self._one is None:
  325             self._one = BasicAlgNumber([[1] + [0] * (self.degree - 1), 1],
  326                               self.polynomial)
  327         return self._one
  328 
  329     one = property(_getOne, None, None, "multiplicative unit.")
  330 
  331     def _getZero(self):
  332         "getter for zero"
  333         if self._zero is None:
  334             self._zero = BasicAlgNumber([[0] + [0] * (self.degree - 1), 1],
  335                               self.polynomial)
  336         return self._zero
  337 
  338     zero = property(_getZero, None, None, "additive unit.")
  339 
  340 
  341 class BasicAlgNumber(object):
  342     """
  343     The class for algebraic number.
  344     """
  345     def __init__(self, valuelist, polynomial, precompute=False):
  346         if len(polynomial) != len(valuelist[0])+1:
  347             raise ValueError
  348         self.value = valuelist
  349         self.coeff = valuelist[0]
  350         self.denom = valuelist[1]
  351         self.degree = len(polynomial) - 1
  352         self.polynomial = polynomial
  353         self.field = NumberField(self.polynomial)
  354         Gcd = gcd.gcd_of_list(self.coeff)
  355         GCD = gcd.gcd(Gcd[0], self.denom)
  356         if GCD != 1:
  357             self.coeff = [i//GCD for i in self.coeff]
  358             self.denom = self.denom//GCD
  359         if precompute:
  360             self.getConj()
  361             self.getApprox()
  362             self.getCharPoly()
  363 
  364     def __repr__(self):
  365         return_str = '%s(%s, %s)' % (self.__class__.__name__, [self.coeff, self.denom], self.polynomial)
  366         return return_str
  367 
  368     def __neg__(self):
  369         coeff = []
  370         for i in range(len(self.coeff)):
  371             coeff.append(-self.coeff[i])
  372         return BasicAlgNumber([coeff, self.denom], self.polynomial)
  373 
  374     def _int_to_algnumber(self, other):
  375         """
  376         other (integer) -> BasicAlgNumber (over self.field)
  377         """
  378         return BasicAlgNumber([[other] + [0] * (self.degree - 1), 1],
  379                               self.polynomial)
  380 
  381     def _rational_to_algnumber(self, other):
  382         """
  383         other (rational) -> BasicAlgNumber (over self.field)
  384         """
  385         return BasicAlgNumber([[other.numerator] + [0] * (self.degree - 1),
  386                                other.denominator], self.polynomial)
  387 
  388     def __add__(self, other):
  389         if isinstance(other, BasicAlgNumber):
  390             d = self.denom*other.denom
  391             coeff = []
  392             for i in range(len(self.coeff)):
  393                 coeff.append(
  394                     other.denom*self.coeff[i] + self.denom*other.coeff[i])
  395             return BasicAlgNumber([coeff, d], self.polynomial)
  396         elif isinstance(other, (int, long)):
  397             return self + self._int_to_algnumber(other)
  398         elif isinstance(other, rational.Rational):
  399             return self + self._rational_to_algnumber(other)
  400         else:
  401             return NotImplemented
  402 
  403     __radd__ = __add__
  404 
  405     def __sub__(self, other):
  406         return self.__add__(-other)
  407 
  408     def __mul__(self, other):
  409         if isinstance(other, BasicAlgNumber):
  410             d = self.denom*other.denom
  411             f = zpoly(self.polynomial)
  412             g = zpoly(self.coeff)
  413             h = zpoly(other.coeff)
  414             j = (g * h).pseudo_mod(f)
  415             jcoeff = [j[i] for i in range(self.degree)]
  416             return BasicAlgNumber([jcoeff, d], self.polynomial)
  417         elif isinstance(other, (int, long)):
  418             Coeff = [i*other for i in self.coeff]
  419             return BasicAlgNumber([Coeff, self.denom], self.polynomial)
  420         elif isinstance(other, rational.Rational):
  421             GCD = gcd.gcd(other.numerator, self.denom)
  422             denom = self.denom * other.denominator
  423             mul = other.numerator
  424             if GCD != 1:
  425                 denom //= GCD
  426                 mul //= GCD
  427             Coeff = [self.coeff[j] * mul for j in range(self.degree)]
  428             return BasicAlgNumber([Coeff, denom], self.polynomial)
  429         else:
  430             return NotImplemented
  431 
  432     __rmul__ = __mul__
  433 
  434     def __pow__(self, exponent, mod=None):
  435         d = self.denom**exponent
  436         f = zpoly(self.polynomial)
  437         g = zpoly(self.coeff)
  438         if mod is None:
  439             j = pow(g, exponent, f)
  440         else:
  441             j = pow(g, exponent, f) % mod
  442             # what does mod exactly means?
  443         jcoeff = [j[i] for i in range(self.degree)]
  444         return BasicAlgNumber([jcoeff, d], self.polynomial)
  445 
  446     def inverse(self):
  447         f = qpoly(self.polynomial)
  448         g = qpoly(self.coeff)
  449         quot = f.extgcd(g)
  450 
  451         icoeff = [self.denom * quot[1][i] for i in range(self.degree)]
  452         #print icoeff
  453         denom_list = []
  454         for i in range(self.degree):
  455             if icoeff[i] == 0:
  456                 denom_list.append(1)
  457             else:
  458                 denom_list.append(icoeff[i].denominator)
  459         new_denom = 1
  460         for i in range(self.degree):
  461             new_denom = gcd.lcm(new_denom, denom_list[i])
  462         icoeff = [icoeff[i].numerator * new_denom // icoeff[i].denominator for i in range(self.degree)]
  463         return BasicAlgNumber([icoeff, new_denom], self.polynomial)
  464 
  465     def __truediv__(self, other):
  466         f = zpoly(self.polynomial)
  467         g = zpoly(self.coeff)
  468         t = BasicAlgNumber([other.coeff, other.denom], self.polynomial)
  469         k = t.inverse()
  470         h = zpoly(k.coeff)
  471         d = self.denom * k.denom
  472         j = (g * h).monic_mod(f)
  473         jcoeff = [j[i] for i in range(self.degree)]
  474         return BasicAlgNumber([jcoeff, d], self.polynomial)
  475 
  476     __div__ = __truediv__
  477     __floordiv__ = __truediv__
  478 
  479     def getConj(self):
  480         """
  481         Return (approximate) solutions of self.polynomial.
  482         We can discriminate the conjugate field of self by these values.
  483         """
  484         if not hasattr(self, "conj"):
  485             self.conj = NumberField(self.polynomial).getConj()
  486         return self.conj
  487 
  488     def getApprox(self):
  489         """
  490         Return the approximations of all conjugates of self.
  491         """
  492         if not hasattr(self, "approx"):
  493             APP = []
  494             for i in range(self.degree):
  495                 Approx = 0
  496                 for j in range(self.degree):
  497                     Approx += self.coeff[j]*(self.getConj()[i]**j)
  498                 APP.append(Approx)
  499             self.approx = APP
  500         return self.approx
  501 
  502     def getCharPoly(self):
  503         """
  504         Return the characteristic polynomial of self
  505         by compute products of (x-self^{(i)}).
  506         """
  507         if not hasattr(self, "charpoly"):
  508             Conj = self.getApprox()
  509             P = uniutil.polynomial({0:-Conj[0], 1:1}, ring.getRing(Conj[0]))
  510             for i in range(1, self.degree):
  511                 P *= uniutil.polynomial({0:-Conj[i], 1:1},
  512                     ring.getRing(Conj[i]))
  513             charcoeff = []
  514             for i in range(self.degree + 1):
  515                 if hasattr(P[i], "real"):
  516                     charcoeff.append(int(math.floor(P[i].real + 0.5)))
  517                 else:
  518                     charcoeff.append(int(math.floor(P[i] + 0.5)))
  519             self.charpoly = charcoeff
  520         return self.charpoly
  521 
  522     def getRing(self):
  523         """
  524         Return the algebraic number field contained self.
  525         """
  526         return NumberField(self.polynomial)
  527 
  528     def trace(self):
  529         """
  530         Compute the trace of self in K.
  531         """
  532         denom = self.denom
  533         n = len(self.polynomial) - 1
  534 
  535         tlist = [n]
  536         for k in range(1, n):
  537             s = 0
  538             for i in range(1, k):
  539                 s += tlist[k - i] * self.polynomial[n - i]
  540             tlist.append(-k * self.polynomial[n - k] - s)
  541 
  542         t = 0
  543         for j in range(len(tlist)):
  544             t += tlist[j] * self.coeff[j]
  545 
  546         if denom == 1:
  547             return t
  548         elif t % denom == 0:
  549             return t // denom
  550         else:
  551             return rational.Rational(t, denom)
  552 
  553     def norm(self):
  554         """
  555         Compute the norm of self in K.
  556         """
  557         f = zpoly(self.polynomial)
  558         g = zpoly(self.coeff)
  559         R = f.resultant(g)
  560 
  561         if self.denom == 1:
  562             return int(R)
  563         else:
  564             denom = self.denom**self.degree
  565             if isinstance(R, int):
  566                 if R % denom == 0:
  567                     return R // denom
  568                 else:
  569                     return rational.Rational(R, denom)
  570             else:
  571                 return R / denom
  572 
  573     def isAlgInteger(self):
  574         """
  575         Determine whether self is an algebraic integer or not.
  576         """
  577         Norm = self.norm()
  578         if isinstance(Norm, int):
  579             return True
  580         else:
  581             return False
  582 
  583     def ch_matrix(self):
  584         """
  585         Change style to MatAlgNumber.
  586         """
  587         list = []
  588         if self.denom == 1:
  589             list = self.coeff
  590         else:
  591             for i in range(self.degree):
  592                 list.append(rational.Rational(self.coeff[i], self.denom))
  593         return MatAlgNumber(list, self.polynomial)
  594 
  595 class MatAlgNumber(object):
  596     """
  597     The class for algebraic number represented by matrix.
  598     """
  599     def __init__(self, coefficient, polynomial):
  600         """
  601         """
  602         self.coeff = coefficient
  603         self.degree = len(coefficient)
  604         List = []
  605         for i in range(self.degree):
  606             stbasis = [0] * self.degree
  607             stbasis[i] = 1
  608             List.append(stbasis)
  609         List.append([- polynomial[i] for i in range(self.degree)])
  610         for l in range(self.degree - 2):
  611             basis1 = []
  612             basis = []
  613             for j in range(self.degree):
  614                 if j == 0:
  615                     basis1.append(0)
  616                 if j < self.degree - 1:
  617                     basis1.append(List[l + self.degree][j])
  618                 elif j == self.degree - 1:
  619                     basis2 = [List[l + self.degree][j] * - polynomial[k] for k in range(self.degree)]
  620             for i in range(self.degree):
  621                 basis.append(basis1[i] + basis2[i])
  622             List.append(basis)
  623         Matrix = []
  624         flag = 0
  625         for i in range(self.degree):
  626             basis3 = []
  627             for j in range(self.degree):
  628                 basis3.append([self.coeff[j] * k for k in List[j+flag]])
  629             for l in range(self.degree):
  630                 t = 0
  631                 for m in range(self.degree):
  632                     t += basis3[m][l]
  633                 Matrix.append(t)
  634             flag += 1
  635 
  636         self.matrix = matrix.createMatrix(self.degree, self.degree, Matrix)
  637         self.polynomial = polynomial
  638         self.field = NumberField(self.polynomial)
  639 
  640     def __repr__(self):
  641         return_str = '%s(%s, %s)' % (self.__class__.__name__, self.matrix.__repr__(), self.polynomial)
  642         return return_str
  643 
  644     def __neg__(self):
  645         mat = - self.matrix
  646         coeff = []
  647         for i in range(mat.row):
  648             coeff.append(mat[0][i])
  649         return MatAlgNumber(coeff, self.polynomial)
  650 
  651     def __add__(self, other):
  652         if isinstance(other, (int, long, rational.Rational)):
  653             other = MatAlgNumber(
  654                 [other] + [0] * (self.degree -1), self.polynomial)
  655         elif not isinstance(other, MatAlgNumber):
  656             return NotImplemented
  657         mat = self.matrix + other.matrix
  658         coeff = []
  659         for i in range(mat.row):
  660             coeff.append(mat[i+1][1])
  661         return MatAlgNumber(coeff, self.polynomial)
  662 
  663     __radd__ = __add__
  664 
  665     def __sub__(self, other):
  666         if isinstance(other, (int, long, rational.Rational)):
  667             other = MatAlgNumber(
  668                 [other] + [0] * (self.degree -1), self.polynomial)
  669         elif not isinstance(other, MatAlgNumber):
  670             return NotImplemented
  671         mat = self.matrix - other.matrix
  672         coeff = []
  673         for i in range(mat.row):
  674             coeff.append(mat[i+1][1])
  675         return MatAlgNumber(coeff, self.polynomial)
  676 
  677     def __mul__(self, other):
  678         if isinstance(other, MatAlgNumber):
  679             mat = self.matrix * other.matrix
  680         elif isinstance(other, (int, long, rational.Rational)):
  681             mat = other * self.matrix
  682         else:
  683             return NotImplemented
  684         coeff = []
  685         for i in range(mat.row):
  686             coeff.append(mat[i+1][1])
  687         return MatAlgNumber(coeff, self.polynomial)
  688 
  689     __rmul__ = __mul__
  690 
  691     def __truediv__(self, other):
  692         if isinstance(other, (int, long, rational.Rational)):
  693             return self * (rational.Rational(other) ** -1)
  694         elif not isinstance(other, MatAlgNumber):
  695             return NotImplemented
  696         other_mat = other.matrix.copy()
  697         other_mat.toFieldMatrix()
  698         mat = other_mat.inverse(self.matrix)
  699         coeff = []
  700         for i in range(mat.row):
  701             coeff.append(mat[i+1][1])
  702         return MatAlgNumber(coeff, self.polynomial)
  703 
  704     __div__ = __truediv__
  705     __floordiv__ = __truediv__
  706 
  707     def __pow__(self, other):
  708         mat = self.matrix ** other
  709         coeff = []
  710         for i in range(mat.row):
  711             coeff.append(mat[i+1][1])
  712         return MatAlgNumber(coeff, self.polynomial)
  713 
  714     def inverse(self):
  715         (self.matrix).toFieldMatrix()
  716         inv = (self.matrix).inverse()
  717         coeff = []
  718         for i in range(inv.row):
  719             coeff.append(inv[i+1][1])
  720         return MatAlgNumber(coeff, self.polynomial)
  721 
  722     def norm(self):
  723         return (self.matrix).determinant()
  724 
  725     def trace(self):
  726         return (self.matrix).trace()
  727 
  728     def getRing(self):
  729         """
  730         Return the algebraic number field contained self.
  731         """
  732         return NumberField(self.polynomial)
  733 
  734     def ch_basic(self):
  735         denom = 1
  736         for i in range(self.degree):
  737             if not isinstance(self.coeff[i], int):
  738                 denom *= gcd.lcm(denom, (self.coeff[i]).denominator)
  739         coeff = []
  740         for i in range(self.degree):
  741             if isinstance(self.coeff[i], int):
  742                 coeff.append(self.coeff[i] * denom)
  743             else:
  744                 coeff.append(int((self.coeff[i]).numerator * denom / (self.coeff[i]).denominator))
  745         return BasicAlgNumber([coeff, denom], self.polynomial)
  746 
  747 def prime_decomp(p, polynomial):
  748     """
  749     Return prime decomposition of (p) over Q[x]/(polynomial).
  750     """
  751     field = NumberField(polynomial)
  752     field_disc = field.field_discriminant()
  753     if (field.disc() // field_disc) % p != 0:
  754         return _easy_prime_decomp(p, polynomial)
  755     else:
  756         return _main_prime_decomp(p, polynomial)
  757 
  758 def _easy_prime_decomp(p, polynomial):
  759     """
  760     prime decomposition by factoring polynomial
  761     """
  762     # import nzmath.module as module
  763     f = fppoly(polynomial, p)
  764     d = f.degree()
  765     factor = f.factor()
  766     p_alg = BasicAlgNumber([[p] + [0] * (d - 1), 1], polynomial)
  767     if len(factor) == 1 and factor[0][1] == 1:
  768         return [(p_alg , 1 )]
  769     else:
  770         dlist = []
  771         for i in range(len(factor)):
  772             basis_list = []
  773             for j in range(d):
  774                 if factor[i][0][j] == 0:
  775                     basis_list.append(0)
  776                 else:
  777                     basis_list.append(factor[i][0][j].toInteger())
  778             dlist.append([BasicAlgNumber([basis_list, 1], polynomial),
  779                           factor[i][1]])
  780         return [( (p_alg, dlist_ele[0]), dlist_ele[1] ) for dlist_ele in dlist]
  781 
  782 def _main_prime_decomp(p, polynomial):
  783     """
  784     main step of prime decomposition
  785     """
  786     # import nzmath.module as module
  787     raise NotImplementedError
  788 
  789 def changetype(a, polynomial=[0, 1]):
  790     """
  791     Change 'a' to be an element of field K defined polynomial
  792     """
  793     if isinstance(a, (int, long)):
  794         coeff = [a] + [0] * (len(polynomial) - 2)
  795         return BasicAlgNumber([coeff, 1], polynomial)
  796     elif isinstance(a, rational.Rational):
  797         coeff = [a.numerator] + [0] * (len(polynomial) - 2)
  798         return BasicAlgNumber([coeff, a.denominator], polynomial)
  799     else:
  800         deg = len(a) - 1
  801         coeff = [0, 1] + [0] * (deg - 2)
  802         if a[-1] != 1:
  803             polynomial = [a[-2], 1]
  804             mul = a[-1]
  805             for j in range(deg - 2, -1, -1):
  806                 polynomial.insert(0, a[j] * mul)
  807                 mul *= a[-1]
  808             return BasicAlgNumber([coeff, a[-1]], polynomial)
  809         else:
  810             return BasicAlgNumber([coeff, 1], a)
  811 
  812 def disc(A):
  813     """
  814     Compute the discriminant of a_i, where A=[a_1,...,a_n]
  815     """
  816     n = A[0].degree
  817     list = []
  818     for i in range(n):
  819         for j in range(n):
  820             s = A[i] * A[j]
  821             list.append(s.trace())
  822     M = matrix.createMatrix(n, n, list)
  823     return M.determinant()
  824 
  825 def qpoly(coeffs):
  826     """
  827     Return a rational coefficient polynomial constructed from given
  828     coeffs.  The coeffs is a list of coefficients in ascending order.
  829     """
  830     terms = [(i, rational.Rational(c)) for (i, c) in enumerate(coeffs)]
  831     return uniutil.polynomial(terms, rational.theRationalField)
  832 
  833 def zpoly(coeffs):
  834     """
  835     Return an integer coefficient polynomial constructed from given
  836     coeffs.  The coeffs is a list of coefficients in ascending order.
  837     """
  838     return uniutil.polynomial(enumerate(coeffs), rational.theIntegerRing)
  839 
  840 def fppoly(coeffs, p):
  841     """
  842     Return a Z_p coefficient polynomial constructed from given
  843     coeffs.  The coeffs is a list of coefficients in ascending order.
  844     """
  845     return uniutil.polynomial(enumerate(coeffs), finitefield.FinitePrimeField(p))