Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. #include <stdlib.h>
  2. #include <time.h>
  3. #include <errno.h>
  4.  
  5. #include <math.h>
  6. #ifndef M_E
  7.  #define M_E     2.7182818284590452354
  8. #endif
  9. #ifndef M_PI
  10.  #define M_PI    3.14159265358979323846
  11. #endif
  12.  
  13. /*
  14.  * following random generators are mutual exclusive.
  15.  */
  16. #define __USE_MERSENNE_TWISTER_RANDOM_GENERATOR
  17. /*#define __USE_WICHMANN_HILL_RANDOM_GENERATOR*/
  18. /*#define __USE_POSIX_RANDOM_GENERATOR*/
  19.  
  20.  
  21. #if defined(__USE_MERSENNE_TWISTER_RANDOM_GENERATOR)
  22.  
  23. /***********************************************************
  24.  * following code is borrowed from Python's _randommodule.c
  25.  ***********************************************************/
  26. /* Random objects */
  27.  
  28. /* ------------------------------------------------------------------
  29.    The code in this module was based on a download from:
  30.       http://www.math.keio.ac.jp/~matumoto/MT2002/emt19937ar.html
  31.  
  32.    It was modified in 2002 by Raymond Hettinger as follows:
  33.  
  34.     * the principal computational lines untouched except for tabbing.
  35.  
  36.     * renamed genrand_res53() to random_random() and wrapped
  37.       in python calling/return code.
  38.  
  39.     * genrand_int32() and the helper functions, init_genrand()
  40.       and init_by_array(), were declared static, wrapped in
  41.       Python calling/return code.  also, their global data
  42.       references were replaced with structure references.
  43.  
  44.     * unused functions from the original were deleted.
  45.       new, original C python code was added to implement the
  46.       Random() interface.
  47.  
  48.    The following are the verbatim comments from the original code:
  49.  
  50.    A C-program for MT19937, with initialization improved 2002/1/26.
  51.    Coded by Takuji Nishimura and Makoto Matsumoto.
  52.  
  53.    Before using, initialize the state by using init_genrand(seed)
  54.    or init_by_array(init_key, key_length).
  55.  
  56.    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
  57.    All rights reserved.
  58.  
  59.    Redistribution and use in source and binary forms, with or without
  60.    modification, are permitted provided that the following conditions
  61.    are met:
  62.  
  63.      1. Redistributions of source code must retain the above copyright
  64.     notice, this list of conditions and the following disclaimer.
  65.  
  66.      2. Redistributions in binary form must reproduce the above copyright
  67.     notice, this list of conditions and the following disclaimer in the
  68.     documentation and/or other materials provided with the distribution.
  69.  
  70.      3. The names of its contributors may not be used to endorse or promote
  71.     products derived from this software without specific prior written
  72.     permission.
  73.  
  74.    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  75.    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  76.    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  77.    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
  78.    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  79.    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  80.    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  81.    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  82.    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  83.    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  84.    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  85.  
  86.  
  87.    Any feedback is very welcome.
  88.    http://www.math.keio.ac.jp/matumoto/emt.html
  89.    email: matumoto@math.keio.ac.jp
  90. */
  91.  
  92. /* Period parameters -- These are all magic.  Don't change. */
  93. #define N 624
  94. #define M 397
  95. #define MATRIX_A 0x9908b0dfUL   /* constant vector a */
  96. #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
  97. #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
  98.  
  99. typedef struct {
  100.     unsigned long state[N];
  101.     int index;
  102.     int has_seed;
  103. } RandomObject;
  104.  
  105. /* generates a random number on [0,0xffffffff]-interval */
  106. static unsigned long
  107. genrand_int32(RandomObject *self)
  108. {
  109.     unsigned long y;
  110.     static unsigned long mag01[2]={0x0UL, MATRIX_A};
  111.     /* mag01[x] = x * MATRIX_A  for x=0,1 */
  112.     unsigned long *mt;
  113.  
  114.     mt = self->state;
  115.     if (self->index >= N) { /* generate N words at one time */
  116.         int kk;
  117.  
  118.         for (kk=0;kk<N-M;kk++) {
  119.             y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
  120.             mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
  121.         }
  122.         for (;kk<N-1;kk++) {
  123.             y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
  124.             mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
  125.         }
  126.         y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
  127.         mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
  128.  
  129.         self->index = 0;
  130.     }
  131.  
  132.     y = mt[self->index++];
  133.     y ^= (y >> 11);
  134.     y ^= (y << 7) & 0x9d2c5680UL;
  135.     y ^= (y << 15) & 0xefc60000UL;
  136.     y ^= (y >> 18);
  137.     return y;
  138. }
  139.  
  140. /* initializes mt[N] with a seed */
  141. static void
  142. init_genrand(RandomObject *self, unsigned long s)
  143. {
  144.     int mti;
  145.     unsigned long *mt;
  146.  
  147.     mt = self->state;
  148.     mt[0]= s & 0xffffffffUL;
  149.     for (mti=1; mti<N; mti++) {
  150.         mt[mti] =
  151.         (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
  152.         /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
  153.         /* In the previous versions, MSBs of the seed affect   */
  154.         /* only MSBs of the array mt[].                */
  155.         /* 2002/01/09 modified by Makoto Matsumoto         */
  156.         mt[mti] &= 0xffffffffUL;
  157.         /* for >32 bit machines */
  158.     }
  159.     self->index = mti;
  160.     return;
  161. }
  162.  
  163. /* initialize by an array with array-length */
  164. /* init_key is the array for initializing keys */
  165. /* key_length is its length */
  166. static void
  167. init_by_array(RandomObject *self, unsigned long init_key[], unsigned long key_length)
  168. {
  169.     unsigned int i, j, k;   /* was signed in the original code. RDH 12/16/2002 */
  170.     unsigned long *mt;
  171.  
  172.     mt = self->state;
  173.     init_genrand(self, 19650218UL);
  174.     i=1; j=0;
  175.     k = (N>key_length ? N : key_length);
  176.     for (; k; k--) {
  177.         mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
  178.              + init_key[j] + j; /* non linear */
  179.         mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
  180.         i++; j++;
  181.         if (i>=N) { mt[0] = mt[N-1]; i=1; }
  182.         if (j>=key_length) j=0;
  183.     }
  184.     for (k=N-1; k; k--) {
  185.         mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
  186.              - i; /* non linear */
  187.         mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
  188.         i++;
  189.         if (i>=N) { mt[0] = mt[N-1]; i=1; }
  190.     }
  191.  
  192.     mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
  193.     return;
  194. }
  195.  
  196. /*----------------------end of Mersenne Twister Algorithm----------------*/
  197.  
  198. /*************************************************************************
  199.  * following are tinypy related stuffs.
  200.  *************************************************************************/
  201.  
  202. static RandomObject _gRandom;       /* random object */
  203.  
  204. static tp_obj random_seed(TP)
  205. {
  206.     tp_obj arg = TP_DEFAULT(tp_None);
  207.  
  208.     if (arg.type == TP_NONE) {
  209.         time_t now;
  210.  
  211.         (void)time(&now);
  212.         init_genrand(&_gRandom, (unsigned long)now);
  213.         _gRandom.has_seed = 1;
  214.     } else if (arg.type == TP_NUMBER) {
  215.         init_genrand(&_gRandom, (unsigned long)arg.number.val);
  216.         _gRandom.has_seed = 1;
  217.     } else if (arg.type == TP_STRING) {
  218.         unsigned long seed;
  219.        
  220.         seed = (unsigned long)tp_hash(tp, arg);
  221.         init_genrand(&_gRandom, seed);
  222.         _gRandom.has_seed = 1;
  223.     } else {
  224.         tp_raise(tp_None,tp_printf(tp, "%s", "invalid argument for seed()"));
  225.     }
  226.    
  227.     return (tp_None);
  228. }
  229.  
  230. static tp_obj random_getstate(TP)
  231. {
  232.     tp_obj state_list = tp_list(tp);
  233.     int i;
  234.  
  235.     for (i = 0; i < N; i++) {
  236.         _tp_list_append(tp, state_list.list.val, tp_number(_gRandom.state[i]));
  237.     }
  238.     _tp_list_append(tp, state_list.list.val, tp_number(_gRandom.index));
  239.  
  240.     return (state_list);
  241. }
  242.  
  243. /*
  244.  * @state_list  must contain exactly N+1(625) integer.
  245.  */
  246. static tp_obj random_setstate(TP)
  247. {
  248.     tp_obj state_list = TP_OBJ();
  249.     tp_obj state_elem;
  250.     tp_obj len;
  251.     int i;
  252.  
  253.     len = tp_len(tp, state_list);
  254.     if (len.number.val != N+1) {
  255.         tp_raise(tp_None,tp_printf(tp, "%s: state vector's size invalid(should be %d)",
  256.                 __func__, N+1));
  257.     }
  258.  
  259.     for (i = 0; i < N; i++) {
  260.         state_elem = tp_get(tp, state_list, tp_number(i));
  261.         _gRandom.state[i] = (unsigned long)state_elem.number.val;
  262.     }
  263.     state_elem = tp_get(tp, state_list, tp_number(i));
  264.     _gRandom.index = (int)state_elem.number.val;
  265.  
  266.     return (tp_None);
  267. }
  268.  
  269. /*
  270.  * Jumpahead should be a fast way advance the generator n-steps ahead, but
  271.  * lacking a formula for that, the next best is to use n and the existing
  272.  * state to create a new state far away from the original.
  273.  *
  274.  * The generator uses constant spaced additive feedback, so shuffling the
  275.  * state elements ought to produce a state which would not be encountered
  276.  * (in the near term) by calls to random().  Shuffling is normally
  277.  * implemented by swapping the ith element with another element ranging
  278.  * from 0 to i inclusive.  That allows the element to have the possibility
  279.  * of not being moved.  Since the goal is to produce a new, different
  280.  * state, the swap element is ranged from 0 to i-1 inclusive.  This assures
  281.  * that each element gets moved at least once.
  282.  *
  283.  * To make sure that consecutive calls to jumpahead(n) produce different
  284.  * states (even in the rare case of involutory shuffles), i+1 is added to
  285.  * each element at position i.  Successive calls are then guaranteed to
  286.  * have changing (growing) values as well as shuffled positions.
  287.  *
  288.  * Finally, the self->index value is set to N so that the generator itself
  289.  * kicks in on the next call to random().  This assures that all results
  290.  * have been through the generator and do not just reflect alterations to
  291.  * the underlying state.
  292.  */
  293. static tp_obj random_jumpahead(TP)
  294. {
  295.     long n = (long)TP_NUM();
  296.     long i, j;
  297.     unsigned long *mt;
  298.     unsigned long tmp;
  299.  
  300.     mt = _gRandom.state;
  301.  
  302.     for (i = N-1; i > 1; i--) {
  303.         j = n % i;
  304.         if (j == -1L) {
  305.             tp_raise(tp_None,tp_printf(tp, "error: %s: j = %ld", __func__, j));
  306.         }
  307.         tmp   = mt[i];
  308.         mt[i] = mt[j];
  309.         mt[j] = tmp;
  310.     }
  311.  
  312.     for (i = 0; i < N; i++)
  313.         mt[i] += i + 1;
  314.  
  315.     _gRandom.index = N;
  316.  
  317.     return (tp_None);
  318. }
  319.  
  320. /* random_random is the function named genrand_res53 in the original code;
  321.  * generates a random number on [0,1) with 53-bit resolution; note that
  322.  * 9007199254740992 == 2**53; I assume they're spelling "/2**53" as
  323.  * multiply-by-reciprocal in the (likely vain) hope that the compiler will
  324.  * optimize the division away at compile-time.  67108864 is 2**26.  In
  325.  * effect, a contains 27 random bits shifted left 26, and b fills in the
  326.  * lower 26 bits of the 53-bit numerator.
  327.  * The orginal code credited Isaku Wada for this algorithm, 2002/01/09.
  328.  */
  329. static tp_obj random_random(TP)
  330. {
  331.     RandomObject *self = &_gRandom;
  332.     unsigned long a, b;
  333.  
  334.     if (! self->has_seed)
  335.         random_seed(tp);
  336.  
  337.     a = genrand_int32(self)>>5;
  338.     b = genrand_int32(self)>>6;
  339.    
  340.     return tp_number((a*67108864.0+b)*(1.0/9007199254740992.0));
  341. }
  342.  
  343.  
  344. #elif defined(__USE_WICHMANN_HILL_RANDOM_GENERATOR)
  345.  
  346. /*************************************************
  347.  * divmod(x, y)
  348.  * a helper function borrowed from Python
  349.  *************************************************/
  350. /* Compute Python divmod(x, y), returning the quotient and storing the
  351.  * remainder into *r.  The quotient is the floor of x/y, and that's
  352.  * the real point of this.  C will probably truncate instead (C99
  353.  * requires truncation; C89 left it implementation-defined).
  354.  * Simplification:  we *require* that y > 0 here.  That's appropriate
  355.  * for all the uses made of it.  This simplifies the code and makes
  356.  * the overflow case impossible (divmod(LONG_MIN, -1) is the only
  357.  * overflow case).
  358.  */
  359. #include <assert.h>
  360.  
  361. static int
  362. divmod(int x, int y, int *r)
  363. {
  364.     int quo;
  365.  
  366.     assert(y > 0);
  367.     quo = x / y;
  368.     *r = x - quo * y;
  369.     if (*r < 0) {
  370.         --quo;
  371.         *r += y;
  372.     }
  373.     assert(0 <= *r && *r < y);
  374.     return quo;
  375. }
  376.  
  377.  
  378. typedef struct WH_RandomObject {
  379.     struct seed {
  380.         unsigned long   x;
  381.         unsigned long   y;
  382.         unsigned long   z;
  383.     } seed;
  384.     int has_seed;
  385. } WH_RandomObject;
  386.  
  387. static WH_RandomObject _gWhRandom;
  388.  
  389. static tp_obj random_seed(TP)
  390. {
  391.     long a;
  392.     int x, y, z;
  393.     tp_obj arg = TP_DEFAULT(tp_None);
  394.  
  395.     if (arg.type == TP_NONE) {
  396.         time_t now;
  397.  
  398.         (void)time(&now);
  399.         a = (long)now * 256;
  400.     } else if (arg.type == TP_NUMBER) {
  401.         a = (long)arg.number.val;
  402.     } else {
  403.         tp_raise(tp_None,tp_printf(tp, "%s", "invalid argument for seed()"));
  404.     }
  405.    
  406.     a = divmod(a, 30268, &x);
  407.     a = divmod(a, 30306, &y);
  408.     a = divmod(a, 30322, &z);
  409.     _gWhRandom.seed.x = (int)x + 1;
  410.     _gWhRandom.seed.y = (int)y + 1;
  411.     _gWhRandom.seed.z = (int)z + 1;
  412.     _gWhRandom.has_seed = 1;
  413.  
  414.     return (tp_None);
  415. }
  416.  
  417. /*
  418.  * following comments are from Python's random.py
  419.  *---------------------------------------
  420.  * Wichman-Hill random number generator.
  421.  *
  422.  * Wichmann, B. A. & Hill, I. D. (1982)
  423.  * Algorithm AS 183:
  424.  * An efficient and portable pseudo-random number generator
  425.  * Applied Statistics 31 (1982) 188-190
  426.  *
  427.  * see also:
  428.  *        Correction to Algorithm AS 183
  429.  *        Applied Statistics 33 (1984) 123
  430.  *
  431.  *        McLeod, A. I. (1985)
  432.  *        A remark on Algorithm AS 183
  433.  *        Applied Statistics 34 (1985),198-200
  434.  * This part is thread-unsafe:
  435.  * BEGIN CRITICAL SECTION
  436.  */
  437. static tp_obj random_random(TP)
  438. {
  439.     long x, y, z;
  440.     double r;
  441.  
  442.     if (! _gWhRandom.has_seed)
  443.         random_seed(tp);
  444.  
  445.     x = _gWhRandom.seed.x;
  446.     y = _gWhRandom.seed.y;
  447.     z = _gWhRandom.seed.z;
  448.  
  449.     x = (171 * x) % 30269;
  450.     y = (172 * y) % 30307;
  451.     z = (170 * z) % 30323;
  452.  
  453.     _gWhRandom.seed.x = x;
  454.     _gWhRandom.seed.y = y;
  455.     _gWhRandom.seed.z = z;
  456.  
  457.     /*
  458.      * Note:  on a platform using IEEE-754 double arithmetic, this can
  459.      * never return 0.0 (asserted by Tim; proof too long for a comment).
  460.      */
  461.     errno = 0;
  462.     r = fmod(((double)x/30269.0+(double)y/30307.0+(double)z/30323.0), 1.0);
  463.     if (errno == EDOM)
  464.         tp_raise(tp_None,tp_printf(tp, "%s", "fmod(): denominator can't be zero"));
  465.  
  466.     return tp_number(r);
  467. }
  468.  
  469. /*
  470.  * for compatibility
  471.  */
  472. static tp_obj random_setstate(TP)
  473. {
  474.     return (tp_None);
  475. }
  476.  
  477. /*
  478.  * for compatibility
  479.  */
  480. static tp_obj random_getstate(TP)
  481. {
  482.     return (tp_None);
  483. }
  484.  
  485. /*
  486.  * FIXME: risk of overflow.
  487.  * following comments are from Python's random.py
  488.  * --------------------------------------------
  489.  * Act as if n calls to random() were made, but quickly.
  490.  *
  491.  * n is an int, greater than or equal to 0.
  492.  *
  493.  * Example use:  If you have 2 threads and know that each will
  494.  * consume no more than a million random numbers, create two Random
  495.  * objects r1 and r2, then do
  496.  *      r2.setstate(r1.getstate())
  497.  *      r2.jumpahead(1000000)
  498.  * Then r1 and r2 will use guaranteed-disjoint segments of the full
  499.  * period.
  500.  */
  501. static tp_obj random_jumpahead(TP)
  502. {
  503.     int n = (int)TP_NUM();
  504.     long x, y, z;
  505.  
  506.     if (n < 0)
  507.         tp_raise(tp_None,tp_printf(tp, "%s: n = %d invalid, should >= 0", __func__, n));
  508.  
  509.     x = _gWhRandom.seed.x;
  510.     y = _gWhRandom.seed.y;
  511.     z = _gWhRandom.seed.z;
  512.  
  513.     x = (int)(x * ((long)pow(171, n) % 30269)) % 30269;
  514.     y = (int)(y * ((long)pow(172, n) % 30307)) % 30307;
  515.     z = (int)(z * ((long)pow(170, n) % 30323)) % 30323;
  516.  
  517.     _gWhRandom.seed.x = x;
  518.     _gWhRandom.seed.y = y;
  519.     _gWhRandom.seed.z = z;
  520.  
  521.     return (tp_None);
  522. }
  523.  
  524.  
  525. #elif defined(__USE_POSIX_RANDOM_GENERATOR)
  526.  
  527. /*
  528.  * judge whether seeded
  529.  */
  530. static int has_seed = 0;
  531.  
  532. static tp_obj random_seed(TP)
  533. {
  534.     tp_obj arg = TP_DEFAULT(tp_None);
  535.  
  536.     if (arg.type == TP_NONE) {
  537.         time_t now;
  538.  
  539.         (void)time(&now);
  540.         srandom((unsigned int)now);
  541.         has_seed = 1;
  542.     } else if (arg.type == TP_NUMBER) {
  543.         srandom((unsigned long)arg.number.val);
  544.         has_seed = 1;
  545.     } else {
  546.         tp_raise(tp_None,tp_printf(tp, "%s", "invalid argument for seed()"));
  547.     }
  548.    
  549.     return (tp_None);
  550. }
  551.  
  552. /*
  553.  * random()
  554.  *
  555.  * generate successive pseudo random number ranging from [0.0, 1.0).
  556.  * usually RAND_MAX is huge number, thus the periods of the success-
  557.  * ive random number is very long, about 16*((2**31)-1).
  558.  * NOTE: if seed() not called before random(), random() will
  559.  * automatically call seed() with current time.
  560.  */
  561. tp_obj random_random(TP)
  562. {
  563.     double r = 0.0;
  564.  
  565.     if (! has_seed)
  566.         random_seed(tp);
  567.  
  568.     r = (tp_num)random()/(tp_num)RAND_MAX;
  569.  
  570.     return (tp_number(r));
  571. }
  572.  
  573. /*
  574.  * setstate(state)
  575.  *
  576.  * for compatibility.
  577.  */
  578. tp_obj random_setstate(TP)
  579. {
  580.     return (tp_None);
  581. }
  582.  
  583. /*
  584.  * getstate()
  585.  *
  586.  * for compatibility.
  587.  */
  588. tp_obj random_getstate(TP)
  589. {
  590.     return (tp_None);
  591. }
  592.  
  593. /*
  594.  * jumpahead()
  595.  *
  596.  * for compatibility.
  597.  */
  598. tp_obj random_jumpahead(TP)
  599. {
  600.     return (tp_None);
  601. }
  602.  
  603. #else
  604.  
  605. #error no underlying random generator is specified
  606.  
  607. #endif
  608.  
  609. /************************************************************
  610.  * some usual distributions
  611.  ************************************************************/
  612.  
  613. /*
  614.  * return real number in range [a, b)
  615.  * a and b can be negtive, but a must be less than b.
  616.  */
  617. tp_obj random_uniform(TP)
  618. {
  619.     double a = TP_NUM();
  620.     double b = TP_NUM();
  621.     double r = 0.0;
  622.     tp_obj rvo;         /* random variable object */
  623.  
  624.     if (a >= b)
  625.         tp_raise(tp_None,tp_printf(tp, "%s: a(%f) must be less than b(%f)", a, b));
  626.  
  627.     rvo = random_random(tp);
  628.     r = a + (b - a) * rvo.number.val;
  629.    
  630.     return (tp_number(r));
  631. }
  632.  
  633. /*
  634.  * Normal distribution
  635.  * @mu      mean
  636.  * @sigma   standard deviation
  637.  *-----------------------------
  638.  * Uses Kinderman and Monahan method. Reference: Kinderman,
  639.  * A.J. and Monahan, J.F., "Computer generation of random
  640.  * variables using the ratio of uniform deviates", ACM Trans
  641.  * Math Software, 3, (1977), pp257-260.
  642.  */
  643. tp_obj random_normalvariate(TP)
  644. {
  645.     double mu = TP_NUM();
  646.     double sigma = TP_NUM();
  647.     double NV_MAGICCONST;
  648.     double u1, u2;
  649.     double z, zz;
  650.     double r = 0.0;
  651.     tp_obj rvo;         /* random variable object */
  652.  
  653.     NV_MAGICCONST = 4.0 * exp(-0.5) / sqrt(2.0);
  654.     while (1) {
  655.         rvo = random_random(tp);
  656.         u1  = rvo.number.val;
  657.         rvo = random_random(tp);
  658.         u2  = 1.0 - rvo.number.val;
  659.         z   = NV_MAGICCONST * (u1 - 0.5) / u2;
  660.         zz  = z * z / 4.0;
  661.         if (zz <= - log(u2))
  662.             break;
  663.     }
  664.  
  665.     r = mu + z * sigma;
  666.  
  667.     return (tp_number(r));
  668. }
  669.  
  670. /*
  671.  * Log normal distribution
  672.  *
  673.  * If take natural logarithm on log normal distribution, normal
  674.  * distribution with mean mu and standard deviation sigma will
  675.  * return.
  676.  * @mu      mean, can be any value
  677.  * @sigma   standard deviation, must be > 0.
  678.  */
  679. tp_obj random_lognormvariate(TP)
  680. {
  681.     double mu = TP_NUM();
  682.     double sigma = TP_NUM();
  683.     tp_obj params;
  684.     tp_obj normvar;     /* normal distribution variate */
  685.     double r = 0.0;
  686.  
  687.     /*
  688.      * call random_normalvariate() actually
  689.      */
  690.     params = tp_params_v(tp, 2, tp_number(mu), tp_number(sigma));
  691.     normvar = tp_ez_call(tp, "random", "normalvariate", params);
  692.     r = exp(normvar.number.val);
  693.  
  694.     return (tp_number(r));
  695. }
  696.  
  697. /*
  698.  * Exponential distribution
  699.  *
  700.  * @lambda      reciprocal of mean.
  701.  * return value range (0, +inf)
  702.  */
  703. tp_obj random_expovariate(TP)
  704. {
  705.     double lambda = TP_NUM();
  706.     double u, r;
  707.     tp_obj rvo;
  708.  
  709.     do {
  710.         rvo = random_random(tp);
  711.         u = rvo.number.val;
  712.     } while (u <= 0.0000001);
  713.  
  714.     r = -log(u) / lambda;
  715.  
  716.     return (tp_number(r));
  717. }
  718.  
  719. /*
  720.  * Circular data distribution.
  721.  *
  722.  * mu is the mean angle, expressed in radians between 0 and 2*pi, and
  723.  * kappa is the concentration parameter, which must be greater than or
  724.  * equal to zero.  If kappa is equal to zero, this distribution reduces
  725.  * to a uniform random angle over the range 0 to 2*pi.
  726.  *
  727.  * mu:    mean angle (in radians between 0 and 2*pi)
  728.  * kappa: concentration parameter kappa (>= 0)
  729.  * if kappa = 0 generate uniform random angle
  730.  *
  731.  * Based upon an algorithm published in: Fisher, N.I.,
  732.  * "Statistical Analysis of Circular Data", Cambridge
  733.  * University Press, 1993.
  734.  *
  735.  * Thanks to Magnus Kessler for a correction to the
  736.  * implementation of step 4.
  737.  */
  738. tp_obj random_vonmisesvariate(TP)
  739. {
  740.     double mu = TP_NUM();
  741.     double kappa = TP_NUM();
  742.     tp_obj rvo;
  743.     double a, b, c, r;
  744.     double u1, u2, u3, z, f;
  745.     double theta;
  746.     double TWOPI = 2.0 * M_PI;
  747.  
  748.     if (kappa <= 1e-6) {
  749.         rvo = random_random(tp);
  750.         theta = TWOPI * rvo.number.val;
  751.         return (tp_number(theta));
  752.     }
  753.  
  754.     a = 1.0 + sqrt(1.0 + 4.0 * kappa * kappa);
  755.     b = (a - sqrt(2.0 * a))/(2.0 * kappa);
  756.     r = (1.0 + b * b)/(2.0 * b);
  757.  
  758.     while (1) {
  759.         rvo = random_random(tp);
  760.         u1 = rvo.number.val;
  761.  
  762.         z = cos(M_PI * u1);
  763.         f = (1.0 + r * z)/(r + z);
  764.         c = kappa * (r - f);
  765.  
  766.         rvo = random_random(tp);
  767.         u2 = rvo.number.val;
  768.  
  769.         if ((u2 < (c * (2.0 - c))) ||
  770.             (u2 <= (c * exp(1.0 - c))))
  771.             break;
  772.     }
  773.  
  774.     rvo = random_random(tp);
  775.     u3 = rvo.number.val;
  776.     if (u3 > 0.5)
  777.         theta = fmod(mu, TWOPI) + acos(f);
  778.     else
  779.         theta = fmod(mu, TWOPI) - acos(f);
  780.  
  781.     return (tp_number(theta));
  782. }
  783.  
  784. /*
  785.  * Gamma distribution.  Not the gamma function!
  786.  *
  787.  * Conditions on the parameters are alpha > 0 and beta > 0.
  788.  */
  789. tp_obj random_gammavariate(TP)
  790. {
  791.     double alpha = TP_NUM();
  792.     double beta  = TP_NUM();
  793.     tp_obj rvo;
  794.     double res;
  795.     double LOG4 = log(4.0);
  796.     double SG_MAGICCONST = 1.0 + log(4.5);
  797.  
  798.     /*
  799.      * alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2
  800.      * Warning: a few older sources define the gamma distribution in terms
  801.      * of alpha > -1.0
  802.      */
  803.     if ((alpha <= 0.0) || (beta <= 0.0))
  804.         tp_raise(tp_None,tp_printf(tp, "%s: alpha(%f) and beta(%f) must be > 0.0",
  805.                 __func__, alpha, beta));
  806.  
  807.     if (alpha > 1.0) {
  808.  
  809.         /*
  810.          * Uses R.C.H. Cheng, "The generation of Gamma
  811.          * variables with non-integral shape parameters",
  812.          * Applied Statistics, (1977), 26, No. 1, p71-74
  813.          */
  814.  
  815.         double ainv;
  816.         double bbb, ccc;
  817.         double u1, u2;
  818.         double v, x, z, r;
  819.  
  820.         ainv = sqrt(2.0 * alpha - 1.0);
  821.         bbb = alpha - LOG4;
  822.         ccc = alpha + ainv;
  823.  
  824.         while (1) {
  825.             rvo = random_random(tp);
  826.             u1 = rvo.number.val;
  827.             if (! ((1e-7 < u1) && (u1 < 0.9999999)))
  828.                 continue;
  829.             rvo = random_random(tp);
  830.             u2 = 1.0 - rvo.number.val;
  831.             v = log(u1 / (1.0 - u1)) / ainv;
  832.             x = alpha * exp(v);
  833.             z = u1 * u1 * u2;
  834.             r = bbb + ccc * v - x;
  835.             if ((r + SG_MAGICCONST - 4.5 * z >= 0.0) ||
  836.                 (r >= log(z))) {
  837.                 res = x * beta;
  838.                 return (tp_number(res));
  839.             }
  840.         }
  841.     }
  842.     else if (alpha == 1.0) {
  843.  
  844.         /*
  845.          * expovariate(1)
  846.          */
  847.  
  848.         double u;
  849.  
  850.         do {
  851.             rvo = random_random(tp);
  852.             u = rvo.number.val;
  853.         } while (u <= 1e-7);
  854.  
  855.         res = - log(u) * beta;
  856.         return (tp_number(res));
  857.     } else {
  858.  
  859.         /*
  860.          * alpha is between 0 and 1 (exclusive)
  861.          *
  862.          * Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle
  863.          */
  864.  
  865.         double b, p, u, u1, x;
  866.  
  867.         while (1) {
  868.             rvo = random_random(tp);
  869.             u = rvo.number.val;
  870.             b = (M_E + alpha) / M_E;
  871.             p = b * u;
  872.             if (p <= 1.0)
  873.                 /*FIXME: x = p ** (1.0/alpha)*/
  874.                 x = pow(p, 1.0/alpha);
  875.             else
  876.                 x = - log((b - p) / alpha);
  877.             rvo = random_random(tp);
  878.             u1 = rvo.number.val;
  879.             if (p > 1.0) {
  880.                 /*FIXME: if u1 <= x ** (alpha - 1.0):*/
  881.                 if (u1 <= pow(x, alpha - 1.0))
  882.                     break;
  883.             }
  884.             else if (u1 <= exp(-x))
  885.                 break;
  886.         }
  887.  
  888.         res = x * beta;
  889.         return (tp_number(res));
  890.     }
  891. }
  892.  
  893. /*
  894.  * Beta distribution.
  895.  *
  896.  *        Conditions on the parameters are alpha > 0 and beta > 0.
  897.  *        Returned values range between 0 and 1.
  898.  *        
  899.  *        # This version due to Janne Sinkkonen, and matches all the std
  900.  *        # texts (e.g., Knuth Vol 2 Ed 3 pg 134 "the beta distribution").
  901.  *        
  902.  * See also:
  903.  * http://sourceforge.net/bugs/?func=detailbug&bug_id=130030&group_id=5470
  904.  * for Ivan Frohne's insightful analysis of why the original implementation:
  905.  *
  906.  *    def betavariate(self, alpha, beta):
  907.  *        # Discrete Event Simulation in C, pp 87-88.
  908.  *
  909.  *        y = self.expovariate(alpha)
  910.  *        z = self.expovariate(1.0/beta)
  911.  *        return z/(y+z)
  912.  *
  913.  * was dead wrong, and how it probably got that way.
  914.  *
  915.  */
  916. tp_obj random_betavariate(TP)
  917. {
  918.     double alpha = TP_NUM();
  919.     double beta  = TP_NUM();
  920.     double t;
  921.     double r = 0.0;    
  922.     tp_obj y;
  923.     tp_obj params;
  924.  
  925.     params = tp_params_v(tp, 2, tp_number(alpha), tp_number(1.0));
  926.     y = tp_ez_call(tp, "random", "gammavariate", params);
  927.     if (y.number.val == 0) {
  928.         return (y);
  929.     } else {
  930.         params = tp_params_v(tp, 2, tp_number(beta), tp_number(1.0));
  931.         t = y.number.val;
  932.         y = tp_ez_call(tp, "random", "gammavariate", params);
  933.         r = t / (t + y.number.val);
  934.         return (tp_number(r));
  935.     }
  936. }
  937.  
  938. /*
  939.  * Pareto distribution.  alpha is the shape parameter.
  940.  * # Jain, pg. 495
  941.  */
  942. tp_obj random_paretovariate(TP)
  943. {
  944.     double alpha = TP_NUM();
  945.     double u;
  946.     double r;
  947.     tp_obj rvo;
  948.    
  949.     rvo = random_random(tp);
  950.     u = 1.0 - rvo.number.val;
  951.     r = 1.0 / pow(u, 1.0/alpha);
  952.    
  953.     return (tp_number(r));
  954. }
  955.  
  956. /*
  957.  * Weibull distribution.
  958.  *
  959.  * alpha is the scale parameter and beta is the shape parameter.
  960.  *
  961.  * Jain, pg. 499; bug fix courtesy Bill Arms
  962.  */
  963. tp_obj random_weibullvariate(TP)
  964. {
  965.     double alpha = TP_NUM();
  966.     double beta  = TP_NUM();
  967.     double u;
  968.     double r;
  969.     tp_obj rvo;
  970.    
  971.     rvo = random_random(tp);
  972.     u = 1.0 - rvo.number.val;
  973.     r = alpha * pow(-log(u), 1.0/beta);
  974.     return (tp_number(r));
  975. }
  976.  
  977. /*
  978.  * randomly select an element from range ([start,] stop[, step])
  979.  *
  980.  * 'stop' must be larger than 'start', both can be negative;
  981.  * 'step' must be integer larger than zero.
  982.  */
  983. tp_obj random_randrange(TP)
  984. {
  985.     tp_obj start = TP_OBJ();
  986.     tp_obj stop = TP_DEFAULT(tp_None);
  987.     tp_obj step = TP_DEFAULT(tp_number(1));
  988.     tp_obj rvo = random_random(tp);
  989.     int istart = (int)start.number.val;
  990.     int istep = (int)step.number.val;
  991.     int istop;
  992.     int iwidth;
  993.     double res;
  994.    
  995.     if (stop.type == TP_NONE) {
  996.         /*
  997.                         * if only one argument, then start just means stop
  998.                         */
  999.         istop = istart;
  1000.         res = (rvo.number.val * istop);
  1001.         return (tp_number(res));
  1002.     } else if (stop.type == TP_NUMBER) {
  1003.         istop = (int)stop.number.val;
  1004.         iwidth = istop - istart;
  1005.         if (iwidth < 0)
  1006.             tp_raise(tp_None,tp_printf(tp, "%s", "stop must be > start"));
  1007.         if (istep <= 0)
  1008.             tp_raise(tp_None,tp_printf(tp, "%s", "step must be integer larger than 0"));
  1009.            
  1010.         if (istep == 1) {
  1011.             res = (int)(istart + (int)(rvo.number.val * iwidth));
  1012.             return (tp_number(res));
  1013.         } else {
  1014.             int n = (iwidth + istep - 1) / istep;
  1015.             res = (int)(istart + istep * (int)(n * rvo.number.val));
  1016.             return (tp_number(res));
  1017.         }
  1018.     } else {
  1019.         tp_raise(tp_None,tp_printf(tp, "%s", "wrong type of stop"));
  1020.     }
  1021. }
  1022.  
  1023. /*
  1024.  * return random integer between [a, b]
  1025.  */
  1026. tp_obj random_randint(TP)
  1027. {
  1028.     double a = TP_NUM();
  1029.     double b = TP_NUM();
  1030.     tp_obj r;
  1031.     tp_obj params;
  1032.  
  1033.     params = tp_params_v(tp, 2, tp_number(a), tp_number(b + 1));
  1034.     r = tp_ez_call(tp, "random", "randrange", params);
  1035.     return (r);
  1036. }
  1037.  
  1038. /*
  1039.  * return a random element of sequence 'seq'. 'seq' mustn't be empty.
  1040.  */
  1041. tp_obj random_choice(TP)
  1042. {
  1043.     tp_obj seq = TP_OBJ();
  1044.     tp_obj len;
  1045.     tp_obj rvo;
  1046.     tp_obj r;
  1047.     int i;
  1048.    
  1049.     len = tp_len(tp, seq);
  1050.     if (len.number.val <= 0)
  1051.         tp_raise(tp_None,tp_printf(tp, "%s", "seq mustn't be empty"));
  1052.    
  1053.     rvo = random_random(tp);
  1054.     i = (int)(len.number.val * rvo.number.val);
  1055.     r = tp_get(tp, seq, tp_number(i));
  1056.    
  1057.     return (r);
  1058. }
  1059.  
  1060. /*
  1061.  * shuffle sequence 'seq' in place, return None
  1062.  */
  1063. tp_obj random_shuffle(TP)
  1064. {
  1065.     tp_obj seq = TP_OBJ();
  1066.     tp_obj elmi;
  1067.     tp_obj elmj;
  1068.     tp_obj params;
  1069.     tp_obj rvo;
  1070.     tp_obj len = tp_len(tp, seq);
  1071.     int i, j;
  1072.    
  1073.     if (len.number.val <= 0)
  1074.         return (tp_None);
  1075.    
  1076.     for (i = len.number.val - 1; i > len.number.val / 2; i--) {
  1077.         /*
  1078.                        * randomly exchange elment i and elment j, element i from the behind end of 'seq', while
  1079.                        * element j from the front end of 'seq'.
  1080.                        */
  1081.         params = tp_params_v(tp, 2, tp_number(0), tp_number(len.number.val / 2));
  1082.         rvo = tp_ez_call(tp, "random", "randint", params);
  1083.         j = (int)rvo.number.val;
  1084.         elmi = tp_get(tp, seq, tp_number(i));
  1085.         elmj = tp_get(tp, seq, tp_number(j));        
  1086.        
  1087.         tp_set(tp, seq, tp_number(i), elmj);
  1088.         tp_set(tp, seq, tp_number(j), elmi);
  1089.     }
  1090.    
  1091.     for (i = len.number.val / 2; i >= 0; i--) {
  1092.         /*
  1093.                        * randomly exchange elment i and elment j, element i from the front end of 'seq', while
  1094.                        * element j from the behind end of 'seq'.
  1095.                        */
  1096.         params = tp_params_v(tp, 2, tp_number(len.number.val / 2), tp_number(len.number.val - 1));
  1097.         rvo = tp_ez_call(tp, "random", "randint", params);
  1098.         j = (int)rvo.number.val;
  1099.         elmi = tp_get(tp, seq, tp_number(i));
  1100.         elmj = tp_get(tp, seq, tp_number(j));
  1101.        
  1102.         tp_set(tp, seq, tp_number(i), elmj);
  1103.         tp_set(tp, seq, tp_number(j), elmi);
  1104.     }
  1105.    
  1106.     return (tp_None);
  1107. }
  1108.