home *** CD-ROM | disk | FTP | other *** search
/ PC Media 2 / PC MEDIA CD02.iso / share / prog / realasm1 / ftsqrt.asm < prev    next >
Encoding:
Assembly Source File  |  1993-07-18  |  4.0 KB  |  100 lines
  1. .286
  2. ;================================================
  3. ; Returns square root of real
  4. ;
  5. ;------------------------------------------------
  6. ; Newton-Raphson approximation method.
  7. ;
  8. ; x[i+1] = 0.5( c/x[i] + x[i] )
  9. ;
  10. ; if c = 25 and seed x[i] = 4
  11. ; 1st iteration = 0.5( 25/4     + 4    ) = 5.125
  12. ; 2nd iteration = 0.5( 25/5.125 + 5.125) = 5.0015
  13. ;
  14. ; Convergence ONLY IF suitable seed is chosen.
  15. ;------------------------------------------------
  16. cseg          segment word public 'code'
  17.               assume  cs:cseg,ss:cseg
  18.               assume  ds:cseg,es:cseg
  19.  
  20.               include math.inc
  21.  
  22. ftsqrt        proc    near real:NPR10
  23.               local   value:REAL10, seed:REAL10, square:REAL10
  24.  
  25.               pusha
  26.               mov     si, real                  ;
  27.               and     word ptr [si]+8, 7fffh    ; positive only
  28.               lea     di, seed                  ;
  29.  
  30.               invoke  movx, di, si, 5           ; seed = real
  31.               mov     dx, word ptr [di]+6       ;
  32.               and     dx, 8000h                 ; DX bit15 = explicit 1^
  33.               mov     ax, word ptr [di]+8       ;
  34.               mov     bx, ax                    ;
  35.               and     bx, 8000h                 ; BX = sign
  36.               and     ax, 7fffh                 ; remove sign
  37.               sub     ax, 3fffh                 ; remove bias
  38.  
  39.               .IF (SWORD PTR ax > 0)            ; if exponent > 0
  40.                  .IF (ax & 1)                   ; if exponent is ODD
  41.                     inc     ax                  ;
  42.                     shr     ax, 1               ;
  43.                     invoke  rshx, di, 4         ;
  44.                     or      word ptr [di]+6, dx ;
  45.                  .ELSE                          ; if exponent is EVEN
  46.                     shr     ax, 1               ;
  47.                  .ENDIF
  48.  
  49.               .ELSEIF (SWORD PTR ax < 0)        ; if exponent < 0
  50.                  neg     ax                     ;
  51.                  .IF (ax & 1)                   ; if exponent is ODD
  52.                     inc     ax                  ;
  53.                     shr     ax, 1               ;
  54.                     invoke  rshx, di, 4         ;
  55.                     or      word ptr [di]+6, dx ;
  56.                  .ELSE                          ; if exponent is EVEN
  57.                     shr     ax, 1               ;
  58.                  .ENDIF
  59.                  neg      ax                    ;
  60.  
  61.               .ELSEIF (ax == 0)                 ;
  62.                  invoke   rshx, di, 4           ;
  63.                  or       word ptr [di]+6, dx   ;
  64.  
  65.               .ENDIF                            ;
  66.               add     ax, 3fffh                 ; add bias
  67.               or      ax, bx                    ; add sign
  68.               mov     word ptr [di]+8, ax       ;
  69.               invoke  ftnormal, di              ;
  70. ;------------------------------------------------
  71.               mov     cx, 16                    ; max 16 iterations
  72.  
  73.               .WHILE (cx)
  74.                  invoke  movx, addr value, si, 5  ; value = real
  75.                  invoke  ftdiv, addr value, di    ; value /= seed
  76.                  invoke  ftadd, addr value, di    ; value += seed
  77.                  invoke  ftdivi, addr value, +2   ; value /= 2
  78.                  invoke  movx, di, addr value, 5  ; seed = value
  79.  
  80.                  invoke  movx, addr square, di, 5 ; square = seed
  81.                  invoke  ftpower, addr square, 2  ; square = square**2
  82.                  invoke  ftsub, addr square, si   ; square -= real
  83.  
  84.                  mov     ax, word ptr square[8] ;
  85.                  or      ax, word ptr square[6] ;
  86.                  or      ax, word ptr square[4] ;
  87.                  or      ax, word ptr square[2] ;
  88.                  .BREAK .IF (ax == 0)           ; exit if 48-bit accuracy
  89.                  dec     cx                     ;
  90.               .ENDW
  91.  
  92.               invoke  movx, si, di, 5           ; return result
  93.               popa
  94.  
  95.               ret
  96. ftsqrt        endp
  97.  
  98. cseg          ends
  99.               end
  100.