home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / realrand.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  5.2 KB  |  124 lines
  1. # Funktionen für Zufallszahlen
  2.  
  3. # Zufallszahlengenerator nach [Knuth: The Art of Computer Programming, Vol. II,
  4. # Seminumerical Algorithms, 3.3.4., Table 1, Line 30], nach C. Haynes:
  5. # X eine 64-Bit-Zahl. Iteration X := (a*X+c) mod m
  6. # mit m=2^64, a=6364136223846793005, c=1.
  7.  
  8. # random_L(randomstate) liefert eine neue Zufallszahl.
  9. # > randomstate: ein Random-State, wird verändert
  10. # < ergebnis: eine 32-Bit-Zufallszahl
  11.   local uint32 random_L (object randomstate);
  12.   local uint32 random_L(randomstate)
  13.     var reg4 object randomstate;
  14.     { var reg5 object seed = # letzte Zahl, ein Simple-Bit-Vektor mit 64 Bits
  15.         The_Random_state(randomstate)->random_state_seed;
  16.       var reg1 uintD* seedMSDptr = (uintD*)(&TheSbvector(seed)->data[0]);
  17.       # Multiplikator a=6364136223846793005 = 0x5851F42D4C957F2D :
  18.       local var uintD multiplier[64/intDsize] =
  19.         { D(0x58,0x51,0xF4,0x2D,_EMA_) D(0x4C,0x95,0x7F,0x2D,_EMA_) } ;
  20.       var uintD product[128/intDsize]; # Produkt
  21.       # multiplizieren:
  22.       mulu_2loop_down(&seedMSDptr[64/intDsize],64/intDsize,
  23.                       &multiplier[64/intDsize],64/intDsize,
  24.                       &product[128/intDsize]
  25.                      );
  26.       # letzte 64 Bits holen:
  27.      {var reg3 uint32 seed_hi = get_32_Dptr(&product[64/intDsize]);
  28.       var reg2 uint32 seed_lo = get_32_Dptr(&product[96/intDsize]);
  29.       seed_lo += 1; if (seed_lo==0) { seed_hi += 1; } # um 1 erhöhen
  30.       # seed neu füllen:
  31.       set_32_Dptr(seedMSDptr,seed_hi); set_32_Dptr(&seedMSDptr[32/intDsize],seed_lo);
  32.       # mittlere 32 Bits als Ergebnis:
  33.       return highlow32(low16(seed_hi),high16(seed_lo));
  34.     }}
  35.  
  36. # random_UDS(randomstate,MSDptr,len) füllt die UDS MSDptr/len/..
  37. # mit len Zufallsdigits.
  38. # > randomstate: ein Random-State, wird verändert
  39. # > MSDptr/len/..: wo die Zufallsdigits abgelegt werden sollen
  40. # > len: gewünschte Anzahl von Zufallsdigits
  41.   local void random_UDS (object randomstate, uintD* MSDptr, uintC len);
  42.   local void random_UDS(randomstate,ptr,len)
  43.     var reg4 object randomstate;
  44.     var reg1 uintD* ptr;
  45.     var reg5 uintC len;
  46.     { var reg3 uintC count;
  47.       dotimesC(count,floor(len,32/intDsize),
  48.         { var reg2 uint32 next = random_L(randomstate); # weitere 32/intDsize Digits besorgen
  49.           set_32_Dptr(ptr,next); ptr += 32/intDsize;
  50.         });
  51.       len = len % (32/intDsize); # Anzahl noch fehlender Digits
  52.       if (len>0)
  53.         { var reg2 uint32 next =  random_L(randomstate); # weitere 32/intDsize Digits besorgen
  54.           set_max32_Dptr(intDsize*len,ptr,next);
  55.         }
  56.     }
  57.  
  58. # I_random_I(randomstate,n) liefert zu einem Integer n>0 ein zufälliges
  59. # Integer x mit 0 <= x < n.
  60. # > randomstate: ein Random-State, wird verändert
  61. # kann GC auslösen
  62.   local object I_random_I (object randomstate, object n);
  63.   local object I_random_I(randomstate,n)
  64.     var reg7 object randomstate;
  65.     var reg2 object n;
  66.     { SAVE_NUM_STACK # num_stack retten
  67.       var reg5 uintD* n_MSDptr;
  68.       var reg3 uintC n_len;
  69.       var reg6 uintD* n_LSDptr;
  70.       I_to_NDS_nocopy(n, n_MSDptr=,n_len=,n_LSDptr=); # Digit sequence >0 zu n
  71.      {var reg4 uintD* MSDptr;
  72.       var reg2 uintC len = n_len + ceiling(16,intDsize); # 16 Bits mehr
  73.       if ((intCsize < 32) && (len < n_len)) { BN_ueberlauf(); }
  74.       # neue UDS mit len Zufallsdigits bilden:
  75.       num_stack_need(len,MSDptr=,_EMA_);
  76.       random_UDS(randomstate,MSDptr,len);
  77.       # und durch n dividieren:
  78.       {var DS q;
  79.        var DS r;
  80.        UDS_divide(MSDptr,len,&MSDptr[(uintP)len], n_MSDptr,n_len,n_LSDptr, &q,&r);
  81.        RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  82.        # Rest in Integer umwandeln:
  83.        return NUDS_to_I(r.MSDptr,r.len);
  84.     }}}
  85.  
  86. # F_random_F(randomstate,n) liefert zu einem Float n>0 ein zufälliges
  87. # Float x mit 0 <= x < n.
  88. # > randomstate: ein Random-State, wird verändert
  89. # kann GC auslösen
  90.   local object F_random_F (object randomstate, object n);
  91.   local object F_random_F(randomstate,n)
  92.     var reg7 object randomstate;
  93.     var reg5 object n;
  94.     {  pushSTACK(n);
  95.      { var reg2 uintL d = F_float_digits(n); # d = (float-digits n) > 0
  96.        # Bilde neue UDS mit d Zufallsbits:
  97.        SAVE_NUM_STACK # num_stack retten
  98.        var reg3 uintL len = ceiling(d,intDsize);
  99.        var reg4 uintD* MSDptr;
  100.        num_stack_need_1(len,MSDptr=,_EMA_);
  101.        random_UDS(randomstate,MSDptr,len); # len (>0) Zufallsdigits
  102.        # von intDsize*ceiling(d/intDsize) auf d Bits herunterschneiden:
  103.        { var reg9 uintL dr = d % intDsize; if (dr>0) { MSDptr[0] &= (bit(dr)-1); } }
  104.        # in Integer umwandeln:
  105.       {var reg1 object mant = UDS_to_I(MSDptr,len);
  106.        RESTORE_NUM_STACK # num_stack zurück
  107.        # Bilde  Zufalls-Float zwischen 0 und 1
  108.        #        = (scale-float (float Zufalls-Integer,d_Bits n) (- d)) :
  109.        mant = I_F_float_F(mant,STACK_0); # in Float vom Typ von n umwandeln
  110.        pushSTACK(mant);
  111.        {var reg4 object minus_d = L_to_I(-d); # (- d)
  112.         mant = popSTACK();
  113.         mant = F_I_scale_float_F(mant,minus_d);
  114.        }
  115.        # Multipliziere es mit n :
  116.        mant = F_F_mal_F(mant,STACK_0);
  117.        # mant ist ein Zufalls-Float >=0, <=n.
  118.        if (eql(mant,popSTACK())) # mit n vergleichen
  119.          # falls (durch Rundung) mant=n, durch 0 ersetzen:
  120.          { mant = I_F_float_F(Fixnum_0,mant); }
  121.        return mant;
  122.     }}}
  123.  
  124.