home *** CD-ROM | disk | FTP | other *** search
/ Garbo / Garbo.cdr / pc / assemblr / tcrnd.arc / RND.ASM < prev    next >
Encoding:
Assembly Source File  |  1988-02-06  |  14.8 KB  |  527 lines

  1. page 60,132
  2.  
  3. ; Copyright (c) 1987, 1988 by Mark Adler  Pasadena, CA
  4. ; This program may be used and freely copied by anyone except that it
  5. ; may not be sold by itself or as part of any package without the
  6. ; permission of the author.
  7.  
  8. ; RND.ASM by Mark Adler  7 Dec 1987
  9. ; updated by Mark Adler  6 Feb 1988  to change randomize() to rndmize()
  10. ;                                    to avoid conflict with Turbo C's
  11. ;                                    randomize() macro.
  12.  
  13. ; Assemble using Microsoft Macro Assembler 5.0 or later.  Specify memory
  14. ; model as follows:
  15. ;
  16. ; C>masm/mx/d__compact__ rnd;
  17. ;
  18. ; where "compact" can be replaced by tiny, small, medium, large, or
  19. ; huge.  Note that two underscores precede "compact", and two
  20. ; underscores follow it.  The file 'tc.asi' must be accessible.  It may
  21. ; be convenient to add the resulting RND.OBJ file to the C library.
  22. ; For example:
  23. ;
  24. ; C>masm/mx/d__compact__ rnd;
  25. ; C>lib cc+rnd;
  26.  
  27.  title Random Number Generators for Turbo C
  28.  
  29.  include tc.asi
  30.  
  31. comment #
  32.  
  33. by Mark Adler --- algorithms shamelessly copied from Knuth's "Art
  34.                   of Computer Programming", vol 2, 2nd ed.
  35.  
  36. This module provides two different random number generators, both
  37. augmented by shuffling.  The two methods are linear conguential and
  38. additive generation.  Both run at approximately the same speed.  The
  39. reason for having two is to allow switching between them to check the
  40. sensitivity of the simulation to the random number generation method.
  41. The routines are:
  42.  
  43. void setseed(s)  - Set the seed to the 32 bit value 's' and initialize
  44. long s;            the tables for the additive generator and shuffling.
  45.                    Also sets the method to linear congruential.
  46.  
  47. long seed()      - Return the current linear congruential seed.
  48.  
  49. long ticks()     - Return the current tick count (IBM PC compatibles).
  50.  
  51. long rndmize()   - Same as setseed(ticks()) and returns the seed used.
  52.  
  53. void rndpick(f)  - Select the random number generator method to be used.
  54. int f;             f=0 picks linear congruential, f=1 picks additive.
  55.  
  56. long lrnd()      - Return a 32 bit random number.
  57.  
  58. int rnd(n)       - Return a random number between 0 and n, inclusive.
  59. int n;             (n is a 16 bit integer greater than 0.)
  60.  
  61. double drnd()    - Return a random floating point number in [0..1].
  62.  
  63. void shuffle(a, n) - Shuffle the (16 bit) integer array 'a' with 'n'
  64. int a[], n;          entries.  (Unrelated to the shuffling used in the
  65.                      random number generation.)
  66.  
  67. The expected use is to call either rndmize() or setseed() to set the
  68. seed and initialize the tables, then call rndpick() to select the
  69. generation method, and then any combination of lrnd()'s, rnd()'s,
  70. drnd()'s, and shuffle()'s.  rndpick() should not be used in the middle
  71. of the simulation to try to get a "more random" sequence.  The same
  72. method should be used throughout.
  73.  
  74. The drnd() routine assumes the presence of an 80X87, but does not check
  75. for it.
  76.  
  77. These routines should be easily customizable for use with other
  78. compilers.
  79.  
  80. #
  81.  
  82.  .8087
  83.  
  84.  
  85. @header
  86.  
  87.  
  88. @dseg
  89.  even
  90. rmul dd 1664525         ;Good multiplier (Knuth Vol. 2, 2nd Ed,
  91.                         ;                 pg. 102, line 26).
  92. scale dw -31            ;Scale to make fraction out of long.
  93. @endd
  94.  
  95.  
  96. @bseg
  97.  even
  98.         ; data for linear conguential method with shuffling.
  99. rseed dw ?,?            ;32 bit random number seed.
  100. rvals dw 128 dup(?)     ;Room for 64 old seeds.
  101.         ; (note: the rvals and radd tables must be adjacent for setup.)
  102.         ; data for additive generator method.
  103. radd dw 110 dup(?)      ;Room for 55 integers for additive sequence.
  104. raj dw ?                ;j and k for additive generator.
  105. rak dw ?
  106.         ; flag for method choice.
  107. rtyp db ?
  108.  
  109. @endb
  110.  
  111.  
  112. @cseg
  113.  
  114. shfval macro
  115.  ;;
  116.  ;; shfval takes the 32 bit random number in DX:AX and uses part of it
  117.  ;; to pick a previous random value from a set of 64, swaps the new one
  118.  ;; with the old one and returns the old one in DX:AX.  This should make
  119.  ;; a random sequence "randomer" (Knuth Vol.2, 2nd Ed., pg 32, algorithm
  120.  ;; B.)  This macro uses the registers AX, BX, and DX.
  121.  ;;
  122.   mov BL,DH
  123.   and BX,0FCh                   ;;Pick off high 6 bits of DX:AX.
  124.   xchg AX,dgroup:rvals[BX]      ;;Exchange new seed with old one.
  125.   xchg DX,dgroup:rvals[BX+2]
  126.  ;;
  127. endm
  128.  
  129.  
  130. lincon macro
  131.  ;;
  132.  ;; lincon uses the linear congruential method to generate a 32 bit
  133.  ;; random number.  The number is left in DX:AX and updates rseed.
  134.  ;; This macro assumes the direction flag is cleared and that ES and DS
  135.  ;; are equal.  It uses the registers AX, BX, CX, DX, SI, and DI.
  136.  ;;
  137.         ;; multiply seed by multiplier modulo 2^32.
  138.   mov SI,offset dgroup:rmul     ;;Point SI to multiplier.
  139.   mov DI,offset dgroup:rseed    ;;Point DI to seed.
  140.   lodsw                 ;;Get low word of multiplier.
  141.   mov BX,AX             ;;Save that.
  142.   mul word ptr [DI+2]   ;;Multiply by high word of seed.
  143.   mov CX,AX             ;;Save low word of that for high word of result.
  144.   lodsw                 ;;Get high word of multiplier.
  145.   mul word ptr [DI]     ;;Multiply by low word of seed.
  146.   add CX,AX             ;;Add low word into high word of result.
  147.   xchg AX,BX            ;;Get low word of multiplier.
  148.   mul word ptr [DI]     ;;Multiply by low word of seed.
  149.   add DX,CX             ;;Add other high word components.
  150.         ;; add constant modulo 2^32.
  151.   add AX,1              ;;Constant = 1 (relatively prime to
  152.   adc DX,0              ;; everything).
  153.         ;; store seed back.
  154.   stosw                 ;;Store low word.
  155.   mov [DI],DX           ;;Store high word.
  156.  ;;
  157. endm
  158.  
  159.  
  160. addgen macro
  161.  local jcyc, kcyc, done
  162.  ;;
  163.  ;; addgen uses the additive generator of Mitchell and Moore (Knuth
  164.  ;; Vol.2, 2nd Ed., pg. 27.) to generate a 32 but random number.  The
  165.  ;; number is left in DX:AX.  This macro uses the registers AX, BX, DX,
  166.  ;; SI, and DI.
  167.  ;;
  168.         ;; get random 32 bit number in DX:AX.
  169.   mov SI,dgroup:raj     ;;Get j.
  170.   mov DI,dgroup:rak     ;;Get k.
  171.   mov BX,offset dgroup:radd     ;;Point to circular list of 55 integers.
  172.   mov AX,[BX+SI]        ;; r = (radd[k] += radd[j]);
  173.   add AX,[BX+DI]
  174.   mov [BX+DI],AX
  175.   inc BX
  176.   inc BX
  177.   mov DX,[BX+SI]
  178.   adc DX,[BX+DI]
  179.   mov [BX+DI],DX
  180.   sub SI,4
  181.   jb jcyc               ;; j = j ? j - 1 : 54;
  182.   sub DI,4
  183.   jb kcyc               ;; k = k ? k - 1 : 54;
  184.   mov dgroup:raj,SI     ;;Update j and k.
  185.   mov dgroup:rak,DI
  186.   jmp short done
  187.         ;; j cycles up to top.
  188.  jcyc:
  189.    sub DI,4             ;;If j bad, then k is OK.
  190.    mov dgroup:raj,54*4  ;;Cycle j.
  191.    mov dgroup:rak,DI    ;;Set k.
  192.    jmp short done
  193.         ;; k cycles up to top.
  194.  kcyc:
  195.    mov dgroup:raj,SI    ;;If k bad, then j is OK.
  196.    mov dgroup:rak,54*4  ;;Cycle k.
  197.         ;; now everything is updated properly.
  198.  done:
  199.  ;;
  200. endm
  201.  
  202.  
  203. choice macro
  204.  local lin,shf
  205.  ;;
  206.  ;; Use either the linear congruential or the additive generator
  207.  ;; depending on the value of rtyp.  Follow with a shuffle.
  208.  ;;
  209.   cmp dgroup:rtyp,0     ;;See if zero.
  210.   je lin                ;;If so, use linear congruential.
  211.    addgen               ;;Else, use additive generator.
  212.    jmp short shf
  213.  lin:
  214.    lincon               ;;Use linear congruential.
  215.  shf:
  216.   shfval                ;;Shuffle with old values.
  217.  ;;
  218. endm
  219.  
  220.  
  221. zerton macro n
  222.  local ismax
  223.  ;;
  224.  ;; zerton takes a random 32 bit number in DX:AX and makes a random
  225.  ;; number in the range 0..n (inclusive).  n may be anything allowed as
  226.  ;; the source of 'mov CX,'.  The result is left in AX and n+1 is left
  227.  ;; in CX.  This macro uses the registers AX, BX, CX, and DX.
  228.  ;;
  229.         ;; get random number in [0..n] (as per Knuth's recommendation).
  230.   ifdifi <n>,<cx>
  231.    mov CX,n             ;;Get n.
  232.   endif
  233.   inc CX                ;;Use n+1 to multiply "fraction" in [0..1).
  234.   jz ismax              ;;If n is 65535, then return high part of seed.
  235.    mov BX,DX            ;;Save high part of seed.
  236.    mul CX               ;;Multiply low part of fraction by n+1.
  237.    xchg AX,BX           ;;Get high part of fraction.
  238.    mov BX,DX            ;;Save high part of product.
  239.    mul CX               ;;Multiply high part of fraction by n+1.
  240.    add AX,BX            ;;Add component from low part of fraction.
  241.    adc DX,0
  242.  ismax:
  243.   xchg AX,DX            ;;Truncate to integer.
  244.  ;;
  245. endm
  246.  
  247.  
  248. @proc shuffle
  249.  a equ [BP+@a]
  250.  n equ [BP+@a+@d]
  251.  ;
  252.  ; void shuffle(a, n)
  253.  ; int a[], n;
  254.  ;
  255.  ; Shuffle n items (algorithm from Knuth).  n must be greater than zero
  256.  ; and 'int a[n]' must fit in the segment.
  257.  ;
  258.   @enter
  259.   push SI
  260.   push DI
  261.         ; set up to use string instructions (for lincon).
  262.   mov AX,DS             ;Get DS.
  263.   mov ES,AX             ;Set ES to same as DS.
  264.   cld                   ;Autoincrement.
  265.         ; do shuffle.
  266.   mov CX,n              ; get n, j = n - 1.
  267.  shflp:
  268.    loop swp             ;Do n-1 swaps.
  269.    jmp shfend
  270.  swp:
  271.    push CX              ;Save j.
  272.    choice               ;Make 32 bit random number.
  273.    pop CX               ;Restore j.
  274.    zerton CX            ;k = random number in 0..j.
  275.    dec CX               ;Restore j.
  276.    cmp AX,CX            ;See if swap needed (k != j).
  277.    je shno              ;If not, skip swap.
  278.     mov SI,AX           ;Put offsets in SI, DI.
  279.     shl SI,1
  280.     mov DI,CX
  281.     shl DI,1
  282.     if __LDATA__
  283.      lds BX,a           ;Point to a[].
  284.     else
  285.      mov BX,a           ;Point to a[].
  286.     endif
  287.     mov AX,[BX+SI]      ;Get a[k].
  288.     xchg AX,[BX+DI]     ;Exchange with a[j].
  289.     mov [BX+SI],AX      ;Set new a[k].
  290.     if __LDATA__
  291.      mov AX,ES          ;Restore DS.
  292.      mov DS,AX
  293.     endif
  294.   shno:
  295.    jmp shflp
  296.         ; done.
  297.  shfend:
  298.   pop DI
  299.   pop SI
  300.   @leave
  301.  ;
  302. @endp shuffle
  303.  
  304.  
  305. @proc lrnd
  306.  ;
  307.  ; long lrnd()
  308.  ;
  309.  ; Return random 32 bit number using linear congruential method.
  310.  ;
  311.   @enter
  312.   push SI
  313.   push DI
  314.         ; set up to use string instructions (for lincon).
  315.   mov AX,DS             ;Set ES to same as DS.
  316.   mov ES,AX
  317.   cld                   ;Autoincrement.
  318.         ; get a new random number in DX:AX.
  319.   choice                ;Make 32 bit random number.
  320.         ; return seed.
  321.   pop DI
  322.   pop SI
  323.   @leave
  324.  ;
  325. @endp lrnd
  326.  
  327.  
  328. @proc rnd
  329.  n equ [BP+@a]
  330.  ;
  331.  ; int rnd(n)
  332.  ; unsigned n;
  333.  ;
  334.  ; Return random number in [0..n] using linear congruential method.
  335.  ;
  336.   @enter
  337.   push SI
  338.   push DI
  339.         ; set up to use string instructions (for lincon).
  340.   mov AX,DS             ;Set ES to same as DS.
  341.   mov ES,AX
  342.   cld                   ;Autoincrement.
  343.         ; get a new random number in DX:AX.
  344.   choice                ;Make 32 bit random number.
  345.         ; convert to random number in 0..n.
  346.   zerton n
  347.         ; return it.
  348.   pop DI
  349.   pop SI
  350.   @leave
  351.  ;
  352. @endp rnd
  353.  
  354.  
  355. @proc drnd
  356.  ;
  357.  ; double drnd()
  358.  ;
  359.  ; Return random number in [0..1].  This routine assumes a numeric
  360.  ; coprocessor is present, but does not check for it.  The fraction
  361.  ; generated has 31 random bits plus the possibility of 1.  The values
  362.  ; of 0 and 1 have half the probability of any other fraction, as it
  363.  ; should be.  The floating point value is left on the numeric processor
  364.  ; stack.  Two levels of stack are used.
  365.  ;
  366.   @enter
  367.   push SI
  368.   push DI
  369.         ; set up to use string instructions (for lincon).
  370.   mov AX,DS             ;Set ES to same as DS.
  371.   mov ES,AX
  372.   cld                   ;Autoincrement.
  373.         ; do NP operation concurrently with random number generation.
  374.   fild dgroup:scale     ;Get scale to make fraction.
  375.         ; get a new random number in DX:AX.
  376.   choice                ;Make 32 bit random number.
  377.         ; make seed into floating point, leave on NP stack.
  378.   push DX               ;Put in memory.
  379.   push AX
  380.   fild dword ptr [BP-8] ;Load seed into stack.
  381.   fabs                  ;Take absolute value.
  382.   fscale                ;Scale to fraction.
  383.   fstp st(1)            ;Pull scale out of stack.
  384.   add SP,4              ;Trash 32 bit value.
  385.         ; done.
  386.   pop DI
  387.   pop SI
  388.   @leave
  389.  ;
  390. @endp drnd
  391.  
  392.  
  393. @proc seed
  394.  ;
  395.  ; long seed()
  396.  ;
  397.  ; Return current seed.
  398.  ;
  399.   @enter
  400.   mov AX,dgroup:rseed   ;Get low word.
  401.   mov DX,dgroup:rseed+2 ;Get high word.
  402.   @leave
  403.  ;
  404. @endp seed
  405.  
  406.  
  407. @proc ticks
  408.  ;
  409.  ; long ticks()
  410.  ;
  411.  ; Get tick count (assumes IBM PC like).
  412.  ;
  413.   @enter
  414.   mov AH,0              ;Read ticks.
  415.   int 1Ah               ;Time of day interrupt.
  416.   mov AX,DX             ;Get low word in AX.
  417.   mov DX,CX             ;Get high word in CX.
  418.   @leave
  419.  ;
  420. @endp ticks
  421.  
  422.  
  423. setup proc near
  424.  ;
  425.  ; Setup radd, raj, and rak for addgen and rvals for shfval.  See Knuth
  426.  ; for why 23 and 54 are the desired values for j and k.
  427.  ; Note: one of the requirements is that not all of the integers in
  428.  ; radd[] can be even.  This is satisfied since lincon returns
  429.  ; alternating even and odd numbers.
  430.  ;
  431.   push SI
  432.   push DI
  433.         ; set up to use string instructions (for lincon).
  434.   mov AX,DS             ;Set ES to same as DS.
  435.   mov ES,AX
  436.   cld                   ;Autoincrement.
  437.         ; fill rvals[] and radd[] tables.
  438.   mov DI,offset dgroup:rvals
  439.   mov CX,64+55          ;Set rvals[0..63] and radd[0..54] (longs).
  440.  fill:
  441.    push DI              ;Save registers.
  442.    push CX
  443.    lincon               ;Get 32 bit random number in DX:AX.
  444.    pop CX               ;Restore registers.
  445.    pop DI
  446.    stosw                ;Store number.
  447.    xchg AX,DX
  448.    stosw
  449.    loop fill
  450.         ; set j and k.
  451.   mov dgroup:raj,23*4   ; j = 23;
  452.   mov dgroup:rak,54*4   ; k = 54;
  453.         ; set rtyp to use linear congruential.
  454.   mov dgroup:rtyp,0
  455.         ; done.
  456.   pop DI
  457.   pop SI
  458.   ret
  459.  ;
  460. setup endp
  461.  
  462.  
  463. @proc rndmize
  464.  ;
  465.  ; long rndmize()
  466.  ;
  467.  ; Randomize seed and set up for addgen.  This sets the seed to the
  468.  ; current tick count, generates random numbers for addgen and shfval,
  469.  ; sets j and k for addgen, and returns the tick count used.  The
  470.  ; returned tick count could be saved and used as the argument for
  471.  ; setseed() to exactly repeat the pseudo-random sequence.
  472.  ;
  473.   @enter
  474.   call _ticks
  475.   mov dgroup:rseed,AX   ;Set low word.
  476.   mov dgroup:rseed+2,DX ;Set high word.
  477.   push DX               ;Save seed.
  478.   push AX
  479.   call setup            ;Setup tables.
  480.   pop AX                ;Return seed.
  481.   pop DX
  482.   @leave
  483.  ;
  484. @endp rndmize
  485.  
  486.  
  487. @proc setseed
  488.  s equ [BP+@a]
  489.  ;
  490.  ; void setseed(s)
  491.  ; long s;
  492.  ;
  493.  ; Set seed to argument.
  494.  ;
  495.   @enter
  496.   mov AX,s              ;Get low word.
  497.   mov dgroup:rseed,AX   ;Set low word.
  498.   mov AX,s+2            ;Get high word.
  499.   mov dgroup:rseed+2,AX ;Set high word.
  500.   call setup            ;Set up tables.
  501.   @leave
  502.  ;
  503. @endp setseed
  504.  
  505.  
  506. @proc rndpick
  507.  f equ [BP+@a]
  508.  ;
  509.  ; void rndpick(f)
  510.  ; int f;
  511.  ;
  512.  ; Set the method to be used for random number generation.  f=0 selects
  513.  ; linear conguential (the default after rndmize() or setseed() called)
  514.  ; and f=1 selects additive generator.
  515.  ;
  516.   @enter
  517.   mov AL,f
  518.   mov dgroup:rtyp,AL
  519.   @leave
  520.  ;
  521. @endp rndpick
  522.  
  523.  
  524. @endc
  525.  
  526.  end
  527.