find-awc-swb-prngs.sh 4.15 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
#!/bin/sh

# 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/.

read -r -d '' python_plot_code <<EOF
import matplotlib.pyplot as pp
import numpy as np
import sys

if len(sys.argv) < 5:
	msg = '5 command line arguments needed, {:d} given'.format(len(sys.argv))
	print(msg, file=sys.stderr)
	sys.exit(1)

type = sys.argv[1]
w = int(sys.argv[2])
r = int(sys.argv[3])
s = int(sys.argv[4])
filename = "luscher-trajectory-{:s}-{:d}-{:d}-{:d}.png".format(type, w, r, s)

if w < 16 or r <= s or s < 2:
	print('Error: w={:d} r={:d} s={:d}'.format(w, r, s), file=sys.stderr)
	sys.exit(1)

xs = np.loadtxt(sys.stdin)

if len(xs.shape) != 1:
	sys.exit(1)

if min(xs) < 0 or max(xs) >= 1:
	print('Error: min={:e} max={:e}'.format(min(xs),max(xs)), file=sys.stderr)
	sys.exit(1)

nnz, = np.nonzero(xs > 0.45)
k = -1 if nnz.size == 0 else min(nnz)
40
k_str = '>{:2d}'.format(xs.size) if k == -1 else '{:3d}'.format(k)
41

42
print('{:4s} {:d} {:2d} {:2d}  {:s}'.format(type, w, r, s, k_str))
43 44 45 46 47 48 49 50 51 52

pp.plot(xs)
pp.grid(which="both")

pp.ylim(top=1.0)
pp.ylim(bottom=2.0**-(w+1))
pp.yscale('log', basey=2)

pp.xlabel("t")
pp.ylabel("d(x(t),y(t))")
53 54
if type == 'swc':
	pp.title("{:s}(2^{:d}, {:d}, {:d})".format('SWB', w, s, r))
55
else:
56
	pp.title("{:s}(2^{:d}, {:d}, {:d})".format(type.upper(), w, r, s))
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74

pp.savefig(filename)
EOF

set -e
set -u


is_modulus_prime()
{
	local type="${1:?}"
	local w="${2:?}"
	local r="${3:?}"
	local s="${4:?}"

	if [ "$type" = awc ]
	then
		local modulus="$(python -c "b = 2**$w; m=b**$r+b**$s-1; print(m)")"
75 76
	elif [ "$type" = swc ]
	then
77
		local modulus="$(python -c "b = 2**$w; m=b**$r-b**$s+1; print(m)")"
78 79 80 81 82 83
	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
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
	fi

	openssl prime "$modulus" | fgrep --quiet 'is prime'
}

configure_cxx_code()
{
	local type="${1:?}"
	local w="${2:?}"
	local r="${3:?}"
	local s="${4:?}"

	sed -e "s/@WORD_SIZE@/$w/" \
		-e "s/@LONG_LAG@/$r/" \
		-e "s/@SHORT_LAG@/$s/" \
		"../compute-luscher-trajectory-$type.cpp.in"
}


103
make_luscher_plot()
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
{
	local type="${1:?}"
	local w="${2:?}"
	local r="${3:?}"
	local s="${4:?}"
	local cxx_filename="compute-luscher-trajectory-$type-$w-$r-$s.cpp"
	local binary="compute-luscher-trajectory-$type-$w-$r-$s"

	if ! $(is_modulus_prime $type $w $r $s)
	then
		>&2 echo "modulus of $type w=$w r=$r s=$s is not prime"
		return
	fi

	configure_cxx_code $type $w $r $s >"$cxx_filename"
	c++ -Wextra -Wall -std=c++11 -pedantic "$cxx_filename" -o "$binary"

121 122 123 124 125 126 127 128 129
	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"
130 131 132 133
}



134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
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"

169
	echo 'type  b  r  s    t log2period        log10period'
170 171 172 173 174 175 176 177 178 179 180
	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"
}


181 182 183 184
cwd="$(pwd)"
mkdir -p -- "luscher-trajectory-plots"
cd -- "luscher-trajectory-plots"

185

186 187 188
# set the python matplotlib back-end
export MPLBACKEND="cairo"

189 190 191 192
find_awc_swb_prngs awc
find_awc_swb_prngs cawc
find_awc_swb_prngs swb
find_awc_swb_prngs swc