home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / pascal / library / dos / math / ultra101 / ultra.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-26  |  7.6 KB  |  222 lines
  1. /*
  2.    ULTRA IN C
  3.  
  4.    This is the ULTRA random number generator written entirely
  5.    in C. Efficiency has always been sacrificed for readability.
  6.  
  7.    Hopefully this may serve as a model for an assembler version
  8.    of this routine. The programmer should avoid simply duplicating
  9.    and should use the features usually present in any assembler
  10.    to increase the speed of this routine.
  11.  
  12.    Especially the subroutine SWB should be replaced by the one
  13.    machine instruction (usually called subtract-with-borrow) that
  14.    is available in almost every hardware.
  15.  
  16.    For people not familiar with 8086 assembler, it may help to
  17.    consult this when reading the assembler code. The output of
  18.    this program should exactly match the output of ULTRADEM.C
  19.    linked with ULTRA.C
  20. */
  21.  
  22. #define N 37            /* The number of 32 bit words in the table */
  23.  
  24. unsigned long swbx[N],
  25.     x[N] =
  26.     {2514488732,  431628406, 1236830750, 4003720720,  823319423,   94684520,
  27.      1527338240,  679520694, 2422692153, 3355680443,  728196457,  701184857,
  28.      1109123469, 3203485719, 1941719477, 2470545560, 1105091547, 2688956238,
  29.      3621637922, 1273782189, 4288988815,  764409422, 2813964156, 2010860282,
  30.       683656203, 3648603015,   73133971, 3968408432, 3402270263, 3808241722,
  31.      2695095159, 2490074496, 1477949423, 3995273805,  760366478, 2466020322,
  32.      2727416986};      /* seed array for the SWB generator. */
  33.  
  34. char    bits,      /* Eight bits to be fed out one at a time by i1bit.   */
  35.     nbits=0,  /* Number of fresh bits remaining in bits.            */
  36.     flags=0;  /* The state of the carry flag for the SWB generator. */
  37. int    swbn=0;      /* Number of fresh bytes remaining in swbneg.        */
  38. long unsigned congx = 0x2C6CE807;
  39.  
  40. /* rinit initializes the constants and creates a seed array one bit at
  41.    a time by taking the leading bit of the xor of a shift register
  42.    and a congruential sequence. The same congruential generator continues
  43.    to be used as a mixing generator for the Subtract-with-borrow generator
  44.    to produce the `ultra' random numbers
  45.  
  46.    Since this is called just once, speed doesn't matter much and it might
  47.    be fine to leave this subroutine coded just as it is.
  48. */
  49.  
  50. void rinit(long unsigned congy,long unsigned shrgx)
  51. { short i,j;
  52.   unsigned long  tidbits;
  53.   congx=congy;
  54.   for (i=0;i<N;i++)
  55.   { bits = 0;
  56.     for (j=32;j>0;j--)
  57.     { congx = congx * 69069;
  58.       shrgx = shrgx ^ (shrgx >> 15);
  59.       shrgx = shrgx ^ (shrgx << 17);
  60.       tidbits = (tidbits>>1) | (0x80000000 & (congx^shrgx));
  61.     }
  62.     x[i] = tidbits;
  63.   }
  64.   swbn = 0;
  65.   nbits = 0;
  66.   flags = 0;
  67. }
  68.  
  69. /* SWB is the subtract-with-borrow operation which should be one line
  70.    in assembler code. In order to do a 32 bit s-w-b, we do two 16 bit
  71.    s-w-b's in 32 bit variables and look for the borrow bit in their
  72.    17th bit.
  73.  
  74.    In a proper impelementation, this would subroutine would not exist,
  75.    and would be be replaced by hardware s-w-b operation in the SWBFill
  76.    subroutine.
  77. */
  78.  
  79. long SWB(long x, long y, char *c)
  80. { long hi,lo;
  81.   lo = (x & 0x0000FFFF) - (y & 0x0000FFFF) - *c;
  82.   *c = ((lo&0x00010000) != 0);   
  83.   lo = lo & 0x0000FFFF;
  84.   x  = x>>16;
  85.   y  = y>>16;
  86.   hi = (x & 0x0000FFFF) - (y & 0x0000FFFF) - *c;
  87.   *c = ((hi&0x00010000) != 0);
  88.   return (hi << 16) | lo;
  89. }
  90.  
  91. /* This is the heart of the system and should be written in assembler
  92.    to be as fast as possible. It may be good to unravel the loop and
  93.    simply write 37 consecutive SWB operations!
  94. */
  95.  
  96. void SWBFill(void)
  97. { short i;
  98.   for (i=0;  i<24; i++) x[i] = SWB(x[i+13],x[i],&flags);
  99.   for (i=24; i<N;  i++) x[i] = SWB(x[i-24],x[i],&flags);
  100.   for (i=0;  i<N;  i++) swbx[i] = x[i] ^ (congx = congx * 69069);
  101. }
  102.  
  103. /* The checkfill routine is written out of laziness. In a proper
  104.    implementation this should not ba a subroutine call. Instead
  105.    each routine should do it's own check internally.
  106.  
  107.    The routine checks to see if there are k bytes left in the array.
  108.    If there aren't it refills the array and sets the appropriate
  109.    counts and pointers. Once there are k bytes available, it returns
  110.    a pointer to the first free byte, and decrements the free byte
  111.    count by k. 
  112.  
  113.    By converting the pointer type, this can be used to implement
  114.    all the integer types of random numbers.
  115. */
  116.  
  117. char *swbp;
  118. char *CheckFill(char k)
  119. { if ( swbn-k < 0 )     /* If k bytes not available    */
  120.   { SWBFill();        /*    Refill the bit array    */
  121.     swbn=4*N;        /*    Reset counter        */
  122.     swbp=(char *) swbx; /*    Reset pointer        */
  123.   };
  124.   swbn = swbn - k;    /* decrement the count by k    */
  125.   swbp = swbp + k;    /* increment the pointer by k    */
  126.   return swbp - k;    /* return the old pointer value    */
  127. }
  128.  
  129. /* The following routines should be coded in assembler or
  130.    perhaps as macros in C */
  131.  
  132. long  i32bit(void) { return *(long  *) CheckFill(4); }
  133. long  i31bit(void) { return *(long  *) CheckFill(4) & 0x7FFFFFFF; }
  134. short i16bit(void) { return *(short *) CheckFill(2); }
  135. short i15bit(void) { return *(short *) CheckFill(2) & 0x7FFF; }
  136. char  i8bit(void)  { return           *CheckFill(1); }
  137. char  i7bit(void)  { return           *CheckFill(1) & 0x7F; }
  138.  
  139. /* The uni/vni and the 1-bit routines are perhaps best done in assembler
  140.    because they involve a lot of minute bit/byte twiddling 
  141.  
  142.    Often I have used 24-bit shifts etc to get access to a byte. Obviously
  143.    in a byte addressable architecture this could be done with great ease.
  144.  
  145.    I freely call i32bit and i31bit to create the uni's and vni's. It
  146.    would probably be faster to duplicate the code of those routines,
  147.    and avoid the subroutine call overhead. This is especially true
  148.    because uni and vni are probably the most often called form of any
  149.    random number generator.
  150. */
  151.  
  152. char i1bit(void)
  153. { char temp;
  154.   if (nbits==0)
  155.   { bits = i8bit();
  156.     nbits = 7;
  157.   } else nbits--;
  158.   temp = bits & 1;
  159.   bits = bits >>1;
  160.   return temp;
  161. }
  162.  
  163. #define two2neg31  ( (2.0/0x10000) / 0x10000 )
  164. #define two2neg32  ( (1.0/0x10000) / 0x10000 )
  165.  
  166. float vni(void)
  167. { float scale;
  168.   long temp; unsigned long leadbyte;
  169.   scale = two2neg31;
  170.   temp = i32bit();
  171.   leadbyte =   temp & 0xFF000000;
  172.   while ((leadbyte == 0xFF000000) | (leadbyte == 0))
  173.   {
  174.     /* This following two lines are executed once every 128 calls,
  175.        and so it doesn't need to be crafted too efficiently */
  176.     leadbyte = (unsigned long) i8bit() << 24;
  177.     scale = scale/128;
  178.   }
  179.   temp = (temp & 0x00FFFFFF) | leadbyte;
  180.   return temp * scale;
  181. }
  182.  
  183. float uni(void)
  184. { float scale;
  185.   long temp,leadbyte;
  186.   scale = two2neg31;
  187.   temp = i31bit();
  188.   leadbyte = temp & 0xFF000000;
  189.   while (leadbyte == 0)
  190.   { /* This following two lines are executed once every 128 calls,
  191.        and so it doesn't need to be crafted too efficiently */
  192.     leadbyte = ((long) i7bit()) << 24;
  193.     scale = scale/128;
  194.   }
  195.   temp = (temp & 0x00FFFFFF) | leadbyte;
  196.   return temp*scale;
  197. }
  198.  
  199. /* The duni and dvni routines are written to take eight consecutive
  200.    bytes and simply float them. No care is taken if the leading bits
  201.    are insignificant becuase with such a high precision, a few leading
  202.    bits shouldn't matter.
  203.  
  204.    It is worth exploring whether it is faster to float a number using
  205.    a floating point processor, or just by counting the number of leading
  206.    bits that are zero, and inserting the exponent by hand (i.e. with
  207.    the main cpu).
  208. */
  209.  
  210. double dvni(void)
  211. { signed long *temp;
  212.   unsigned long lo;
  213.   temp = (signed long *) CheckFill(8);
  214.   lo = (unsigned long) *(temp++);
  215.   return ((*temp) + lo * two2neg32) * two2neg31; }
  216.  
  217. double duni(void)
  218. { unsigned long *temp, lo;
  219.   temp = (unsigned long *) CheckFill(8);
  220.   lo = *(temp++);
  221.   return ( ((*temp)&0x7FFFFFFF) + lo * two2neg32) * two2neg31; }
  222.