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

Text File | 1992-03-26 | 12.0 KB | 432 lines

; 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.
;

cong  macro     ; Congruential generator ==============================
;
; argument in dx:ax
; result   in dx:ax = low32bits of 69069 * arg
; clobbered   bx
;
; We need the lower 32 bits of 69069 * dx:ax
; Since we only have 16 bit multiplys with 32 bit answers,
; we use long multiplication. Let (a,b) represent a*2^16 + b.
;
;            69069 =                      (  1, 3533)
;            dx:ax =                 x    ( hi, lo  )
;                                 -------------------
;                     3533 * lo =         ( a1, a2  )
;                     3533 * hi =     ( b1, b2 )
;                        1 * lo =     (  0, lo )
;                                 -------------------
;            low 32 bits of answer = (a1+b2+lo, a2  )
;
      mov  bx,ax
      mov  ax,3533
      mul  dx                     ; dx:ax = b1,b2 = 3533*hi
      xchg bx,ax
      add  bx,ax                  ; bx = b2+lo;     ax = lo;
      mov  dx,3533
      mul  dx                     ; dx:ax = a1,a2 = 3533*lo
      add  dx,bx                  ; dx = a1+b2+lo
      endm

shreg macro     ; Shift register sequence =============================
;
;  argument is dx:ax
;  result   in dx:ax
;  clobbered   bx
;
;  consider the ax register as one sign bit `as' and the rest `a..r'
;  similarly decompose dx. The operation we want to do is:
;
;          x = x xor (x shr 15);
;          x = x xor (x shl 17);
;
;  Which can be shown in terms of the decomposed registers as:
;
;     x        = ds, d..r, as, a..r
;     x shr 15 =       ds, d..r  as
;                ------------------
;     x        = ds, d..r, as, a..r
;     x shl 17 = a..r   0
;                ------------------
;
;  The above operation is equivalent to:
;
;           ax = ax xor [d..r,as]
;           dx = dx xor [a..r,ds]
;
      mov  bh,ah
      shl  bh,1                 ; mov as into carry
      mov  bx,dx
      rcl  bx,1                 ; bx = [d..r,as]
      pushf                     ; save carry flag (=ds)
      xor  ax,bx                ; ax = ax xor bx
      popf                      ; restore the carry
      mov  bx,ax
      rcl  bx,1                 ; bx = [a..r,ds]
      xor  dx,bx                ; dx = dx xor bx
      endm

RinitProcStart         ; INITIALIZE SEED ARRAY =======================
;
; 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
;   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]
;   bl contains the bits until eight of them are gathered
;   es:di point to where x[i]
;   dx:ax are used for computations
;
    mov ch,N*4
    mov di,offset x     ; do loop for x[0], ... , x[4*n] (bytes)
nextbyte:
    mov bl,0        ; 8 bits of bl are made bit by bit
    mov cl,8
nextbit:    push bx                     ; compute the cong and the sh-reg
        mov  ax,conglo              ; generators in place
        mov  dx,conghi
        cong
        mov  conglo,ax
        mov  conghi,dx
        mov  ax,shrglo
        mov  dx,shrghi
        shreg
        mov  shrglo,ax
        mov  shrghi,dx
        pop  bx

        xor  dh,byte ptr conghi+1   ; xor the top byte of arg1 and arg2
        rcl  dh,1           ; move the sign bit into the carry
        rcr  bl,1           ; shift it into bl
    dec cl
    jnz nextbit
    mov al,bl
    stosb               ; 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
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*2             ; will do first loop 48 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 8086, the instructions `lodsw', `stosw' and `loop'
    ; all change registers automatically as noted in parentheses

    lodsw                   ;    ax = y[i+26]     ( si = si+2 )
    sbb ax,[di]             ;    ax = ax-y[i]-carry
    stosw                   ;    y[i] = ax        ( di = di+2 )
    loop loop1              ; loop for i=0..47    ( cx = cx-1 )

    mov cx,13*2             ; will do next loop 26 times
    mov si,offset(x)        ; set up ds:si -> y[0]
                ; es:di is already set up
loop2:
    lodsw                   ;    ax = y[i-48]     ( si = si+2 )
    sbb ax,[di]             ;    ax = ax-y[i]-carry
    stosw                   ;    x[i] = ax        ( di = di+2 )
    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  ax,conglo
    mov  dx,conghi
loopc:  cong
    mov  bx,ax
    lodsw
    xor  ax,bx
    stosw
    lodsw
    xor  ax,dx
    stosw
    mov  ax,bx
    loop loopc
    mov  conglo,ax
    mov  conghi,dx

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
    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  ax,word ptr swbneg[bx-4]
    mov  shrglo,ax
    mov  al,byte ptr swbneg[bx-2]
    mov  byte ptr shrghi,al          ; and carry along the three lsb's.
ugetleadbyte:
    fisub seven                      ; decrease exponent by seven
    call i8bit                       ; and get 7 new leading bytes.
    and  al,07Fh
    jz   ugetleadbyte
    mov  byte ptr shrghi[1],al
    fild shrgx
    jmp  uscale_n_exit
uni endp

vni proc
    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  ax,word ptr swbneg[bx-4]
    mov  shrglo,ax
    mov  al,byte ptr swbneg[bx-2]
    mov  byte ptr shrghi,al          ; and carry along the three lsb's.
vgetleadbyte:
    fisub seven                      ; decrease exponent by seven
    call i8bit                       ; and get 7 new leading bytes.
    cbw
    cmp  al,ah
    je   vgetleadbyte
    mov  byte ptr shrghi[1],al
    fild 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