"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 
   12 class ReducedQuadraticForm(object):
   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)
  129         self.elements.append(newestl)
  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)
  254     ret_list.append(ReducedQuadraticForm(unit, unit))
  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 
  277 class ReducedQuadraticFormForBSGS(ReducedQuadraticForm):
  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):
  427                 self.add_to_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