home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 8 / FreshFishVol8-CD2.bin / bbs / comm / amitcp-3.0ß2.lha / AmiTCP / src / devs / agnet / lrandom.c < prev    next >
C/C++ Source or Header  |  1993-10-08  |  7KB  |  240 lines
  1. RCS_ID_C="$Id: lrandom.c,v 3.1 93/10/07 19:24:23 ppessi Exp $";
  2. /*
  3.  * random.c --- generate random integers
  4.  *
  5.  * Author: ppessi <Pekka.Pessi@hut.fi>
  6.  *
  7.  * Copyright (c) 1993 OHT-AmiTCP/IP Group,
  8.  *                    Helsinki University of Technology, Finland.
  9.  *                    All rights reserved.
  10.  *
  11.  * Created      : Wed Feb 24 08:24:13 1993 ppessi
  12.  * Last modified: Thu Oct  7 18:20:20 1993 ppessi
  13.  *
  14.  * $Log:    lrandom.c,v $
  15.  * Revision 3.1  93/10/07  19:24:23  ppessi
  16.  * Release 2.1 version
  17.  * 
  18.  * Revision 2.0  93/05/14  16:50:37  ppessi
  19.  * Release version.
  20.  * 
  21.  */
  22.  
  23. #include <exec/types.h>
  24. #include <exec/io.h>
  25. #include <devices/timer.h>
  26.  
  27. #include <clib/timer_protos.h>
  28. #include <pragmas/timer_pragmas.h>
  29.  
  30. #ifndef TEST
  31. #include "agnet.h"
  32. #include "bases.h"
  33. #else
  34. extern struct Device *TimerBase ;
  35. #endif
  36.  
  37. /****** random.module/LRandom ******************************************
  38. *
  39. *   NAME    
  40. *    LRandom -- Random number generator
  41. *
  42. *   SYNOPSIS
  43. *        LRandom()
  44. *
  45. *    LONG InitLRandom(VOID)
  46. *
  47. *   FUNCTION
  48. *      Generates a random long integer. The random number generator
  49. *      must have been initialized before calling LRandom() or strange
  50. *      things will happen.
  51. *
  52. *      LRandom() generates fairly good random numbers, all bits in the
  53. *      long word seem to be useful.
  54. *
  55. *   BUGS
  56. *      None known.
  57. *
  58. *   SEE ALSO
  59. *    InitLRandom, timer.device/GetSysTime()
  60. *
  61. ******************************************************************************
  62. */
  63. static union random_seed {
  64.   struct timeval time;
  65.   ULONG          longs[2];
  66.   UWORD          words[4];
  67. } seed;
  68.  
  69. static const UWORD random_offset[4] = { 0x7823, 0xab34, 0x93b4, 0x7673 };
  70.  
  71. static const UWORD random_eor[4] = { 0xc97d, 0x6988, 0x32e9, 0x8487 };
  72.  
  73. ULONG LRandom(VOID)
  74. {
  75.   ULONG carry = 0;
  76.   int  i;
  77.  
  78.   for (i = 3; i >= 0; i--) {
  79.     carry += seed.words[i]*34639L + random_offset[i];
  80.     seed.words[i] = carry ^ random_eor[i];
  81.     carry >>= 16;
  82.   }
  83.  
  84.   return (ULONG)(seed.words[3] << 24) + (seed.words[2]<<16) + 
  85.     (seed.words[1] << 8) + seed.words[0];
  86. }
  87.  
  88. /****** random.module/InitLRandom ******************************************
  89. *
  90. *   NAME    
  91. *    InitLRandom -- initialize LRandom generator
  92. *
  93. *   SYNOPSIS
  94. *        InitLRandom()
  95. *
  96. *    VOID InitLRandom(VOID)
  97. *
  98. *   FUNCTION
  99. *      Initializes the random number generator with the current time.
  100. *
  101. *   NOTES
  102. *      This function uses GetSysTime() to get the initial seed for the
  103. *      random number generator. Thus it needs the V36 timer.device.
  104. *
  105. *   BUGS
  106. *      None known.
  107. *
  108. *   SEE ALSO
  109. *    LRandom, timer.device/GetSysTime()
  110. *
  111. ******************************************************************************
  112. */
  113. VOID InitLRandom(VOID)
  114. {
  115.   GetSysTime(&seed.time);
  116.   (VOID) LRandom();  (VOID) LRandom();  (VOID) LRandom();  (VOID) LRandom();
  117.   (VOID) LRandom();  (VOID) LRandom();  (VOID) LRandom();  (VOID) LRandom();
  118. }
  119.  
  120. /****** random.module/RandomDev ******************************************
  121. *
  122. *   NAME    
  123. *    RandomDev -- Binomial distribution generator
  124. *
  125. *   SYNOPSIS
  126. *        RandomDev(mean, deviation)
  127. *
  128. *    ULONG RandomDev(ULONG mean, ULONG deviation)
  129. *
  130. *   FUNCTION
  131. *      RandomDev generates a random number in a binomial (n=16)
  132. *      distribution, which have specified mean and standard deviation.
  133. *      Whole distribution should be in the domain of the ULONG,
  134. *      otherwise the mean and/or the deviation of generated numbers
  135. *      may have unexpected values.
  136. *
  137. *   BUGS
  138. *      There are only a 16 distinct values for the deviation.
  139. *
  140. *   SEE ALSO
  141. *    InitLRandom, timer.device/GetSysTime()
  142. *
  143. ******************************************************************************
  144. */
  145. static const UWORD sdev[257] = 
  146. { 65535, 63519, 61565, 59671, 57835, 56056, 54332, 52661, 
  147.  51041, 49471, 47950, 46476, 45048, 43664, 42322, 41022, 
  148.  39762, 38542, 37359, 36212, 35102, 34025, 32982, 31972, 
  149.  30992, 30044, 29124, 28233, 27370, 26534, 25723, 24938, 
  150.  24177, 23440, 22725, 22033, 21363, 20713, 20083, 19473, 
  151.  18882, 18310, 17755, 17217, 16696, 16192, 15702, 15229, 
  152.  14770, 14325, 13894, 13476, 13071, 12679, 12299, 11931, 
  153.  11574, 11229, 10894, 10569, 10255, 9950, 9655, 9369, 
  154.  9091, 8823, 8562, 8310, 8065, 7828, 7599, 7376, 
  155.  7160, 6951, 6749, 6552, 6362, 6177, 5999, 5825, 
  156.  5657, 5494, 5337, 5184, 5035, 4891, 4752, 4617, 
  157.  4486, 4359, 4236, 4116, 4000, 3888, 3779, 3674, 
  158.  3571, 3472, 3376, 3282, 3192, 3104, 3018, 2936, 
  159.  2856, 2778, 2702, 2629, 2558, 2489, 2422, 2357, 
  160.  2294, 2233, 2174, 2116, 2061, 2006, 1954, 1903, 
  161.  1853, 1805, 1758, 1713, 1668, 1626, 1584, 1544, 
  162.  1504, 1466, 1429, 1393, 1358, 1324, 1291, 1259, 
  163.  1228, 1197, 1168, 1139, 1111, 1084, 1058, 1032, 
  164.  1007,  983,  959,  936,  914,  892,  871,  850, 
  165.   830,  811,  792,  773,  755,  738,  720,  704, 
  166.   687,  671,  656,  641,  626,  612,  598,  584, 
  167.   571,  558,  545,  533,  521,  509,  497,  486, 
  168.   475,  464,  454,  443,  433,  423,  414,  404, 
  169.   395,  386,  377,  368,  360,  352,  343,  335, 
  170.   328,  320,  312,  305,  298,  291,  284,  277, 
  171.   270,  263,  257,  250,  244,  238,  232,  226, 
  172.   220,  214,  208,  203,  197,  192,  186,  181, 
  173.   176,  170,  165,  160,  155,  150,  145,  141, 
  174.   136,  131,  126,  122,  117,  113,  108,  104, 
  175.    99,   95,   90,   86,   82,   78,   73,   69, 
  176.    65,   61,   57,   53,   48,   44,   40,   36, 
  177.    32,   28,   24,   20,   16,   12,    8,    4 };
  178.  
  179. const static BYTE ones[16] = 
  180. { 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4 };
  181.  
  182. ULONG RandomDev(ULONG mean, ULONG deviation)
  183. {
  184.   WORD m, k, actual, done;
  185.   ULONG relative;
  186.   ULONG x = LRandom();
  187.   WORD offset;
  188.  
  189.   /* Define the multiplier (actually mpr/256) */
  190.   relative = (deviation<<8)/mean;
  191.   /* find suitable multiplier */
  192.   m = 128; k=64;
  193.   while (k) {
  194.     if (sdev[m] > relative) 
  195.       m += k;
  196.     else 
  197.       m -= k;
  198.     k >>=1;
  199.   }
  200.  
  201.   /*
  202.    * Count '1' bits in a lower word of x
  203.    */
  204.   actual = ones[x >> 12 & 0xf] + ones[x >> 8 & 0xf] 
  205.     + ones[x >> 4 & 0xf] + ones[x & 0xf];
  206.   /* 
  207.    * delay * (m/256)**(16-actual) * ((512-m)/256)**actual
  208.    */
  209.   x = mean; offset = 0;
  210.   /* scale x */ 
  211.   while (x < 0x4000) {
  212.     x <<= 8; offset += 8;
  213.   }
  214.   while (x < 0x40000) {
  215.     x <<= 4; offset += 4;
  216.   }
  217.   while (x < 0x400000) {
  218.     x <<= 1; offset += 1;
  219.   }
  220.  
  221.   done = actual > 16 - actual ? 16 - actual : actual;
  222.   for (k = done; k > 0; k--) {
  223.     x = (x * (512 - m) >> 8) * m >> 8;
  224.   }
  225.   for (k = done; k < 16 - done; k++) {
  226.     if (k < actual) {
  227.       x = x * (512 - m) >> 8;
  228.       /* multiplication may overflow without this scaling */
  229.       while (x >= 0x800000) { x >>= 1; offset--; }
  230.     }
  231.     else if (k >= actual) {
  232.       /* we do not miss the lost accuracy, do we? */ 
  233.       x = x * m >> 8;
  234.     }
  235.   }
  236.   return x >> offset;
  237. }
  238.  
  239.  
  240.