Commit 00197987 authored by Christoph Conrads's avatar Christoph Conrads

Fully automated AWC, SWB parameter search

parent f2f37eab
#!/bin/bash
#!/usr/bin/python2
# Copyright 2019 Christoph Conrads
#
......@@ -6,95 +6,276 @@
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
from __future__ import print_function
import collections
import fractions
import operator
import math
import primefac
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
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_64_26_4 = [3, 5, 17, 23, 89, 257, 353, 397, 641, 683, 1409, 2113,
65537, 229153, 5304641, 274177, 119782433, 43872038849, 1258753,
6700417, 2931542417, 67280421310721, 441995541378330835457,
60299259845689822028046342401,
275509565477848842604777623828011666349761,
3210843755324367119258027752661239735297]
ps_64_26_4 = collections.Counter(factors_64_26_4)
ps_64_26_4[2] = 64*4
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])
prime_factors[('awc', 32, 33, 13)] = collections.Counter([
2, 65537, 3, 5, 17, 257, 331, 46337, 259723, 812864267, 612985394553226439])
prime_factors[('awc', 32, 33, 27)] = collections.Counter([
2, 3, 3, 5, 5, 7, 13, 17, 97, 193, 241, 257, 673, 937, 45139, 9918151937,
65537, 8985577, 22253377])
prime_factors[('awc', 32, 39, 11)] = collections.Counter([
2, 3, 5, 5, 17, 109, 257, 913676655973])
prime_factors[('swb',16,21,16)] = collections.Counter([
2, 2545177, 82666588531051373,
332659968319342960548379646340356862657525511747393830774592419703526959320019])
prime_factors[('swb',16,25,16)] = collections.Counter([
2, 16283821, 3618769732453, 363017649779303260447])
prime_factors[('swb',16,27,19)] = collections.Counter([
2,5084251061237, 249186967913293, 72542073395875053228967309])
prime_factors[('swb',32,19,16)] = collections.Counter([2, 2713, 20507, 3744619])
prime_factors[('swb',32,34,9)] = collections.Counter([2, 569, 641, 1433, 6700417])
prime_factors[('swb',32,34,28)] = collections.Counter([2, 883, 118438939])
prime_factors[('swc',16,24,5)] = collections.Counter([3, 5, 17, 229, 257, 457,
27361, 1217, 148961, 174763, 524287, 525313, 24517014940753,
69394460463940481, 11699557817717358904481])
prime_factors[('swc',16,24,5)][2] = 16*5
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,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',64,13,7)] = collections.Counter([3, 3, 5, 7, 13, 17, 97,
193, 241, 257, 641, 673, 769, 274177, 6700417, 22253377, 65537,
18446744069414584321, 67280421310721, 442499826945303593556473164314770689])
prime_factors[('swc',64,13,7)][2] = 64*7
def main():
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
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_awc_32_8_3 = [2, 87956635234305, 37217928793913440210506431,
17685937523735152413434380411781231297]
ps_awc_32_8_3 = collections.Counter(factors_awc_32_8_3)
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_32_8_3 = [2, 87956635234305, 37217928793913440210506431,
17685937523735152413434380411781231297]
ps_awc_32_8_3 = collections.Counter(factors_awc_32_8_3)
factors_64_26_4 = [3, 5, 17, 23, 89, 257, 353, 397, 641, 683, 1409, 2113,
65537, 229153, 5304641, 274177, 119782433, 43872038849, 1258753,
6700417, 2931542417, 67280421310721, 441995541378330835457,
60299259845689822028046342401,
275509565477848842604777623828011666349761,
3210843755324367119258027752661239735297]
ps_64_26_4 = collections.Counter(factors_64_26_4)
ps_64_26_4[2] = 64*4
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_8_3 = [2, 87956635234305, 37217928793913440210506431,
17685937523735152413434380411781231297]
ps_awc_32_8_3 = collections.Counter(factors_awc_32_8_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)
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,))
key = (kind, w, r, s)
e = 64
r = 26
s = 4
b = 2**e
if key in prime_factors:
ps = prime_factors[key]
n = m
m = b**r - b**s + 1
for p in ps:
multiplicity = ps[p]
print('b=2^{:d} r={:d} s={:d}'.format(e, r, s))
assert n % p**multiplicity == 0
d = m-1
n = n // p**multiplicity
ps = ps_64_26_4
q = reduce(operator.mul, factors_64_26_4, 1)
if n == 1 or n > 10**30:
return ps
print(((m-1) / b**s) % q)
print(((m-1) / b**s) // q)
qs = f(n)
factors = collections.Counter(qs) + ps
print('Found {:d} prime factors'.format(len(ps)))
elif kind == 'swc':
b = 2**w
np1 = (2**(w*(r-s)/2) + 1)
nm1 = (2**(w*(r-s)/2) - 1)
#assert reduce(operator.mul, ps, 1) == m-1
print((m-1) % (q * b**s))
print((m-1) // (q * b**s))
assert m % b**s == 0
assert m % np1 == 0
assert m % nm1 == 0
assert b**s * np1 * nm1 == m
assert np1 % 2 == 1
assert nm1 % 2 == 1
fp1 = f(np1)
fm1 = f(nm1)
factors = collections.Counter(fp1) + collections.Counter(fm1)
assert 2 not in factors
factors[2] = s*w
elif w < 33:
factors = f(m)
else:
return None
for f in factors:
assert m % f == 0
return factors
# @param[in] ps factors of k. ideally these are prime or you will get a lower
# bound on the order of a mod m.
def compute_order_mod_m(a, k, m, ps):
assert isinstance(a, int) or isinstance(a, long)
assert isinstance(k, int) or isinstance(k, long)
assert isinstance(m, int) or isinstance(m, long)
assert a > 0
assert k > 0
assert m > a
assert isinstance(ps, dict)
if fractions.gcd(k, m) != 1:
raise ArithmeticError('gcd(k,m) != 1')
k_0 = k
product = 1
# https://math.stackexchange.com/questions/1025578/is-there-a-better-way-of-finding-the-order-of-a-number-modulo-n
for p in sorted(ps.keys()):
for n in xrange(ps[p]):
r = pow(b, d/p, m)
assert k % p == 0
num = ps[p]
product = product * p**num
for n in xrange(num):
remainder = pow(a, k/p, m)
if remainder == 1:
k = k // p
assert product <= k_0
assert k_0 % product == 0
max_order = k
min_order = k
if k_0 != product:
div = k_0 // product
min_order = k // div
return [min_order, max_order]
def main():
if len(sys.argv) < 5:
sys.exit(1)
kind = sys.argv[1]
w = int(sys.argv[2])
r = int(sys.argv[3])
s = int(sys.argv[4])
b = 2**w
compute_modulus = \
(lambda b, r, s: b**r + b**s - 1) if kind == 'awc' else \
(lambda b, r, s: b**r + b**s + 1) if kind == 'cawc' else \
(lambda b, r, s: b**r - b**s - 1) if kind == 'swb' else \
(lambda b, r, s: b**r - b**s + 1) if kind == 'swc' else \
sys.exit(3)
m = compute_modulus(b, r, s)
#print('b=2^{:d} r={:d} s={:d}'.format(w, 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_64_26_4 if kind == 'swc' and w == 64 and r == 26 and s == 4 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)
#print('Found {:d} prime factors'.format(len(ps)))
assert ps is None or (m-1) % reduce(operator.mul, ps, 1) == 0
if r == 1:
d = d // p
min_period, max_period = \
(1, m-1) if ps is None else compute_order_mod_m(b, m-1, m, ps)
p10_min = math.log10(min_period)
p10_max = math.log10(max_period)
period10 = math.log10(d)
period2 = math.log(d) / math.log(2)
p2_min = math.log(min_period) / math.log(2)
p2_max = math.log(max_period) / math.log(2)
print('Period length 2^{:.2f} 10^{:.2f}'.format(period2, period10))
#print('Period length 2^{:.2f} 2^{:.2f}'.format(p2_min, p2_max))
#print('Period length 10^{:.2f} 10^{:.2f}'.format(p10_min, p10_max))
msg = '{:7.2f} {:7.2f} {:6.2f} {:6.2f}'
print(msg.format(p2_min, p2_max, p10_min, p10_max))
if __name__ == '__main__':
......
......@@ -37,7 +37,7 @@ if min(xs) < 0 or max(xs) >= 1:
nnz, = np.nonzero(xs > 0.45)
k = -1 if nnz.size == 0 else min(nnz)
k_str = '{:2d}+'.format(xs.size) if k == -1 else '{:2d}'.format(k)
k_str = '>{:2d}'.format(xs.size) if k == -1 else '{:3d}'.format(k)
print('{:4s} {:d} {:2d} {:2d} {:s}'.format(type, w, r, s, k_str))
......@@ -118,8 +118,15 @@ make_luscher_plot()
configure_cxx_code $type $w $r $s >"$cxx_filename"
c++ -Wextra -Wall -std=c++11 -pedantic "$cxx_filename" -o "$binary"
"./$binary" | awk '{print $2}' | \
python3 -c "$python_plot_code" $type $w $r $s
local stats="$("./$binary" | awk '{print $2}' | \
python3 -c "$python_plot_code" $type $w $r $s)"
echo -n "$stats" ' '
local period="$( \
python2 "$cwd/compute-awc-swb-period-length.py" $type $w $r $s)"
echo "$period"
}
......@@ -159,7 +166,7 @@ find_awc_swb_prngs()
echo "Found $num_generators $type parameter combinations"
echo 'type b r s t'
echo 'type b r s t log2period log10period'
while read gen
do
local b="$(awk '{print $2}' <<<"$gen")"
......
......@@ -24,7 +24,7 @@ compute_modulus = \
sys.exit(3)
for e in [16, 32, 64]:
for e in [16,32,64]:
b = 2**e
r_max = \
30 if e == 16 else \
......@@ -32,7 +32,7 @@ for e in [16, 32, 64]:
50 if e == 64 else -1
for r in range(3, r_max):
s_max = r-1 if r < 10 else r-2
s_max = r-1 if r < 10 else r-2 if r < 20 else r-3
for s in range(2, s_max):
assert r > s
......
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