home *** CD-ROM | disk | FTP | other *** search
/ Virtual Reality Homebrewer's Handbook / vr.iso / vr386 / sqrti.asm < prev    next >
Assembly Source File  |  1996-03-19  |  7KB  |  441 lines
  1.     TITLE    SQRTI - INTEGER SQRT AND MAGNITUDE
  2.  
  3.     COMMENT $
  4.  
  5. /* Routines for integer sqrt and magnitude for REND386 */
  6.  
  7. // All code by Dave Stampe, last updated 23/12/93
  8. // These routines are (c) 1993 by Dave Stampe
  9.  
  10. /*
  11.  This code is part of the VR-386 project, created by Dave Stampe.
  12.  VR-386 is a desendent of REND386, created by Dave Stampe and
  13.  Bernie Roehl.  Almost all the code has been rewritten by Dave
  14.  Stampre for VR-386.
  15.  
  16.  Copyright (c) 1994 by Dave Stampe:
  17.  May be freely used to write software for release into the public domain
  18.  or for educational use; all commercial endeavours MUST contact Dave Stampe
  19.  (dstampe@psych.toronto.edu) for permission to incorporate any part of
  20.  this software or source code into their products!  Usually there is no
  21.  charge for under 50-100 items for low-cost or shareware products, and terms
  22.  are reasonable.  Any royalties are used for development, so equipment is
  23.  often acceptable payment.
  24.  
  25.  ATTRIBUTION:  If you use any part of this source code or the libraries
  26.  in your projects, you must give attribution to VR-386 and Dave Stampe,
  27.  and any other authors in your documentation, source code, and at startup
  28.  of your program.  Let's keep the freeware ball rolling!
  29.  
  30.  DEVELOPMENT: VR-386 is a effort to develop the process started by
  31.  REND386, improving programmer access by rewriting the code and supplying
  32.  a standard API.  If you write improvements, add new functions rather
  33.  than rewriting current functions.  This will make it possible to
  34.  include you improved code in the next API release.  YOU can help advance
  35.  VR-386.  Comments on the API are welcome.
  36.  
  37.  CONTACT: dstampe@psych.toronto.edu
  38. */
  39.  
  40.         $
  41.  
  42.     .MODEL large
  43.  
  44.     .CODE INTMATH
  45.  
  46.                ; one iteration of 32->16 sqrt
  47. SQRT32  MACRO
  48.     LOCAL skip
  49.     shld    edx,eax,2  ;; get 2 bits of input to error
  50.     shl    eax,2
  51.     add    ebx,ebx    ;; estimate*2
  52.     mov    ecx,ebx    ;; temp = est*2
  53.     add    ecx,ecx
  54.     cmp    edx,ecx    ;; error>2*est?
  55.     jle    skip
  56.     inc    ebx        ;; yes, update for new bit added
  57.     inc    ecx
  58.     sub    edx,ecx
  59. skip:
  60.     ENDM
  61.  
  62.  
  63. ;long squareroot16(long arg)
  64.  
  65. ; takes root of 32 bit number to 16 bit result
  66. ; about 220 clocks worst case:
  67. ; 3 us on 486/66 and 10 us on 386/25
  68.  
  69.  
  70. larg    equ    DWORD PTR [bp+8]
  71.  
  72.     PUBLIC    _squareroot32
  73.  
  74. _squareroot32 proc    far
  75.  
  76.     .386
  77.     push    ebp
  78.     mov    ebp,esp
  79.  
  80.     push    esi
  81.     push    ecx
  82.  
  83.     xor    edx,edx
  84.     xor    ebx,ebx
  85.  
  86.     mov    eax,DWORD PTR larg
  87.     test    eax,0FFFF0000h        ; can we cut it in half?
  88.     jne    hasupper
  89.  
  90.     shl    eax,16        ; yes, so prescale
  91.     test    eax,0FF000000h        ; half again?
  92.     jne    do16
  93.  
  94.     shl    eax,8        ; yes, prescale
  95.     jmp    do8        ; do 8 loops
  96.  
  97. hasupper:
  98.     test    eax,0FF000000h        ; half again?
  99.     jne    do32
  100.     shl    eax,8
  101.     jmp    do24
  102. do32:
  103.     SQRT32
  104.     SQRT32
  105.     SQRT32
  106.     SQRT32
  107. do24:
  108.     SQRT32
  109.     SQRT32
  110.     SQRT32
  111.     SQRT32
  112. do16:
  113.     SQRT32
  114.     SQRT32
  115.     SQRT32
  116.     SQRT32
  117. do8:
  118.     SQRT32
  119.     SQRT32
  120.     SQRT32
  121.     SQRT32
  122.  
  123.     mov    eax,ebx
  124.     shld    edx,eax,16       ; returns in both eax and dx:ax
  125.  
  126.     pop    ecx
  127.     pop    esi
  128.  
  129.     mov    esp,ebp
  130.     pop    ebp
  131.     ret
  132.  
  133. _squareroot32  endp
  134.  
  135.  
  136.                 ; one iteration of 62->31 sqrt
  137. SQRT64  MACRO
  138.     LOCAL skip
  139.     shld    edi,edx,2
  140.     shld    edx,eax,2  ; get 2 bits of input to error
  141.     shl    eax,2
  142.     add    ebx,ebx    ; estimate*2
  143.     mov    ecx,ebx    ; temp = est*2
  144.     add    ecx,ecx
  145.     cmp    edi,ecx    ; error>2*est?
  146.     jle    skip
  147.     inc    ebx        ; yes, update for new bit added
  148.     inc    ecx
  149.     sub    edi,ecx
  150. skip:
  151.     ENDM
  152.  
  153.  
  154.  
  155. ;long squareroot62(long hiarg, long loarg)
  156.  
  157. ; takes root of 62 bit number to 31 bit result
  158. ; about 500 clocks worst case:
  159. ; 8 us on 486/66 and 20 us on 386/25
  160.  
  161.  
  162. hiarg    equ    DWORD PTR [bp+8]
  163. loarg    equ    DWORD PTR [bp+12]
  164.  
  165.     PUBLIC    _squareroot62
  166.  
  167. _squareroot62 proc    far
  168.  
  169.     .386
  170.     push    ebp
  171.     mov    ebp,esp
  172.  
  173.     mov    edx,DWORD PTR hiarg
  174.     mov    eax,DWORD PTR loarg
  175.     or    edx,edx
  176.     jne    dohigh
  177.  
  178.     push    eax        ; can use short root!
  179.     call    _squareroot32
  180.     sub    esp,4
  181.     mov    esp,ebp
  182.     pop    ebp
  183.     ret
  184.  
  185. dohigh:
  186.     push    ecx        ; have to do 2 dwords
  187.     push    esi
  188.     push    edi
  189.  
  190.     xor    edi,edi
  191.     xor    ebx,ebx
  192.  
  193.     shld    edx,eax,2    ; prescale for 62 bits in 64 bit word
  194.     shl    eax,2
  195.  
  196.     test    edx,0FFFF0000h        ; can we cut it in half?
  197.     jne    hashigh
  198.  
  199.     shld    edx,eax,16      ; yes, so prescale
  200.     shl    eax,16
  201.     test    edx,0FF000000h    ; half again?
  202.     jne    do48            ; no, do 48 loops
  203.  
  204.     shld    edx,eax,8
  205.     shl    eax,8        ; yes, prescale
  206.     jmp    do40        ; do 40 loops
  207.  
  208. hashigh:
  209.     test    edx,0FF000000h        ; half again?
  210.     jne    do64
  211.     shld    edx,eax,8
  212.     shl    eax,8
  213.     jmp    do56
  214.  
  215. do64:
  216.     SQRT64
  217.     SQRT64
  218.     SQRT64
  219.     SQRT64
  220. do56:
  221.     SQRT64
  222.     SQRT64
  223.     SQRT64
  224.     SQRT64
  225. do48:
  226.     SQRT64
  227.     SQRT64
  228.     SQRT64
  229.     SQRT64
  230. do40:
  231.     SQRT64
  232.     SQRT64
  233.     SQRT64
  234.     SQRT64
  235.  
  236.     SQRT64
  237.     SQRT64
  238.     SQRT64
  239.     SQRT64
  240.  
  241.     SQRT64
  242.     SQRT64
  243.     SQRT64
  244.     SQRT64
  245.  
  246.     SQRT64
  247.     SQRT64
  248.     SQRT64
  249.     SQRT64
  250.  
  251.     SQRT64
  252.     SQRT64
  253.     SQRT64
  254.         ; one missing because of prescale
  255.  
  256.     mov    eax,ebx
  257.  
  258.     pop    edi
  259.     pop    esi
  260.     pop    ecx
  261.  
  262.     shld    edx,eax,16    ; returns in both eax and dx:ax
  263.  
  264.     mov    esp,ebp
  265.     pop    ebp
  266.     ret
  267.  
  268. _squareroot62  endp
  269.  
  270.  
  271. ;long magnitude32(long x, long y, long z)
  272.  
  273. ; computes overall magnitude of vector
  274. ; no scaling or shortcuts: does 3x32-bit multiplies
  275. ; time: worst case of 650 clocks, best case of 200
  276. ; 3 to 10 us on 486/66, 8 to 25 us on 386/25
  277.  
  278.  
  279. x    equ    DWORD PTR [bp+8]
  280. y    equ    DWORD PTR [bp+12]
  281. z    equ    DWORD PTR [bp+16]
  282.  
  283.     PUBLIC    _magnitude32
  284.  
  285. _magnitude32 proc    far
  286.  
  287.     .386
  288.     push    ebp
  289.     mov    ebp,esp
  290.     push    ecx
  291.  
  292.     mov    eax,x           ; sum of squares
  293.     imul    x
  294.     mov    ebx,eax
  295.     mov    ecx,edx
  296.  
  297.     mov    eax,y
  298.     imul    y
  299.     add    ebx,eax
  300.     adc    ecx,edx
  301.  
  302.     mov    eax,z
  303.     imul    z
  304.     add    ebx,eax
  305.     adc    ecx,edx
  306.  
  307.     push    ebx        ; square root
  308.     push    ecx
  309.     call    _squareroot62
  310.     add    esp,8
  311.  
  312.     pop    ecx
  313.     mov    esp,ebp
  314.     pop    ebp
  315.     ret
  316.  
  317. _magnitude32  endp
  318.  
  319.  
  320. ;long magnitude16(int x, int y, int z)
  321.  
  322. ; computes overall magnitude of vector
  323. ; no scaling or shortcuts: does 3x16-bit multiplies
  324. ; time: worst case of 300 clocks, best case of 150
  325. ; 2 to 5 us on 486/66, 6 to 12 us on 386/25
  326.  
  327.  
  328. x    equ    WORD PTR [bp+8]
  329. y    equ    WORD PTR [bp+10]
  330. z    equ    WORD PTR [bp+12]
  331.  
  332.     PUBLIC    _magnitude16
  333.  
  334. _magnitude16 proc    far
  335.  
  336.     .386
  337.     push    ebp
  338.     mov    ebp,esp
  339.     push    ecx
  340.  
  341.     mov    ax,x           ; sum of squares
  342.     imul    x
  343.     mov    bx,ax
  344.     mov    cx,dx
  345.  
  346.     mov    ax,y
  347.     imul    y
  348.     add    bx,ax
  349.     adc    cx,dx
  350.  
  351.     mov    ax,z
  352.     imul    z
  353.     add    bx,ax
  354.     adc    cx,dx
  355.  
  356.     push    cx        ; square root
  357.     push    bx
  358.     call    _squareroot32
  359.     add    esp,4
  360.  
  361.     pop    ecx
  362.     mov    esp,ebp
  363.     pop    ebp
  364.     ret
  365.  
  366. _magnitude16  endp
  367.  
  368.  
  369.  
  370. ;void set_vector_length32(long length, long *xp, long *yp, long *zp)
  371.  
  372. ; sets overall magnitude of vector
  373.  
  374.  
  375. length    equ    DWORD PTR [bp+8]
  376. xp    equ    DWORD PTR [bp+12]
  377. yp    equ    DWORD PTR [bp+16]
  378. zp    equ    DWORD PTR [bp+20]
  379.  
  380.  
  381.     PUBLIC    _set_vector_magnitude32
  382.  
  383. _set_vector_magnitude32 proc    far
  384.  
  385.     .386
  386.     push    ebp
  387.     mov    ebp,esp
  388.     sub    esp,20
  389.  
  390.     push    ecx
  391.     push    edi
  392.     push    esi
  393.  
  394.     les    bx,xp              ; compute magnitude
  395.     push    DWORD PTR es:[bx]
  396.     les    bx,yp
  397.     push    DWORD PTR es:[bx]
  398.     les    bx,zp
  399.     push    DWORD PTR es:[bx]
  400.  
  401.     call    _magnitude32
  402.     add    esp,12
  403.     mov    esi,eax
  404.     or    eax,eax
  405.     je    zero_magnitude
  406.  
  407.     les    bx,xp              ; scale each part
  408.     mov    eax, es:[bx]
  409.     imul    length
  410.     idiv    esi
  411.     les    bx,xp
  412.     mov    es:[bx],eax
  413.  
  414.     les    bx,yp
  415.     mov    eax, es:[bx]
  416.     imul    length
  417.     idiv    esi
  418.     les    bx,yp
  419.     mov    es:[bx],eax
  420.  
  421.     les    bx,zp
  422.     mov    eax, es:[bx]
  423.     imul    length
  424.     idiv    esi
  425.     les    bx,zp
  426.     mov    es:[bx],eax
  427.  
  428. zero_magnitude:
  429.     pop    esi
  430.     pop    edi
  431.     pop    ecx
  432.  
  433.     mov    esp,ebp
  434.     pop    ebp
  435.     ret
  436.  
  437. _set_vector_magnitude32  endp
  438.  
  439.  
  440.  
  441.     end