Power-Programmierung/CD1.mdf/pascal/library/dos/math/ultra101/ult32cod.inc

Text File | 1992-03-26 | 9.6 KB | 363 lines

.386

; This file is included in after a header file which defines
; language dependent elements. (see ultra.doc for details)
;
; Specifically it expects the following items to be defined:
;
;       RinitProcStart  FillProcStart   EnterProcedure
;       RinitProcEnd    FillProcEnd     ExitProcedure
;       DwordFn         WordFn          ByteFn
;       RealFn          DoubleFn
;
;    The DATA and STACK segments.
;

RinitProcStart         ; INITIALIZE SEED ARRAY ========================
;
; On entry:
;    eax and ebx contain congruential and shift-register seeds
;
; On exit:
;    The x array is set using the McGill Super-Duper generator,
;    which is a mix of a congruential and shift-register sequence.
;    Flags and counters are reset.
;
; Registers Clobbered:  AX,BX,CX,DX,SI,DI.
;
; Algorithm:
;    Starting from lsm of x[0] to the msb of x[N-1], each bit is
;    formed using the sign bit of the xor of a congruential generator
;    with seed ConX and a shift register sequence with seed ShrX.
;
; Register usage
;   eax contains congruential seed
;   ebx contains the shift-register seed
;   cx contains count for two loops
;      ch counts the outer loop (for each byte of the array x)
;      cl counts for each bit of x[i]
;   es:di point to where x[i] is
;   edx is a temporary register
;
    mov ch,N*4
    mov di,offset x     ; do loop for x[0], ... , x[4*n] (bytes)
nextbyte:
      mov cl,8
nextbit:
    mov  edx,dword ptr 69069
    mul  edx

    mov  edx,ebx
    shr  edx,15
    xor  ebx,edx
    mov  edx,ebx
    shl  edx,17
    xor  ebx,edx

    mov  edx,eax
    xor  edx,ebx
    rcl  edx,1

    rcr  byte ptr [di],1   ; shift it into the answer
    dec  cl
    jnz  nextbit
    inc  di             ; Store the answer in x[i]
    dec ch
    jnz nextbyte
;
; reset all counters, flags etc. and return.
;
    mov swbn,0
    mov nbits,0
    mov flags,0
    mov congx,eax
RinitProcEnd

