home *** CD-ROM | disk | FTP | other *** search
/ C!T ROM 2 / ctrom_ii_b.zip / ctrom_ii_b / PROGRAM / PASCAL / TPL60N19 / ARISOURC / F48FSRT.ASM < prev    next >
Assembly Source File  |  1993-01-24  |  10KB  |  200 lines
  1.  
  2. ; *******************************************************
  3. ; *                                                     *
  4. ; *     Turbo Pascal Runtime Library Version 6.0        *
  5. ; *     Real Square Root                                *
  6. ; *                                                     *
  7. ; *     Copyright (C) 1989-1992 Norbert Juffa           *
  8. ; *                                                     *
  9. ; *******************************************************
  10.  
  11.              TITLE   F48FSRT
  12.  
  13.              INCLUDE SE.ASM
  14.  
  15.  
  16. CODE         SEGMENT BYTE PUBLIC
  17.  
  18.              ASSUME  CS:CODE
  19.  
  20. ; Externals
  21.  
  22.              EXTRN   HaltError:NEAR
  23.  
  24. ; Publics
  25.  
  26.              PUBLIC  RSqrt
  27.  
  28. ;-------------------------------------------------------------------------------
  29. ; RSqrt computes the square root of it argument. For a negative argument,
  30. ; runtime error 207 is invoked. The routine uses a estimate-and-correct method
  31. ; similar to that used in the division routine. This algorithm is based on the
  32. ; longhand method taught in school. Since there can be no tie case in rounding
  33. ; when computing the square root, no guard and sticky flags are needed to round
  34. ; to nearest or even in compliance with the IEEE floating point standard.
  35. ;
  36. ; INPUT:     DX:BX:AX  argument
  37. ;
  38. ; OUTPUT:    DX:BX:AX  square root of argument
  39. ;
  40. ; DESTROYS:  AX,BX,CX,DX,SI,DI,Flags
  41. ;-------------------------------------------------------------------------------
  42.  
  43. ROOT_EXT     PROC    FAR
  44. $zero_sqrt:  RET                       ; return argument
  45. $sqrt_err:   MOV     AX, 0CFh          ; load error code
  46.              JMP     HaltError         ; execute error handler
  47. ROOT_EXT     ENDP
  48.  
  49.              ALIGN   4
  50.  
  51. RSqrt        PROC    FAR
  52.              OR      AL, AL            ; argument zero ?
  53.              JZ      $zero_sqrt        ; yes, return zero
  54.              OR      DH, DH            ; argument negative ?
  55.              JS      $sqrt_err         ; yes, error
  56.              PUSH    BP                ; save TURBO-Pascal frame pointer
  57.              SHR     AL, 1             ; divide biased exponent by 2
  58.              SBB     CL, CL            ; CL = 0, if expo even, else CL = -1
  59.              ADC     AL, 40h           ; make new biased exponent
  60.              PUSH    AX                ; save new exponent
  61.              XOR     AL, AL            ; clear LSB of a3
  62.              OR      DH, 80h           ; set hidden bit
  63.              NEG     CL                ; CL = 0, if expo even, else CL = 1
  64.              SHR     DX, CL            ; divide argument
  65.              RCR     BX, CL            ;  by 2 if
  66.              RCR     AX, CL            ;   expo odd
  67.              XCHG    AX, CX            ; argument in DX:BX:CX
  68.              MOV     SI, DX            ; save a1
  69.              MOV     DI, BX            ; save a2
  70.              MOV     BP, CX            ; save a3
  71.              MOV     BH, DH            ; get MSB of a1
  72.              STC                       ; generate 4 bit approximation
  73.              RCR     BH, 1             ;  in hi-nibble of BH
  74.              NEG     AL                ; adjust approximation (subtract 10)
  75.              AND     AL, 10            ;  for arguments
  76.              SUB     BH, AL            ;   between 4000 0000 00 and 7FFF FFFF FF
  77.              MOV     AL, 0FFh          ; default quotient = FFh
  78.              CMP     DH, BH            ; division overflow ?
  79.              JNB     $no_div0          ; yes, use default quotient
  80.              MOV     AX, DX            ; get a1
  81.              DIV     BH                ; compute a1 / approximation
  82. $no_div0:    ADD     BH, AL            ; add result to approximation
  83.              RCR     BH, 1             ;  and divide by 2 yields 8 bit approx.
  84.              MOV     BL, 0FFh          ; BX >= Trunc(Sqrt(a1))
  85.              MOV     AX, 0FFFFh        ; default quotient
  86.              CMP     DX, BX            ; division overflow ?
  87.              JNB     $no_div1          ; yes, use default quotient
  88.              MOV     AX, DI            ; get a2
  89.              DIV     BX                ; compute a1:a2 / approximation
  90. $no_div1:    ADD     AX, BX            ; add result to approximation
  91.              RCR     AX, 1             ;  and divide by 2 yields 16 bit approx.
  92.              MOV     BX, AX            ; save q1
  93.              MUL     AX                ; DX:AX = q1*q1
  94.              SUB     DI, AX            ; compute
  95.              SBB     SI, DX            ;  remainder (SI:DI)
  96.              JNC     $quot1_ok         ; remainder must be positive
  97.  
  98.              ALIGN   4
  99.  
  100. $corr_quot1: DEC     BX                ; try next smaller quotient q1
  101.              STC                       ; adjust
  102.              ADC     DI, BX            ;  remainder
  103.              ADC     SI, 0             ;   to new
  104.              ADD     DI, BX            ;     value
  105.              ADC     SI, 0             ;      of quotient
  106.              JS      $corr_quot1       ; until remainder becomes positive
  107. $quot1_ok:   XCHG    AX, CX            ; concat rest of argument to remainder
  108.              MOV     DX, DI            ; get remainder lo-word
  109.              SHR     SI, 1             ; divide remainder (from here SI=0 !!!)
  110.              RCR     DX, 1             ;  by two so it fits into 32 bits
  111.              RCR     AX, 1             ; DX:AX = remainder / 2
  112.              MOV     CX, 0FFFFh        ; load default quotient
  113.              CMP     DX, BX            ; division overflow ?
  114.              JNB     $no_div2          ; yes, skip division and use default quot
  115.              DIV     BX                ; AX = q2
  116.              XCHG    AX, CX            ; BX:CX = q1:q2
  117. $no_div2:    MOV     AX, CX            ; get q2
  118.              MUL     BX                ; q1*q2
  119.              SUB     BP, AX            ; subtract product
  120.              SBB     DI, DX            ;  from
  121.              SUB     BP, AX            ;   old
  122.              SBB     DI, DX            ;    remainder
  123.              MOV     AX, CX            ; get q2
  124.              MUL     AX                ; q2*q2
  125.              NEG     AX                ; compute
  126.              SBB     BP, DX            ;   new
  127.              SBB     DI, SI            ;    remainder (DI:BP:AX)
  128.              JNS     $quot2_ok         ; remainder must be positive
  129. $corr_quot2: DEC     CX                ; try next smaller value for q2
  130.              STC                       ; adjust
  131.              ADC     AX, CX            ;  value
  132.              ADC     BP, BX            ;   of remainder
  133.              ADC     DI, SI            ;    according
  134.              ADD     AX, CX            ;     to changed
  135.              ADC     BP, BX            ;      value
  136.              ADC     DI, SI            ;       of q2
  137.              JS      $corr_quot2       ; until remainder positive
  138. $quot2_ok:   MOV     SI, AX            ; DI:BP:SI = remainder, BX:CX = quot
  139.              MOV     DX, BP            ; divide
  140.              SHR     DI, 1             ;  remainder by two and
  141.              RCR     DX, 1             ;   stuff 32 most significant bits
  142.              RCR     AX, 1             ;    into DX:AX
  143.              MOV     DI, 0FFFFh        ; load default quotient q3
  144.              CMP     DX, BX            ; division overflow ?
  145.              JNB     $no_div3          ; yes, use default quotient and skip div.
  146.              DIV     BX                ; AX = q3
  147.              XCHG    AX, DI            ; BX:CX:DI = q1:q2:q3, BP:SI:0:0 = rem
  148. $no_div3:    MOV     AX, DI            ; get q3
  149.              MUL     BX                ; q1*q3
  150.              SUB     SI, AX            ; subtract
  151.              SBB     BP, DX            ;  2*q1*q3
  152.              SUB     SI, AX            ;   from
  153.              SBB     BP, DX            ;    remainder
  154.              MOV     AX, DI            ; get q3
  155.              MUL     CX                ; q2*q3
  156.              PUSH    BX                ; save q1
  157.              XOR     BX, BX            ; BP:SI:BX:0 = remainder
  158.              SUB     BX, AX            ; subtract
  159.              SBB     SI, DX            ;  2*q2*q3
  160.              SBB     BP, 0             ;   from
  161.              SUB     BX, AX            ;    remainder
  162.              SBB     SI, DX            ;     "
  163.              SBB     BP, 0             ;      "
  164.              MOV     AX, DI            ; get q3
  165.              MUL     AX                ; q3*q3
  166.              NEG     AX                ; subtract
  167.              SBB     BX, DX            ;  q3*q3
  168.              SBB     SI, 0             ;   from remainder
  169.              SBB     BP, 0             ; BP:SI:BX:AX = remainder
  170.              POP     DX                ; DX:CX:DI = q1:q2:q3
  171.              JNS     $quot3_ok         ; remainder must be positive
  172. $corr_quot3: DEC     DI                ; try next smaller value for q3
  173.              STC                       ; adjust
  174.              ADC     AX, DI            ;  value
  175.              ADC     BX, CX            ;   of
  176.              ADC     SI, DX            ;    remainder
  177.              ADC     BP, 0             ;     according
  178.              ADD     AX, DI            ;      to
  179.              ADC     BX, CX            ;       changed
  180.              ADC     SI, DX            ;        value of
  181.              ADC     BP, 0             ;         quotient
  182.              JS      $corr_quot3       ; until remainder positive
  183. $quot3_ok:   XCHG    AX, DI            ; get result
  184.              MOV     BX, CX            ;  into DX:BX:AX
  185.              ADD     AX, 80h           ; round result
  186.              ADC     BX, 0             ;  (no tie case and
  187.              ADC     DX, 0             ;   no mantissa overflow possible)
  188.              POP     CX                ; get saved exponent
  189.              MOV     AL, CL            ; stuff exponent into result
  190.              AND     DH, 7Fh           ; clear sign bit
  191.              POP     BP                ; restore TURBO-Pascal frame pointer
  192.              RET                       ; done
  193. RSqrt        ENDP
  194.  
  195.              ALIGN   4
  196.  
  197. CODE         ENDS
  198.  
  199.              END
  200.