Commit 86248e74 authored by Christoph Conrads's avatar Christoph Conrads

Use SymPy instead of primefac

parent 922ae864
#!/usr/bin/python2
#!/usr/bin/python3
# Copyright 2019 Christoph Conrads
#
......@@ -12,36 +12,14 @@ import collections
import fractions
import operator
import math
import primefac
import sympy.ntheory
import sys
prime_factors = {}
factors_32_17_3 = [3, 5, 17, 29, 43, 113, 127, 257, 449, 641, 5153, 2689,
65537, 6700417, 15790321, 54410972897, 183076097,
167773885276849215533569, 358429848460993,
37414057161322375957408148834323969]
ps_32_17_3 = collections.Counter(factors_32_17_3)
ps_32_17_3[2] = 32*3
prime_factors[('awc',16,16,6)] = collections.Counter([2, 65537, 1342091265,
37217928793913440210506431, 17685937523735152413434380411781231297])
factors_32_26_14 = [3, 3, 5, 7, 13, 17, 97, 193, 241, 257, 641, 673, 769,
65537, 274177, 6700417, 22253377, 18446744069414584321,
67280421310721, 442499826945303593556473164314770689]
ps_32_26_14 = collections.Counter(factors_32_26_14)
ps_32_26_14[2] = 32*14
factors_32_37_24 = [3, 5, 17, 53, 157, 257, 1613, 2731, 8191, 928513,
858001, 65537, 308761441, 18558466369, 23877647873, 21316654212673,
715668470267111297, 78919881726271091143763623681]
ps_32_37_24 = collections.Counter(factors_32_37_24)
ps_32_37_24[2] = 32*24
factors_awc_16_16_6 = [2, 65537, 1342091265, 37217928793913440210506431,
17685937523735152413434380411781231297]
ps_awc_16_16_6 = collections.Counter(factors_awc_16_16_6)
prime_factors[('awc',16,16,6)] = ps_awc_16_16_6
prime_factors[('awc',16,18,13)] = collections.Counter([
2, 3, 5, 17, 29, 257, 8419, 24606269473977504808342979,
631584071400516786177158057197336443993342659104877])
......@@ -63,7 +41,6 @@ prime_factors[('swb',8,17,7)] = collections.Counter([2, 491, 8985274326667063,
9872714575566177404203])
prime_factors[('swb',16,21,16)] = collections.Counter([
2, 2545177, 82666588531051373,
332659968319342960548379646340356862657525511747393830774592419703526959320019])
......@@ -112,18 +89,47 @@ prime_factors[('swc',24,28,8)] = collections.Counter([
11531540326642020196410750234673421368541704174300765223823234633011775])
prime_factors[('swc',24,28,8)][2] = 24*8
prime_factors[('awc',32,8,3)] = collections.Counter([2, 87956635234305,
37217928793913440210506431, 17685937523735152413434380411781231297])
prime_factors[('awc',32,16,3)] = collections.Counter([2, 3, 5, 11, 17, 29, 257,
82387, 65537, 156687811973560733, 2696785382316285340445273,
140551974055502473133117533007407559678958241229948687231152346988198330311395265590556890646801])
prime_factors[('swc',32,17,3)] = collections.Counter([3, 5, 17, 29, 43, 113,
127, 257, 449, 641, 5153, 2689, 65537, 6700417, 15790321, 54410972897,
183076097, 167773885276849215533569, 358429848460993,
37414057161322375957408148834323969])
prime_factors[('swc',32,17,3)][2] = 32*3
prime_factors[('swc',32,21,6)] = collections.Counter([ 394783681, 46908728641,
4278255361, 44479210368001, 414721, 4562284561, 61681, 65537, 22253377, 3,
3, 5, 5, 7, 11, 13, 17, 31, 41, 61, 97, 151, 193, 241, 257, 331, 673, 1321,
23041, 14768784307009061644318236958041601])
prime_factors[('swc',32,21,6)][2] = 32*6
prime_factors[('swc',32,26,14)] = collections.Counter([3, 3, 5, 7, 13, 17, 97,
193, 241, 257, 641, 673, 769, 65537, 274177, 6700417, 22253377,
18446744069414584321, 67280421310721, 442499826945303593556473164314770689])
prime_factors[('swc',32,26,14)][2] = 32*14
prime_factors[('awc',32,29,17)] = collections.Counter([2, 3, 5, 17, 23, 257,
65537, 7165729, 16261579, 296174737, 9286471429,
35834095934305889246278160558607808213536611909101400963631213603354576647919834164787160509080753251812869864745539853004298715922465713710357103697213980386700856246837790064184164935087482021935440410696100869139134472551164881019937])
prime_factors[('swc',32,34,19)] = collections.Counter([ 4278255361, 22253377,
394783681, 65537, 46908728641, 44479210368001,
14768784307009061644318236958041601, 3, 3, 5, 5, 7, 11, 13, 17, 31, 41, 61,
97, 151, 193, 241, 257, 331, 673, 1321, 23041, 414721, 61681, 4562284561])
prime_factors[('swc',32,34,19)][2] = 32*19
prime_factors[('swc',32,37,24)] = collections.Counter([3, 5, 17, 53, 157, 257,
1613, 2731, 8191, 928513, 858001, 65537, 308761441, 18558466369,
23877647873, 21316654212673, 715668470267111297,
78919881726271091143763623681])
prime_factors[('swc',32,37,24)][2] = 32*24
prime_factors[('swc',64,13,7)] = collections.Counter([3, 3, 5, 7, 13, 17, 97,
193, 241, 257, 641, 673, 769, 274177, 6700417, 22253377, 65537,
......@@ -143,24 +149,8 @@ prime_factors[('swc',64,26,4)] = collections.Counter([3, 5, 17, 23, 89, 257, 353
prime_factors[('swc',64,26,4)][2] = 64*4
factors_awc_32_8_3 = [2, 87956635234305, 37217928793913440210506431,
17685937523735152413434380411781231297]
ps_awc_32_8_3 = collections.Counter(factors_awc_32_8_3)
factors_awc_32_16_3 = [2, 3, 5, 11, 17, 29, 257, 82387, 65537,
156687811973560733, 2696785382316285340445273,
140551974055502473133117533007407559678958241229948687231152346988198330311395265590556890646801]
ps_awc_32_16_3 = collections.Counter(factors_awc_32_16_3)
factors_awc_32_29_17 = [2, 3, 5, 17, 23, 257, 65537, 7165729, 16261579,
296174737, 9286471429,
35834095934305889246278160558607808213536611909101400963631213603354576647919834164787160509080753251812869864745539853004298715922465713710357103697213980386700856246837790064184164935087482021935440410696100869139134472551164881019937]
ps_awc_32_29_17 = collections.Counter(factors_awc_32_29_17)
def factorize(kind, w, r, s, m):
f = lambda n: primefac.factorint(n, methods=(primefac.pollardRho_brent,))
f = lambda n: sympy.ntheory.factorint(n, limit=10*1000)
key = (kind, w, r, s)
if key in prime_factors:
......@@ -174,16 +164,13 @@ def factorize(kind, w, r, s, m):
n = n // p**multiplicity
if n == 1 or n > 10**30:
return ps
qs = f(n)
factors = collections.Counter(qs) + ps
elif kind == 'swc' and w*(r-s)/2 < 200:
elif kind == 'swc':
b = 2**w
np1 = (2**(w*(r-s)/2) + 1)
nm1 = (2**(w*(r-s)/2) - 1)
np1 = (2**(w*(r-s)//2) + 1)
nm1 = (2**(w*(r-s)//2) - 1)
assert m % b**s == 0
assert m % np1 == 0
......@@ -199,10 +186,8 @@ def factorize(kind, w, r, s, m):
assert 2 not in factors
factors[2] = s*w
elif w*r < 300:
factors = f(m)
else:
return None
factors = f(m)
for f in factors:
assert m % f == 0
......@@ -236,11 +221,14 @@ def compute_order_mod_m(a, k, m, ps):
for p in sorted(ps.keys()):
assert k % p == 0
if not sympy.ntheory.isprime(p):
continue
num = ps[p]
product = product * p**num
for n in xrange(num):
remainder = pow(a, k/p, m)
for n in range(num):
remainder = pow(a, k//p, m)
if remainder == 1:
k = k // p
......@@ -279,17 +267,7 @@ def main():
m = compute_modulus(b, r, s)
key = (kind, w, r, s)
ps = \
ps_32_17_3 if kind == 'swc' and w == 32 and r == 17 and s == 3 else \
ps_32_26_14 if kind == 'swc' and w == 32 and r == 26 and s == 14 else \
ps_32_37_24 if kind == 'swc' and w == 32 and r == 37 and s == 24 else \
ps_awc_32_8_3 if kind== 'awc' and w == 32 and r == 8 and s == 3 else \
ps_awc_32_16_3 if kind=='awc' and w == 32 and r == 16 and s == 3 else \
ps_awc_32_29_17 if kind=='awc' and w== 32 and r == 29 and s == 17 else \
factorize(kind, w, r, s, m-1)
assert ps is None or (m-1) % reduce(operator.mul, ps, 1) == 0
ps = factorize(kind, w, r, s, m-1)
min_period, max_period = \
(1, m-1) if ps is None else compute_order_mod_m(b, m-1, m, ps)
......
......@@ -72,7 +72,7 @@ make_luscher_plot()
echo -n "$stats" ' '
local period="$( \
python2 "$cwd/compute-awc-swb-period-length.py" $type $w $r $s)"
python3 "$cwd/compute-awc-swb-period-length.py" $type $w $r $s)"
echo "$period"
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment