"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
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.
674         """
675         Module.__init__(self, pair_mat_repr, number_field, base, ishnf)
676
677     issubideal = Module.issubmodule
678     issuperideal = Module.issupermodule
679
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
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
```