1#!/usr/bin/env python
2#
3# Copyright 2013 The Rust Project Developers. See the COPYRIGHT
4# file at the top-level directory of this distribution and at
5# https://rust-lang.org/COPYRIGHT.
6#
7# Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
8# https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
9# <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
10# option. This file may not be copied, modified, or distributed
11# except according to those terms.
12
13# This creates the tables used for distributions implemented using the
14# ziggurat algorithm in `rand::distributions;`. They are
15# (basically) the tables as used in the ZIGNOR variant (Doornik 2005).
16# They are changed rarely, so the generated file should be checked in
17# to git.
18#
19# It creates 3 tables: X as in the paper, F which is f(x_i), and
20# F_DIFF which is f(x_i) - f(x_{i-1}). The latter two are just cached
21# values which is not done in that paper (but is done in other
22# variants). Note that the adZigR table is unnecessary because of
23# algebra.
24#
25# It is designed to be compatible with Python 2 and 3.
26
27from math import exp, sqrt, log, floor
28import random
29
30# The order should match the return value of `tables`
31TABLE_NAMES = ['X', 'F']
32
33# The actual length of the table is 1 more, to stop
34# index-out-of-bounds errors. This should match the bitwise operation
35# to find `i` in `zigurrat` in `libstd/rand/mod.rs`. Also the *_R and
36# *_V constants below depend on this value.
37TABLE_LEN = 256
38
39# equivalent to `zigNorInit` in Doornik2005, but generalised to any
40# distribution. r = dR, v = dV, f = probability density function,
41# f_inv = inverse of f
42def tables(r, v, f, f_inv):
43    # compute the x_i
44    xvec = [0]*(TABLE_LEN+1)
45
46    xvec[0] = v / f(r)
47    xvec[1] = r
48
49    for i in range(2, TABLE_LEN):
50        last = xvec[i-1]
51        xvec[i] = f_inv(v / last + f(last))
52
53    # cache the f's
54    fvec = [0]*(TABLE_LEN+1)
55    for i in range(TABLE_LEN+1):
56        fvec[i] = f(xvec[i])
57
58    return xvec, fvec
59
60# Distributions
61# N(0, 1)
62def norm_f(x):
63    return exp(-x*x/2.0)
64def norm_f_inv(y):
65    return sqrt(-2.0*log(y))
66
67NORM_R = 3.6541528853610088
68NORM_V = 0.00492867323399
69
70NORM = tables(NORM_R, NORM_V,
71              norm_f, norm_f_inv)
72
73# Exp(1)
74def exp_f(x):
75    return exp(-x)
76def exp_f_inv(y):
77    return -log(y)
78
79EXP_R = 7.69711747013104972
80EXP_V = 0.0039496598225815571993
81
82EXP = tables(EXP_R, EXP_V,
83             exp_f, exp_f_inv)
84
85
86# Output the tables/constants/types
87
88def render_static(name, type, value):
89    # no space or
90    return 'pub static %s: %s =%s;\n' % (name, type, value)
91
92# static `name`: [`type`, .. `len(values)`] =
93#     [values[0], ..., values[3],
94#      values[4], ..., values[7],
95#      ... ];
96def render_table(name, values):
97    rows = []
98    # 4 values on each row
99    for i in range(0, len(values), 4):
100        row = values[i:i+4]
101        rows.append(', '.join('%.18f' % f for f in row))
102
103    rendered = '\n    [%s]' % ',\n     '.join(rows)
104    return render_static(name, '[f64, .. %d]' % len(values), rendered)
105
106
107with open('ziggurat_tables.rs', 'w') as f:
108    f.write('''// Copyright 2013 The Rust Project Developers. See the COPYRIGHT
109// file at the top-level directory of this distribution and at
110// https://rust-lang.org/COPYRIGHT.
111//
112// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
113// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
114// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
115// option. This file may not be copied, modified, or distributed
116// except according to those terms.
117
118// Tables for distributions which are sampled using the ziggurat
119// algorithm. Autogenerated by `ziggurat_tables.py`.
120
121pub type ZigTable = &\'static [f64, .. %d];
122'''  % (TABLE_LEN + 1))
123    for name, tables, r in [('NORM', NORM, NORM_R),
124                            ('EXP', EXP, EXP_R)]:
125        f.write(render_static('ZIG_%s_R' % name, 'f64', ' %.18f' % r))
126        for (tabname, table) in zip(TABLE_NAMES, tables):
127            f.write(render_table('ZIG_%s_%s' % (name, tabname), table))
128