Skip to content

Instantly share code, notes, and snippets.

@nightuser
Created February 22, 2015 19:12
Show Gist options
  • Save nightuser/4c16f64e0fd72e8686ee to your computer and use it in GitHub Desktop.
Save nightuser/4c16f64e0fd72e8686ee to your computer and use it in GitHub Desktop.
Simple Lehmer random number generator and Pi calculus. No NumPy, pure Python only.
#!/usr/bin/env python3
# -*- encoding: utf-8 -*-
from functools import reduce
from math import sqrt
from random import random
class MyRandomizer:
def __init__(self, m, n):
assert m < n
self.m = m
self.n = n
def random(self, count):
seed = int(random() * (self.n - 1)) + 1
return reduce(lambda x, _: x + [(x[-1] * self.m) % self.n],
range(count - 1),
[seed])
def test1d(func, randomizer, counts, tries, n):
for count in counts:
for _ in range(tries):
xs = [(2 * x - n) / n for x in randomizer.random(count)]
ys = map(func, xs)
area = (2 / count) * sum(ys)
print('count = %s\tarea = %s' % (count, area))
def test2d(func, randomizer, counts, tries, n):
for count in counts:
for _ in range(tries):
points = ((2 * x - n) / n for x in randomizer.random(2 * count))
points = zip(*[iter(points)] * 2)
points_inside = sum(func(x, y) for x, y in points)
area = 4 * points_inside / count
print('count = %s\tarea = %s' % (count, area))
def main():
n = 2 ** 32
tries = 3
counts = 4
badrandomizer = MyRandomizer(3, n)
goodrandomizer = MyRandomizer(1664525, n)
counts_list = [10 ** n for n in range(1, counts + 1)]
# 1-d
def half_circle(x): return sqrt(1 - x ** 2)
print('# 1D')
print()
print('Bad randomizer (M = 3)')
test1d(half_circle, badrandomizer, counts_list, tries, n)
print()
print('Good randomizer (M = 2^16+3)')
test1d(half_circle, goodrandomizer, counts_list, tries, n)
print()
# 2-d
def circle_indicator(x, y): return x ** 2 + y ** 2 <= 1
print('# 2D')
print()
print('Bad randomizer (M = 3)')
test2d(circle_indicator, badrandomizer, counts_list, tries, n)
print()
print('Good randomizer (M = 2^16+3)')
test2d(circle_indicator, goodrandomizer, counts_list, tries, n)
if __name__ == '__main__':
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment