home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / System / Ultra / UltraLib.c < prev   
Encoding:
C/C++ Source or Header  |  1992-12-08  |  15.4 KB  |  590 lines  |  [TEXT/KAHL]
  1. /**************************    ULTRA (v1.0)   **************************
  2.  
  3. This is a Macintosh/Motorola implementation of the Ultra pseudo-random
  4. number generator --  the state-of-the-art as of 11/92.  It requires a
  5. 68020/68881 (or better) chip set.
  6.  
  7. References --
  8.  
  9.     Creator:    Prof. George Marsaglia and Co-workers
  10.                 Dept. of Statistics
  11.                 Florida State University
  12.                 
  13.     Code by:    Dr. Michael P. McLaughlin
  14.                 MITRE Corp.
  15.                 MS W337
  16.                 7525 Colshire Dr.
  17.                 McLean, VA  22102
  18.                 
  19.                 E-mail: mpmcl@mitre.org
  20.                 
  21.     Please acknowledge the work of others by leaving this header intact.
  22.     
  23. ************************************************************************
  24.     
  25.             ---------------  General Remarks  ---------------
  26.             
  27. This code was written in a Think C <TM> environment but should be readily
  28. portable.  Assembly language was used to reduce execution time.  However,
  29. to ensure portablilty, C was used in a few places for functions returning
  30. a float or double.  Things to watch out for include the following:
  31.  
  32.     Compiler options for 68020 and 68881 are REQUIRED.
  33.     
  34.     sizeof(double) is compiler-option dependent and is here unspecified.
  35.     Longs and floats are 4 bytes.
  36.     Integers are 2 bytes.
  37.     
  38.     Functions returning bytes or bits are pre-cast to integer.
  39.     
  40.     Registers overwritten --> D0,D1,D2,A0,A1,FP0,FP1,FP2,FP3
  41.     
  42.     WARNING!  Be sure to compile and run UltraDemo in your desired
  43.     environment.  Check that your output is identical to that supplied
  44.     with this package.
  45.     
  46.     Users are cautioned that this code is highly optimized.  Even seemingly
  47.     innocuous changes, such as rearranging the functions below, may cause
  48.     significant performance degradation.
  49.      
  50.                 ---------------  Timings  ---------------
  51.  
  52. The execution times for the functions supplied are somewhat variable and
  53. are increased by factors such as instruction-cache overflow and program
  54. segmentation.  The timings below are averages of 1-10 million calls using
  55. a standard MacIIci (68030/68882) with extensions off and the calling module
  56. and UltraLib in the same segment. ("Native Floating Point" was not used.)
  57. To facilitate comparisons, the calling overhead has been subtracted out. 
  58. Your timings may differ.
  59.  
  60.               Function                    Microseconds per Call     
  61.  
  62.             Ultra_long32                         4.7
  63.             Ultra_long31                         4.9
  64.         
  65.             Ultra_int16                              2.9
  66.             Ultra_int15                              3.1
  67.  
  68.             Ultra_int8                              2.2
  69.             Ultra_int8u                             2.2
  70.             Ultra_int7                              2.3
  71.  
  72.             Ultra_int1                              1.6
  73.  
  74.             Ultra_uni                            14.3
  75.             Ultra_vni                            16.2
  76.  
  77.             Ultra_duni                            31.8
  78.             Ultra_dvni                            29.6
  79.  
  80.             Ultra_norm                            50.5
  81.             Ultra_expo                            39.5
  82.  
  83.             Ultra_SaveStart                        48.0
  84.             Ultra_RecallStart                    49.0
  85.  
  86. ***********************************************************************/
  87.  
  88.  
  89.  
  90.  
  91. #include <math.h>
  92.  
  93. struct {                /* to restart Ultra at a known beginning */
  94.         float          Ugauss;
  95.         unsigned long UFSR[37],USWB[37],Ubits,Useed1,Useed2;
  96.         int              Ubrw,Ubyt,Ubit,Udum;
  97.         char          *Uptr;
  98.     } Ultra_Remember;    /* 324 bytes */
  99.  
  100. double          Ultra_2n63,Ultra_2n31,Ultra_2n7;
  101. float          Ultra_gauss;        /* remaining Ultra_norm */
  102. unsigned long Ultra_FSR[37],    /* feedback-shift-register output */
  103.               Ultra_SWB[37],    /* lagged-Fibonacci output */
  104.               Ultra_bits,        /* bits for Ultra_int1 */
  105.               Ultra_seed1,        /* seeds MUST be initialized with */
  106.               Ultra_seed2;        /*    non-zero values */
  107. int              Ultra_brw,        /* borrow bit */
  108.               Ultra_byt,        /* # bytes left in Ultra_FSR[37] */
  109.               Ultra_bit,        /* # bits left in Ultra_bits */
  110.               Ultra_dum;        /* dummy, to give multiple of 4 */
  111. char          *Ultra_ptr;        /* running pointer to Ultra_FSR */
  112.  
  113. /*    Ultra_Refill() is the core of Ultra (see Marsaglia and Zaman, 1991,
  114.     "A New Class of Random Number Generators", Annals of Applied
  115.     Probability, vol. 1(3), 426-480). It refills Ultra_SWB[37] and uses
  116.     a feedback-shift algorithm to overlay a second PRNG, thus creating
  117.     Ultra_FSR[37] which supplies all of the Ultra-random bytes. */
  118.     
  119. void Ultra_Refill()        /* Unrolling the loops is actually slower! */
  120. {
  121. asm    {
  122.         LEA        Ultra_SWB,A1        ; Ultra_SWB[0]
  123.         LEA        52(A1),A0            ; Ultra_SWB[13]
  124.         MOVEQ    #0,D0                ; restore extend bit
  125.         SUB        Ultra_brw,D0
  126.         MOVEQ    #23,D2                ; 24 of these
  127. @cont1    MOVE.L    (A0)+,D0
  128.         MOVE.L    (A1),D1
  129.         SUBX.L    D1,D0
  130.         MOVE.L    D0,(A1)+
  131.         DBRA    D2,@cont1
  132.         LEA        Ultra_SWB,A0
  133.         MOVEQ    #12,D2                ; 13 of these
  134. @cont2    MOVE.L    (A0)+,D0
  135.         MOVE.L    (A1),D1
  136.         SUBX.L    D1,D0                ; "subtract with borrow"
  137.         MOVE.L    D0,(A1)+
  138.         DBRA    D2,@cont2
  139.         MOVEQ    #0,D0
  140.         MOVE.L    D0,D1
  141.         ADDX    D1,D0                ; get borrow bit
  142.         MOVE    D0,Ultra_brw        ;    and save it
  143.         LEA        Ultra_SWB,A0
  144.         LEA        Ultra_FSR,A1
  145.         MOVE.L    A1,Ultra_ptr        ; reinitialize running pointer
  146.         MOVE.L    Ultra_seed1,D0
  147.         MOVE.L    #69069,D1            ; overlay feedback-shift PRNG
  148.         MOVEQ    #36,D2                ; 37 of these
  149. @cont3    MOVE.L    (A0)+,(A1)
  150.         MULU.L    D1,D0
  151.         EOR.L    D0,(A1)+
  152.         DBRA    D2,@cont3
  153.         MOVE.L    D0,Ultra_seed1        ; save global for next time
  154.         MOVE    #148,Ultra_byt        ; 4*37 bytes now available
  155. }
  156. }
  157.  
  158. /*    Ultra_int7() returns a two-byte integer, U[0,127]. */
  159.  
  160. int Ultra_int7()
  161. {
  162. asm {
  163.         SUBQ    #1,Ultra_byt        ; one byte left?
  164.         BMI.S    @cont2                ; branch if not
  165. @cont1    MOVE.L    Ultra_ptr,A0
  166.         MOVEQ    #0,D0                ; clear long for Ultra_uni()
  167.         MOVE.B    (A0),D0
  168.         BCLR.L    #7,D0                ; clear sign bit
  169.         ADDQ.L    #1,Ultra_ptr
  170.         RTS
  171. @cont2    JSR        Ultra_Refill
  172.         SUBQ    #1,Ultra_byt
  173.         BRA        @cont1
  174. }    
  175. }
  176.  
  177. /*    Ultra_long32() returns a four-byte long, U[-2147483648,2147483647]. It may,
  178.     of course, be cast to unsigned long. */
  179.  
  180. long Ultra_long32()
  181. {
  182. asm {
  183.         SUBQ    #4,Ultra_byt        ; four bytes left?
  184.         BMI.S    @cont2                ; branch if not
  185. @cont1    MOVE.L    Ultra_ptr,A0
  186.         MOVE.L    (A0),D0
  187.         ADDQ.L    #4,Ultra_ptr
  188.         RTS
  189. @cont2    JSR        Ultra_Refill
  190.         SUBQ    #4,Ultra_byt
  191.         BRA        @cont1
  192. }
  193. }
  194.  
  195. /*    Ultra_uni() returns a four-byte float, U(0,1), with a mantissa of at least
  196.     25-bit precision.  Ultra_uni() achieves this precision by continually
  197.     examining the leading byte of the mantissa.  If this byte is zero, it is
  198.     replaced with a new random byte and the decimal point readjusted.  This
  199.     procedure also ensures that Ultra_uni() never returns zero. */
  200.  
  201. float Ultra_uni()
  202. {
  203.     float result;        /* inserted for portability */
  204.     
  205. asm    {
  206.         FMOVE.X    Ultra_2n31,FP0        ; modulus
  207.         SUBQ    #4,Ultra_byt        ; modified Ultra_long31
  208.         BMI.S    @cont3
  209. @cont1    MOVE.L    Ultra_ptr,A0
  210.         MOVE.L    (A0),D1
  211.         BCLR.L    #31,D1
  212.         ADDQ.L    #4,Ultra_ptr
  213.         MOVE.L    D1,D2
  214.         AND.L    #0xFF000000,D2        ; test leading byte
  215.         BEQ.S    @cont4                ; branch if zero
  216. @cont2    FMUL.L    D1,FP0                ; normalize
  217.         FMOVE.S    FP0,result
  218. }
  219.  
  220.     return(result);
  221.  
  222. asm {
  223. @cont3    JSR        Ultra_Refill
  224.         SUBQ    #4,Ultra_byt
  225.         BRA        @cont1
  226. @cont4    JSR        Ultra_int7            ; get another leading byte
  227.         FMUL.X    Ultra_2n7,FP0        ;    and re-scale
  228.         TST.B    D0                    ; make sure new byte is non-zero
  229.         BEQ        @cont4                ;    else, try again
  230.         ROR.L    #8,D0                ; low -> high byte    
  231.         OR.L    D0,D1                ; replace high byte of longword
  232.         BRA        @cont2
  233. }
  234. }
  235.  
  236. /*    Ultra_vni() returns a four-byte float, U(-1,1), with the same features as
  237.     described above for Ultra_uni(). */
  238.     
  239. float Ultra_vni()
  240. {
  241.     float result;        /* inserted for portability */
  242.     
  243. asm    {
  244.         FMOVE.X    Ultra_2n31,FP0        ; modulus
  245.         SUBQ    #4,Ultra_byt        ; modified Ultra_long32
  246.         BMI.S    @cont4
  247. @cont1    MOVE.L    Ultra_ptr,A0
  248.         MOVE.L    (A0),D1
  249.         BPL.S    @cont2
  250.         FNEG.X    FP0                    ; else modulus = -modulus
  251.         NEG.L    D1                    ;    and D1 = abs(D1)
  252. @cont2    ADDQ.L    #4,Ultra_ptr        ; continue as in Ultra_uni()
  253.         MOVE.L    D1,D2
  254.         AND.L    #0xFF000000,D2        ; test leading byte
  255.         BEQ.S    @cont5                ; branch if zero
  256. @cont3    FMUL.L    D1,FP0                ; normalize
  257.         FMOVE.S    FP0,result
  258. }
  259.  
  260.     return(result);
  261.  
  262. asm {
  263. @cont4    JSR        Ultra_Refill
  264.         SUBQ    #4,Ultra_byt
  265.         BRA        @cont1
  266. @cont5    JSR        Ultra_int7            ; get another leading byte
  267.         FMUL.X    Ultra_2n7,FP0        ;    and re-scale
  268.         TST.B    D0                    ; make sure new byte is non-zero
  269.         BEQ        @cont5                ;    else, try again
  270.         ROR.L    #8,D0                ; low -> high byte    
  271.         OR.L    D0,D1                ; replace high byte of longword
  272.         BRA        @cont3
  273. }
  274. }
  275.  
  276. /*    Ultra_norm() returns a four-byte float, N(mu,sigma).  The normal variates
  277.     returned are exact, not approximate.  Ultra_norm() uses Ultra_vni() so
  278.     there is no possibility of a result exactly equal to mu.  Note that mu and
  279.     sigma must also be floats, not doubles. */
  280.     
  281. float Ultra_norm(float mu, float sigma)
  282. {
  283.     float result;        /* inserted for portability */
  284.     
  285. asm {
  286.         MOVE.L    Ultra_gauss,result    ; is there one left?
  287.         BEQ.S    @cont1
  288.         CLR.L    Ultra_gauss
  289. }
  290.  
  291.     return(sigma*result + mu);
  292.     
  293. asm {
  294. @cont1    BSR.S    @vni                ; "Polar" method
  295.         FMOVE.X    FP0,FP1
  296.         BSR.S    @vni
  297.         FMOVE.X    FP0,FP2
  298.         FMOVE.X    FP1,FP3
  299.         FMUL.X    FP2,FP2
  300.         FMUL.X    FP3,FP3
  301.         FADD.X    FP3,FP2
  302.         FCMP.X    #1,FP2
  303.         FBGE    @cont1                ; Prob(branch) = 1 - π/4
  304.         FMOVE.X    FP2,FP3
  305.         FLOGN.X    FP2
  306.         FADD.X    FP2,FP2
  307.         FNEG.X    FP2
  308.         FDIV.X    FP3,FP2
  309.         FSQRT    FP2
  310.         FMUL.X    FP2,FP0                ; --> standard normals
  311.         FMUL.X    FP2,FP1
  312.         FMOVE.S    FP0,Ultra_gauss        ; save the second one
  313.         FMUL.S    sigma,FP1
  314.         FADD.S    mu,FP1
  315.         FMOVE.S    FP1,result            ;    and return the first
  316. }
  317.  
  318.     return(result);
  319.     
  320. asm {
  321. @vni    FMOVE.X    Ultra_2n31,FP0        ; modified Ultra_vni
  322.         SUBQ    #4,Ultra_byt        ; modified Ultra_long32
  323.         BMI.S    @vni4
  324. @vni1    MOVE.L    Ultra_ptr,A0
  325.         MOVE.L    (A0),D1
  326.         BPL.S    @vni2
  327.         FNEG.X    FP0                    ; else modulus = -modulus
  328.         NEG.L    D1                    ;    and D1 = abs(D1)
  329. @vni2    ADDQ.L    #4,Ultra_ptr        ; continue as in Ultra_uni()
  330.         MOVE.L    D1,D2
  331.         AND.L    #0xFF000000,D2        ; test leading byte
  332.         BEQ.S    @vni5                ; branch if zero
  333. @vni3    FMUL.L    D1,FP0                ; normalize
  334.         RTS
  335. @vni4    JSR        Ultra_Refill
  336.         SUBQ    #4,Ultra_byt
  337.         BRA        @vni1
  338. @vni5    JSR        Ultra_int7            ; get another leading byte
  339.         FMUL.X    Ultra_2n7,FP0        ;    and re-scale
  340.         TST.B    D0                    ; make sure new byte is non-zero
  341.         BEQ        @vni5                ;    else, try again
  342.         ROR.L    #8,D0                ; low -> high byte    
  343.         OR.L    D0,D1                ; replace high byte of longword
  344.         BRA        @vni3
  345. }
  346. }
  347.  
  348. /*    Ultra_int1() returns a two-byte integer, U[0,1].  It uses Ultra_long32()
  349.     and returns the bits one at a time.  Unlike most PRNGs, all the bits are
  350.     equally random. */
  351.  
  352. int Ultra_int1()
  353. {
  354. asm {
  355.         MOVE    Ultra_bit,D1        ; one bit left?
  356.         BPL.S    @cont1                ; faster than BMI.S!?
  357.         JSR        Ultra_long32        ;    else, get 32 more
  358.         MOVE.L    D0,Ultra_bits
  359.         MOVEQ    #31,D1
  360.         MOVE    D1,Ultra_bit        ; bit number to test
  361. @cont1    MOVE.L    Ultra_bits,D2        ; to see all 32 bits with BTST
  362.         MOVEQ    #0,D0
  363.         BTST.L    D1,D2                ; examine bit #D1
  364.         BEQ.S    @cont2
  365.         MOVEQ    #1,D0
  366. @cont2    SUBQ    #1,Ultra_bit
  367.         RTS
  368. }    
  369. }
  370.  
  371. /*    Ultra_long31() returns a four-byte integer, U[0,2147483647]. */
  372.  
  373. long Ultra_long31()
  374. {
  375. asm {
  376.         SUBQ    #4,Ultra_byt        ; four bytes left?
  377.         BMI.S    @cont2                ; branch if not
  378. @cont1    MOVE.L    Ultra_ptr,A0
  379.         MOVE.L    (A0),D0
  380.         BCLR.L    #31,D0                ; clear sign bit
  381.         ADDQ.L    #4,Ultra_ptr
  382.         RTS
  383. @cont2    JSR        Ultra_Refill
  384.         SUBQ    #4,Ultra_byt
  385.         BRA        @cont1
  386. }    
  387. }
  388.  
  389. /*    Ultra_int16() returns a two-byte integer, U[-32768,32767]. */
  390.  
  391. int Ultra_int16()
  392. {
  393. asm {
  394.         SUBQ    #2,Ultra_byt        ; two bytes left?
  395.         BMI.S    @cont2                ; branch if not
  396. @cont1    MOVE.L    Ultra_ptr,A0
  397.         MOVE    (A0),D0
  398.         ADDQ.L    #2,Ultra_ptr
  399.         RTS
  400. @cont2    JSR        Ultra_Refill
  401.         SUBQ    #2,Ultra_byt
  402.         BRA        @cont1
  403. }    
  404. }
  405.  
  406. /*    Ultra_int15() returns a two-byte integer, U[0,32767]. */
  407.  
  408. int Ultra_int15()
  409. {
  410. asm {
  411.         SUBQ    #2,Ultra_byt        ; two bytes left?
  412.         BMI.S    @cont2                ; branch if not
  413. @cont1    MOVE.L    Ultra_ptr,A0
  414.         MOVE    (A0),D0
  415.         BCLR.L    #15,D0                ; clear sign bit
  416.         ADDQ.L    #2,Ultra_ptr
  417.         RTS
  418. @cont2    JSR        Ultra_Refill
  419.         SUBQ    #2,Ultra_byt
  420.         BRA        @cont1
  421. }    
  422. }
  423.  
  424. /*    Ultra_int8() returns a two-byte integer, U[-128,127].  It gets a random
  425.     byte and casts it to integer.  This operation extends the sign bit.
  426.     Consequently, you may NOT cast this function to unsigned integer or long
  427.     (see Ultra_int8u() below). */
  428.  
  429. int Ultra_int8()
  430. {
  431. asm {
  432.         SUBQ    #1,Ultra_byt        ; one byte left?
  433.         BMI.S    @cont2                ; branch if not
  434. @cont1    MOVE.L    Ultra_ptr,A0
  435.         MOVE.B    (A0),D0
  436.         EXT.W    D0                    ; make D0 an int
  437.         ADDQ.L    #1,Ultra_ptr
  438.         RTS
  439. @cont2    JSR        Ultra_Refill
  440.         SUBQ    #1,Ultra_byt
  441.         BRA        @cont1
  442. }    
  443. }
  444.  
  445. /*    Ultra_int8u() returns a two-byte integer, U[0,255].  It proceeds as in
  446.     Ultra_int8() but clears the high byte instead of extending the sign bit. */
  447.  
  448. int Ultra_int8u()        /* unsigned byte cast to int */
  449. {
  450. asm {
  451.         SUBQ    #1,Ultra_byt        ; one byte left?
  452.         BMI.S    @cont2                ; branch if not
  453. @cont1    MOVE.L    Ultra_ptr,A0
  454.         MOVEQ    #0,D0                ; clear high byte
  455.         MOVE.B    (A0),D0
  456.         ADDQ.L    #1,Ultra_ptr
  457.         RTS
  458. @cont2    JSR        Ultra_Refill
  459.         SUBQ    #1,Ultra_byt
  460.         BRA        @cont1
  461. }    
  462. }
  463.  
  464. /*    Ultra_duni() and Ultra_dvni return double-precision U[0,1) and U(-1,1).  In
  465.     both cases, zero IS a possibility.  These functions are intended for those
  466.     occasions when seven significant figures are not enough.  If you need type
  467.     double, but not 63 bits of precision, then it is MUCH faster to use Ultra_uni()
  468.     or Ultra_vni() and cast -- implicitly or explicitly.  Because C compilers
  469.     usually support several different double formats, each with its own set-up
  470.     requirements, these functions were written in C. */
  471.  
  472. double Ultra_duni()
  473. {
  474.     return(Ultra_long31()*Ultra_2n31 + ((unsigned long) Ultra_long32())*Ultra_2n63);
  475. }
  476.  
  477. double Ultra_dvni()
  478. {
  479.     return(Ultra_long32()*Ultra_2n31 + ((unsigned long) Ultra_long32())*Ultra_2n63);
  480. }
  481.  
  482. /*    Ultra_expo() returns a four-byte float, Exponential(mu).  The parameter, mu,
  483.     is both the mean and variance of this distribution.  It must be a float
  484.     greater than zero. */
  485.  
  486. float Ultra_expo(float mu)
  487. {
  488.     float result;        /* inserted for portability */
  489.     
  490. asm {
  491.         FMOVE.X    Ultra_2n31,FP0        ; modified Ultra_uni
  492.         SUBQ    #4,Ultra_byt        ; modified Ultra_long31
  493.         BMI.S    @cont3
  494. @cont1    MOVE.L    Ultra_ptr,A0
  495.         MOVE.L    (A0),D1
  496.         BCLR.L    #31,D1
  497.         ADDQ.L    #4,Ultra_ptr
  498.         MOVE.L    D1,D2
  499.         AND.L    #0xFF000000,D2        ; test leading byte
  500.         BEQ.S    @cont4                ; branch if zero
  501. @cont2    FMUL.L    D1,FP0                ; normalize
  502.         FLOGN.X    FP0
  503.         FNEG.X    FP0
  504.         FMUL.S    mu,FP0
  505.         FMOVE.S    FP0,result
  506. }
  507.  
  508.     return(result);
  509.     
  510. asm {
  511. @cont3    JSR        Ultra_Refill
  512.         SUBQ    #4,Ultra_byt
  513.         BRA        @cont1
  514. @cont4    JSR        Ultra_int7            ; get another leading byte
  515.         FMUL.X    Ultra_2n7,FP0        ;    and re-scale
  516.         TST.B    D0                    ; make sure new byte is non-zero
  517.         BEQ        @cont4                ;    else, try again
  518.         ROR.L    #8,D0                ; low -> high byte    
  519.         OR.L    D0,D1                ; replace high byte of longword
  520.         BRA        @cont2
  521. }
  522. }
  523.  
  524. /*    Ultra_SaveStart() and Ultra_RecallStart() save and restore, respectively,
  525.     the global variables of Ultra, from Ultra_gauss through Ultra_ptr.  If you
  526.     think it may be necessary to recall a sequence of random numbers exactly,
  527.     first call Ultra_SaveStart().  To recover the sequence later, call
  528.     Ultra_RecallStart().  If you wish to terminate the program and still recover
  529.     the same random numbers, you must save the array Ultra_Remember[] to a file
  530.     and read it back upon restart.  Ultra_Remember[] is declared in the header
  531.     file as an array of unsigned longs even though only some of its variables
  532.     are actually of this type. */
  533.  
  534. void Ultra_SaveStart()
  535. {
  536. asm {
  537.         LEA        Ultra_gauss,A0
  538.         LEA        Ultra_Remember,A1
  539.         MOVEQ    #80,D0                ; 81 longs
  540. @cont    MOVE.L    (A0)+,(A1)+
  541.         DBRA    D0,@cont    
  542. }
  543. }
  544.  
  545. void Ultra_RecallStart()
  546. {
  547. asm {
  548.         LEA        Ultra_gauss,A0
  549.         LEA        Ultra_Remember,A1
  550.         MOVEQ    #80,D0                ; 81 longs
  551. @cont    MOVE.L    (A1)+,(A0)+
  552.         DBRA    D0,@cont    
  553. }
  554. }
  555.  
  556. /*    Ultra_Init() computes a few global constants, initializes others and fills
  557.     in the initial Ultra_SWB array using the user-supplied seeds which must be
  558.     non-zero unsigned longs.  It terminates by calling Ultra_SaveStart() so that
  559.     you may recover the whole sequence of random numbers by calling Ultra_Init(),
  560.     with the same seeds, then calling Ultra_RecallStart(). */
  561.  
  562. void Ultra_Init()
  563. {
  564.     unsigned long tempbits,ul;
  565.     int i,j;
  566.     
  567.     if ((Ultra_seed1 == 0) || (Ultra_seed2 == 0)) {
  568.         puts("\nYou forgot to supply seeds for Ultra!!");
  569.         exit(1);
  570.     }
  571.      for (i=0;i<37;i++) {
  572.         tempbits = 0;
  573.         for (j=32;j>0;j--) {
  574.             Ultra_seed1 *= 69069;
  575.             Ultra_seed2 ^= (Ultra_seed2 >> 15);
  576.             Ultra_seed2 ^= (Ultra_seed2 << 17);
  577.             ul = Ultra_seed1 ^ Ultra_seed2;
  578.             tempbits = (tempbits >> 1) | (0x80000000 & ul);
  579.         }
  580.         Ultra_SWB[i] = tempbits;
  581.     }
  582.     Ultra_2n31 = ((2.0/65536)/65536);
  583.     Ultra_2n63 = 0.5*Ultra_2n31*Ultra_2n31;
  584.     Ultra_2n7 = 1.0/128;
  585.     Ultra_gauss = 0.0;
  586.     Ultra_byt = Ultra_brw = 0;
  587.     Ultra_bit = -1;
  588.     Ultra_SaveStart();
  589. }
  590.