diff --git a/compute-awc-swb-period-length.py b/compute-awc-swb-period-length.py index e1721bc0788ab821a8385832973fa08dcc765f7a..f742448683e70d677ea761ac79968185cdd9d653 100644 --- a/compute-awc-swb-period-length.py +++ b/compute-awc-swb-period-length.py @@ -1,4 +1,4 @@ -#!/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__': diff --git a/find-awc-swb-prngs.sh b/find-awc-swb-prngs.sh index f4fa0bcd2f72c45fa8145c728c1f2ab8d1c88ef2..d2cf71e3c62cbaece05b6e4aecc9168cd786dd59 100644 --- a/find-awc-swb-prngs.sh +++ b/find-awc-swb-prngs.sh @@ -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")" diff --git a/generate-awc-swb-moduli.py b/generate-awc-swb-moduli.py index 2df8695cf06194d09af3e107a3121c652475ffa9..ccff414e6e14f280e04d69c3b5e7ea1d44e9166d 100644 --- a/generate-awc-swb-moduli.py +++ b/generate-awc-swb-moduli.py @@ -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