"Fossies" - the Fresh Open Source Software Archive

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

```    1 from __future__ import division
2 import math
3 import copy
4 import warnings
5 import nzmath.gcd as gcd
6 import nzmath.arith1 as arith1
7 import nzmath.prime as prime
8 import nzmath.factor.misc as misc
9 import nzmath.bigrandom as bigrandom
10
11
13     """
14     The class is for reduced quadratic form.
15     """
16     def __init__(self, element, unit):
17         self.element = element # form = (a_1, a_2, a_3)
18         self.unit = unit
19
20     def __repr__(self):
21         return_str = 'ReducedQuadraticForm(%d, %d, %d)' % self.element
22         return return_str
23
24     def __str__(self):
25         return_str = '[%d, %d, %d]' % self.element
26         return return_str
27
28     def __mul__(self, other):
29         if not isinstance(other, ReducedQuadraticForm):
30             return NotImplemented
31         composition = compositePDF(self.element, other.element)
32         return self.__class__(composition, self.unit)
33
34     def __pow__(self, exp):
35         if not isinstance(exp, (int, long)):
36             raise TypeError("powering index must be an integer.")
37         powered_form = powPDF(self.element, exp, self.unit)
38         return self.__class__(powered_form, self.unit)
39
40     def __truediv__(self, other):
41         invel = other.inverse()
42         composition = compositePDF(self.element, invel.element)
43         return self.__class__(composition, self.unit)
44
45     def inverse(self):
46         """
47         Return inverse of the form.
48         """
49         if self.element == self.unit:
50             return copy.deepcopy(self)
51         else:
52             f_1, f_2, f_3 = self.element
53             return self.__class__(reducePDF((f_1, -f_2, f_3)), self.unit)
54
55     def __eq__(self, other):
56         if self.__class__ != other.__class__:
57             return False
58         if self.element == other.element and self.unit == other.unit:
59             return True
60         else:
61             return False
62
63     def __hash__(self):
64         val = sum([hash(ele) for ele in self.element])
65         val += sum([hash(ele) for ele in self.unit])
66         return val
67
68     def __ne__(self, other):
69         return not self.__eq__(other)
70
71     def disc(self):
72         """
73         Return discriminant of the form.
74         """
75         return self.element[1]**2 - 4*self.element[0]*self.element[2]
76
77     def repOfModule(self):
78         ld = self.disc()
79         a_m2 = 2*self.element[0]
80         rb = -self.element[1]
81
82         return_str = '%s + root(%s) / %s' % (rb, ld, a_m2)
83         return return_str
84
85     def __iter__(self):
86         """
87         Return iterator for form (a, b, c).
88         """
89         return iter(self.element)
90
91     def __getitem__(self, index):
92         """
93         Return each component of form (a, b, c),
94         """
95         return self.element[index]
96
97
98 class ClassGroup(object):
99     """
100     The class is for class group.
101     """
102     def __init__(self, disc, classnum, elements=None):
103         # element is an element of some class (for example ReducedQuadraticForm
104         self.disc = disc
105         self.rootoftree = []
106         self.rootornot = 0
107         if not elements:
108             self.elements = []
109         else:
110             self.elements = copy.deepcopy(elements)
111         self.classnum = classnum
112         self.expunit = unit_form(disc)
113
114     def __str__(self):
115         return_str = "class of ClassGroup:\n"
116         return_str = return_str + 'disc is %s\n' % self.disc
117         return_str = return_str + 'rootoftree is %s' % self.rootoftree
118         return return_str
119
120     def __repr__(self):
121         return "ClassGroup(%d, %d, %s)" % (self.disc, self.classnum, self.elements)
122
123     def inserels(self, newlist):
124         for newel in newlist:
125             self.inserel(newel)
126
127     def inserel(self, newel):
128         newestl = copy.deepcopy(newel)
130
131     def inststree(self, newlist):
132         for newel in newlist:
133             self.insttree(newel)
134
135     def insttree(self, newel0):
136         newel = copy.deepcopy(newel0)
137         disc = newel.element[1]**2 - 4*newel.element[0]*newel.element[2]
138         if disc != self.disc:
139             raise ValueError("this value is not an element of the discriminant")
140         if self.rootornot == 0:
141             self.rootoftree = [newel, [], []]
142             self.rootornot = 1
143             return True
144         else:
145             curntpnt = self.rootoftree
146         while curntpnt != []:
147             if newel.element == curntpnt[0].element:
148                 return True
149             elif newel.element < curntpnt[0].element:
150                 curntpnt = curntpnt[1]
151             else:
152                 curntpnt = curntpnt[2]
153
154         curntpnt.append(newel)
155         self.elements.append(newel)
156         curntpnt.append([])
157         curntpnt.append([])
158
159     def search(self, tarel):
160         curntpnt = self.rootoftree
161         while (curntpnt != []):
162             if tarel.element == curntpnt[0].element:
163                 return curntpnt[0]
164             elif tarel.element < curntpnt[0].element:
165                 curntpnt = curntpnt[1]
166             else:
167                 curntpnt = curntpnt[2]
168         return False
169
170     def retel(self):
171         copyofroot = copy.deepcopy(self.rootoftree)
172         tpa = []
173         while copyofroot:
174             curntpnt = copyofroot
175             while True:
176                 if curntpnt[1]:
177                     curntpnt = curntpnt[1]
178                 elif curntpnt[2]:
179                     curntpnt = curntpnt[2]
180                 else:
181                     tpa.append(curntpnt[0])
182                     del curntpnt[:3]
183                     break
184         return tpa
185
186 def class_formula(disc, uprbd):
187     """
188     Return the approximation of class number 'h' with the given discriminant.
189     h = sqrt(|D|)/pi (1 - (D/p)(1/p))^{-1} where p is less than ubound.
190     """
191     ht = math.sqrt(abs(disc)) / math.pi
192     ml = number_unit(disc) / 2
193
194     for factor in prime.generator_eratosthenes(uprbd):
195         ml = ml * (1 - (kronecker(disc, factor) / factor))**(-1)
196     return int(ht * ml + 0.5)
197
198 def class_number(disc, limit_of_disc=100000000):
199     """
200     Return class number with the given discriminant by counting reduced forms.
201     Not only fundamental discriminant.
202     """
203     if disc & 3 not in (0, 1):
204         raise ValueError("a discriminant must be 0 or 1 mod 4")
205     if disc >= 0:
206         raise ValueError("a discriminant must be negative")
207     if -disc >= limit_of_disc:
208         warnings.warn("the discriminant seems to have too big absolute value")
209
210     h = 1
211     b = disc & 1
212     c_b = int(math.sqrt(-disc / 3))
213
214     while b <= c_b:
215         q = (b**2 - disc) >> 2
216         a = b
217         if a <= 1:
218             a = 2
219         while a**2 <= q:
220             if q % a == 0 and gcd.gcd_of_list([a, b, q//a])[0] == 1:
221                 if a == b or a**2 == q or b == 0:
222                     h += 1
223                 else:
224                     h += 2
225             a += 1
226         b += 2
227     return h
228
229 def class_group(disc, limit_of_disc=100000000):
230     """
231     Return the class number and the class group with the given discriminant
232     by counting reduced forms. Not only fundamental discriminant.
233     """
234     if disc & 3 not in (0, 1):
235         raise ValueError("a discriminant must be 0 or 1 mod 4")
236     if disc >= 0:
237         raise ValueError("a discriminant must be negative")
238     if -disc >= limit_of_disc:
239         warnings.warn("the discriminant seems to have too big absolute value")
240
241     h = 1
242     b = disc & 1
243     c_b = int(math.sqrt(-disc / 3))
244
245     ret_list = []
246     f_a = 1
247     if disc & 3 == 0:
248         f_b = 0
249         f_c = -(disc >> 2)
250     else:
251         f_b = 1
252         f_c = -((disc - 1) >> 2)
253     unit = (f_a, f_b, f_c)
255
256     while b <= c_b:
257         q = (b**2 - disc) >> 2
258         a = b
259         if a <= 1:
260             a = 2
261         while a**2 <= q:
262             if q % a == 0 and gcd.gcd_of_list([a, b, q//a])[0] == 1:
263                 f_c = (b**2 - disc) // (4 * a)
264                 if a == b or a**2 == q or b == 0:
265                     h += 1
266                     ret_list.append(ReducedQuadraticForm((a, b, f_c), unit))
267                 else:
268                     h += 2
269                     ret_list.append(ReducedQuadraticForm((a, b, f_c), unit))
270                     ret_list.append(ReducedQuadraticForm((a, -b, f_c), unit))
271             a += 1
272         b += 2
273
274     return h, ret_list
275
276
278     """
279     The class is for reduced quadratic form for *bsgs functions.
280     """
281     def __init__(self, element, unit):
282         ReducedQuadraticForm.__init__(self, element, unit)
283         self.ind = -1
284         self.alpha = []
285         self.beta = []
286         self.s_parent = 0
287         self.g_parent = 0
288
289
290 def class_number_bsgs(disc):
291     """
292     Return the class number with the given discriminant.
293     """
294     if disc & 3 not in (0, 1):
295         raise ValueError("a discriminant must be 0 or 1 mod 4")
296
297     if disc >= 0:
298         raise ValueError("a discriminant must be negative")
299
300     lx = max(arith1.floorpowerroot(abs(disc), 5), 500 * (math.log(abs(disc)))**2)
301     uprbd = int(class_formula(disc, int(lx)) * 3 / 2)
302     lwrbd = (uprbd >> 1) + 1
303     bounds = [lwrbd, uprbd]
304
305     # get the unit
306     element = RetNext(disc)
307     ut = element.unit()
308
309     # append the unit to subset of G
310     sossp = ClassGroup(disc, 0, [])
311     sogsp = ClassGroup(disc, 0, [])
312     sossp.insttree(ut)
313     sogsp.insttree(ut)
314
315     h = 1 # order
316     finished = False
317     while not finished:
318         mstp1 = bounds[1] - bounds[0]
319         if mstp1 <= 1:
320             q = 1
321         else:
322             q = arith1.floorsqrt(mstp1)
323             if misc.primePowerTest(mstp1)[1] != 2:
324                 q += 1
325         # get next element
326         nt = element.retnext()
327         # x is the set of elements of G
328         x = [ut, nt ** h]
329         if q > 2:
330             x.extend([0] * (q - 2))
331         # compute small steps
332         if x[1] == ut:
333             # compute the order of nt
334             n = trorder(h, sossp, sogsp, nt, disc)
335         else:
336             n = trbabysp(q, x, bounds, sossp, sogsp, ut, h, nt, disc)
337
338         # finished?
339         finished, h, sossp, sogsp = isfinished_trbsgs(lwrbd, bounds, h, n, sossp, sogsp, nt, disc)
340
341     return h
342
343 def class_group_bsgs(disc, classnum, qin):
344     """
345     Return the construction of the class group with the given discriminant.
346     """
347     if disc & 3 not in (0, 1):
348         raise ValueError("a discriminant must be 0 or 1 mod 4")
349     if disc >= 0:
350         raise ValueError("a discriminant must be negative")
351
352     matla = []
353     lstofelg = []
354     lpt = []
355
356     # compute bounds
357     qpt = qin[0] ** qin[1]
358     uprbd = qpt + 1
359     lwrbd = (uprbd >> 1) + 1
360     if lwrbd > uprbd:
361         raise TypeError("lower bound needs to be less than upper bound")
362     if lwrbd <= (uprbd / 2):
363         raise TypeError("upper bound / 2 needs to be more than lower bound")
364     bounds = [lwrbd, uprbd]
365
366     # get the unit
367     uto = unit_form(disc)
368     ut = ReducedQuadraticFormForBSGS(uto, uto)
369
370     # append the unit to subset of G
371     sossp = ClassGroup(disc, classnum, []) # a subset of G
372     sogsp = ClassGroup(disc, classnum, []) # a subset of G
373     utwi = copy.deepcopy(ut)
374     utwi.alpha.append([0, ut, 1])
375     utwi.beta.append([0, ut, 1])
376     sossp.insttree(utwi)
377     sogsp.insttree(utwi)
378     finished = False
379
380     # take a new element of the group.
381     indofg = 1
382     while not finished:
383         # get next element
384         while True:
385             nt = rand_generator(disc, classnum, qin)
386             if nt not in lstofelg:
387                 lstofelg.append(nt)
388                 break
389
390         # compute small steps
391         assert nt != ut, "random element may not be equal to unit"
392         n, tmp_ss, tmp_gs = babyspcv(bounds, sossp, sogsp, utwi, nt, disc, classnum)
393         matla.append(setind(n, indofg, tmp_ss, tmp_gs))
394         finished, sossp, sogsp = isfinished_bsgscv(n, sossp, sogsp, nt, lpt, qpt, disc, classnum, indofg)
395         indofg = indofg + 1
396     return lstofelg, matla
397
398 ##############################################################
399 # following functions are sub functions for above functions. #
400 ##############################################################
401
402 class RetNext:
403     """
404     Next element producer.
405     """
406     def __init__(self, disc):
407         self.disc = disc
408         self.utroot = unit_form(disc)
409         self._unit = ReducedQuadraticForm(self.utroot, self.utroot)
410         self.cnt = 1
411         self.previous = []
412         self.elhash = range(int(math.sqrt(abs(disc) // 3)) + 2)
413
414     def unit(self):
415         """
416         Return reduced quadratic unit form.
417         """
418         return self._unit
419
420     def retnext(self):
421         """
422         Return the next random element.
423         """
424         while True:
425             next_form = ReducedQuadraticForm(self._randomele1(), self.utroot)
426             if not self.has_in_hash(next_form):
428                 return next_form
429
430     def _randomele1(self):
431         """
432         Return a reduced random form with the given discriminant.
433         """
434         while True:
435             nextp = prime.nextPrime(self.cnt)
436             self.cnt += 1
437             if kronecker(self.disc, nextp) == 1:
438                 nxtfm = sqroot(self.disc, nextp)
439                 if not self.previous or nxtfm != self.previous:
440                     self.previous = nxtfm
441                     return nxtfm
442
443     def add_to_hash(self, form):
444         """
445         Add a form to hash
446         """
447         key = form.element[0]
448         if isinstance(self.elhash[key], (int, long)):
449             self.elhash[key] = [form]
450         else:
451             self.elhash[key].append(form)
452
453     def has_in_hash(self, form):
454         """
455         check hash
456         """
457         key = form.element[0]
458         if isinstance(self.elhash[key], (int, long)):
459             return False
460         return form in self.elhash[key]
461
462
463 def disc(f):
464     """
465     Return the discriminant of the given quadratic form 'f'.
466     f = [a, b, c]
467     """
468     if len(f) != 3:
469         raise TypeError("form must be composed of 3 integers")
470     for i in f:
471         if not isinstance(i, (int, long)):
472             raise TypeError("all components must be integers")
473     return (f[1]*f[1] - 4*f[0]*f[2])
474
475 def reducePDF(f):
476     """
477     Return the reduced form of the given positive definite form 'f'.
478     f = (f_1, f_2, f_3)
479     """
480     f_1, f_2, f_3 = f
481     if f_1 < 0:
482         raise ValueError("f_1 must be positive in quadratic form f=(f_1,f_2,f_3).")
483     if (f_2**2 - 4*f_1*f_3) >= 0:
484         raise ValueError("discriminant (D= f_2^2 - 4*f_1*f_3) must be negative.")
485     if -f_1 < f_2 <= f_1:
486         if f_1 > f_3:
487             f_2 = -f_2
488             f_1, f_3 = f_3, f_1
489         else:
490             if (f_1 == f_3) and (f_2 < 0):
491                 f_2 = -f_2
492             return (f_1, f_2, f_3)
493     while 1:
494         q = f_2 // (2*f_1)
495         r = f_2 - q*(2*f_1)
496         if r > f_1:
497             r -= 2*f_1
498             q += 1
499         f_3 -= ((f_2 + r)>>1)*q
500         f_2 = r
501         if f_1 > f_3:
502             f_2 = -f_2
503             f_1, f_3 = f_3, f_1
504             continue
505         else:
506             if (f_1 == f_3) and (f_2 < 0):
507                 f_2 = -f_2
508             return (f_1, f_2, f_3)
509
510 def sqrPDF(f):
511     """
512     Return the square of the given quadratic form 'f'.
513     """
514     # compute disc and etc
515     D = disc(f)
516     sogsp = arith1.floorpowerroot(int(abs(D / 4)), 4)
517     (u, v, d_1) = euclid_exd(f[1], f[0])
518
519     la = f[0] // d_1
520     lb = f[1] // d_1
521     lc = (-f[2] * u) % la
522     c_1 = la - lc
523     if c_1 < lc:
524         lc = -c_1
525
526     # partial reduction
527     v_2, v_3, z, d, v = parteucl(la, lc, sogsp)
528
529     if z == 0:
530         g = (lb * v_3 + f[2]) // d
531         a_2 = d ** 2
532         c_2 = v_3 ** 2
533         b_2 = f[1] + (d + v_3)**2 - a_2 - c_2
534         c_2 = c_2 + g * d_1
535         return reducePDF((a_2, b_2, c_2))
536
537     e = (f[2] * v + lb * d) // la
538     g = (e * v_2 - lb) // v
539     b_2 = e * v_2 + v * g
540     if d_1 > 1:
541         b_2 = d_1 * b_2
542         v = d_1 * v
543         v_2 = d_1 * v_2
544
545     a_2 = d ** 2
546     c_2 = v_3 ** 2
547     b_2 = b_2 + (d + v_3) ** 2 - a_2 - c_2
548     a_2 = a_2 + e * v
549     c_2 = c_2 + g * v_2
550     return reducePDF((a_2, b_2, c_2))
551
552 def powPDF(f, exp, ut=None):
553     """
554     Return the powering 'exp' of the given quadratic form 'f'.
555     """
556     if ut is None:
557         D = disc(f)
558         ut = unit_form(D)
559
560     if exp == 0:
561         return ut
562     elif exp == 1:
563         return f
564     elif f == ut:
565         return ut
566     if exp < 0:
567         lexp = -exp
568         sz = (f[0], -f[1], f[2])
569     else:
570         lexp = exp
571         sz = f
572     # Right-Left Binary algorithm
573     sy = ut
574     while True:
575         if (lexp & 1) == 1:
576             sy = compositePDF(sz, sy)
577         lexp >>= 1
578         if lexp == 0:
579             return sy
580         else:
581             sz = sqrPDF(sz)
582
583 def compositePDF(f_1, f_2):
584     """
585     Return the reduced form of composition of the given forms 'f_1' and 'f_2'.
586     'f_1' and 'f_2' are quadratic forms with same disc.
587     """
588     if gcd.gcd_of_list(f_1)[0] != 1:
589         raise ValueError("coefficients of a quadratic form must be relatively prime")
590     if gcd.gcd_of_list(f_2)[0] != 1:
591         raise ValueError("coefficients of a quadratic form must be relatively prime")
592     if disc(f_1) != disc(f_2):
593         raise ValueError("two quadratic forms must have same discriminant")
594
595     if f_1[0] > f_2[0]:
596         f_1, f_2 = f_2, f_1
597
598     s = (f_1[1] + f_2[1]) >> 1
599     n = f_2[1] - s
600
601     if f_2[0] % f_1[0] == 0:
602         y_1 = 0
603         d = f_1[0]
604     else:
605         u, v, d = euclid_exd(f_2[0], f_1[0])
606         y_1 = u
607
608     if s % d == 0:
609         y_2 = -1
610         x_2 = 0
611         d_1 = d
612     else:
613         u, v, d_1 = euclid_exd(s, d)
614         x_2 = u
615         y_2 = -v
616
617     v_1 = f_1[0] // d_1
618     v_2 = f_2[0] // d_1
619     r = (y_1*y_2*n - x_2*f_2[2]) % v_1
620
621     b_3 = f_2[1] + 2*v_2*r
622     a_3 = v_1*v_2
623     c_3 = (f_2[2]*d_1 + r*(f_2[1] + v_2*r)) // v_1
624
625     return reducePDF((a_3, b_3, c_3))
626
627 def unit_form(disc):
628     """
629     Return generated quadratic form with the given discriminant.
630     """
631     if disc & 3 == 0:
632         a = 1
633         b = 0
634         c = disc // -4
635     elif disc & 3 == 1:
636         a = 1
637         b = 1
638         c = (disc - 1) // -4
639     else:
640         raise ValueError("discriminant is not 0 or 1 mod 4.")
641     return (a, b, c)
642
643 def kronecker(a, b):
644     """
645     Compute the Kronecker symbol (a/b) using algo 1.4.10 in Cohen's book.
646     """
647     tab2 = (0, 1, 0, -1, 0, -1, 0, 1)
648     if b == 0:
649         if abs(a) != 1:
650             return 0
651         if abs(a) == 1:
652             return 1
653     if a & 1 == 0 and b & 1 == 0:
654         return 0
655
656     v, b = arith1.vp(b, 2)
657     if v & 1 == 0:
658         k = 1
659     else:
660         k = tab2[a & 7]
661     if b < 0:
662         b = -b
663         if a < 0:
664             k = -k
665     while a:
666         v, a = arith1.vp(a, 2)
667         if v & 1:
668             k *= tab2[b & 7]
669         if a & b & 2:
670             # both a and b are 3 mod 4
671             k = -k
672         r = abs(a)
673         a = b % r
674         b = r
675     if b > 1:
676         # a and be are not coprime
677         return 0
678     return k
679
680 def number_unit(disc):
681     """
682     Return the number of units with the given discriminant.
683     """
684     if disc < -4:
685         return 2
686     elif disc == -4:
687         return 4
688     elif disc == -3:
689         return 6
690     else:
691         raise ValueError
692
693 def crt(inlist):
694     """
695     Chinese Remainder Theorem, Algo. 1.3.11 of Cohen's Book.
696     """
697     k = len(inlist)
698     if k < 2:
699         raise ValueError("nothing to be done for one element")
700     ccj = [None] # gabage to simplify loop index
701     ellist = list(inlist)
702     ellist.sort()
703     modulus = 1
704     for j in range(1, k):
705         modulus *= ellist[j - 1][1]
706         d, tpl = gcd.gcd_of_list([modulus % ellist[j][1], ellist[j][1]])
707         if d > 1:
708             raise ValueError("moduli are not pairwise coprime")
709         ccj.append(tpl[0])
710
711     yj = [ellist[0][0] % ellist[0][1]]
712     for j in range(1, k):
713         ctp = yj[-1]
714         for i in range(j - 2, -1, -1):
715             ctp = yj[i] + ellist[i][1] * ctp
716         yj.append(((ellist[j][0] - ctp) * ccj[j]) % ellist[j][1])
717     outp = yj.pop()
718     for j in range(k - 2, -1, -1):
719         outp = yj.pop() + (ellist[j][1]) * outp
720     return outp
721
722 def rand_generator(disc, classnum, qin):
723     """
724     Return the reduced random quadratic form with given discriminant and order t,
725     where t = classnum / a ** b and qin = [a, b].
726     """
727     q = qin[0]**qin[1]
728     unit = unit_form(disc)
729     while True:
730         elgt1 = randomele(disc, unit)
731         elg1 = elgt1 ** (classnum // q)
732         if elg1.element == elg1.unit:
733             continue
734         return elg1
735
736 def sqroot(disc, p):
737     """
738     Return a reduced quadratic form with the given discriminant.
739     'disc' is a quadratic residue mod 'p'.
740     """
741     if p == 2: # if 8 | disc => (disc / 8) = 0, 8 not | disc but 4 | disc => 2
742         if (disc & 7) == 0:
743             bp = disc
744         elif (disc & 3) == 0: # 4 - 4 * odd % 8 => 0
745             bp = 2
746         elif (disc & 7) == 1: # disc is odd and disc % 8 is 1
747             bp = disc
748         else: # disc is odd and disc & 3 is 1 => impossible (-5 / 2) = -1
749             raise ValueError("disc is odd and disc & 3 is 1 => impossible (-5 / 2) = -1")
750     else:
751         bpf1 = arith1.modsqrt(disc, p)
752         bpf2 = disc
753         bp = crt([(bpf1, p), (bpf2, 4)])
754     if bp > p:
755         bp = 2 * p - bp
756
757     fpt = reducePDF((p, bp, ((bp ** 2) - disc) // (4 * p)))
758     return fpt
759
760 def randomele(disc, unit):
761     """
762     Return a reduced random form with the given discriminant and the given unit.
763     Also random element is not unit.
764     """
765     limit = int(math.sqrt(-disc / 3))
766     while True:
767         a = bigrandom.randrange(1, limit + 1)
768         ind = 0
769         while ind < 2*a:
770             b = bigrandom.randrange(a)
771             if bigrandom.randrange(2):
772                 b = -b
773             tp = disc - b**2
774             if tp % (-4 * a) == 0:
775                 c = tp // (-4 * a)
776                 if gcd.gcd_of_list([a, b, c])[0] != 1:
777                     continue
778                 red = reducePDF((a, b, c))
779                 if red != unit:
780                     return ReducedQuadraticForm(red, unit)
781             ind += 1
782
783
784 def isfundamental(disc):
785     """
786     Determine whether the given discriminant is fundamental or not.
787     """
788     if disc == 1:
789         return False
790     spt = misc.squarePart(abs(disc))
791     if (disc & 3) == 1:
792         return spt == 1
793     elif (disc & 3) == 0:
794         spt //= 2
795         if spt != 1:
796             return False
797         discof = (disc >> 2) & 3
798         return discof == 2 or discof == 3
799     return False
800
801 def euclid_exd(a, b):
802     """
803     Return a tuple (u, v, d); they are the greatest common divisor d
804     of two integers a and b and u, v such that d = a * u + b * v.
805     """
806     if not isinstance(a, (int, long)) or not isinstance(b, (int, long)):
807         raise TypeError
808     u = 1
809     d = a
810     if b == 0:
811         v = 0
812         return (u, v, d)
813     else:
814         v_1 = 0
815         v_3 = b
816
817         while 1:
818             if v_3 == 0:
819                 v = (d - a*u) // b
820                 return (u, v, d)
821             q = d // v_3
822             t_3 = d % v_3
823             t_1 = u - q*v_1
824             u = v_1
825             d = v_3
826             v_1 = t_1
827             v_3 = t_3
828
829 def parteucl(a, b, sogsp):
830     """
831     Do extended partial Euclidean algorithm on 'a' and 'b'.
832     """
833     v = 0
834     d = a
835     v_2 = 1
836     v_3 = b
837     z = 0
838
839     while 1:
840         if abs(v_3) > sogsp:
841             # euclidean step
842             q = d // v_3
843             t_3 = d % v_3
844             t_2 = v - (q * v_2)
845             v = v_2
846             d = v_3
847             v_2 = t_2
848             v_3 = t_3
849             z = z + 1
850         else:
851             if z & 1:
852                 v_2 = -v_2
853                 v_3 = -v_3
854             return (v_2, v_3, z, d, v)
855
856 def isfinished_bsgscv(n, sossp, sogsp, nt, lpt, qpt, disc, classnum, indofg):
857     """
858     Determine whether the bsgs algorithm is finished or not yet.
859     This is a submodule called by the bsgs module.
860     """
861     lpt.append(n)
862     sumn = 1
863     for nid in lpt:
864         sumn = sumn * nid
865     if sumn == qpt:
866         return True, sossp, sogsp
867     elif sumn > qpt:
868         raise ValueError
869
870     if n == 1:
871         return False, sossp, sogsp
872     else:
873         tpsq = misc.primePowerTest(n)
874         if (tpsq[1] != 0) and ((tpsq[1] & 1) == 0):
875             q = arith1.floorsqrt(n)
876         else:
877             q = arith1.floorsqrt(n) + 1
878
879     ss = sossp.retel()
880     new_sossp = ClassGroup(disc, classnum, [])
881     tnt = copy.deepcopy(nt)
882     for i in range(q):
883         base = tnt ** i
884         for ssi in ss:
885             newel = base * ssi
886             if new_sossp.search(newel) is False:
887                 newel.alpha = ssi.alpha[:]
888                 lenal = len(newel.alpha)
889                 sfind = indofg - lenal
890                 for sit in range(sfind):
891                     newel.alpha.append([lenal + sit, 0, 0])
892                 newel.alpha.append([indofg, tnt, i])
893                 new_sossp.insttree(newel) # multiple of two elements of G
894
895     y = nt ** q
896     ltl = sogsp.retel()
897     new_sogsp = ClassGroup(disc, classnum, [])
898     for i in range(q + 1):
899         base = y ** (-i)
900         for eol in ltl:
901             newel2 = base * eol
902             if new_sogsp.search(newel2) is False:
903                 newel2.beta = eol.beta[:]
904                 lenbt = len(newel2.beta)
905                 gfind = indofg - lenbt
906                 for git in range(gfind):
907                     newel2.beta.append([lenbt + git, 0, 0])
908                 newel2.beta.append([indofg, tnt, q * (-i)])
909                 new_sogsp.insttree(newel2) # multiple of two elements of G
910
911     return False, new_sossp, new_sogsp
912
913 def ordercv(n, sossp, sogsp, nt, disc, classnum, tmp_ss, tmp_gs):
914     """
915     n: int
916     """
917     factor_n = misc.FactoredInteger(n)
918     while True:
919         c_s1 = ClassGroup(disc, classnum, []) # a subset of G
920         lstp_ls = sossp.retel()
921         sogsptp = sogsp.retel()
922         for p, e in factor_n:
923             base = nt ** (int(factor_n) // p)
924             for ttp_ls in lstp_ls:
925                 tmp_c_s1 = base * ttp_ls
926                 tmp_c_s1.s_parent = ttp_ls
927                 c_s1.insttree(tmp_c_s1)
928             for tmp_ell in sogsptp:
929                 rete = c_s1.search(tmp_ell)
930                 if rete != False:
931                     factor_n //= p
932                     tmp_ss = rete.s_parent
933                     tmp_gs = tmp_ell
934                     break
935             else:
936                 continue
937             break
938         else:
939             break
940     return int(factor_n), tmp_ss, tmp_gs
941
942 def giantspcv(q, sz, y, c_s1, bounds, sogsp):
943     """
944     giant step called from babyspcv.
945
946     q: int
947     sz, y: element
948     """
949     n = bounds[0]
950     # sz == x[1] ** n
951     # y  == x[1] ** q
952     while 1:
953         for tpw in sogsp.retel():
954             sz1 = sz * tpw
955             sz1.g_parent = tpw
956             rete = c_s1.search(sz1)
957             if rete is not False:
958                 return n - rete.ind, rete.s_parent, sz1.g_parent
959         # continue (sp. 5)
960         sz = sz * y
961         n = n + q
962         if n - q + 1 > bounds[1]:
963             raise ValueError("the order is larger than upper bound")
964
965 def babyspcv(bounds, sossp, sogsp, utwi, nt, disc, classnum):
966     """
967     Compute small steps
968     """
969     mstp1 = bounds[1] - bounds[0]
970     if (mstp1 == 0) or (mstp1 == 1):
971         q = 1
972     else:
973         tppm = misc.primePowerTest(mstp1)
974         q = arith1.floorsqrt(mstp1)
975         if (tppm[1] == 0) or (tppm[1] & 1):
976             q += 1
977
978     n_should_be_set = True
979     # initialize
980     c_s1 = ClassGroup(disc, classnum, []) # a subset of G
981
982     # extracting i = 0 case of main loop
983     for ttr in sossp.retel():
984         tmpx = ttr
985         tmpx.s_parent = ttr # tmpx belongs ttr in the set of smallstep
986         # index of the element
987         tmpx.ind = 0
988         c_s1.insttree(tmpx)
989
990     # main loop
991     x_i = nt
992     for i in range(1, q):
993         for ttr in sossp.retel():
994             tmpx = x_i * ttr
995             tmpx.s_parent = ttr # tmpx belongs ttr in the set of smallstep
996             if n_should_be_set and tmpx == utwi:
997                 n = i
998                 tmp_ss = tmpx.s_parent
999                 tmp_gs = utwi
1000                 n_should_be_set = False
1001             # index of the element
1002             tmpx.ind = i
1003             c_s1.insttree(tmpx)
1004         x_i = nt * x_i
1005     assert x_i == nt ** q
1006
1007     if n_should_be_set:
1008         sz = nt ** bounds[0]
1009         n, tmp_ss, tmp_gs = giantspcv(q, sz, x_i, c_s1, bounds, sogsp)
1010     return ordercv(n, sossp, sogsp, nt, disc, classnum, tmp_ss, tmp_gs)
1011
1012 def trbabysp(q, x, bounds, sossp, sogsp, ut, h, nt, disc):
1013     """
1014     Compute small steps.
1015
1016     q, h: int
1017     ut: unit element
1018     nt: element
1019     """
1020     c_s1 = ClassGroup(disc, 0, []) # a subset of G
1021     n_should_be_set = True
1022
1023     # extracting i = 0 case simplifies the main loop
1024     for tmpx in sossp.retel():
1025         tmpx.ind = 0
1026         c_s1.insttree(tmpx)
1027
1028     # main loop
1029     x_i = x[1]
1030     for i in range(1, q):
1031         for ttr in sossp.retel():
1032             tmpx = x_i * ttr
1033             if n_should_be_set and tmpx == ut:
1034                 n = i
1035                 n_should_be_set = False
1036             tmpx.ind = i
1037             c_s1.insttree(tmpx)
1038             # sort ( if you want to sort it with your estimate,
1039             # you have to implement '__ge__' method of the class with your way.)
1040         x_i = x[1] * x_i
1041
1042     if n_should_be_set:
1043         n = trgiantsp(q, x[1] ** bounds[0], x_i, c_s1, bounds, sogsp)
1044     return trorder(n * h, sossp, sogsp, nt, disc)
1045
1046 def trgiantsp(stride_index, pivot, stride, c_s1, bounds, sogsp):
1047     """
1048     Compute giant steps.
1049
1050     stride_index: int
1051     pivot, stride: element
1052     """
1053     pivot_index = bounds[0]
1054     # pivot == x[1] ** pivot_index
1055     # stride = x[1] ** stride_index
1056     while True:
1057         for tpw in sogsp.retel():
1058             rete = c_s1.search(pivot * tpw)
1059             if rete is not False:
1060                 return pivot_index - rete.ind
1061         pivot, pivot_index = pivot * stride, pivot_index + stride_index
1062         if pivot_index - stride_index + 1 > bounds[1]:
1063             raise ValueError("the order is larger than upper bound")
1064
1065 def trorder(n, sossp, sogsp, nt, disc):
1066     """
1067     Compute the order.
1068
1069     n: int
1070     nt: element
1071     """
1072     factor_n = misc.FactoredInteger(n)
1073     while True:
1074         c_s1 = ClassGroup(disc, 0, [])
1075         lstp_ls = sossp.retel()
1076         sogsptp = sogsp.retel()
1077         for p, e in factor_n:
1078             # initialize c_s1
1079             base = nt ** (int(factor_n) // p)
1080             for ttp_ls in lstp_ls:
1081                 c_s1.insttree(base * ttp_ls)
1082             # search in c_s1
1083             for tmp_ell in sogsptp:
1084                 rete = c_s1.search(tmp_ell)
1085                 if rete != False:
1086                     factor_n //= p
1087                     break
1088             else:
1089                 continue
1090             break
1091         else:
1092             return int(factor_n)
1093
1094 def isfinished_trbsgs(lwrbd, bounds, h, n, sossp, sogsp, nt, disc):
1095     """
1096     Determine whether bsgs is finished or not yet.
1097     This is a submodule called by the bsgs module.
1098
1099     lwrbd, h, n: int
1100     nt: element
1101     """
1102     h *= n
1103     if h >= lwrbd:
1104         result = True
1105     elif n == 1:
1106         result = False
1107     else:
1108         bounds[0] = (bounds[0] + n - 1) // n # ceil of lower bound // n
1109         bounds[1] = bounds[1] // n # floor of upper bound // n
1110         q = arith1.floorsqrt(n)
1111         if misc.primePowerTest(n)[1] != 2:
1112             q = arith1.floorsqrt(n) + 1 # ceil of sqrt
1113         sossp, sogsp = _update_subgrps(q, nt, sossp, sogsp, disc)
1114         result = False
1115
1116     return result, h, sossp, sogsp
1117
1118 def _update_subgrps(q, element, sossp, sogsp, discriminant):
1119     """
1120     update sossp and sogsp
1121     """
1122     new_sossp = ClassGroup(discriminant, 0, [])
1123     new_sogsp = ClassGroup(discriminant, 0, [])
1124
1125     unit = element ** 0
1126     ithpow = unit
1127     for i in range(q):
1128         for ssi in sossp.retel():
1129             newel = ithpow * ssi
1130             if new_sossp.search(newel) is False:
1131                 new_sossp.insttree(newel)
1132         ithpow = ithpow * element
1133     assert ithpow == element ** q
1134
1135     qithpow = unit
1136     for i in range(q + 1):
1137         for eol in sogsp.retel():
1138             newel = qithpow * eol
1139             if new_sogsp.search(newel) is False:
1140                 new_sogsp.insttree(newel)
1141         qithpow *= ithpow
1142
1143     return new_sossp, new_sogsp
1144
1145 def setind(n, indofg, tmp_ss, tmp_gs):
1146     """
1147     """
1148     lgtinlst = indofg
1149     if lgtinlst == 1:
1150         return [n]
1151     tmp_mt = [n]
1152     for idofel in range(lgtinlst):
1153         if idofel == 0:
1154             continue
1155         try:
1156             if type(tmp_ss.alpha[idofel][1]) != int:
1157                 ioind = tmp_ss.alpha[idofel][2]
1158             else:
1159                 ioind = 0
1160         except IndexError:
1161             ioind = 0
1162         except:
1163             raise ValueError
1164         try:
1165             if type(tmp_gs.beta[idofel][1]) != int:
1166                 joind = tmp_gs.beta[idofel][2]
1167             else:
1168                 joind = 0
1169         except IndexError:
1170             joind = 0
1171         except:
1172             raise ValueError
1173         tmp_mt.append(ioind - joind)
1174     return tmp_mt
```