"Fossies" - the Fresh Open Source Software Archive

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

    1 """
    2 module, ideal etc. for number field
    3 """
    4 from __future__ import division
    5 
    6 from nzmath.algfield import BasicAlgNumber, MatAlgNumber
    7 import nzmath.arith1 as arith1
    8 import nzmath.gcd as gcd
    9 import nzmath.vector as vector
   10 import nzmath.matrix as matrix
   11 import nzmath.rational as rational
   12 import nzmath.ring as ring
   13 import nzmath.prime as prime
   14 
   15 
   16 class Submodule(matrix.RingMatrix):
   17     """
   18     Submodule is a class for submodules (typically, Z^n).
   19     we assume that coeff_ring is PID.
   20     (i.e. hermite normal form (HNF) for matrices exists.)
   21     """
   22     
   23     def __init__(self, row, column, compo=0, coeff_ring=0, ishnf=None):
   24         """
   25         Submodule(row, column [,components, coeff_ring, ishnf])
   26         """
   27         if isinstance(compo, bool):
   28             ishnf = compo
   29             compo = 0
   30         elif isinstance(coeff_ring, bool):
   31             ishnf = coeff_ring
   32             coeff_ring = 0
   33         self._initialize(row, column, compo, coeff_ring)
   34         self.ishnf = ishnf
   35 
   36     @classmethod
   37     def fromMatrix(cls, mat, ishnf=None):
   38         """
   39         A constructor class method, which creates Submodule from a
   40         Matrix instance.
   41         """
   42         compo = []
   43         for row in mat.compo:
   44             compo += row
   45         return cls(mat.row, mat.column, compo, mat.coeff_ring, ishnf)
   46 
   47     def getGenerators(self):
   48         """
   49         return (current) generator of self (i.e. essentially, self.compo)
   50         """
   51         return [self[j] for j in range(1, self.column+1)]
   52 
   53     def isSubmodule(self, other):
   54         """ 
   55         Check self is in other as submodule
   56         """
   57         try:
   58             return self.sumOfSubmodules(other).isEqual(other)
   59         except:
   60             return False
   61 
   62     def isEqual(self, other):
   63         """
   64         determine self == other as module
   65         """
   66         if self.ishnf:
   67             if other.ishnf:
   68                 return self == other
   69             else:
   70                 other_copy = Submodule.fromMatrix(other)
   71                 other_copy.toHNF()
   72                 return self == other_copy
   73         else:
   74             if other.ishnf:
   75                 self_copy = Submodule.fromMatrix(self)
   76                 self_copy.toHNF()
   77                 return self_copy == other
   78             else:
   79                 self_copy = Submodule.fromMatrix(self)
   80                 self_copy.toHNF()
   81                 other_copy = Submodule.fromMatrix(other)
   82                 other_copy.toHNF()
   83                 return self_copy == other_copy
   84 
   85     def isContains(self, other):
   86         """
   87         determine whether other is in self or not.
   88         if you want other to be represented as self, use represent_element.
   89         """
   90         self_copy = Submodule.fromMatrix(self)
   91         return not isinstance(self_copy.represent_element(other), bool)
   92 
   93     def toHNF(self):
   94         """
   95         Change matrix to HNF(hermite normal form)
   96         Note that HNF do not always give basis of self
   97         (i.e. HNF may be redundant)
   98         """
   99         if not self.ishnf:
  100             hnf = self.hermiteNormalForm(True)
  101             if hnf == None: #zero module
  102                 hnf = matrix.zeroMatrix(self.row, 1, self.coeff_ring)
  103             self.compo = hnf.compo
  104             self.column = hnf.column
  105             self.ishnf = True
  106 
  107     def sumOfSubmodules(self, other):
  108         """
  109         Return module which is sum of self and other.
  110         """
  111         if self.row != other.row:
  112             raise matrix.MatrixSizeError()
  113         N = self.copy()
  114         N.extendColumn(other)
  115         module = Submodule.fromMatrix(N)
  116         module.toHNF()
  117         return module
  118 
  119     def intersectionOfSubmodules(self, other):
  120         """
  121         Return module which is intersection of self and other.
  122         """
  123         if self.row != other.row:
  124             raise matrix.MatrixSizeError()
  125         M1 = self.copy()
  126         M1.extendColumn(other)
  127         N = M1.kernelAsModule()
  128         if N == None: # zero kernel
  129             N = matrix.zeroMatrix(self.row, 1, self.coeff_ring)
  130         N1 = N.getBlock(1, 1, self.column, N.column) # N.column is dim(ker(M1))
  131         module = Submodule.fromMatrix(self * N1)
  132         module.toHNF()
  133         return module
  134 
  135     def represent_element(self, other):
  136         """
  137         represent other as a linear combination with hnf generators
  138         if other not in self, return False
  139         this method execute self.toHNF()
  140         """
  141         self.toHNF()
  142         return self._HNF_solve(other)
  143 
  144     def linear_combination(self, coeff):
  145         """
  146         return a vector corresponding a linear combination of (current) basis
  147         with given Z-coefficients.
  148         """
  149         sol = coeff[0] * self[1]
  150         for i in range(2, self.column + 1):
  151             sol += coeff[i - 1] * self[i]
  152         return sol
  153 
  154     def _HNF_solve(self, other):
  155         """
  156         return X over coeff_ring s.t. self * X = other,
  157         where self is HNF, other is vector
  158         """
  159         if self.row != len(other):
  160             raise vector.VectorSizeError()
  161         n = self.column
  162         zero = self.coeff_ring.zero
  163         X = vector.Vector([zero] * n)
  164         # search non-zero entry
  165         for i in range(1, self.row + 1)[::-1]:
  166             if self[i, n] != zero:
  167                 break
  168             if other[i] != zero:
  169                 return False
  170         else:
  171             return X
  172         # solve HNF system (only integral solution)
  173         for j in range(1, n + 1)[::-1]:
  174             sum_i = other[i]
  175             for k in range(j + 1, n + 1):
  176                 sum_i -= X[k] * self[i, k]
  177             try:
  178                 X[j] = ring.exact_division(sum_i, self[i, j])
  179             except ValueError: # division is not exact
  180                 return False
  181             ii = i
  182             i -= 1
  183             for i in range(1, ii)[::-1]:
  184                 if j != 1 and self[i, j - 1] != zero:
  185                     break
  186                 sum_i = X[j] * self[i, j]
  187                 for k in range(j + 1, n + 1):
  188                     sum_i += X[k] * self[i, k]
  189                 if sum_i != other[i]:
  190                     return False
  191             if i == 0:
  192                 break
  193         return X
  194 
  195 class Module(object):
  196     """
  197     for computing module with HNF over a number field.
  198     A module is a finitely generated sub Z-module.
  199     (Note that we do not assume rank of a module is deg(number_field).)
  200     """
  201     def __init__(self, pair_mat_repr, number_field, base=None, ishnf=False):
  202         """
  203         pair_mat_repr: generators represented with respect to base_module
  204                        1: [[deg(number_field))-integral tuple/vector], denominator]
  205                        2: [integral matrix whose the number of row is deg(number_field), 
  206                               denominator]
  207                        3: rational matrix whose the number of row is deg(number_field)
  208         number_field:  an instance of NumberField
  209         base:          a base module represented with respect to Z[theta]
  210                        1: [deg(number_field)-rational tuple/vector], which represents basis
  211                        2: deg(number_field)-square non-singular rational matrix
  212         ishnf:         if hnf_repr is already HNF-form, set True
  213         """
  214         self.number_field = number_field
  215         size = self.number_field.degree
  216         if isinstance(base, bool):
  217             ishnf = base
  218             base = None
  219         if base == None: # power basis
  220             self.base = matrix.unitMatrix(size, rational.theIntegerRing)
  221             self.base.toFieldMatrix()
  222         elif isinstance(base, list): # list repr
  223             self.base = matrix.FieldSquareMatrix(size, base)
  224         else: # matrix repr
  225             if base.row != size:
  226                 raise TypeError("The number of row of base is incorrect.")
  227             self.base = base.copy()
  228             self.base.toFieldMatrix()
  229         if not isinstance(pair_mat_repr, list): #rational matrix repr
  230             self.mat_repr, self.denominator = _toIntegerMatrix(
  231                 pair_mat_repr)
  232         else:
  233             self.denominator = pair_mat_repr[1]
  234             if isinstance(pair_mat_repr[0], list): # list repr
  235                 column = len(pair_mat_repr[0])
  236                 self.mat_repr = Submodule(
  237                     size, column, pair_mat_repr[0], ishnf)
  238             else: # integral matrix repr
  239                 self.mat_repr = Submodule.fromMatrix(
  240                     pair_mat_repr[0].copy(), ishnf)
  241 
  242     def _simplify(self):
  243         """
  244         simplify self.denominator
  245         """
  246         mat_repr, numer = _toIntegerMatrix(self.mat_repr, 2)
  247         denom_gcd = gcd.gcd(numer, self.denominator)
  248         if denom_gcd != 1:
  249             self.denominator = ring.exact_division(self.denominator, denom_gcd)
  250             self.mat_repr = Submodule.fromMatrix(
  251                 ring.exact_division(numer, denom_gcd) * mat_repr)
  252 
  253     def toHNF(self):
  254         """
  255         transform mat_repr to hnf form
  256         """
  257         self.mat_repr.toHNF()
  258 
  259     def __repr__(self):
  260         """
  261         formal representation
  262         """
  263         return_str = '%s([%s, %s], %s, %s)' % ( self.__class__.__name__, 
  264           repr(self.mat_repr), repr(self.denominator), 
  265           repr(self.number_field), repr(self.base) )
  266         return return_str
  267 
  268     def __str__(self):
  269         """
  270         simple representation
  271         """
  272         return str(
  273             (self.mat_repr, self.denominator)) + "\n over \n" + str(
  274             (self.base, self.number_field))
  275 
  276     def __eq__(self, other):
  277         """
  278         Check whether self == other or not (as module)
  279         """
  280         #if not(other in self.number_field):
  281         #    return False
  282         changed_self = self.change_base_module(other.base)
  283         return (changed_self.denominator == other.denominator
  284             ) and changed_self.mat_repr.isEqual(other.mat_repr)
  285 
  286     def __ne__(self, other):
  287         return self != other
  288 
  289     def __contains__(self, other):
  290         """
  291         Check whether other is an element of self or not
  292         """
  293         #if not(other in self.number_field):
  294         #    return False
  295         if isinstance(other, (tuple, vector.Vector)): #vector repr
  296             other_int, other_denom = _toIntegerMatrix(vector.Vector(other))
  297         elif isinstance(other, list):
  298             if isinstance(other[0], (tuple, vector.Vector)):
  299                 other_int = vector.Vector(other[0])
  300                 other_denom = other[1]
  301             else:
  302                 raise ValueError("input [vector, denominator]!")
  303         else: # field element
  304             if not isinstance(other, BasicAlgNumber):
  305                 other_copy = other.ch_basic()
  306             else:
  307                 other_copy = other
  308             other_base_repr = self.base.inverse(
  309                 vector.Vector(other_copy.coeff))
  310             other_int, other_int_denom = _toIntegerMatrix(other_base_repr)
  311             other_denom = other_int_denom * other_copy.denom
  312         if isinstance(other_int, matrix.Matrix):
  313             other_int = other_int[1]
  314         try:
  315             numer = ring.exact_division(self.denominator, other_denom)
  316         except ValueError: # division is not exact
  317             return False
  318         if numer != 1:
  319             other_vect = numer * other_int
  320         else:
  321             other_vect = other_int
  322         return not isinstance(
  323             self.mat_repr.represent_element(other_vect), bool)
  324 
  325     def __add__(self, other):
  326         """
  327         return self + other (as module)
  328         """
  329         # if self.number_field != other.number_field:
  330         #     raise ValueError("based on different number field")
  331         new_self = self.change_base_module(other.base)
  332         new_denominator = gcd.lcm(new_self.denominator, other.denominator)
  333         if new_self.denominator != new_denominator:
  334             multiply_denom = ring.exact_division(new_denominator,
  335                 new_self.denominator)
  336             new_self_mat_repr = multiply_denom * new_self.mat_repr
  337         else:
  338             new_self_mat_repr = new_self.mat_repr
  339         if other.denominator != new_denominator:
  340             multiply_denom = ring.exact_division(new_denominator,
  341                 other.denominator)
  342             new_other_mat_repr = multiply_denom * other.mat_repr
  343         else:
  344             new_other_mat_repr = other.mat_repr
  345         new_mat_repr = Submodule.fromMatrix(
  346             new_self_mat_repr).sumOfSubmodules(new_other_mat_repr)
  347         new_module = self.__class__(
  348             [new_mat_repr, new_denominator], self.number_field, other.base)
  349         new_module._simplify()
  350         return new_module
  351 
  352     def __mul__(self, other):
  353         """
  354         return self * other (as ideal computation)
  355         """
  356         if isinstance(other, 
  357             (int, long, rational.Integer, rational.Rational)):
  358             return self._rational_mul(other)
  359         #try:
  360         #     use __contains__ function?
  361         #    if other in self.number_field:
  362         if isinstance(other, (BasicAlgNumber, MatAlgNumber)):
  363             return self._scalar_mul(other)
  364         #except: # other is a module
  365         return self._module_mul(other)
  366 
  367     __rmul__ = __mul__
  368 
  369     def _rational_mul(self, other):
  370         """
  371         return other * self, assuming other is a rational element
  372         """
  373         if rational.isIntegerObject(other):
  374             other_numerator = rational.Integer(other)
  375             other_denominator = rational.Integer(1)
  376         else:
  377             other_numerator = other.numerator
  378             other_denominator = other.denominator
  379         denom_gcd = gcd.gcd(self.denominator, other_numerator)
  380         if denom_gcd != 1:
  381             new_denominator = ring.exact_division(
  382                 self.denominator, denom_gcd) * other_denominator
  383             multiply_num = other_numerator.exact_division(denom_gcd)
  384         else:
  385             new_denominator = self.denominator * other_denominator
  386             multiply_num = other_numerator
  387         new_module = self.__class__(
  388                          [self.mat_repr * multiply_num, new_denominator],
  389                          self.number_field, self.base, self.mat_repr.ishnf)
  390         new_module._simplify()
  391         return new_module
  392 
  393     def _scalar_mul(self, other):
  394         """
  395         return other * self, assuming other is a scalar
  396         (i.e. an element of self.number_field)
  397         """
  398         # use other is an element of higher or lower degree number field ?
  399         #if other.getRing() != self.number_field:
  400         #    raise NotImplementedError(
  401         #        "other is not a element of number field")
  402         if not isinstance(other, BasicAlgNumber):
  403             try:
  404                 other = other.ch_basic()
  405             except:
  406                 raise NotImplementedError(
  407                     "other is not an element of a number field")
  408         # represent other with respect to self.base
  409         other_repr, pseudo_other_denom = _toIntegerMatrix(
  410             self.base.inverse(vector.Vector(other.coeff)))
  411         other_repr = other_repr[1]
  412         # compute mul using self.base's multiplication
  413         base_mul = self._base_multiplication()
  414         n = self.number_field.degree
  415         mat_repr = []
  416         for k in range(1, self.mat_repr.column +1):
  417             mat_repr_ele = vector.Vector([0] * n)
  418             for i1 in range(1, n + 1):
  419                 for i2 in range(1, n + 1):
  420                     mat_repr_ele += self.mat_repr[i1, k] * other_repr[
  421                         i2] * _symmetric_element(i1, i2, base_mul)
  422             mat_repr.append(mat_repr_ele)
  423         mat_repr, denom, numer = _toIntegerMatrix(
  424             matrix.FieldMatrix(n, len(mat_repr), mat_repr), 1)
  425         denom = self.denominator * pseudo_other_denom * other.denom * denom
  426         gcds = gcd.gcd(denom, numer)
  427         if gcds != 1:
  428             denom = ring.exact_division(denom, gcds)
  429             numer = ring.exact_division(numer, gcds)
  430         if numer != 1:
  431             mat_repr = numer * mat_repr
  432         else:
  433             mat_repr = mat_repr
  434         return self.__class__([mat_repr, denom], self.number_field, self.base)
  435 
  436     def _module_mul(self, other):
  437         """
  438         return self * other as the multiplication of modules
  439         """
  440         #if self.number_field != other.number_field:
  441         #    raise NotImplementedError
  442         if self.base != other.base:
  443             new_self = self.change_base_module(other.base)
  444         else:
  445             new_self = self.copy()
  446         base_mul = other._base_multiplication()
  447         n = self.number_field.degree
  448         new_mat_repr = []
  449         for k1 in range(1, new_self.mat_repr.column + 1):
  450             for k2 in range(1, other.mat_repr.column + 1):
  451                 new_mat_repr_ele = vector.Vector([0] * n)
  452                 for i1 in range(1, n + 1):
  453                     for i2 in range(1, n + 1):
  454                         new_mat_repr_ele += new_self.mat_repr[
  455                             i1, k1] * other.mat_repr[i2, k2
  456                             ] * _symmetric_element(i1, i2, base_mul)
  457                 new_mat_repr.append(new_mat_repr_ele)
  458         new_mat_repr, denom, numer = _toIntegerMatrix(
  459             matrix.FieldMatrix(n, len(new_mat_repr), new_mat_repr), 1)
  460         denom = new_self.denominator * other.denominator * denom
  461         gcds = gcd.gcd(denom, numer)
  462         if gcds != 1:
  463             denom = ring.exact_division(denom, gcds)
  464             numer = ring.exact_division(numer, gcds)
  465         if numer != 1:
  466             mat_repr = numer * new_mat_repr
  467         else:
  468             mat_repr = new_mat_repr
  469         sol = self.__class__(
  470             [mat_repr, denom], self.number_field, other.base)
  471         sol.toHNF()
  472         return sol
  473 
  474     def __pow__(self, other):
  475         """
  476         self ** other (based on ideal multiplication)
  477         """
  478         if other <= 0:
  479             raise ValueError, "only support non-negative powering" 
  480         mul_part = self.copy()
  481         index = other
  482         while True:
  483             if index & 1:
  484                 try:
  485                     sol *= mul_part
  486                 except NameError:
  487                     sol = mul_part
  488             index >>= 1
  489             if not index:
  490                 return sol
  491             else:
  492                 mul_part *= mul_part    
  493 
  494     def _base_multiplication(self):
  495         """
  496         return [base[i] * base[j]] (as a numberfield element)
  497         this is a precomputation for computing _module_mul
  498         """
  499         if not hasattr(self, "_base_multiply"):
  500             base_mul = _base_multiplication(self.base, self.number_field)
  501             self._base_multiply = base_mul
  502         return self._base_multiply
  503 
  504     def copy(self):
  505         """
  506         create copy of self
  507         """
  508         new_pair_mat_repr = [[self.mat_repr.copy()[i] for i in range(1, 
  509             self.mat_repr.column + 1)], self.denominator]
  510         new_base_module = [self.base.copy()[i] for i in range(
  511             1, self.base.column + 1)]
  512         return self.__class__(new_pair_mat_repr, self.number_field, 
  513             new_base_module, True)
  514 
  515     def intersect(self, other):
  516         """
  517         return intersection of self and other
  518         """
  519         # if self.number_field != other.number_field:
  520         #     raise ValueError("based on different number field")
  521         if self.base != other.base:
  522             new_self = self.change_base_module(other.base)
  523         else:
  524             new_self = self.copy()
  525         new_denominator = gcd.lcm(new_self.denominator, other.denominator)
  526         if new_self.denominator != new_denominator:
  527             multiply_denom = ring.exact_division(
  528                 new_denominator, new_self.denominator)
  529             new_self_mat_repr = multiply_denom * new_self.mat_repr
  530         else:
  531             new_self_mat_repr = new_self.mat_repr
  532         if other.denominator != new_denominator:
  533             multiply_denom = ring.exact_division(
  534                 new_denominator, other.denominator)
  535             new_other_mat_repr = multiply_denom * other.mat_repr
  536         else:
  537             new_other_mat_repr = other.mat_repr
  538         new_mat_repr = Submodule.fromMatrix(
  539             new_self_mat_repr).intersectionOfSubmodules(new_other_mat_repr)
  540         new_module = self.__class__(
  541             [new_mat_repr, new_denominator], self.number_field, other.base)
  542         new_module._simplify()
  543         return new_module
  544     
  545     def issubmodule(self, other):
  546         """
  547         Check self is submodule of other.
  548         """
  549         if self.__class__.__name__ != other.__class__.__name__:
  550             return False
  551         if (self + other) == other: # as module
  552             return True
  553         else:
  554             return False
  555 
  556     def issupermodule(self, other):
  557         """
  558         Check other is supermodule of self.
  559         """
  560         if self.__class__.__name__ != other.__class__.__name__:
  561             return False
  562         else:
  563             return other.issubmodule(self)
  564 
  565     def represent_element(self, other):
  566         """
  567         represent other as a linear combination with generators of self
  568         if other not in self, return False
  569         Note that we do not assume self.mat_repr is HNF
  570         """
  571         #if other not in self.number_field:
  572         #    return False
  573         theta_repr = (self.base * self.mat_repr)
  574         theta_repr.toFieldMatrix()
  575         pseudo_vect_repr = theta_repr.inverseImage(
  576             vector.Vector(other.coeff))
  577         pseudo_vect_repr = pseudo_vect_repr[1]
  578         gcd_self_other = gcd.gcd(self.denominator, other.denom)
  579         multiply_numerator = self.denominator // gcd_self_other
  580         multiply_denominator = other.denom // gcd_self_other
  581 
  582         def getNumerator_and_Denominator(ele):
  583             """
  584             get the pair of numerator and denominator of ele
  585             """
  586             if rational.isIntegerObject(ele):
  587                 return (ele, 1)
  588             else: # rational
  589                 return (ele.numerator, ele.denominator)
  590 
  591         list_repr = []
  592         for j in range(1, len(pseudo_vect_repr) + 1):
  593             try:
  594                 numer, denom = getNumerator_and_Denominator(
  595                     pseudo_vect_repr[j])
  596                 list_repr_ele = ring.exact_division(
  597                     numer * multiply_numerator, denom * multiply_denominator)
  598                 list_repr.append(list_repr_ele)
  599             except ValueError: # division is not exact
  600                 return False
  601         return list_repr
  602 
  603     def change_base_module(self, other_base):
  604         """
  605         change base_module to other_base_module
  606         (i.e. represent with respect to base_module)
  607         """
  608         if self.base == other_base:
  609             return self
  610         n = self.number_field.degree
  611         if isinstance(other_base, list):
  612             other_base = matrix.FieldSquareMatrix(n,
  613                 [vector.Vector(other_base[i]) for i in range(n)])
  614         else: # matrix repr
  615             other_base = other_base.copy()
  616         tr_base_mat, tr_base_denom = _toIntegerMatrix(
  617             other_base.inverse(self.base))
  618         denom = tr_base_denom * self.denominator
  619         mat_repr, numer = _toIntegerMatrix(
  620             tr_base_mat * self.mat_repr, 2)
  621         d = gcd.gcd(denom, numer)
  622         if d != 1:
  623             denom = ring.exact_division(denom, d)
  624             numer = ring.exact_division(numer, d)
  625         if numer != 1:
  626             mat_repr = numer * mat_repr
  627         return self.__class__([mat_repr, denom], self.number_field, other_base)
  628 
  629     def index(self):
  630         """
  631         return order of a residue group over base_module 
  632         (i.e. [base_module : module] or [module : base_module] ** (-1))
  633         """
  634         if not  self.mat_repr.ishnf:
  635             mat_repr = Submodule.fromMatrix(self.mat_repr.copy())
  636             mat_repr.toHNF()
  637         else:
  638             mat_repr = self.mat_repr
  639         n = self.base.column
  640         if mat_repr.column < n:
  641             raise ValueError("index is infinite")
  642         det = 1
  643         for i in range(1, n + 1):
  644             # for HNF, determinant = product of diagonal elements
  645             det *= mat_repr[i, i]
  646         if self.denominator != 1:
  647             return  rational.Rational(det, self.denominator ** n)
  648         else:
  649             return det
  650 
  651     def smallest_rational(self):
  652         """
  653         return the Z-generator of intersection of self and rational field
  654         """
  655         if not  self.mat_repr.ishnf:
  656             mat_repr = Submodule.fromMatrix(self.mat_repr.copy())
  657             mat_repr.toHNF()
  658         else:
  659             mat_repr = self.mat_repr
  660         theta_repr, pseudo_denom = _toIntegerMatrix(
  661             self.base * mat_repr)
  662         return rational.Rational(theta_repr.hermiteNormalForm()[1, 1],
  663             self.denominator * pseudo_denom)
  664 
  665 
  666 class Ideal(Module):
  667     """
  668     for computing ideal with HNF (as Z-module).
  669     """
  670     def __init__(self, pair_mat_repr, number_field, base=None, ishnf=False):
  671         """
  672         Ideal is subclass of Module.
  673         Please refer to Module.__init__.__doc__
  674         """
  675         Module.__init__(self, pair_mat_repr, number_field, base, ishnf)
  676 
  677     issubideal = Module.issubmodule
  678     issuperideal = Module.issupermodule
  679 
  680     gcd = Module.__add__
  681     lcm = Module.intersect
  682 
  683     def __pow__(self, other):
  684         if other < 0:
  685             return Module.__pow__(self.inverse(), -other)
  686         elif other == 0:
  687             field = self.number_field
  688             n = field.degree
  689             int_basis = field.integer_ring()
  690             return Ideal([matrix.unitMatrix(n), 1], field, int_basis)
  691         else:
  692             return Module.__pow__(self, other)
  693 
  694     __pow__.__doc__ = Module.__pow__.__doc__   
  695 
  696     def inverse(self):
  697         """
  698         Return the inverse ideal of self
  699         """
  700         # Algorithm 4.8.21 in CCANT
  701         field = self.number_field
  702         # step 1 (Note that det(T)T^-1=(adjugate matrix of T))
  703         T = Ideal._precompute_for_different(field)
  704         d = int(T.determinant())
  705         inv_different = Ideal([T.adjugateMatrix(), 1], field,
  706                               field.integer_ring())
  707         # step 2
  708         inv_I = self * inv_different # already base is taken for integer_ring
  709         inv_I.toHNF()
  710         # step 3
  711         inv_mat, denom = _toIntegerMatrix(
  712             (inv_I.mat_repr.transpose() * T).inverse(), 0)
  713         numer = d * inv_I.denominator
  714         gcd_denom = gcd.gcd(denom, numer)
  715         if gcd_denom != 1:
  716             denom = ring.exact_division(denom, gcd_denom)
  717             numer = ring.exact_division(numer, gcd_denom)
  718         return Ideal([numer * inv_mat, denom], field, field.integer_ring())
  719 
  720     @classmethod
  721     def _precompute_for_different(cls, number_field):
  722         """
  723         Return T such that T^-1 represents HNF of inverse of different (codifferent)
  724         """
  725         field = number_field
  726         n = field.degree
  727         int_ring_mul = _base_multiplication(field.integer_ring(), field)
  728         T_lst = []
  729         for i in range(1, n + 1):
  730             each_T_lst = []
  731             for j in range(1, n + 1):
  732                 each_T_lst.append(MatAlgNumber(
  733                     _symmetric_element(
  734                     i, j, int_ring_mul).compo, field.polynomial
  735                     ).trace())
  736             T_lst.append(each_T_lst)
  737         T = matrix.RingSquareMatrix(n, T_lst)
  738         return T        
  739 
  740     def twoElementRepresentation(self):
  741         # Exercise 31 + Proposition 4.7.8 in CCANT
  742         # (need prime ideal decompostion)
  743         raise NotImplementedError
  744         # return Ideal_with_generator([alpha, beta], self.number_field)
  745 
  746     def norm(self):
  747         """
  748         return the norm of self
  749         (Note that Norm(I)=[Z_K : I] for an ideal I and an integral ring Z_K)
  750         """
  751         # use method of returning integral ring for a number field
  752         return self.change_base_module(
  753             self.number_field.integer_ring()).index()
  754 
  755     def isIntegral(self):
  756         """
  757         determine whether self is integral ideal or not
  758         """
  759         mdl = self.change_base_module(self.number_field.integer_ring())
  760         return mdl.denominator == 1
  761 
  762     def isPrime(self):
  763         """
  764         determine whether self is prime ideal or not
  765         """
  766         if not self.isIntegral():
  767             return False
  768         size = self.number_field.degree
  769         nrm = self.norm()
  770         p = arith1.floorpowerroot(nrm, size)
  771         if p ** size != nrm or not prime.primeq(p):
  772             return False
  773         return None #undefined
  774 
  775 
  776 class Ideal_with_generator(object):
  777     """
  778     ideal given as a generator
  779     (i.e. (a_1,...,a_n) = a_1 Z_K + ... + a_n Z_K)
  780     """
  781     def __init__(self, generator):
  782         """
  783         Ideal_with_generator(generator)
  784         generator: list of instances in BasicAlgNumber
  785                    (over same number_field)
  786         """
  787         self.generator = generator
  788         self.number_field = self.generator[0].field
  789 
  790     def __repr__(self):
  791         """
  792          formal representation
  793         """
  794         return '%s(%s)' % ( self.__class__.__name__, repr(self.generator) )
  795 
  796     def __str__(self):
  797         """
  798         simple representation
  799         """
  800         return str(self.generator)
  801 
  802     def __add__(self, other):
  803         """
  804         self + other as ideal
  805         other must be an instance of Ideal_with_generator
  806         """
  807         #if self.number_field != other.number_field:
  808         if self.number_field.polynomial != other.number_field.polynomial:
  809             raise NotImplementedError("use same number field")
  810         new_generator = self.generator[:]
  811         new_generator.extend(other.generator)
  812         return Ideal_with_generator(new_generator)
  813     
  814     def __mul__(self, other):
  815         """
  816         self * other as ideal
  817         other must be an instance of Ideal_with_generator
  818         """
  819         # if self.number_field != other.number_field:
  820         if self.number_field.polynomial != other.number_field.polynomial:
  821             raise NotImplementedError("use same number field")
  822         new_generator = [gen1 * gen2 
  823             for gen1 in self.generator for gen2 in other.generator]
  824         return Ideal_with_generator(new_generator)
  825 
  826     def __pow__(self, other):
  827         """
  828         self ** other (based on ideal multiplication)
  829         """
  830         if other <= 0:
  831             raise ValueError, "only support non-negative powering" 
  832         mul_part = self.copy()
  833         index = other
  834         while True:
  835             if index & 1:
  836                 try:
  837                     sol *= mul_part
  838                 except NameError:
  839                     sol = mul_part
  840             index >>= 1
  841             if not index:
  842                 return sol
  843             else:
  844                 mul_part *= mul_part
  845 
  846     def copy(self):
  847         """
  848         create copy of self
  849         """
  850         new_generator = []
  851         for gen in self.generator:
  852             new_generator.append(
  853                 BasicAlgNumber([gen.coeff[:], gen.denom], gen.polynomial))
  854         return self.__class__(new_generator)
  855 
  856     def to_HNFRepresentation(self):
  857         """
  858         Transform self to the corresponding ideal as HNF representation
  859         """
  860         int_ring = self.number_field.integer_ring()
  861         n = self.number_field.degree
  862         polynomial = self.number_field.polynomial
  863         k = len(self.generator)
  864         # separate coeff and denom for generators and base
  865         gen_denom = reduce( gcd.lcm,
  866                             (self.generator[j].denom for j in range(k)) )
  867         gen_list = [gen_denom * self.generator[j] for j in range(k)]
  868         base_int, base_denom = _toIntegerMatrix(int_ring)
  869         
  870         base_alg = [BasicAlgNumber([base_int[j], 1], 
  871         polynomial) for j in range(1, n + 1)]
  872         repr_mat_list = []
  873         for gen in gen_list:
  874             for base in base_alg:
  875                 new_gen = gen * base
  876                 repr_mat_list.append(
  877                     vector.Vector([new_gen.coeff[j] for j in range(n)]))
  878         mat_repr = matrix.RingMatrix(n, n * k, repr_mat_list)
  879         new_mat_repr, numer = _toIntegerMatrix(mat_repr, 2)
  880         denom = gen_denom * base_denom
  881         denom_gcd = gcd.gcd(numer, denom)
  882         if denom_gcd != 1:
  883             denom = ring.exact_division(denom, denom_gcd)
  884             numer = ring.exact_division(numer, denom_gcd)
  885             if numer != 1:
  886                 new_mat_repr = numer * new_mat_repr
  887         else:
  888             new_mat_repr = mat_repr
  889         return Ideal([new_mat_repr, denom], self.number_field)
  890 
  891     def twoElementRepresentation(self):
  892         # refer to [Poh-Zas] 
  893         # Algorithm 4.7.10 and Exercise 29 in CCANT
  894         """
  895         Reduce the number of generator to only two elements
  896         Warning: If self is not a prime ideal, this method is not efficient
  897         """
  898         k = len(self.generator)
  899         gen_denom = reduce( gcd.lcm, 
  900             (self.generator[j].denom for j in range(k)) )
  901         gen_list = [gen_denom * self.generator[j] for j in range(k)]
  902         int_I = Ideal_with_generator(gen_list)
  903         R = 1
  904         norm_I = int_I.norm()
  905         l_I = int_I.smallest_rational()
  906         if l_I.denominator > 1:
  907             raise ValueError, "input an integral ideal"
  908         else:
  909             l_I = l_I.numerator
  910         while True:
  911             lmda = [R for i in range(k)]
  912             while lmda[0] > 0:
  913                 alpha = lmda[0] * gen_list[0]
  914                 for i in range(1, k):
  915                     alpha += lmda[i] * gen_list[i]
  916                 if gcd.gcd(
  917                   norm_I, ring.exact_division(
  918                   alpha.norm(), norm_I)
  919                   ) == 1 or gcd.gcd(
  920                   norm_I, ring.exact_division(
  921                   (alpha + l_I).norm(), norm_I)) == 1:
  922                     l_I_ori = BasicAlgNumber(
  923                         [[l_I] + [0] * (self.number_field.degree - 1),
  924                         gen_denom], self.number_field.polynomial)
  925                     alpha_ori = BasicAlgNumber([alpha.coeff, 
  926                         gen_denom], self.number_field.polynomial)
  927                     return Ideal_with_generator([l_I_ori, alpha_ori])
  928                 for j in range(k)[::-1]:
  929                     if lmda[j] != -R:
  930                         lmda[j] -= 1
  931                         break
  932                     else:
  933                         lmda[j] = R
  934             R += 1
  935 
  936 mono_method = ["index", "smallest_rational", "inverse", "norm"]
  937 di_method = ["__eq__", "__ne__", "__contains__", "intersect", 
  938     "issubideal", "issuperideal"]
  939 
  940 for new_method in mono_method:
  941     exec("def new_func(myself): " + 
  942         "return myself.to_HNFRepresentation()." + new_method + "()")
  943     exec("Ideal_with_generator." + new_method + " = new_func")
  944 for new_method in di_method:
  945     exec("def new_func(myself, myother): " + 
  946         "return myself.to_HNFRepresentation()." + 
  947         new_method + "(myother.to_HNFRepresentation())")
  948     exec("Ideal_with_generator." + new_method + " = new_func")
  949 
  950 
  951 ######################################
  952 #functions for internal computations
  953 ######################################
  954 def _toIntegerMatrix(mat, option=0):
  955     """
  956     transform a (including integral) rational matrix to 
  957         some integral matrix as the following
  958     [option]
  959     0: return integral-matrix, denominator
  960        (mat = 1/denominator * integral-matrix)
  961     1: return integral-matrix, denominator, numerator
  962        (mat = numerator/denominator * reduced-int-matrix)
  963     2: return integral-matrix, numerator (assuming mat is integral)
  964        (mat = numerator * numerator-reduced-rational-matrix)
  965     """
  966     def getDenominator(ele):
  967         """
  968         get the denominator of ele
  969         """
  970         if rational.isIntegerObject(ele):
  971             return 1
  972         else: # rational
  973             return ele.denominator
  974     def getNumerator(ele):
  975         """
  976         get the numerator of ele
  977         """
  978         if rational.isIntegerObject(ele):
  979             return ele
  980         else: # rational
  981             return ele.numerator
  982     if isinstance(mat, vector.Vector):
  983         mat = mat.toMatrix(True)
  984     if option <= 1:
  985         denom = mat.reduce(
  986             lambda x, y: gcd.lcm(x, getDenominator(y)), 1)
  987         new_mat = mat.map(
  988             lambda x: getNumerator(x) * ring.exact_division(
  989             denom, getDenominator(x)))
  990         if option == 0:
  991             return Submodule.fromMatrix(new_mat), denom
  992     else:
  993         new_mat = mat
  994     numer = new_mat.reduce(
  995         lambda x, y: gcd.gcd(x, getNumerator(y)))
  996     if numer != 0:
  997         new_mat2 = new_mat.map(
  998             lambda x: ring.exact_division(getNumerator(x), numer))
  999     else:
 1000         new_mat2 = new_mat
 1001     if option == 1:
 1002         return Submodule.fromMatrix(new_mat2), denom, numer
 1003     else:
 1004         return Submodule.fromMatrix(new_mat2), numer
 1005 
 1006 def _symmetric_element(i, j, mat):
 1007     """
 1008     get (i, j)-element from lst ordered (1, 1), (1, 2), (2, 2), ...
 1009     we assume that
 1010     1: (j, i)-element is same as (i, j)-element (i.e. symmetric)
 1011     2: index of mat starts 1 (i.e. for Matrix only)
 1012     """
 1013     if i <= j:
 1014         return mat[int((j - 1) * j / 2 + i)]
 1015     else:
 1016         return mat[int((i - 1) * i / 2 + j)]
 1017 
 1018 def _base_multiplication(base, number_field):
 1019     """
 1020     return [base[i] * base[j]] (as a numberfield element)
 1021     this is a precomputation for computing _module_mul
 1022     """
 1023     n = number_field.degree
 1024     polynomial = number_field.polynomial
 1025     base_int, base_denom = _toIntegerMatrix(base)
 1026     base_alg = [BasicAlgNumber([base_int[j], 1], 
 1027         polynomial) for j in range(1, n + 1)]
 1028     base_ij_list = []
 1029     for j in range(n):
 1030         for i in range(j + 1):
 1031             base_ij = base_alg[i] * base_alg[j]
 1032             base_ij_list.append(
 1033                 vector.Vector([base_ij.coeff[k] * 
 1034                 rational.Rational(1, base_denom ** 2) for k in range(n)]))
 1035     base_ij = base.inverse(
 1036         matrix.FieldMatrix(n, len(base_ij_list), base_ij_list))
 1037     return base_ij