FillProcStart           ; SUBTRACT-WITH-BORROW ========================
;
;  On Entry:
;     DS should point to the data segment.
;
;  On Exit:
;     The x array contains all new values computed by the
;     subtract-with-carry generator. The CARRY byte contains the
;     state of the carry flag due to the last subtraction.
;
;  Registers Clobbered: AX,BX,CX,DX,SI,DI
;
;  Algorithm:
;     The following subtractions are performed from right to left.
;     The carry propagates as though these were two long numbers,
;     with each x[i] being a `digit' in base 2^32.
;
;        x[12] ...  x[ 0]      x[36] ...  x[13]
;       -x[36] ... -x[24]     -x[23] ... -x[ 0]
;       ---------------------------------------
;        x[36] ...  x[24]      x[23] ...  x[ 0]
;
;     The x's also could be considered as pairs of base 16 digits,
;     so that x[i] is the pair y[2i+1]y[2i]. This allows us to use
;     only 16 bit subtractions with carry, perfectly suited for all
;     80x86 processors. The same idea could be extended for machines
;     with only eight bit, or even only 1 bit arithmetic.
;
    mov ah,byte ptr flags   ; set carry flag to what it was the
    sahf                    ; last time we exited this routine

    mov cx,24               ; will do first loop 24 times
    mov si,offset(x)+13*4   ; set up ds:si -> x[13]
    mov di,offset(x)        ; set up es:di -> x[0]

loop1:  ; On a 80386, the instructions `lodsd', `stosd' and `loop'
    ; all change registers automatically as noted in parentheses

    lodsd                   ;    ax = x[i+13]     ( si = si+4 )
    sbb eax,[di]            ;    ax = ax-x[i]-carry
    stosd                   ;    y[i] = ax        ( di = di+4 )
    loop loop1              ; loop for i=0..47    ( cx = cx-1 )

    mov cx,13               ; will do next loop 13 times
    mov si,offset(x)        ; set up ds:si -> x[0]
                ; es:di is already set up
loop2:
    lodsd                   ;    ax = x[i-24]     ( si = si+4 )
    sbb eax,[di]            ;    ax = ax-x[i]-carry
    stosd                   ;    x[i] = ax        ( di = di+4 )
    loop loop2              ; loop for i=48..73   ( cx = cx-1 )

    lahf
    mov byte ptr flags,ah   ; save carry flag for next time
;
; XOR the elements of x with a congruential generator and put the
; result in swbx
;
    mov  cx,N
    mov  si,offset(x)
    mov  di,offset(swbx)
    mov  eax,congx
    mov  ebx,69069
loopc:  mul  ebx
    mov  edx,eax
    lodsd
    xor  eax,edx
    stosd
    mov  eax,edx
    loop loopc
    mov  congx,eax
FillProcEnd

; Random Number procedures ============================================
;
; The following procedures all rely on one table SWBX, and a count
; SWBN. The count is negative of the number of bytes remaining.
; The array SWBNEG is defined so that
;   SWBNEG[-4N] is equivalenced to SWBX[0]
;   SWBNEG[ -1] is equivalenced to SWBX[4N-1]
; Thus SWBNEG[SWBN] is the index of the first free byte.
; In FORTRAN, you could think of it as:
;   Dimension SWBX(0:4n-3), SWBNEG(-4n:-1)
;   equivalence(SWBX,SWBNEG)
; This odd setup is so that we can test SWBN against 0, to check
; for exhausting the array, and yet still use the array entries
; in the order that they are stored in the array.
;
; All the routines first increment SWBN by the amount of bytes they
; intend to use. Then they check SWBN. If it has become positive,
; the array is to be refilled, and SWBN set to indicate 4n-nbytes.
;
; If there are enough bytes in the array, they are returned as
; the function value.
;

BitProc MACRO nbits
  local ok
  nbytes=(nbits+1) shr 3
  i&nbits&bit proc
    EnterProcedure
    mov  bx,swbn
  if nbytes eq 1
    inc  bx
  else
    add  bx,nbytes
  endif                         ; -bl is n. of bytes free after this
    jle  ok                 ; if bx>0, nbytes bytes are not remaining
    call fillswb            ; so refill table
    mov  bx,-(N*4-nbytes)   ; and set swbn to  -(N*4-nbytes)
  ok:   mov  swbn,bx
ENDM

EndProc MACRO nbits
    ExitProcedure
  i&nbits&bit endp
ENDM

BitProc 32
    mov  ax,swbneg[bx-4]
    mov  dx,swbneg[bx-2]
    DwordFn
EndProc 32

BitProc 31
    mov  ax,swbneg[bx-4]
    mov  dx,swbneg[bx-2]
    and  dh,7Fh
    DwordFn
EndProc 31

BitProc 16
    mov  ax,swbneg[bx-2]
    WordFn
EndProc 16

BitProc 15
    mov  ax,swbneg[bx-2]
    and  ah,7Fh
    WordFn
EndProc 15

BitProc 8
    mov  al,byte ptr swbneg[bx-1]
    ByteFn
EndProc 8

BitProc 7
    mov  al,byte ptr swbneg[bx-1]
    and  al,7Fh
    ByteFn
EndProc 7

i1bit proc
    EnterProcedure

    dec  nbits
    jge  ok02

    mov  bx,swbn
    inc  bx
    jle  ok01
    call fillswb
    mov  bx,-(N*4-1)
ok01:   mov  swbn,bx

    mov  ah,byte ptr swbneg[bx-1]
    mov  nbits,7
    jmp  short ok03

ok02:   mov  ah,bits

ok03:   mov  al,1
    and  al,ah
    shr  ah,1
    mov  bits,ah

    and  ax,1
    ByteFn
    ExitProcedure
i1bit endp

uni proc
    local shrgx:dword
    EnterProcedure

    fild scale31
    mov  bx,swbn
    add  bx,4
    jle  uchecklead
    call fillswb                     ; refill the table
    mov  bx,-N*4+4
uchecklead:
    mov  swbn,bx             ; set the counter correctly
    and  byte ptr swbneg[bx-1],07Fh  ; if leading byte is zero
    jz   uprecise                    ; then increase precision
                     ; else
    fild dword ptr swbneg[bx-4]      ; load the 32-bit number
uscale_n_exit:
    fscale                           ; scale the integer by the exponent
    fstp st(1)                       ; dump the exponent from the fpp
    RealFn
    ExitProcedure

uprecise:
    mov  eax,dword ptr swbneg[bx-4]  ; put three lsbytes into shrgx
    mov  shrgx,eax                   ; as a temporary storage
ugetleadbyte:
    fisub seven                      ; decrease exponent by seven
    call i8bit                       ; and get 7 new leading bytes.
    and  al,07Fh
    jz   ugetleadbyte
    mov  byte ptr shrgx[3],al
    fild dword ptr shrgx
    jmp  uscale_n_exit
uni endp

vni proc
    local shrgx:dword
    EnterProcedure

    fild scale31
    mov  bx,swbn
    add  bx,4
    jle  vchecklead
    call fillswb                     ; refill the table
    mov  bx,-N*4+4
vchecklead:
    mov  swbn,bx             ; set the counter correctly
    mov  al,byte ptr swbneg[bx-1]    ; if leading byte is 0 or -1
    cbw
    cmp  ah,al
    je   uprecise                    ; then increase precision
                     ; else
    fild dword ptr swbneg[bx-4]      ; load the 32-bit number
vscale_n_exit:
    fscale                           ; scale the integer by the exponent
    fstp st(1)                       ; dump the exponent from the fpp
    RealFn
    ExitProcedure

vprecise:
    mov  eax,dword ptr swbneg[bx-4]  ; put three lsbytes into shrgx
    mov  shrgx,eax                   ; as a temporary storage
vgetleadbyte:
    fisub seven                      ; decrease exponent by seven
    call i8bit                       ; and get 7 new leading bytes.
    cbw
    cmp  al,ah
    je   vgetleadbyte
    mov  byte ptr shrgx[3],al
    fild dword ptr shrgx
    jmp  vscale_n_exit
vni endp

dUVniProc MACRO nbits,name
  local foundit
  name proc
    EnterProcedure

    fild scale63
    mov  bx,swbn
    add  bx,8
    jle  foundit            ; check table if empty
    call fillswb                ; refill table
    mov  bx,-N*4+8
  foundit:
  if     nbits eq 63
    and  byte ptr swbneg[bx-1],07Fh     ; knock off sign bit for duni
  endif
    fild qword ptr swbneg[bx-8] ; load the 64-bit number
    mov  swbn,bx                ; set the counter correctly
    fscale                      ; scale the integer by the exponent
    fstp st(1)                  ; dump the exponent from the fpp
    DoubleFn
    ExitProcedure
  name endp
ENDM

dUVniProc 63,duni
dUVniProc 64,dvni

swbfill proc
    EnterProcedure
    call fillswb
    ExitProcedure
swbfill endp