home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / gnu / g__lib / rng.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-23  |  2.2 KB  |  97 lines
  1. #include <values.h>
  2. #include <assert.h>
  3. #include <builtin.h>
  4. #include <RNG.h>
  5.  
  6. //
  7. //    The scale constant is 2^-31. It is used to scale a 31 bit
  8. //    long to a double.
  9. //
  10.  
  11. const double randomDoubleScaleConstant = 4.656612873077392578125e-10;
  12. const float  randomFloatScaleConstant = 4.656612873077392578125e-10;    
  13.  
  14. unsigned long RNG::asLong()
  15. {
  16.   (*lib_error_handler)("RNG", "Unimplemented");
  17. }
  18.  
  19. void RNG::reset()
  20. {
  21.   (*lib_error_handler)("RNG", "Unimplemented");
  22. }
  23.  
  24.  
  25.  
  26. static char initialized = 0;
  27.  
  28. RNG::RNG()
  29. {
  30.     if (!initialized) {
  31.     if (sizeof(double) != 2 * sizeof(unsigned long)) {
  32.       (*lib_error_handler)("RNG", "sizeof(Double) != 2 Sizeof(long)\n");
  33.     };
  34.     //
  35.     //    The following is a hack that I attribute to
  36.     //    Andres Nowatzyk at CMU. The intent of the loop
  37.     //    is to form the smallest number 0 <= x < 1.0,
  38.     //    which is then used as a mask for two longwords.
  39.     //    this gives us a fast way way to produce double
  40.     //    precision numbers from longwords.
  41.     //
  42.     //    I know that this works for IEEE and VAX floating
  43.     //    point representations.
  44.     //
  45.     //    A further complication is that gnu C will blow
  46.     //    the following loop, unless compiled with -ffloat-store,
  47.     //    because it uses extended representations for some of
  48.     //    of the comparisons. Thus, we have the following hack.
  49.     //    If we could specify #pragma optimize, we wouldn't need this.
  50.     //
  51.  
  52.     PrivateRNGDoubleType t;
  53.     PrivateRNGSingleType s;
  54.  
  55. #if _IEEE == 1
  56.     
  57.     t.d = 1.5;
  58.     if ( t.u[1] == 0 ) {        // sun word order?
  59.         t.u[0] = 0x3fffffff;
  60.         t.u[1] = 0xffffffff;
  61.     }
  62.     else {
  63.         t.u[0] = 0xffffffff;    // encore word order?
  64.         t.u[1] = 0x3fffffff;
  65.     }
  66.  
  67.     s.u = 0x3fffffff;
  68. #else
  69.     double x = 1.0;
  70.     double y = 0.5;
  71.     do {                // find largest fp-number < 2.0
  72.         t.d = x;
  73.         x += y;
  74.         y *= 0.5;
  75.     } while (x != t.d && x < 2.0);
  76.  
  77.     float xx = 1.0;
  78.     float yy = 0.5;
  79.     do {                // find largest fp-number < 2.0
  80.         s.s = xx;
  81.         xx += yy;
  82.         yy *= 0.5;
  83.     } while (xx != s.s && xx < 2.0);
  84. #endif
  85.     // set doubleMantissa to 1 for each doubleMantissa bit
  86.     doubleMantissa.d = 1.0;
  87.     doubleMantissa.u[0] ^= t.u[0];
  88.     doubleMantissa.u[1] ^= t.u[1];
  89.  
  90.     // set singleMantissa to 1 for each singleMantissa bit
  91.     singleMantissa.s = 1.0;
  92.     singleMantissa.u ^= s.u;
  93.  
  94.     initialized = 1;
  95.     }
  96. }
  97.