src/random.c

Go to the documentation of this file.
00001 /*
00002    libit - Library for basic source and channel coding functions
00003    Copyright (C) 2005-2005 Vivien Chappelier, Herve Jegou
00004 
00005    This library is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU Library General Public
00007    License as published by the Free Software Foundation; either
00008    version 2 of the License, or (at your option) any later version.
00009 
00010    This library is distributed in the hope that it will be useful,
00011    but WITHOUT ANY WARRANTY; without even the implied warranty of
00012    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013    Library General Public License for more details.
00014 
00015    You should have received a copy of the GNU Library General Public
00016    License along with this library; if not, write to the Free
00017    Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00018 */
00019 
00020 /*
00021   Random numbers generators
00022 
00023   This Libit code is: 
00024   Copyright (C) 2005-2006 Vivien Chappelier, François Cayre  
00025   
00026   Based on the mt19937 implementation (LGPL)
00027     Copyright (C) 1997, 1999 Makoto Matsumoto and Takuji Nishimura.
00028     http://www.math.keio.ac.jp/matumoto/emt.html
00029   REFERENCE
00030     M. Matsumoto and T. Nishimura,
00031     "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
00032      Pseudo-Random Number Generator",                             
00033     ACM Transactions on Modeling and Computer Simulation,
00034     Vol. 8, No. 1, January 1998, pp 3--30.
00035 
00036   Uses ZIGGURAT method for GAUSSIAN pseudo-random numbers generation 
00037      Copyright (C) 2006 François Cayre
00038 
00039   Our Gaussian PRNG uses 256 ziggurat levels and tries to stick 
00040   as much as possible to the original paper. In particular, 
00041   compared to the GSL implementation of J. Voss, we do not model 
00042   the tail of the Gaussian with an exponential distribution. 
00043   Informal tests on an x86_64 machine do not reveal significant 
00044   speed improvement of one implementation over the other. 
00045 
00046   REFERENCE 
00047      George Marsagli, Wai Wan Tsang, 
00048      "The Ziggurat Method for Generating Random Variables",
00049      Journal of Statistical Software, vol. 5 (2000), no. 8.
00050 
00051 
00052     TODO: 
00053     - Add exponential PRNG based on Ziggurat
00054 
00055     Changelog (06/16/06) 
00056     - Replaced Box-Muller based it_randn() function by Ziggurat method 
00057 
00058     Changelog (06/15/06) 
00059     - Now based on: 
00060     http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/MTARCOK/mt19937ar-cok.c 
00061     - Provides more functions for various real output 
00062 
00063 */
00064 
00065 
00066 #include <math.h>
00067 #ifdef WIN32
00068 #include <windows.h>    /* sigh... */
00069 #else
00070 #include <sys/time.h>
00071 #endif
00072 #include <stdlib.h>
00073 #include <it/vec.h>
00074 #include <it/random.h>
00075 
00076 typedef unsigned int it_uint32_t; /* Solaris and Windows don't have stdint... */
00077 
00078 /* MT19937 related stuff */
00079 
00080 /* MT199737 constants */
00081 #define MSB 0x80000000UL  /* most significant bit */
00082 #define LSB 0x7fffffffUL  /* least significant bits */
00083 #define N          624
00084 #define M          397
00085 #define A   0x9908b0dfUL
00086 #define MIXBITS(u,v) ( ((u) & MSB) | ((v) & LSB) )
00087 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? A : 0UL))
00088 
00089 
00090 /* some random default state (generated from a random seed) */
00091 #define MT199737_STATE_INITIALIZER { 0x7036c3e9UL, 0x67b6b695UL, 0x0b2ac651UL, 0xe3a3ebddUL, 0x891b18f9UL, 0x43f0a465UL, 0x4471c5e1UL, 0x5c08e22dUL, 0x08c66709UL, 0x8124f735UL, 0x425ca671UL, 0x7f16057dUL, 0xe866be19UL, 0x10067f05UL, 0xa8abf801UL, 0xc54ea5cdUL, 0xa4332e29UL, 0x447d0bd5UL, 0x3d914a91UL, 0xf8b3131dUL, 0x527bc739UL, 0x1d756da5UL, 0xfb3f2e21UL, 0x95109d6dUL, 0xd3b99949UL, 0x91b17475UL, 0x227932b1UL, 0xf55194bdUL, 0x2b9eb459UL, 0xb197f045UL, 0xde23e841UL, 0x9dcd490dUL, 0x13262869UL, 0x6e04b115UL, 0x07d4ded1UL, 0xf3980a5dUL, 0xd3a40579UL, 0xd41886e5UL, 0x9d62a661UL, 0xb0d328adUL, 0x7ad55b89UL, 0x8e0941b5UL, 0x7774cef1UL, 0x65fcf3fdUL, 0x77f03a99UL, 0x78f1b185UL, 0xd113e881UL, 0x5840bc4dUL, 0xb1b3b2a9UL, 0x1fa1a655UL, 0x30398311UL, 0x0cc6d19dUL, 0x2577d3b9UL, 0xea6df025UL, 0x3f602ea1UL, 0xd10483edUL, 0x1f3dadc9UL, 0xd4005ef5UL, 0x28137b31UL, 0x900c233dUL, 0x1abf50d9UL, 0x7327c2c5UL, 0xfe7ff8c1UL, 0x44dcff8dUL, 0x5d7fcce9UL, 0x29a7eb95UL, 0xce033751UL, 0x59b368ddUL, 0x59db31f9UL, 0x4809a965UL, 0xd6bbc6e1UL, 0x4458af2dUL, 0xeb169009UL, 0xf66acc35UL, 0x8c193771UL, 0xaf73227dUL, 0x7a6ff719UL, 0xb24e2405UL, 0xe46c1901UL, 0x80d612cdUL, 0x1d2e7729UL, 0x316b80d5UL, 0x8975fb91UL, 0x0cd1d01dUL, 0xbbb22039UL, 0x797fb2a5UL, 0x79f96f21UL, 0xc683aa6dUL, 0x51840249UL, 0xfd1c8975UL, 0xac4a03b1UL, 0xbd25f1bdUL, 0x56662d59UL, 0x8d78d545UL, 0x41dc4941UL, 0x365ff60dUL, 0x6063b169UL, 0xf1406615UL, 0xdbd5cfd1UL, 0xb596075dUL, 0x0ee09e79UL, 0xf0640be5UL, 0xa09d2761UL, 0xc03975adUL, 0x4eaa0489UL, 0xa4e996b5UL, 0x8269dff1UL, 0xaf1890fdUL, 0x0705f399UL, 0xe0bbd685UL, 0x56d48981UL, 0xdcaea94dUL, 0x3fc37ba9UL, 0x787a9b55UL, 0x4f66b411UL, 0x80740e9dUL, 0xd04aacb9UL, 0x434ab525UL, 0x632aefa1UL, 0x872e10edUL, 0xa7ac96c9UL, 0x9fa5f3f5UL, 0x393ccc31UL, 0xb83f003dUL, 0xbdb349d9UL, 0x4d2b27c5UL, 0x2458d9c1UL, 0x77f62c8dUL, 0xbcf1d5e9UL, 0x6b6e2095UL, 0xbf6ca851UL, 0x76dfe5ddUL, 0x75d44af9UL, 0x6dc7ae65UL, 0xbb26c7e1UL, 0x9e157c2dUL, 0x2aafb909UL, 0xd425a135UL, 0x6c86c871UL, 0x888d3f7dUL, 0xc4d23019UL, 0x78dac905UL, 0xac6d3a01UL, 0xd96a7fcdUL, 0x0292c029UL, 0x436ef5d5UL, 0x982bac91UL, 0xbf4d8d1dUL, 0xae617939UL, 0x106ef7a5UL, 0xc314b021UL, 0xf4a3b76dUL, 0xeed76b49UL, 0x9e3c9e75UL, 0x690bd4b1UL, 0x8cf74ebdUL, 0xbfc6a659UL, 0x4edeba45UL, 0x3215aa41UL, 0xdf3fa30dUL, 0xa44a3a69UL, 0x8ed11b15UL, 0x16e7c0d1UL, 0xdd31045dUL, 0xa1d63779UL, 0xb0d490e5UL, 0xf678a861UL, 0x278cc2adUL, 0x9447ad89UL, 0x0ebeebb5UL, 0x6c8ff0f1UL, 0x2f712dfdUL, 0xeaf4ac99UL, 0x3f4afb85UL, 0x79562a81UL, 0xb4a9964dUL, 0xdebc44a9UL, 0x30e89055UL, 0x89e4e511UL, 0xf0fe4b9dUL, 0x311685b9UL, 0xf98c7a25UL, 0x71d6b0a1UL, 0xc0849dedUL, 0x84247fc9UL, 0x2b8088f5UL, 0xe5d71d31UL, 0x16eedd3dUL, 0x5bc042d9UL, 0x7f338cc5UL, 0x0732bac1UL, 0x11dc598dUL, 0xd78cdee9UL, 0xa2095595UL, 0x90671951UL, 0xf82962ddUL, 0x360663f9UL, 0xfa2ab365UL, 0x32b2c8e1UL, 0x763f492dUL, 0x3091e209UL, 0x2f557635UL, 0xb4a55971UL, 0x67645c7dUL, 0x408d6919UL, 0x48ac6e05UL, 0x61af5b01UL, 0x7c0beccdUL, 0xdd600929UL, 0x2f876ad5UL, 0x5ab25d91UL, 0x0d264a1dUL, 0xc389d239UL, 0x67433ca5UL, 0x5790f121UL, 0x6c70c46dUL, 0x54b3d449UL, 0xca11b375UL, 0x69bea5b1UL, 0x01c5abbdUL, 0x20c01f59UL, 0x1ac99f45UL, 0x4fd00b41UL, 0x856c500dUL, 0xa7d9c369UL, 0x3bb6d015UL, 0xea0ab1d1UL, 0xa769015dUL, 0x6584d079UL, 0xda6a15e5UL, 0x5ff52961UL, 0x73cd0fadUL, 0x34ae5689UL, 0x608940b5UL, 0x86e701f1UL, 0xc406cafdUL, 0x1cbc6599UL, 0xf99f2085UL, 0x1998cb81UL, 0x0d31834dUL, 0x979e0da9UL, 0x7deb8555UL, 0x50b41611UL, 0xdb65889dUL, 0x60db5eb9UL, 0x12333f25UL, 0x6c6371a1UL, 0x4a082aedUL, 0xdda568c9UL, 0x4c901df5UL, 0xbee26e31UL, 0xc91bba3dUL, 0x2de63bd9UL, 0xae40f1c5UL, 0xc80d9bc1UL, 0x7f8f868dUL, 0xf650e7e9UL, 0x42798a95UL, 0xf1f28a51UL, 0x9a8fdfddUL, 0xf3717cf9UL, 0x3232b865UL, 0x7e5fc9e1UL, 0xd9d6162dUL, 0x65bd0b09UL, 0x1cfa4b35UL, 0x3574ea71UL, 0xa8f8797dUL, 0x66a1a219UL, 0x06c31305UL, 0x65327c01UL, 0x15ba59cdUL, 0x36965229UL, 0xaab4dfd5UL, 0xc20a0e91UL, 0xf35c071dUL, 0x942b2b39UL, 0x02fc81a5UL, 0xb86e3221UL, 0x7aead16dUL, 0x2c193d49UL, 0xd59bc875UL, 0xbf6276b1UL, 0xb89108bdUL, 0x32529859UL, 0x16398445UL, 0x3c0b6c41UL, 0x15e5fd0dUL, 0x34124c69UL, 0xecf18515UL, 0x863ea2d1UL, 0x513dfe5dUL, 0x32ec6979UL, 0x32249ae5UL, 0x9e12aa61UL, 0x31fa5cadUL, 0x18ddff89UL, 0x2f4895b5UL, 0x226f12f1UL, 0x49d967fdUL, 0x955d1e99UL, 0x74b84585UL, 0x189c6c81UL, 0x1346704dUL, 0x7368d6a9UL, 0x94837a55UL, 0x14d44711UL, 0xbca9c59dUL, 0x789937b9UL, 0x923f0425UL, 0x53d132a1UL, 0xf0b8b7edUL, 0xdd2f51c9UL, 0xd7d4b2f5UL, 0x555ebf31UL, 0xebc5973dUL, 0x6d2534d9UL, 0x7f5356c5UL, 0x87e97cc1UL, 0x2e0fb38dUL, 0x623df0e9UL, 0xc1bebf95UL, 0x950efb51UL, 0x1b135cddUL, 0x071595f9UL, 0x5adfbd65UL, 0xdf2dcae1UL, 0xd5d9e32dUL, 0x33313409UL, 0xb2142035UL, 0xbff57b71UL, 0xaa49967dUL, 0xb00edb19UL, 0x981eb805UL, 0x17f69d01UL, 0x5375c6cdUL, 0x97359b29UL, 0x69f754d5UL, 0xbf32bf91UL, 0x6eeec41dUL, 0xb9458439UL, 0x689ac6a5UL, 0x66ac7321UL, 0x6d11de6dUL, 0x1e07a649UL, 0x15dadd75UL, 0x7af747b1UL, 0x4e5965bdUL, 0xad7e1159UL, 0x662e6945UL, 0x97c7cd41UL, 0x7dacaa0dUL, 0x11f3d569UL, 0x97813a15UL, 0x1c8393d1UL, 0x17affb5dUL, 0xe30d0279UL, 0x7d041fe5UL, 0x71d12b61UL, 0xef14a9adUL, 0x29d6a889UL, 0x0ffceab5UL, 0x902823f1UL, 0x9de904fdUL, 0x4dd6d799UL, 0x15966a85UL, 0x57610d81UL, 0xf3e85d4dUL, 0x7b1c9fa9UL, 0xa9b06f55UL, 0x47457811UL, 0x11cb029dUL, 0x915010b9UL, 0x7eafc925UL, 0x291ff3a1UL, 0x819644edUL, 0xabc23ac9UL, 0xa24e47f5UL, 0x3a4c1031UL, 0x9bec743dUL, 0x527d2dd9UL, 0x976abbc5UL, 0x67c65dc1UL, 0x8a5ce08dUL, 0x6453f9e9UL, 0x94d8f495UL, 0x2abc6c51UL, 0x36b3d9ddUL, 0xc9f2aef9UL, 0xb931c265UL, 0x961ccbe1UL, 0x774ab02dUL, 0x01ee5d09UL, 0x03a2f535UL, 0x25270c71UL, 0xc857b37dUL, 0x95d51419UL, 0xe1bf5d05UL, 0xdafbbe01UL, 0xe23e33cdUL, 0x883de429UL, 0x224ec9d5UL, 0x432c7091UL, 0x7cde811dUL, 0xcbd8dd39UL, 0x1d1e0ba5UL, 0xe34bb421UL, 0x8fe5eb6dUL, 0xd37f0f49UL, 0xdfcef275UL, 0xad7d18b1UL, 0x601ec2bdUL, 0x4b428a59UL, 0x2fa84e45UL, 0x04052e41UL, 0xa9c0570dUL, 0x0a7e5e69UL, 0x3065ef15UL, 0xddd984d1UL, 0x37bef85dUL, 0x4ee69b79UL, 0x8008a4e5UL, 0x9c30ac61UL, 0x381bf6adUL, 0x50985189UL, 0x97a63fb5UL, 0x211234f1UL, 0x9d35a1fdUL, 0x3f299099UL, 0x41398f85UL, 0xb6e6ae81UL, 0xdc174a4dUL, 0xb7b968a9UL, 0xf2726455UL, 0x5907a911UL, 0x57c93f9dUL, 0xc3ffe9b9UL, 0xdc858e25UL, 0xed4fb4a1UL, 0xc9a0d1edUL, 0x725e23c9UL, 0x80fcdcf5UL, 0xfeaa6131UL, 0xf690513dUL, 0x16ee26d9UL, 0x9b8720c5UL, 0x88a43ec1UL, 0x01770d8dUL, 0x459302e9UL, 0x30c82995UL, 0x63fadd51UL, 0xaa7156ddUL, 0x9508c7f9UL, 0x9228c765UL, 0xe42ccce1UL, 0xcb287d2dUL, 0x3af48609UL, 0x26a6ca35UL, 0x36099d71UL, 0x6022d07dUL, 0x90f44d19UL, 0xc8a50205UL, 0x0f41df01UL, 0x6f13a0cdUL, 0x92af2d29UL, 0x88bb3ed5UL, 0x3ef72191UL, 0x1a2b3e1dUL, 0x64e53639UL, 0xa58650a5UL, 0xaf4bf521UL, 0x3066f86dUL, 0xf57f7849UL, 0x88780775UL, 0x67f3e9b1UL, 0x8ae11fbdUL, 0xc4a00359UL, 0x97a73345UL, 0x21c38f41UL, 0x8721040dUL, 0xe6b1e769UL, 0xac9fa415UL, 0xfb4075d1UL, 0xee6af55dUL, 0x4f793479UL, 0x003229e5UL, 0xde312d61UL, 0x9a1043adUL, 0x7622fa89UL, 0x5b4494b5UL, 0x262d45f1UL, 0x24bf3efdUL, 0x62554999UL, 0x5ca1b485UL, 0x182d4f81UL, 0xf8d3374dUL, 0x323f31a9UL, 0xa3c95955UL, 0xbb1ada11UL, 0x0ba47c9dUL, 0x29a8c2b9UL, 0xb0c05325UL, 0xa16075a1UL, 0x95d85eedUL, 0x5a030cc9UL, 0x48e071f5UL, 0x3379b231UL, 0x18b12e3dUL, 0xf3781fd9UL, 0x30a885c5UL, 0x0b831fc1UL, 0x005e3a8dUL, 0x4efb0be9UL, 0x0a8c5e95UL, 0xf1ca4e51UL, 0x334bd3ddUL, 0xc157e0f9UL, 0x2ac4cc65UL, 0x0a5dcde1UL, 0xde734a2dUL, 0x4743af09UL, 0x301f9f35UL, 0xc39d2e71UL, 0xceaaed7dUL, 0x1a6c8619UL, 0x31cfa705UL, 0x15c90001UL, 0xa6f60dcdUL, 0x3f897629UL, 0x523cb3d5UL, 0xa392d291UL, 0x43d4fb1dUL, 0x1d6a8f39UL, 0x86d395a5UL, 0x4bad3621UL, 0x9b95056dUL, 0x2d08e149UL, 0x64d61c75UL, 0xbb5bbab1UL, 0x6ba07cbdUL, 0xd2967c59UL, 0xc32b1845UL, 0x9202f041UL, 0x02ceb10dUL, 0x6f8e7069UL, 0x012e5915UL, 0xa5b866d1UL, 0x78b3f25dUL, 0xbdc4cd79UL, 0xc280aee5UL, 0xf8d2ae61UL, 0xa1f190adUL, 0x8376a389UL, 0xefd7e9b5UL, 0xf07956f1UL, 0x1185dbfdUL, 0xb05a0299UL, 0xccced985UL, 0x5c34f081UL, 0x771c244dUL, 0xf3adfaa9UL, 0xf2b54e55UL, 0xde7f0b11UL, 0xaa5cb99dUL, 0xdb4a9bb9UL, 0x00601825UL, 0x465236a1UL, 0xb33cebedUL, 0x8bb0f5c9UL, 0xcef906f5UL, 0x69ba0331UL, 0x1f4f0b3dUL, 0x211b18d9UL, 0xfbceeac5UL, 0x116300c1UL, 0xf412678dUL, 0xc98c14e9UL, 0x97259395UL, 0x852abf51UL, 0x8e4350ddUL, 0xa7dff9f9UL, 0xc805d165UL, 0x49afcee1UL, 0xbe2b172dUL, 0x8fdbd809UL, 0x350d7435UL, 0x9ee1bf71UL, 0x70f00a7dUL, 0xab3dbf19UL, 0x023f4c05UL, 0x4f912101UL, 0x36e57acdUL, 0x17ccbf29UL, 0x33d328d5UL, 0x61ff8391UL, 0xf6dbb81dUL, 0x8e68e839UL, 0x4605daa5UL, 0x396f7721UL, 0x1e70126dUL, 0x231b4a49UL, 0xc9e93175UL, 0xb8b48bb1UL, 0x9f5cd9bdUL, 0x2e25f559UL, 0xd733fd45UL, 0xf5c35141UL, 0x09c95e0dUL, 0x6e13f969UL, 0x23120e15UL, 0x0e4157d1UL, 0x1399ef5dUL, 0x72c96679UL, 0x8bf433e5UL, 0xad152f61UL, 0xdcbfddadUL, 0x61934c89UL, 0xea603eb5UL, 0xd0f667f1UL, 0x408978fdUL, 0x2237bb99UL, 0xf6c0fe85UL, 0x63fd9181UL, 0x83f2114dUL }
00092 
00093 static it_uint32_t state[N] = MT199737_STATE_INITIALIZER; /* state vector  */
00094 static int left = 1;
00095 static int initf = 0;
00096 static it_uint32_t *next;
00097 
00098 static int idx = 0;   /* current index */
00099 
00100 /* Ziggurat related stuff */
00101 
00102 
00103 #define ZIGLEVELS 256
00104 
00105 #define ZIGR 3.6554204190269
00106 
00107 #define ZIGRINV 0.2735663440503
00108 
00109 #define ZIGV 0.0049287609634869
00110 
00111 static const unsigned long int zkn[ZIGLEVELS] = {
00112   3995918566u, 3230394625u, 3660926751u, 3844040550u,
00113   3945094684u, 4009015195u, 4053037110u, 4085175722u,
00114   4109658338u, 4128923088u, 4144473829u, 4157288020u,
00115   4168028249u, 4177159436u, 4185017354u, 4191850516u,
00116   4197846767u, 4203150717u, 4207875493u, 4212110858u,
00117   4215928942u, 4219388367u, 4222537267u, 4225415529u,
00118   4228056489u, 4230488221u, 4232734536u, 4234815761u,
00119   4236749356u, 4238550398u, 4240231977u, 4241805513u,
00120   4243281010u, 4244667275u, 4245972086u, 4247202343u,
00121   4248364183u, 4249463089u, 4250503969u, 4251491232u,
00122   4252428852u, 4253320414u, 4254169169u, 4254978064u,
00123   4255749782u, 4256486769u, 4257191258u, 4257865296u,
00124   4258510757u, 4259129366u, 4259722710u, 4260292251u,
00125   4260839340u, 4261365226u, 4261871068u, 4262357938u,
00126   4262826834u, 4263278682u, 4263714347u, 4264134635u,
00127   4264540296u, 4264932033u, 4265310504u, 4265676324u,
00128   4266030071u, 4266372285u, 4266703477u, 4267024124u,
00129   4267334677u, 4267635559u, 4267927170u, 4268209888u,
00130   4268484069u, 4268750048u, 4269008145u, 4269258660u,
00131   4269501878u, 4269738070u, 4269967490u, 4270190383u,
00132   4270406979u, 4270617497u, 4270822145u, 4271021120u,
00133   4271214612u, 4271402799u, 4271585851u, 4271763932u,
00134   4271937197u, 4272105792u, 4272269859u, 4272429531u,
00135   4272584939u, 4272736203u, 4272883440u, 4273026764u,
00136   4273166279u, 4273302089u, 4273434291u, 4273562977u,
00137   4273688238u, 4273810159u, 4273928821u, 4274044303u,
00138   4274156678u, 4274266019u, 4274372394u, 4274475869u,
00139   4274576505u, 4274674363u, 4274769501u, 4274861972u,
00140   4274951829u, 4275039123u, 4275123900u, 4275206207u,
00141   4275286087u, 4275363581u, 4275438730u, 4275511570u,
00142   4275582138u, 4275650468u, 4275716591u, 4275780540u,
00143   4275842342u, 4275902026u, 4275959617u, 4276015139u,
00144   4276068617u, 4276120070u, 4276169520u, 4276216984u,
00145   4276262481u, 4276306026u, 4276347634u, 4276387317u,
00146   4276425089u, 4276460960u, 4276494939u, 4276527035u,
00147   4276557254u, 4276585601u, 4276612082u, 4276636699u,
00148   4276659454u, 4276680347u, 4276699378u, 4276716545u,
00149   4276731843u, 4276745269u, 4276756815u, 4276766475u,
00150   4276774239u, 4276780097u, 4276784038u, 4276786047u,
00151   4276786109u, 4276784209u, 4276780328u, 4276774446u,
00152   4276766542u, 4276756592u, 4276744570u, 4276730451u,
00153   4276714204u, 4276695798u, 4276675201u, 4276652375u,
00154   4276627284u, 4276599887u, 4276570140u, 4276537998u,
00155   4276503413u, 4276466332u, 4276426702u, 4276384463u,
00156   4276339555u, 4276291913u, 4276241467u, 4276188144u,
00157   4276131868u, 4276072556u, 4276010120u, 4275944470u,
00158   4275875508u, 4275803129u, 4275727225u, 4275647679u,
00159   4275564367u, 4275477159u, 4275385915u, 4275290486u,
00160   4275190717u, 4275086438u, 4274977471u, 4274863626u,
00161   4274744700u, 4274620477u, 4274490724u, 4274355195u,
00162   4274213623u, 4274065726u, 4273911198u, 4273749712u,
00163   4273580915u, 4273404430u, 4273219848u, 4273026727u,
00164   4272824590u, 4272612923u, 4272391163u, 4272158705u,
00165   4271914885u, 4271658984u, 4271390211u, 4271107707u,
00166   4270810524u, 4270497623u, 4270167859u, 4269819968u,
00167   4269452547u, 4269064041u, 4268652718u, 4268216642u,
00168   4267753642u, 4267261279u, 4266736796u, 4266177069u,
00169   4265578544u, 4264937154u, 4264248226u, 4263506362u,
00170   4262705290u, 4261837680u, 4260894904u, 4259866734u,
00171   4258740950u, 4257502823u, 4256134423u, 4254613696u,
00172   4252913182u, 4250998234u, 4248824459u, 4246333994u,
00173   4243449892u, 4240067428u, 4236040118u, 4231156253u,
00174   4225097395u, 4217360117u, 4207095967u, 4192747394u,
00175   4171088088u, 4134068812u, 4053578179u, 0
00176 };
00177 
00178 static const double zfn[ZIGLEVELS] = {
00179   1, 0.97710143059782, 0.95987861830851, 0.94519830763172,
00180   0.93205927637968, 0.91999056461098, 0.90872536850199, 0.89809472692918,
00181   0.88798334868776, 0.87830823196574, 0.86900715699489, 0.86003198695311,
00182   0.85134452455213, 0.84291382270764, 0.8347143689422, 0.82672481886199,
00183   0.81892708786231, 0.81130568410951, 0.80384720853967, 0.79653997325589,
00184   0.7893737056314, 0.78233931561091, 0.77542871038842, 0.76863464512957,
00185   0.76195060148557, 0.75537068779605, 0.74888955640602, 0.74250233462364,
00186   0.73620456665171, 0.7299921644229, 0.72386136571625, 0.71780869827208,
00187   0.71183094888197, 0.70592513663134, 0.70008848962889, 0.69431842467991,
00188   0.6886125294583, 0.68296854680955, 0.67738436087973, 0.67185798481587,
00189   0.6663875498242, 0.66097129540631, 0.65560756062102, 0.65029477624233,
00190   0.64503145770312, 0.63981619872978, 0.63464766558632, 0.6295245918578,
00191   0.62444577371207, 0.619410065587, 0.61441637625712, 0.60946366523946,
00192   0.60455093950318, 0.59967725045204, 0.59484169115244, 0.59004339378279,
00193   0.58528152728292, 0.58055529518449, 0.57586393360571, 0.57120670939507,
00194   0.56658291841091, 0.56199188392465, 0.55743295513706, 0.55290550579771,
00195   0.54840893291904, 0.54394265557723, 0.53950611379255, 0.53509876748314,
00196   0.53072009548603, 0.52636959464053, 0.52204677892881, 0.51775117866965,
00197   0.51348233976116, 0.50923982296903, 0.50502320325688, 0.50083206915568,
00198   0.49666602216963, 0.49252467621576, 0.48840765709505, 0.48431460199287,
00199   0.48024515900682, 0.47619898670008, 0.47217575367868, 0.46817513819109,
00200   0.46419682774868, 0.46024051876584, 0.45630591621844, 0.45239273331948,
00201   0.44850069121096, 0.444629518671, 0.4407789518352, 0.43694873393151,
00202   0.43313861502772, 0.42934835179101, 0.42557770725868, 0.42182645061953,
00203   0.41809435700533, 0.41438120729177, 0.41068678790836, 0.40701089065686,
00204   0.40335331253773, 0.39971385558425, 0.39609232670375, 0.39248853752587,
00205   0.38890230425711, 0.38533344754175, 0.38178179232853, 0.37824716774291,
00206   0.37472940696465, 0.37122834711051, 0.3677438291216, 0.3642756976555,
00207   0.36082380098257, 0.35738799088659, 0.35396812256921, 0.35056405455838,
00208   0.34717564862021, 0.3438027696745, 0.34044528571341, 0.33710306772341,
00209   0.33377598961025, 0.33046392812682, 0.3271667628038, 0.32388437588306,
00210   0.32061665225349, 0.31736347938947, 0.31412474729155, 0.31090034842955,
00211   0.30769017768774, 0.30449413231222, 0.30131211186028, 0.29814401815175,
00212   0.29498975522227, 0.29184922927829, 0.28872234865397, 0.2856090237697,
00213   0.2825091670923, 0.27942269309683, 0.27634951822993, 0.27328956087471,
00214   0.27024274131705, 0.26720898171339, 0.26418820605984, 0.2611803401627,
00215   0.25818531161026, 0.25520304974587, 0.25223348564236, 0.24927655207758,
00216   0.2463321835112, 0.24340031606275, 0.24048088749069, 0.23757383717278,
00217   0.2346791060875, 0.23179663679657, 0.22892637342863, 0.22606826166403,
00218   0.22322224872061, 0.22038828334068, 0.21756631577903, 0.21475629779197,
00219   0.21195818262754, 0.20917192501673, 0.20639748116582, 0.20363480874977,
00220   0.2008838669068, 0.19814461623394, 0.19541701878386, 0.19270103806273,
00221   0.18999663902929, 0.18730378809511, 0.18462245312601, 0.18195260344482,
00222   0.17929420983526, 0.17664724454729, 0.17401168130367, 0.171387495308,
00223   0.16877466325412, 0.16617316333711, 0.16358297526567, 0.16100408027628,
00224   0.15843646114894, 0.1558801022247, 0.15333498942496, 0.15080111027283,
00225   0.14827845391637, 0.14576701115404, 0.14326677446239, 0.14077773802608,
00226   0.13829989777048, 0.13583325139685, 0.13337779842039, 0.13093354021134,
00227   0.12850048003918, 0.12607862312034, 0.12366797666949, 0.12126854995483,
00228   0.11888035435751, 0.11650340343556, 0.11413771299272, 0.11178330115245,
00229   0.10944018843764, 0.10710839785639, 0.10478795499442, 0.1024788881147,
00230   0.10018122826488, 0.097895009393179, 0.095620268473705, 0.093357045641809,
00231   0.091105384340662, 0.088865331480041, 0.086636937608583, 0.084420257100895,
00232   0.082215348361065, 0.080022274044367, 0.077841101299147, 0.075671902031198,
00233   0.073514753193218, 0.071369737102363, 0.069236941789315, 0.06711646138285,
00234   0.065008396534479, 0.062912854888512, 0.060829951603757, 0.058759809934132,
00235   0.056702561876794, 0.054658348897885, 0.052627322747973, 0.050609646381549,
00236   0.048605494997894, 0.046615057224181, 0.044638536466257, 0.042676152458222,
00237   0.040728143049248, 0.038794766275447, 0.036876302776783, 0.034973058635058,
00238   0.033085368730228, 0.031213600740953, 0.029358159954282, 0.027519495103399,
00239   0.02569810552843, 0.023894550064357, 0.022109458219765, 0.020343544449317,
00240   0.018597626690655, 0.016872650919093, 0.015169724428837, 0.013490162179794,
00241   0.011835553465957, 0.010207861685007, 0.008609581203028, 0.0070440001618338,
00242   0.0055156798328771, 0.0040314406569058, 0.0026028041937557,
00243     0.0012544610762767
00244 };
00245 
00246 static const double zwn[ZIGLEVELS] = {
00247   9.1478851695481e-10, 5.0115208832011e-11, 6.6630615736164e-11,
00248     7.8170456525769e-11,
00249   8.7340273839945e-11, 9.5086087855792e-11, 1.0186831870851e-10,
00250     1.079489492412e-10,
00251   1.1349259815484e-10, 1.186101026495e-10, 1.2337999544575e-10,
00252     1.2786014996961e-10,
00253   1.3209456738668e-10, 1.3611756277377e-10, 1.3995646787861e-10,
00254     1.4363344318803e-10,
00255   1.4716673191116e-10, 1.5057155148405e-10, 1.5386074229221e-10,
00256     1.5704524939719e-10,
00257   1.6013448669183e-10, 1.6313661656548e-10, 1.6605876773433e-10,
00258     1.6890720707322e-10,
00259   1.716874767228e-10, 1.7440450463239e-10, 1.7706269453364e-10,
00260     1.7966599981084e-10,
00261   1.8221798463526e-10, 1.8472187493298e-10, 1.8718060116695e-10,
00262     1.8959683447525e-10,
00263   1.919730173773e-10, 1.943113900078e-10, 1.9661401264493e-10,
00264     1.9888278514931e-10,
00265   2.0111946381326e-10, 2.0332567602737e-10, 2.0550293309828e-10,
00266     2.0765264149312e-10,
00267   2.097761127391e-10, 2.118745721685e-10, 2.1394916666864e-10,
00268     2.1600097157111e-10,
00269   2.1803099679343e-10, 2.2004019232975e-10, 2.2202945317215e-10,
00270     2.2399962373315e-10,
00271   2.2595150182919e-10, 2.2788584227718e-10, 2.2980336014879e-10,
00272     2.3170473372137e-10,
00273   2.3359060715914e-10, 2.3546159295419e-10, 2.3731827415288e-10,
00274     2.3916120639031e-10,
00275   2.4099091975259e-10, 2.4280792048441e-10, 2.4461269255729e-10,
00276     2.4640569911216e-10,
00277   2.4818738378843e-10, 2.4995817195017e-10, 2.5171847181915e-10,
00278     2.5346867552316e-10,
00279   2.5520916006735e-10, 2.5694028823533e-10, 2.5866240942634e-10,
00280     2.6037586043383e-10,
00281   2.620809661706e-10, 2.6377804034486e-10, 2.6546738609138e-10,
00282     2.6714929656126e-10,
00283   2.6882405547385e-10, 2.7049193763359e-10, 2.7215320941476e-10,
00284     2.738081292165e-10,
00285   2.7545694789034e-10, 2.7709990914247e-10, 2.7873724991253e-10,
00286     2.803692007307e-10,
00287   2.8199598605465e-10, 2.8361782458787e-10, 2.852349295807e-10,
00288     2.8684750911524e-10,
00289   2.8845576637534e-10, 2.9005989990275e-10, 2.9166010384025e-10,
00290     2.9325656816283e-10,
00291   2.9484947889764e-10, 2.9643901833348e-10, 2.9802536522063e-10,
00292     2.9960869496149e-10,
00293   3.011891797929e-10, 3.027669889605e-10, 3.0434228888573e-10,
00294     3.0591524332598e-10,
00295   3.074860135284e-10, 3.0905475837765e-10, 3.1062163453817e-10,
00296     3.1218679659121e-10,
00297   3.1375039716712e-10, 3.1531258707309e-10, 3.1687351541674e-10,
00298     3.1843332972583e-10,
00299   3.1999217606444e-10, 3.2155019914576e-10, 3.2310754244186e-10,
00300     3.246643482906e-10,
00301   3.2622075800001e-10, 3.2777691195019e-10, 3.2933294969319e-10,
00302     3.3088901005074e-10,
00303   3.3244523121042e-10, 3.3400175082003e-10, 3.3555870608073e-10,
00304     3.3711623383883e-10,
00305   3.3867447067658e-10, 3.4023355300202e-10, 3.4179361713812e-10,
00306     3.4335479941134e-10,
00307   3.449172362397e-10, 3.4648106422067e-10, 3.4804642021884e-10,
00308     3.4961344145369e-10,
00309   3.511822655875e-10, 3.5275303081355e-10, 3.5432587594486e-10,
00310     3.5590094050346e-10,
00311   3.5747836481053e-10, 3.5905829007737e-10, 3.6064085849751e-10,
00312     3.6222621334e-10,
00313   3.6381449904416e-10, 3.6540586131581e-10, 3.670004472253e-10,
00314     3.6859840530739e-10,
00315   3.7019988566324e-10, 3.718050400646e-10, 3.7341402206059e-10,
00316     3.7502698708696e-10,
00317   3.7664409257835e-10, 3.7826549808356e-10, 3.7989136538412e-10,
00318     3.8152185861644e-10,
00319   3.831571443977e-10, 3.847973919558e-10, 3.8644277326365e-10,
00320     3.8809346317801e-10,
00321   3.897496395833e-10, 3.914114835406e-10, 3.9307917944223e-10,
00322     3.9475291517225e-10,
00323   3.9643288227326e-10, 3.9811927611991e-10, 3.9981229609958e-10,
00324     4.0151214580058e-10,
00325   4.0321903320851e-10, 4.0493317091112e-10, 4.0665477631241e-10,
00326     4.083840718564e-10,
00327   4.1012128526127e-10, 4.1186664976452e-10, 4.1362040437995e-10,
00328     4.1538279416707e-10,
00329   4.1715407051395e-10, 4.1893449143423e-10, 4.2072432187938e-10,
00330     4.225238340672e-10,
00331   4.2433330782756e-10, 4.261530309668e-10, 4.2798329965182e-10,
00332     4.2982441881539e-10,
00333   4.3167670258426e-10, 4.3354047473143e-10, 4.3541606915463e-10,
00334     4.3730383038276e-10,
00335   4.3920411411232e-10, 4.4111728777626e-10, 4.4304373114743e-10,
00336     4.449838369796e-10,
00337   4.4693801168865e-10, 4.4890667607729e-10, 4.5089026610667e-10,
00338     4.5288923371877e-10,
00339   4.5490404771345e-10, 4.5693519468507e-10, 4.5898318002327e-10,
00340     4.6104852898369e-10,
00341   4.6313178783442e-10, 4.6523352508497e-10, 4.6735433280495e-10,
00342     4.6949482804066e-10,
00343   4.7165565433842e-10, 4.7383748338474e-10, 4.7604101677407e-10,
00344     4.7826698791667e-10,
00345   4.805161641e-10, 4.8278934871908e-10, 4.8508738369253e-10,
00346     4.8741115208372e-10,
00347   4.8976158094802e-10, 4.9213964443049e-10, 4.9454636714088e-10,
00348     4.9698282783659e-10,
00349   4.9945016344811e-10, 5.0194957348615e-10, 5.0448232487504e-10,
00350     5.0704975726297e-10,
00351   5.096532888673e-10, 5.1229442292105e-10, 5.1497475479687e-10,
00352     5.176959798963e-10,
00353   5.2045990240549e-10, 5.2326844503496e-10, 5.2612365987966e-10,
00354     5.2902774055865e-10,
00355   5.3198303582032e-10, 5.349920648321e-10, 5.3805753441236e-10,
00356     5.4118235850975e-10,
00357   5.4436968029303e-10, 5.4762289728469e-10, 5.5094569005839e-10,
00358     5.5434205512748e-10,
00359   5.5781634278442e-10, 5.6137330081797e-10, 5.6501812524383e-10,
00360     5.6875651945056e-10,
00361   5.7259476350115e-10, 5.7653979576767e-10, 5.8059930964246e-10,
00362     5.8478186881163e-10,
00363   5.8909704555766e-10, 5.9355558786838e-10, 5.981696229009e-10,
00364     6.0295290676928e-10,
00365   6.0792113397877e-10, 6.1309232454029e-10, 6.1848731352045e-10,
00366     6.2413037753653e-10,
00367   6.3005004712642e-10, 6.3628017568877e-10, 6.4286136930223e-10,
00368     6.4984293499996e-10,
00369   6.5728559197658e-10, 6.652653367448e-10, 6.7387910993869e-10,
00370     6.8325338236121e-10,
00371   6.9355768403019e-10, 7.0502695963406e-10, 7.1800075173536e-10,
00372     7.3299724341921e-10,
00373   7.5086784216768e-10, 7.7316823738581e-10, 8.0326004350637e-10,
00374     8.5109388898754e-10
00375 };
00376 
00377 /* BEGIN MT19937 CODE */
00378 
00379 /* seed the mt19937 random number generator */
00380 void mt19937_srand (it_uint32_t seed)
00381 {
00382   state[0] = seed & 0xffffffffUL;
00383   for (idx = 1; idx < N; idx++) {
00384     state[idx] =
00385       (1812433253UL * (state[idx - 1] ^ (state[idx - 1] >> 30)) + idx);
00386     /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
00387     /* In the previous versions, MSBs of the seed affect   */
00388     /* only MSBs of the array state[].                        */
00389     /* 2002/01/09 modified by Makoto Matsumoto             */
00390     state[idx] &= 0xffffffffUL; /* for >32 bit machines */
00391   }
00392 
00393   left = 1;
00394   initf = 1;
00395 
00396   return;
00397 }
00398 
00399 void mt19937_srand_by_array (it_uint32_t init_key[], it_uint32_t key_length)
00400 {
00401   int  i, j, k;
00402 
00403   mt19937_srand (19650218UL);
00404   i = 1;
00405   j = 0;
00406   k = (N > key_length ? N : key_length);
00407   for (; k; k--) {
00408     state[i] =
00409       (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1664525UL))
00410       + init_key[j] + j;  /* non linear */
00411     state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
00412     i++;
00413     j++;
00414     if (i >= N) {
00415       state[0] = state[N - 1];
00416       i = 1;
00417     }
00418     if (j >= key_length)
00419       j = 0;
00420   }
00421   for (k = N - 1; k; k--) {
00422     state[i] =
00423       (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * 1566083941UL))
00424       - i;      /* non linear */
00425     state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
00426     i++;
00427     if (i >= N) {
00428       state[0] = state[N - 1];
00429       i = 1;
00430     }
00431   }
00432 
00433   state[0] = 0x80000000UL;  /* MSB is 1; assuring non-zero initial array */
00434   left = 1;
00435   initf = 1;
00436 
00437   return;
00438 }
00439 
00440 static void mt19937_next_state (void)
00441 {
00442   it_uint32_t *p = state;
00443   int  j;
00444 
00445   /* if init_genrand() has not been called, */
00446   /* a default initial seed is used         */
00447   if (initf == 0)
00448     mt19937_srand (5489UL);
00449 
00450   left = N;
00451   next = state;
00452 
00453   for (j = N - M + 1; --j; p++)
00454     *p = p[M] ^ TWIST (p[0], p[1]);
00455 
00456   for (j = M; --j; p++)
00457     *p = p[M - N] ^ TWIST (p[0], p[1]);
00458 
00459   *p = p[M - N] ^ TWIST (p[0], state[0]);
00460 
00461   return;
00462 }
00463 
00464 /* generates a random number on [0,0xffffffff]-interval */
00465 unsigned int mt19937_rand_int32 (void)
00466 {
00467   it_uint32_t y;
00468 
00469   if (--left == 0)
00470     mt19937_next_state ();
00471   y = *next++;
00472 
00473   /* Tempering */
00474   y ^= (y >> 11);
00475   y ^= (y << 7) & 0x9d2c5680UL;
00476   y ^= (y << 15) & 0xefc60000UL;
00477   y ^= (y >> 18);
00478 
00479   return y;
00480 }
00481 
00482 /* generates a random number on [0,0x7fffffff]-interval */
00483 int mt19937_rand_int31 (void)
00484 {
00485   it_uint32_t y;
00486 
00487   if (--left == 0)
00488     mt19937_next_state ();
00489   y = *next++;
00490 
00491   /* Tempering */
00492   y ^= (y >> 11);
00493   y ^= (y << 7) & 0x9d2c5680UL;
00494   y ^= (y << 15) & 0xefc60000UL;
00495   y ^= (y >> 18);
00496 
00497   return (long) (y >> 1);
00498 }
00499 
00500 /* generates a random number on [0,1]-real-interval */
00501 double mt19937_rand_real1 (void)
00502 {
00503   it_uint32_t y;
00504 
00505   if (--left == 0)
00506     mt19937_next_state ();
00507   y = *next++;
00508 
00509   /* Tempering */
00510   y ^= (y >> 11);
00511   y ^= (y << 7) & 0x9d2c5680UL;
00512   y ^= (y << 15) & 0xefc60000UL;
00513   y ^= (y >> 18);
00514 
00515   return (double) y *(1.0 / 4294967295.0);
00516   /* divided by 2^32-1 */
00517 }
00518 
00519 /* generates a random number on [0,1)-real-interval */
00520 double mt19937_rand_real2 (void)
00521 {
00522   it_uint32_t y;
00523 
00524   if (--left == 0)
00525     mt19937_next_state ();
00526   y = *next++;
00527 
00528   /* Tempering */
00529   y ^= (y >> 11);
00530   y ^= (y << 7) & 0x9d2c5680UL;
00531   y ^= (y << 15) & 0xefc60000UL;
00532   y ^= (y >> 18);
00533 
00534   return (double) y *(1.0 / 4294967296.0);
00535   /* divided by 2^32 */
00536 }
00537 
00538 /* generates a random number on (0,1)-real-interval */
00539 double mt19937_rand_real3 (void)
00540 {
00541   it_uint32_t y;
00542 
00543   if (--left == 0)
00544     mt19937_next_state ();
00545   y = *next++;
00546 
00547   /* Tempering */
00548   y ^= (y >> 11);
00549   y ^= (y << 7) & 0x9d2c5680UL;
00550   y ^= (y << 15) & 0xefc60000UL;
00551   y ^= (y >> 18);
00552 
00553   return ((double) y + 0.5) * (1.0 / 4294967296.0);
00554   /* divided by 2^32 */
00555 }
00556 
00557 /* generates a random number on [0,1) with 53-bit resolution*/
00558 double mt19937_rand_res53 (void)
00559 {
00560   it_uint32_t a = mt19937_rand_int32 () >> 5, b = mt19937_rand_int32 () >> 6;
00561   return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
00562 }
00563 
00564 /* These real versions are due to Isaku Wada, 2002/01/09 added */
00565 
00566 /* -------------------------------------------------------------------------- */
00567 
00568 /* generate a random number using the mt19937 algorithm */
00569 /* this number has a precision of 32 bits and is stored */
00570 /* in double-precision in the [0, 1) range.             */
00571 static double mt19937_rand ()
00572 {
00573 
00574   return (mt19937_rand_real2 ());
00575 
00576 }
00577 
00578 /* END MT19937 CODE */
00579 
00580 /* create and initialize the random number generator */
00581 void it_randomize (void)
00582 {
00583   unsigned int usec;
00584 
00585 #ifdef WIN32
00586   /* the windows version */
00587   SYSTEMTIME SystemTime;
00588 
00589   GetLocalTime (&SystemTime);
00590   usec = (unsigned int) SystemTime.wMilliseconds;
00591 #else
00592   /* the unix version */
00593   struct timeval tp;
00594 
00595   /* seed the random number generator with the current time */
00596   gettimeofday (&tp, NULL);
00597   usec = tp.tv_usec;
00598 #endif
00599 
00600   it_seed (usec);
00601 }
00602 
00603 void it_seed (int seed)
00604 {
00605   mt19937_srand ((it_uint32_t) seed);
00606 }
00607 
00608 /* generate a random value uniformly distributed in [0,1) */
00609 double it_rand (void)
00610 {
00611   return (mt19937_rand ());
00612 }
00613 
00614 /*
00615   Use this method to generate Gaussian distributed values 
00616   with N(0,1) model. This code closely follows original 
00617   implementation from Marsaglia and Tsang. 
00618  */
00619 double it_randn (void)
00620 {
00621 
00622   unsigned long int i, j;
00623   double x, y;
00624   int  s;
00625 
00626   while (1) {
00627 
00628     j = mt19937_rand_int32 ();
00629 
00630     i = j & 0x000000FF;   /* "Form i from the last 8 bits of j" */
00631 
00632     x = j * zwn[i];
00633 
00634     s = j & 0x00000800 ? 1. : -1.;  /* Get a sign from any bit of j */
00635 
00636     if (j < zkn[i])   /* Should exit here 99% of the time */
00637       return (s * x);
00638 
00639 
00640     if (!i) {     /* "return an x from the tail" */
00641       do {
00642   x = -log (mt19937_rand_real3 ()) * ZIGRINV;
00643   y = -log (mt19937_rand_real3 ());
00644       } while (y + y < x * x);
00645 
00646       return (s * (ZIGR + x));
00647     }
00648 
00649     if (mt19937_rand_real3 () * (zfn[i - 1] - zfn[i]) <
00650   exp (-.5 * x * x) - zfn[i])
00651       return (s * x);
00652   }
00653 
00654 }
00655 
00656 /* generate a random variable from its probability
00657    density function using the acceptance-rejection method.
00658    the pdf is assumed to be zero outside [a, b]
00659    and centered on its maximum value.
00660 */
00661 double it_randpdf (double a, double b, it_function_t pdf, it_args_t args)
00662 {
00663   double x, y;
00664   double max;
00665 
00666   max = pdf (0, args);    /* evaluate the maximum value of the pdf. */
00667   do {
00668     /* draw x and y uniformly in [a,b]x[0,1] */
00669     x = a + it_rand () * (b - a);
00670     y = max * it_rand ();
00671 
00672     /* check if (x,y) is below the pdf curve */
00673     if (y < pdf (x, args))
00674       return (x);
00675   } while (1);
00676 }
00677 
00678 
00679 
00680 int it_rand_memoryless (vec pdf)
00681 {
00682   double x;
00683   int  s, e, h, r;    /* r is the result */
00684   vec  pdfcs = vec_cum_sum (pdf);
00685   vec_ins (pdfcs, 0, 0.0);
00686 
00687   /* generate a uniformly distributed value  */
00688   /* and search the integral of the pdf for  */
00689   /* the abciss corresponding to that random */
00690   /* value.                                  */
00691   x = it_rand ();
00692   s = 0;
00693   e = vec_length (pdfcs) - 1;
00694 
00695   while (e > s + 1) {
00696     /* compare with the median of the range */
00697     /* and narrow down the range. */
00698     h = (s + e) / 2;
00699 
00700     if (x < pdfcs[h])
00701       e = h;
00702     else
00703       s = h;
00704   }
00705 
00706   if (x < pdfcs[e])
00707     r = s;
00708   else
00709     r = e;
00710 
00711   vec_delete (pdfcs);
00712   return r;
00713 }

Hosted by
Copyright (C) 2005-2006 Hervé Jégou
Vivien Chappelier
Francois Cayre
libit logo courtesy of Jonathan Delhumeau