Commit f2f37eab authored by Christoph Conrads's avatar Christoph Conrads

Automate search for AWC, SWB parameters

parent dc59ea34
......@@ -9,7 +9,6 @@
#include <cstddef>
#include <cstdio>
#include <limits>
#define private public
#include <random>
#include <type_traits>
#include <vector>
......@@ -25,8 +24,7 @@ struct add_with_carry_engine
static_assert(r > s, "");
// prototype implementation
static_assert(w == 32u, "");
static_assert(std::is_same<T, std::uint32_t>::value, "");
static_assert(w == std::numeric_limits<T>::digits, "");
using result_type = T;
......@@ -54,20 +52,18 @@ struct add_with_carry_engine
{
auto i = index_;
auto j = index_ >= s ? index_ - s : index_ + r - s;
auto u = xs_[i];
auto v = xs_[j];
assert(carry_ == 0 or carry_ == 1);
auto x64 = std::uint64_t{u} + v + carry_;
auto x = std::uint32_t(x64);
auto c = x64 >> 32;
auto a = T{carry_};
auto b = T(a + xs_[i]);
auto c = T(b + xs_[j]);
xs_[i] = x;
carry_ = c;
xs_[i] = c;
carry_ = b < a or c < b ? 1u : 0u;
index_ = index_ + 1u == r ? 0u : index_ + 1u;
return x;
return c;
}
......@@ -96,7 +92,7 @@ std::size_t distance(T m, const T (&u)[N], const T (&v)[N])
auto p = u[i] >= v[i] ? u[i]-v[i] : v[i]-u[i];
auto q = std::min(p, m-p);
d = std::max(d, q);
d = std::max(T(d), T(q));
}
return d;
......@@ -109,14 +105,23 @@ int main()
constexpr auto s = std::size_t{@SHORT_LAG@};
constexpr auto r = std::size_t{@LONG_LAG@};
using Generator = add_with_carry_engine<std::uint32_t, w, s, r>;
using T = Generator::result_type;
using T = std::conditional<
w == 16, std::uint16_t, std::conditional<
w == 32, std::uint32_t, std::conditional<
w == 64, std::uint64_t, void
>::type>::type>::type;
using Generator = add_with_carry_engine<T, w, s, r>;
constexpr auto num_trajectories = 1000u;
constexpr auto t_max = 40u;
constexpr auto t_max =
w == 16u ? 30u :
w == 32u ? 40u :
w == 64u ? 50u : 0u
;
constexpr auto M = w == std::numeric_limits<T>::digits
? std::numeric_limits<T>::max()
: T{1} << Generator::word_size
: T(T{1} << Generator::word_size)
;
auto seed_gen = std::mt19937();
......
......@@ -9,12 +9,77 @@
#include <cstddef>
#include <cstdio>
#include <limits>
#define private public
#include <random>
#include <type_traits>
#include <vector>
template<typename T, std::size_t w, std::size_t r, std::size_t s>
struct subtract_with_borrow
{
static_assert(std::is_integral<T>::value, "");
static_assert(not std::is_signed<T>::value, "");
static_assert(w <= std::numeric_limits<T>::digits, "");
static_assert(s > 0u, "");
static_assert(r > s, "");
// prototype implementation
static_assert(w == std::numeric_limits<T>::digits, "");
using result_type = T;
static constexpr auto long_lag = r;
static constexpr auto short_lag = s;
static constexpr auto word_size = w;
static constexpr T max() { return std::numeric_limits<T>::max(); }
static constexpr T min() { return std::numeric_limits<T>::min(); }
explicit subtract_with_borrow(unsigned seed=123456u)
{
auto gen = std::minstd_rand(seed);
for(auto i = std::size_t{0}; i < r; ++i)
xs_[i] = gen();
discard(r);
}
T operator() ()
{
auto i = index_;
auto j = index_ >= s ? index_ - s : index_ + r - s;
assert(carry_ == 0 or carry_ == 1);
auto x = xs_[i];
auto y = T(x - xs_[j]);
auto z = T(y - carry_);
xs_[i] = y;
carry_ = y > x or z > y ? 1 : 0;
index_ = index_ + 1u == r ? 0u : index_ + 1u;
return z;
}
void discard(std::size_t n)
{
for(auto i = std::size_t{0}; i < n; ++i)
(*this)();
}
T xs_[r] = { 0 };
std::uint8_t carry_ = 0;
std::size_t index_ = 0;
};
template<typename T, std::size_t N>
std::size_t distance(T m, const T (&u)[N], const T (&v)[N])
{
......@@ -27,7 +92,7 @@ std::size_t distance(T m, const T (&u)[N], const T (&v)[N])
auto p = u[i] >= v[i] ? u[i]-v[i] : v[i]-u[i];
auto q = std::min(p, m-p);
d = std::max(d, q);
d = std::max(T(d), T(q));
}
return d;
......@@ -40,15 +105,23 @@ int main()
constexpr auto s = std::size_t{@SHORT_LAG@};
constexpr auto r = std::size_t{@LONG_LAG@};
using Generator =
std::subtract_with_carry_engine<std::uint64_t, w, s, r>;
using T = Generator::result_type;
using T = std::conditional<
w == 16, std::uint16_t, std::conditional<
w == 32, std::uint32_t, std::conditional<
w == 64, std::uint64_t, void
>::type>::type>::type;
using Generator = subtract_with_borrow<T, w, r, s>;
constexpr auto num_trajectories = 1000u;
constexpr auto t_max = 40u;
constexpr auto t_max =
w == 16u ? 30u :
w == 32u ? 40u :
w == 64u ? 50u : 0u
;
constexpr auto M = w == std::numeric_limits<T>::digits
? std::numeric_limits<T>::max()
: T{1} << Generator::word_size
: T(T{1} << Generator::word_size)
;
auto seed_gen = std::mt19937();
......@@ -60,11 +133,11 @@ int main()
auto gen_1 = Generator(seed);
auto gen_2 = Generator(seed);
gen_2._M_x[it%r] = (gen_1._M_x[it%r] + 1u) % M;
gen_2.xs_[it%r] = (gen_1.xs_[it%r] + 1u) % M;
for(auto t = 0u; t < t_max; ++t)
{
auto d = distance(M, gen_1._M_x, gen_2._M_x);
auto d = distance(M, gen_1.xs_, gen_2.xs_);
ds[t] += 1.0 * d / M;
......
// Copyright 2019 Christoph Conrads
//
// This Source Code Form is subject to the terms of the Mozilla Public
// 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/.
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <limits>
#define private public
#include <random>
#include <type_traits>
#include <vector>
template<typename T, std::size_t N>
std::size_t distance(T m, const T (&u)[N], const T (&v)[N])
{
static_assert(not std::is_signed<T>::value, "");
auto d = T{0};
for(auto i = std::size_t{0}; i < N; ++i)
{
auto p = u[i] >= v[i] ? u[i]-v[i] : v[i]-u[i];
auto q = std::min(p, m-p);
d = std::max(d, q);
}
return d;
}
int main()
{
constexpr auto w = std::size_t{@WORD_SIZE@};
constexpr auto s = std::size_t{@SHORT_LAG@};
constexpr auto r = std::size_t{@LONG_LAG@};
using Generator =
std::subtract_with_carry_engine<std::uint64_t, w, s, r>;
using T = Generator::result_type;
constexpr auto num_trajectories = 1000u;
constexpr auto t_max =
w == 16u ? 30u :
w == 32u ? 40u :
w == 64u ? 50u : 0u
;
constexpr auto M = w == std::numeric_limits<T>::digits
? std::numeric_limits<T>::max()
: T{1} << Generator::word_size
;
auto seed_gen = std::mt19937();
auto ds = std::vector<double>(t_max);
for(auto it = 0u; it < num_trajectories; ++it)
{
auto seed = seed_gen();
auto gen_1 = Generator(seed);
auto gen_2 = Generator(seed);
gen_2._M_x[it%r] = (gen_1._M_x[it%r] + 1u) % M;
for(auto t = 0u; t < t_max; ++t)
{
auto d = distance(M, gen_1._M_x, gen_2._M_x);
ds[t] += 1.0 * d / M;
gen_1.discard(r);
gen_2.discard(r);
}
}
for(auto t = std::size_t{0}; t < ds.size(); ++t)
{
std::printf("%2zu %23.17e\n", t, ds[t] / num_trajectories);
}
}
......@@ -39,7 +39,7 @@ 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)
print('{:s} {:d} {:2d} {:2d} {:s}'.format(type, w, r, s, k_str))
print('{:4s} {:d} {:2d} {:2d} {:s}'.format(type, w, r, s, k_str))
pp.plot(xs)
pp.grid(which="both")
......@@ -50,10 +50,10 @@ pp.yscale('log', basey=2)
pp.xlabel("t")
pp.ylabel("d(x(t),y(t))")
if type == 'awc':
pp.title("{:s}(2^{:d}, {:d}, {:d})".format(type.upper(), w, r, s))
if type == 'swc':
pp.title("{:s}(2^{:d}, {:d}, {:d})".format('SWB', w, s, r))
else:
pp.title("{:s}(2^{:d}, {:d}, {:d})".format(type.upper(), w, s, r))
pp.title("{:s}(2^{:d}, {:d}, {:d})".format(type.upper(), w, r, s))
pp.savefig(filename)
EOF
......@@ -72,8 +72,15 @@ is_modulus_prime()
if [ "$type" = awc ]
then
local modulus="$(python -c "b = 2**$w; m=b**$r+b**$s-1; print(m)")"
else
elif [ "$type" = swc ]
then
local modulus="$(python -c "b = 2**$w; m=b**$r-b**$s+1; print(m)")"
elif [ "$type" = swb ]
then
local modulus="$(python -c "b = 2**$w; m=b**$r-b**$s-1; print(m)")"
else
>&2 echo BUG is_modulus_prime
exit 1
fi
openssl prime "$modulus" | fgrep --quiet 'is prime'
......@@ -93,7 +100,7 @@ configure_cxx_code()
}
make_plot()
make_luscher_plot()
{
local type="${1:?}"
local w="${2:?}"
......@@ -117,44 +124,62 @@ make_plot()
find_prng_parameters()
{
local type="${1:?}"
local candidates_file="$type-candidates.txt"
local results_file="$type-primality-checks.txt"
python3 "$cwd/generate-awc-swb-moduli.py" "$type" >"$candidates_file"
cat "$candidates_file" \
| awk '{print $5}' \
| xargs -n 1 openssl prime >"$results_file"
paste "$candidates_file" "$results_file" \
| fgrep 'is prime' \
| awk '{printf("%6s %2d %2d %2d\n", $1, $2, $3, $4)}'
}
find_awc_swb_prngs()
{
type="${1:?}"
echo "Searching for $type parameter combinations with prime modulus..."
local generators="$(find_prng_parameters "$type")"
local num_generators="$(wc -l <<<"$generators")"
if [ -z "$generators" ]
then
echo "Found no $type parameter combinations"
return
fi
echo "Found $num_generators $type parameter combinations"
echo 'type b r s t'
while read gen
do
local b="$(awk '{print $2}' <<<"$gen")"
local r="$(awk '{print $3}' <<<"$gen")"
local s="$(awk '{print $4}' <<<"$gen")"
make_luscher_plot "$type" $b $r $s
done <<<"$generators"
}
cwd="$(pwd)"
mkdir -p -- "luscher-trajectory-plots"
cd -- "luscher-trajectory-plots"
# set the python matplotlib back-end
export MPLBACKEND="cairo"
make_plot swb 64 9 8
make_plot swb 64 13 7
make_plot swb 64 26 4
make_plot swb 64 30 6
make_plot swb 32 7 3
make_plot swb 32 10 9
make_plot swb 32 16 9
make_plot swb 32 17 3
make_plot swb 32 18 16
make_plot swb 32 21 6
make_plot swb 32 24 19
make_plot swb 32 26 14
make_plot swb 32 31 26
make_plot swb 32 34 19
make_plot swb 32 37 24
make_plot swb 24 24 10
make_plot swb 24 25 11
make_plot swb 24 28 8
make_plot swb 16 4 2
make_plot swb 16 5 3
make_plot swb 16 6 2
make_plot swb 16 11 3
make_plot swb 16 11 5
make_plot swb 16 14 6
make_plot awc 32 8 3
make_plot awc 32 16 3
make_plot awc 32 29 17
make_plot awc 32 33 13
make_plot awc 32 33 27
make_plot awc 32 39 11
find_awc_swb_prngs awc
find_awc_swb_prngs cawc
find_awc_swb_prngs swb
find_awc_swb_prngs swc
#!/usr/bin/python3
import sys
if len(sys.argv) < 2:
msg = 'No PRNG type given, try `awc`, `cawc`, `swb`, or `swc`'
print(msg, file=sys.stderr)
sys.exit(1)
kind = sys.argv[1]
if kind not in ['awc', 'cawc', 'swb', 'swc']:
print('unknown PRNG type {:s}'.format(kind), file=sys.stderr)
sys.exit(2)
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)
for e in [16, 32, 64]:
b = 2**e
r_max = \
30 if e == 16 else \
40 if e == 32 else \
50 if e == 64 else -1
for r in range(3, r_max):
s_max = r-1 if r < 10 else r-2
for s in range(2, s_max):
assert r > s
if r % s == 0:
continue
m = compute_modulus(b, r, s)
print('{:6s} {:2d} {:2d} {:2d} {:d}'.format(kind, e, r, s, m))
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