src/etc/ziggurat_tables.py

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