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# http://rust-lang.org/COPYRIGHT. 6# 7# Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or 8# http://www.apache.org/licenses/LICENSE-2.0> or the MIT license 9# <LICENSE-MIT or http://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// http://rust-lang.org/COPYRIGHT. 111// 112// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or 113// http://www.apache.org/licenses/LICENSE-2.0> or the MIT license 114// <LICENSE-MIT or http://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