Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | Download | RSS feed

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