chiark - git - mdw - catacomb/blob - utils/findpock.sage

   1 #! /usr/local/bin/sage
   2
   3 from sys import argv
   4 from itertools import combinations
   5
   6 sievebits = 32
   7
   8 SEQ = 0
   9 LABEL = {}
  10
  11 def label(p):
  12   global SEQ
  13   try: lab = LABEL[p]
  14   except KeyError:
  15     if p.nbits() < sievebits: lab = LABEL[p] = '%d' % p
  16     else:
  17       lab = LABEL[p] = '$%s.%d' % (name, SEQ)
  18       SEQ += 1
  19   return lab
  20
  21 SMALL = set()
  22 POCK = []
  23 DONE = {}
  24
  25 def pock(p, rr):
  26   r = prod(rr)
  27   ll = map(recurse, rr)
  28   lab = DONE[p] = label(p)
  29   a = 2
  30   while True:
  31     if pow(a, p - 1, p) != 1:
  32       raise ValueError('%d not prime (%d is Fermat witness)' % (p, a))
  33     win = True
  34     for q in rr:
  35       g = gcd(pow(a, (p - 1)/q, p) - 1, p)
  36       if 1 < g < p: raise ValueError('%d not prime (%d divides)' % (p, g))
  37       if g != 1: win = False
  38     if win: break
  39     a += 1
  40   POCK.append('pock %s = %d, %d, [%s]' %
  41               (lab, a, (p - 1)/(2*r), ', '.join(ll)))
  42   return lab
  43
  44 def recurse(p):
  45   try: return DONE[p]
  46   except KeyError: pass
  47
  48   if p.nbits() < sievebits:
  49     lab = DONE[p] = str(p)
  50     SMALL.add(p)
  51     return lab
  52
  53   best, score = None, p
  54   qq = [q for (q, e) in factor((p - 1)/2)]
  55   for n in xrange(1, len(qq) + 1):
  56     for rr in combinations(qq, n):
  57       r = prod(rr)
  58       if r^2 <= p: continue
  59       if r < score: best, score = rr, r
  60
  61   best = list(best); best.sort()
  62   return pock(p, best)
  63
  64 name, p = argv[1], Integer(argv[2])
  65 if len(argv) == 3:
  66   LABEL[p] = name
  67   recurse(p)
  68   for q in sorted(SMALL): print 'small %d = %d' % (q, q)
  69   print
  70 else:
  71   qq = map(Integer, argv[3:])
  72   for i in xrange(len(qq)): LABEL[qq[i]] = DONE[qq[i]] = 'q%d' % i
  73   q = prod(qq)
  74   if not p: p = 2*q + 1
  75   elif p%(2*q) != 1: raise ValueError('incorrect factorization')
  76   if q^2 <= p: raise ValueError('factorization insufficient')
  77   LABEL[p] = name
  78   pock(p, qq)
  79
  80 for line in POCK: print line
  81 print 'check %s, %d, %d' % (name, p.nbits(), p)