home *** CD-ROM | disk | FTP | other *** search
/ Fish 'n' More 2 / fishmore-publicdomainlibraryvol.ii1991xetec.iso / fish / applications / xlispstat / src / src2.lzh / XLisp-Stat / uni.c < prev    next >
C/C++ Source or Header  |  1990-10-02  |  3KB  |  166 lines
  1. /* uni - uniform random number generator - Marsaglia's portable gen.   */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include <time.h>
  8. #include "xlisp.h"
  9. #include "osdef.h"
  10. #ifdef ANSI
  11. #include "xlproto.h"
  12. #include "xlsproto.h"
  13. #else
  14. 3include "xlfun.h"
  15. #include "xlsfun.h"
  16. #endif ANSI
  17. #include "xlsvar.h"
  18.  
  19. #ifdef ANSI
  20. double uni1(void);
  21. void checkstate(LVAL);
  22. int random_state_p(LVAL);
  23. #else
  24. double uni1();
  25. void checkstate();
  26. int random_state_p();
  27. #endif ANSI
  28.  
  29. #define STATELENGTH 22
  30. #define MDIG 32
  31. #define M1 18
  32. #define M2 19
  33. #define I  20
  34. #define J  21
  35.  
  36. /* seed: seed, m[1] - m[17], m1, m2, i, j */
  37.  
  38. LVAL newrandomstate()
  39. {
  40.   LVAL state;
  41.   unsigned long seed;
  42.   long m1, m2, k0, k1, j0, j1;
  43.   int i;
  44.   
  45.   xlsave1(state);
  46.   state = newvector(STATELENGTH);
  47.   
  48.   /* set seed from clock */
  49.   seed = time((unsigned long *) NULL);
  50.   
  51.   m1 = floor(pow(2.0, (MDIG - 2.0)) + .1)
  52.      + floor(pow(2.0, (MDIG - 2.0)) + .1) - 1;
  53.   m2 = floor(pow(2.0, (MDIG / 2.0)) + .1);
  54.  
  55.   setelement(state, 0, cvfixnum((FIXTYPE) seed));
  56.  
  57.   setelement(state, M1, cvfixnum((FIXTYPE) m1));
  58.   setelement(state, M2, cvfixnum((FIXTYPE) m2));
  59.   setelement(state,  I,  cvfixnum((FIXTYPE) 5));
  60.   setelement(state,  J,  cvfixnum((FIXTYPE) 17));
  61.  
  62.   if (seed % 2 == 0) seed--;
  63.   k0 = 9069 % m2;
  64.   k1 = 9069 / m2;
  65.   j0 = seed % m2;
  66.   j1 = seed / m2;
  67.   for (i = 1; i <= 17; i++) {
  68.     seed = j0 * k0;
  69.     j1 = ((seed / m2) + j0 * k1 + j1 * k0) % (m2 / 2);
  70.     j0 = seed % m2;
  71.     setelement(state, i, cvfixnum((FIXTYPE) j0 + m2 * j1));
  72.   }
  73.   xlpop();
  74.   return(state);
  75. }
  76.  
  77. static double uni1()
  78. {
  79.   LVAL state;
  80.   long k, i, j, m1;
  81.  
  82.   state = getvalue(s_random_state);
  83.   checkstate(state);
  84.   
  85.   m1 = getfixnum(getelement(state, M1));
  86.   i =  getfixnum(getelement(state,  I));
  87.   j =  getfixnum(getelement(state,  J));
  88.  
  89.   k =  getfixnum(getelement(state, i)) - getfixnum(getelement(state, j));
  90.   if (k < 0) k = k + m1;
  91.   setelement(state, j, cvfixnum((FIXTYPE) k));
  92.   i = i - 1;
  93.   if (i <= 0) i = 17;
  94.   j = j - 1;
  95.   if (j <= 0) j = 17;
  96.   setelement(state, I, cvfixnum((FIXTYPE) i));
  97.   setelement(state, J, cvfixnum((FIXTYPE) j));
  98.   
  99.   return(((double) k) / m1);
  100. }
  101.  
  102. double uni()
  103. {
  104.   double x;
  105.   do { x = uni1(); } while (x <= 0.0 || x >= 1.0);
  106.   return (x);
  107. }
  108.  
  109. int osrand(n)
  110.   unsigned n;
  111. {
  112.   int k;
  113.   
  114.   do {
  115.     k = n * uni();
  116.   } while (k < 0 || k >= n);
  117.   return(k);
  118. }
  119.  
  120. LVAL xsmake_random_state()
  121. {
  122.   LVAL arg = NIL;
  123.   
  124.   if (moreargs()) arg = xlgetarg();
  125.   xllastarg();
  126.   
  127.   if (arg == NIL) return(copyvector(getvalue(s_random_state)));
  128.   else if (arg == s_true) return(newrandomstate());
  129.   else {
  130.     checkstate(arg);
  131.     return(copyvector(arg));
  132.   }
  133. }
  134.  
  135. LVAL xsrandom_state_p()
  136. {
  137.   LVAL state;
  138.   
  139.   state = xlgetarg();
  140.   return((random_state_p(state)) ? s_true : NIL);
  141. }
  142.  
  143. static random_state_p(state)
  144.     LVAL state;
  145. {
  146.   long i, j, m1;
  147.   
  148.   if (! vectorp(state) || getsize(state) != STATELENGTH)
  149.     return(FALSE);
  150.     
  151.   m1 = getfixnum(getelement(state, M1));
  152.   i =  getfixnum(getelement(state,  I));
  153.   j =  getfixnum(getelement(state,  J));
  154.  
  155.   if (m1 <= 0 || i <= 0 || i > 17 || j <= 0 || j > 17)
  156.     return(FALSE);
  157.  
  158.   return(TRUE);
  159. }
  160.  
  161. static void checkstate(state)
  162.     LVAL state;
  163. {
  164.   if (! random_state_p(state)) xlerror("bad random state", state);
  165. }
  166.