Irwin Hall Gaussian Approximation

irwin_hall_constants.py

from mpmath import mp

mp.prec = 1024

def represent(x, base=2**96):
    return int(mp.nint(x * base))

samples = 20
denom = mp.sqrt(mp.mpf(samples) / mp.mpf(12.0))
scaling = (mp.mpf(10**18) / mp.mpf(1 + 0x1fffffffffffffff)) / denom
shift = (mp.mpf(samples // 2) * mp.mpf(10**18) ) / denom

print("scaling (base 2**96)", represent(scaling))
print("shift", int(shift))
print("area outside cdf", (mp.mpf(1) - mp.ncdf(samples // 2)) * mp.mpf(2))

irwin_hall_ks_test.py

import random
import numpy as np
from scipy.stats import kstest

UINT256_MAX = 0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
UINT127_MAX = 0x7fffffffffffffffffffffffffffffff
WAD = 1000000000000000000

def _add(a, b):
    return (a + b) & UINT256_MAX

def _and(a, b):
    return (a & b) & UINT256_MAX    

def _mul(a, b):
    return (a * b) & UINT256_MAX

def _not(a):
    return (UINT256_MAX ^ a) & UINT256_MAX

def _mulmod(a, b, d):
    return ((a * b) % d) & UINT256_MAX

def _mod(a, d):
    return (a % d) & UINT256_MAX

def _shr(s, x):
    return (x >> s) & UINT256_MAX

def _sar(s, x):
    if x < 0:
        return -_shr(s, -x)
    return _shr(s, x)

def _shl(s, x):
    return (x << s) & UINT256_MAX

def _sub(a, b):
    return a - b

def _div(a, b):
    return (a // b) & UINT256_MAX

def generate():
    n = 0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff43
    a = 0x100000000000000000000000000000051
    m = 0x1fffffffffffffff1fffffffffffffff1fffffffffffffff1fffffffffffffff
    s = 0x1000000000000000100000000000000010000000000000001
    r = random.randint(0, UINT256_MAX)
    r1 = _mulmod(r, a, n)
    r2 = _mulmod(r1, a, n)
    r3 = _mulmod(r2, a, n)
    result = _sub(_sar(96, _mul(26614938895861601847173011183,
        _add(_add(_shr(192, _mul(s, _add(_and(m, r), _and(m, r1)))),
        _shr(192, _mul(s, _add(_and(m, r2), _and(m, r3))))),
        _shr(192, _mul(_and(m, _mulmod(r3, a, n)), s))))),
        7745966692414833770)
    return result


def samples() :
    n = 200000
    # return np.random.standard_normal(n)
    a = []
    for i in range(n):
        a.append(generate() / WAD)
    return a

# Kolmogorov-Smirnov test.
(D, p) = kstest(samples(), 'norm')
print("(D, p)", (D, p))

primitive_roots.py

from sympy.ntheory import factorint, is_primitive_root, isprime
from sympy import primerange
from mpmath import mp
import random 

mp.prec = 1024

def find_primitive_root(start, m):
    factors = factorint(m - 1)  # Factorize m - 1
    for a in primerange(start, m):  # Iterate over potential candidates for 'a'
        if all(pow(a, (m - 1) // p, m) != 1 for p in factors):
            return a  # Found a primitive root
    return None


m = 2 ** 256 - 1
while not isprime(m):
    m -= 1

start = int(mp.sqrt(m))
print("m", hex(m))
a = find_primitive_root(start, m)
print("a", hex(a))

assert isprime(a)