NZMATH  1.2.0
About: NZMATH is a Python based number theory oriented calculation system.
  Fossies Dox: NZMATH-1.2.0.tar.gz  ("inofficial" and yet experimental doxygen-generated source code documentation)  

quad.py
Go to the documentation of this file.
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 
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 
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 
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)
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
nzmath.quad.ReducedQuadraticFormForBSGS.beta
beta
Definition: quad.py:285
nzmath.quad.ReducedQuadraticForm.__init__
def __init__(self, element, unit)
Definition: quad.py:16
nzmath.quad.ClassGroup.inststree
def inststree(self, newlist)
Definition: quad.py:131
nzmath.quad.sqrPDF
def sqrPDF(f)
Definition: quad.py:510
nzmath.bigrandom
Definition: bigrandom.py:1
nzmath.quad.unit_form
def unit_form(disc)
Definition: quad.py:627
nzmath.quad.rand_generator
def rand_generator(disc, classnum, qin)
Definition: quad.py:722
nzmath.bigrange.range
def range(start, stop=None, step=None)
Definition: bigrange.py:19
nzmath.quad.ClassGroup.search
def search(self, tarel)
Definition: quad.py:159
nzmath.quad.ClassGroup.rootoftree
rootoftree
Definition: quad.py:105
nzmath.quad.trgiantsp
def trgiantsp(stride_index, pivot, stride, c_s1, bounds, sogsp)
Definition: quad.py:1046
nzmath.quad.ReducedQuadraticFormForBSGS.__init__
def __init__(self, element, unit)
Definition: quad.py:281
nzmath.quad.class_formula
def class_formula(disc, uprbd)
Definition: quad.py:186
nzmath.quad.ClassGroup.retel
def retel(self)
Definition: quad.py:170
nzmath.factor.misc
Definition: misc.py:1
nzmath.quad.ReducedQuadraticForm.__mul__
def __mul__(self, other)
Definition: quad.py:28
nzmath.quad.ClassGroup.inserel
def inserel(self, newel)
Definition: quad.py:127
nzmath.quad.RetNext._unit
_unit
Definition: quad.py:409
nzmath.quad.RetNext.has_in_hash
def has_in_hash(self, form)
Definition: quad.py:453
nzmath.quad.disc
def disc(f)
Definition: quad.py:463
nzmath.quad.setind
def setind(n, indofg, tmp_ss, tmp_gs)
Definition: quad.py:1145
nzmath.quad.class_group
def class_group(disc, limit_of_disc=100000000)
Definition: quad.py:229
nzmath.quad.ReducedQuadraticForm.__str__
def __str__(self)
Definition: quad.py:24
nzmath.quad.ReducedQuadraticForm.unit
unit
Definition: quad.py:18
nzmath.quad.RetNext.utroot
utroot
Definition: quad.py:408
nzmath.gcd
Definition: gcd.py:1
nzmath.quad.isfundamental
def isfundamental(disc)
Definition: quad.py:784
nzmath.quad.kronecker
def kronecker(a, b)
Definition: quad.py:643
nzmath.quad.ReducedQuadraticFormForBSGS
Definition: quad.py:277
nzmath.quad.ClassGroup.inserels
def inserels(self, newlist)
Definition: quad.py:123
nzmath.quad.class_group_bsgs
def class_group_bsgs(disc, classnum, qin)
Definition: quad.py:343
nzmath.quad.ReducedQuadraticForm.__truediv__
def __truediv__(self, other)
Definition: quad.py:40
nzmath.quad.ClassGroup.__repr__
def __repr__(self)
Definition: quad.py:120
nzmath.quad.ReducedQuadraticFormForBSGS.s_parent
s_parent
Definition: quad.py:286
nzmath.quad.RetNext.add_to_hash
def add_to_hash(self, form)
Definition: quad.py:443
nzmath.quad.ReducedQuadraticForm.element
element
Definition: quad.py:17
nzmath.quad.sqroot
def sqroot(disc, p)
Definition: quad.py:736
nzmath.quad.RetNext.elhash
elhash
Definition: quad.py:412
nzmath.quad.ReducedQuadraticForm
Definition: quad.py:12
nzmath.quad.RetNext.__init__
def __init__(self, disc)
Definition: quad.py:406
nzmath.quad.ReducedQuadraticForm.__repr__
def __repr__(self)
Definition: quad.py:20
nzmath.quad.isfinished_trbsgs
def isfinished_trbsgs(lwrbd, bounds, h, n, sossp, sogsp, nt, disc)
Definition: quad.py:1094
nzmath.quad.RetNext._randomele1
def _randomele1(self)
Definition: quad.py:430
nzmath.quad.ReducedQuadraticFormForBSGS.g_parent
g_parent
Definition: quad.py:287
nzmath.quad.ClassGroup.classnum
classnum
Definition: quad.py:111
nzmath.quad.giantspcv
def giantspcv(q, sz, y, c_s1, bounds, sogsp)
Definition: quad.py:942
nzmath.quad.RetNext
following functions are sub functions for above functions.
Definition: quad.py:402
nzmath.quad.ReducedQuadraticForm.__hash__
def __hash__(self)
Definition: quad.py:63
nzmath.quad.compositePDF
def compositePDF(f_1, f_2)
Definition: quad.py:583
nzmath.quad.RetNext.previous
previous
Definition: quad.py:411
nzmath.quad.ReducedQuadraticFormForBSGS.ind
ind
Definition: quad.py:283
nzmath.quad._update_subgrps
def _update_subgrps(q, element, sossp, sogsp, discriminant)
Definition: quad.py:1118
nzmath.quad.ReducedQuadraticForm.repOfModule
def repOfModule(self)
Definition: quad.py:77
nzmath.quad.RetNext.unit
def unit(self)
Definition: quad.py:414
nzmath.quad.ClassGroup
Definition: quad.py:98
nzmath.quad.randomele
def randomele(disc, unit)
Definition: quad.py:760
nzmath.quad.RetNext.cnt
cnt
Definition: quad.py:410
nzmath.quad.isfinished_bsgscv
def isfinished_bsgscv(n, sossp, sogsp, nt, lpt, qpt, disc, classnum, indofg)
Definition: quad.py:856
nzmath.arith1
Definition: arith1.py:1
nzmath.quad.ClassGroup.expunit
expunit
Definition: quad.py:112
nzmath.quad.powPDF
def powPDF(f, exp, ut=None)
Definition: quad.py:552
nzmath.quad.ReducedQuadraticForm.__getitem__
def __getitem__(self, index)
Definition: quad.py:91
nzmath.quad.ClassGroup.disc
disc
Definition: quad.py:104
nzmath.quad.euclid_exd
def euclid_exd(a, b)
Definition: quad.py:801
nzmath.quad.ClassGroup.__init__
def __init__(self, disc, classnum, elements=None)
Definition: quad.py:102
nzmath.quad.babyspcv
def babyspcv(bounds, sossp, sogsp, utwi, nt, disc, classnum)
Definition: quad.py:965
nzmath.quad.RetNext.retnext
def retnext(self)
Definition: quad.py:420
nzmath.prime
Definition: prime.py:1
nzmath.quad.ReducedQuadraticForm.__pow__
def __pow__(self, exp)
Definition: quad.py:34
nzmath.quad.crt
def crt(inlist)
Definition: quad.py:693
nzmath.quad.trbabysp
def trbabysp(q, x, bounds, sossp, sogsp, ut, h, nt, disc)
Definition: quad.py:1012
nzmath.quad.class_number_bsgs
def class_number_bsgs(disc)
Definition: quad.py:290
nzmath.quad.ReducedQuadraticForm.disc
def disc(self)
Definition: quad.py:71
nzmath.quad.ReducedQuadraticForm.__ne__
def __ne__(self, other)
Definition: quad.py:68
nzmath.quad.parteucl
def parteucl(a, b, sogsp)
Definition: quad.py:829
nzmath.quad.ClassGroup.rootornot
rootornot
Definition: quad.py:106
nzmath.quad.ReducedQuadraticForm.__iter__
def __iter__(self)
Definition: quad.py:85
nzmath.quad.ClassGroup.elements
elements
Definition: quad.py:108
nzmath.quad.reducePDF
def reducePDF(f)
Definition: quad.py:475
nzmath.quad.ClassGroup.__str__
def __str__(self)
Definition: quad.py:114
nzmath.quad.trorder
def trorder(n, sossp, sogsp, nt, disc)
Definition: quad.py:1065
nzmath.quad.class_number
def class_number(disc, limit_of_disc=100000000)
Definition: quad.py:198
nzmath.quad.ordercv
def ordercv(n, sossp, sogsp, nt, disc, classnum, tmp_ss, tmp_gs)
Definition: quad.py:913
nzmath.quad.ClassGroup.insttree
def insttree(self, newel0)
Definition: quad.py:135
nzmath.quad.ReducedQuadraticForm.inverse
def inverse(self)
Definition: quad.py:45
nzmath.quad.number_unit
def number_unit(disc)
Definition: quad.py:680
nzmath.quad.ReducedQuadraticForm.__eq__
def __eq__(self, other)
Definition: quad.py:55
nzmath.quad.ReducedQuadraticFormForBSGS.alpha
alpha
Definition: quad.py:284
nzmath.quad.RetNext.disc
disc
Definition: quad.py:407