Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Copyright 2013 Kwant authors.
2#
3# This file is part of Kwant. It is subject to the license terms in the file
4# LICENSE.rst found in the top-level directory of this distribution and at
5# https://kwant-project.org/license. A list of Kwant authors can be found in
6# the file AUTHORS.rst at the top-level directory of this distribution and at
7# https://kwant-project.org/authors.
9"""Random-access random numbers
11This module provides routines that given some input compute a "random" output
12that depends on the input in a (cryptographically) intractable way.
14This turns out to be very useful when one needs to map some irregular objects to
15random numbers in a deterministic and reproducible way.
17Internally, the md5 hash algorithm is used. The randomness thus generated is
18good enough to pass the "dieharder" battery of tests: see the function `test` of
19this module.
20"""
22from math import pi, log, sqrt, cos
23from hashlib import md5
24from struct import unpack
26__all__ = ['uniform', 'gauss', 'test']
29TWOPI = 2 * pi
30BPF = 53 # Number of bits in a Python float.
31BPF_MASK = 2**53 - 1
32RECIP_BPF = 2**-BPF
35def str_to_bytes(s):
36 """Return bytes if the input is a string, else return the object as-is."""
37 try:
38 return s.encode('utf8')
39 except AttributeError:
40 return s
42def uniform2(input, salt=''):
43 """Return two independent [0,1)-distributed numbers."""
44 input = str_to_bytes(input)
45 salt = str_to_bytes(salt)
46 input = memoryview(input).tobytes() + salt
47 a, b = unpack('qq', md5(input).digest())
48 a &= BPF_MASK
49 b &= BPF_MASK
50 return a * RECIP_BPF, b * RECIP_BPF
53def uniform(input, salt=''):
54 """md5-hash `input` and `salt` and map the result to the [0,1) interval.
56 `input` must be some object that supports the buffer protocol (i.e. a string
57 or a numpy/tinyarray array). `salt` must be a string or a bytes object.
58 """
59 return uniform2(input, salt)[0]
62def gauss(input, salt=''):
63 """md5-hash `input` and `salt` and return the result as a standard normal
64 distributed variable.
66 `input` must be some object that supports the buffer protocol (i.e. a string
67 or a numpy/tinyarray array). `salt` must be a string or a bytes object.
68 """
69 # This uses the Box-Muller transform. Only one of the two results is
70 # computed.
71 a, b = uniform2(input, salt)
72 return cos(a * TWOPI) * sqrt(-2.0 * log(1.0 - b))
75def test(n=20000): # skip coverage
76 """Test the generator with the dieharder suite generating n**2 samples.
78 Executing this function may take a very long time.
79 """
80 import os
81 import tempfile
82 import subprocess
83 from tinyarray import array
84 from struct import pack
86 f = tempfile.NamedTemporaryFile(delete=False)
87 try:
88 for x in range(n):
89 for y in range(n):
90 a = array((x, y))
91 i = int(2**32 * uniform(a))
92 f.write(pack('I', i))
93 f.close()
94 subprocess.call(['dieharder', '-a', '-g', '201', '-f', f.name])
95 finally:
96 os.remove(f.name)