home *** CD-ROM | disk | FTP | other *** search
/ Virtual Reality Homebrewer's Handbook / vr.iso / vr386 / intrig.asm < prev    next >
Assembly Source File  |  1996-03-19  |  8KB  |  406 lines
  1.     TITLE    INTRIG - Triginometric routines in fixed point
  2.  
  3.     COMMENT $
  4.  
  5. // All code and algorithms used are by Dave Stampe.
  6. // These are some of the fastest integer trig routines
  7. // in 32 bit fixed point around, and are accurate to 3 LSB
  8.  
  9. // See INTMATH.H for formats
  10.  
  11. /*
  12.  This code is part of the VR-386 project, created by Dave Stampe.
  13.  VR-386 is a desendent of REND386, created by Dave Stampe and
  14.  Bernie Roehl.  Almost all the code has been rewritten by Dave
  15.  Stampre for VR-386.
  16.  
  17.  Copyright (c) 1994 by Dave Stampe:
  18.  May be freely used to write software for release into the public domain
  19.  or for educational use; all commercial endeavours MUST contact Dave Stampe
  20.  (dstampe@psych.toronto.edu) for permission to incorporate any part of
  21.  this software or source code into their products!  Usually there is no
  22.  charge for under 50-100 items for low-cost or shareware products, and terms
  23.  are reasonable.  Any royalties are used for development, so equipment is
  24.  often acceptable payment.
  25.  
  26.  ATTRIBUTION:  If you use any part of this source code or the libraries
  27.  in your projects, you must give attribution to VR-386 and Dave Stampe,
  28.  and any other authors in your documentation, source code, and at startup
  29.  of your program.  Let's keep the freeware ball rolling!
  30.  
  31.  DEVELOPMENT: VR-386 is a effort to develop the process started by
  32.  REND386, improving programmer access by rewriting the code and supplying
  33.  a standard API.  If you write improvements, add new functions rather
  34.  than rewriting current functions.  This will make it possible to
  35.  include you improved code in the next API release.  YOU can help advance
  36.  VR-386.  Comments on the API are welcome.
  37.  
  38.  CONTACT: dstampe@psych.toronto.edu
  39. */
  40.  
  41.         $
  42.  
  43.     .MODEL large
  44.  
  45.     .CODE INTMATH
  46.  
  47. ;include 3dstruct.inc
  48.  
  49.     .DATA
  50.  
  51. extrn   _sintable    ; arrays of 258 dwords
  52. extrn   _atantable   ; <3.29> args -> <16.16> angle
  53.  
  54.  
  55. ; # define XFSC 536870912   /* 2**29 for shifting xform coeffs to long */
  56.  
  57.     .CODE INTMATH
  58.  
  59. ; /*************** INTEGER SINE, COSINE **************/
  60.  
  61. ;/* return(XFSC * sin(3.14159/180*angle/65536.0));*/
  62. ;long isine(long angle)           /* returns <3.29> sine of <16.16> angle */
  63.  
  64. angle    equ    [bp+8]          ; arguments
  65.  
  66.     PUBLIC    _isine
  67.  
  68. _isine     proc    far
  69.  
  70.     .386
  71.     push    ebp
  72.     mov    ebp,esp
  73.  
  74.     push    esi
  75.     push    edi
  76.     push    ecx
  77.  
  78.     mov    eax,angle
  79.     mov    edx,5B05B05Bh
  80.     imul    edx
  81.     shrd    eax,edx,29       ;/* convert so 90 deg -> 256<<16 */
  82.     adc    eax,0              ;/* 2 upper bits = quadrant      */
  83.     mov    angle,eax          ;/* 16 lsb used to interp.       */
  84.  
  85.     mov    ebx,eax
  86.     mov    esi,eax
  87.     mov    edx,eax
  88.  
  89.     shr    ebx,14             ;/* 8 bit 2.84*degree index */
  90.     and    ebx,03FCh
  91.     test    esi,01000000h
  92.     jz    forward
  93.     not    edx
  94.     xor    ebx,03FCh
  95. forward:
  96.     les    di,DWORD PTR _sintable
  97.     mov    ecx,DWORD PTR es:[bx+di]
  98.     mov    eax,DWORD PTR es:[bx+di+4]  ;/* lookup this, next entry */
  99.  
  100.     sub    eax,ecx
  101.  
  102.     and    edx,0000FFFFh       ;/* compute interp. factor */
  103.     mul     edx
  104.     shrd    eax,edx,16
  105.     adc    ecx,eax             ;/* add in */
  106.  
  107.     test    esi,02000000h
  108.     jz    positive
  109.     neg    ecx
  110. positive:
  111.     mov    eax,ecx        ;return result in dx:ax
  112.     shld    edx,eax,16      ; dlete to rtn in eax
  113.  
  114.     pop    ecx
  115.     pop    edi
  116.     pop    esi
  117.  
  118.     mov    esp,ebp
  119.     pop    ebp
  120.     ret
  121.  
  122. _isine    endp
  123.  
  124.  
  125. ; long icosine(long angle)
  126.  
  127. angle    equ    [bp+8]          ; arguments
  128.  
  129.     PUBLIC    _icosine
  130.  
  131. _icosine    proc    far
  132.  
  133.     .386
  134.     push    ebp
  135.     mov    ebp,esp
  136.  
  137.     mov    eax,angle
  138.     add    eax,005A0000h      ; just ask for sine(90+angle)
  139.     push    eax
  140.     call    far ptr _isine
  141.     add    esp,4
  142.  
  143.     mov    esp,ebp
  144.     pop    ebp
  145.     ret
  146.  
  147. _icosine    endp
  148.  
  149.  
  150. ;/************ INVERSE TRIG FUNCTIONS *************/
  151.  
  152. ;long arcsine(long x)      /* <3.29> args -> <16.16> angle */
  153. ;{                         /* have to use table in reverse */
  154.  
  155. x    equ    [bp+8]          ; arguments
  156.  
  157. sign    equ    [bp-2]        ; locals
  158.  
  159.         PUBLIC    _arcsine
  160.  
  161. _arcsine     proc    far
  162.  
  163.     .386
  164.     push    ebp
  165.     mov    ebp,esp
  166.     sub    esp,14
  167.  
  168.     push    esi
  169.     push    edi
  170.     push    ecx
  171.  
  172.     mov    WORD PTR sign,0
  173.  
  174.     mov    eax,x
  175.     or    eax,eax
  176.     je    ret_eax        ; 0->0 mapping common
  177.     jg    posarg
  178.  
  179.     inc    WORD PTR sign    ; negative argument
  180.     neg    eax
  181. posarg:
  182.     cmp    eax,20000000h    ; 1.0 or greater?
  183.     jl    notone
  184.     mov    eax,005A0000h    ; 90<<16 returned
  185.     jmp    ret_signed
  186. notone:
  187.     les    bx,DWORD PTR _sintable  ; begin binary search
  188.  
  189.     xor    ecx,ecx
  190.     cmp    eax,es:[bx+512]
  191.     jb    b1
  192.     add    ebx,512
  193.     or    ecx,00800000h
  194. b1:
  195.     cmp    eax,es:[bx+256]
  196.     jb    b2
  197.     add    ebx,256
  198.     or    ecx,00400000h
  199. b2:
  200.     cmp    eax,es:[bx+128]
  201.     jb    b3
  202.     add    ebx,128
  203.     or    ecx,00200000h
  204. b3:
  205.     cmp    eax,es:[bx+64]
  206.     jb    b4
  207.     add    ebx,64
  208.     or    ecx,00100000h
  209. b4:
  210.     cmp    eax,es:[bx+32]
  211.     jb    b5
  212.     add    ebx,32
  213.     or    ecx,00080000h
  214. b5:
  215.     cmp    eax,es:[bx+16]
  216.     jb    b6
  217.     add    ebx,16
  218.     or    ecx,00040000h
  219. b6:
  220.     cmp    eax,es:[bx+8]
  221.     jb    b7
  222.     add    ebx,8
  223.     or    ecx,00020000h
  224. b7:
  225.     cmp    eax,es:[bx+4]
  226.     jb    b8
  227.     add    ebx,4
  228.     or    ecx,00010000h
  229. b8:
  230.     sub    eax,es:[bx]    ; OK, got entry lower than x
  231.     je    nointerp    ; compute difference
  232.     mov    esi,es:[bx+4]    ; and interpolate
  233.     sub    esi,es:[bx]
  234.     je    okinterp
  235.     cmp    eax,esi
  236.     jb    okinterp
  237.     add    ecx,00010000h
  238.     jmp    nointerp
  239. okinterp:
  240.     cdq
  241.     shld    edx,eax,16
  242.     shl    eax,16
  243.     idiv    esi
  244.     mov    cx,ax
  245. nointerp:
  246.     mov    eax,05A000000h    ;/* convert to <16.16> angle */
  247.     imul    ecx
  248.     mov    eax,edx
  249. ret_signed:                       ; add sign to result
  250.     test    WORD PTR sign,-1
  251.     je    ret_eax
  252.     neg    eax
  253. ret_eax:
  254.     shld    edx,eax,16    ; eax -> dx:ax
  255.  
  256.     pop    ecx
  257.     pop    edi
  258.     pop    esi
  259.     mov    esp,ebp
  260.     pop    ebp
  261.     ret
  262.  
  263. _arcsine    endp
  264.  
  265.  
  266. ;long arccosine(long x)
  267.  
  268. x    equ    [bp+8]          ; arguments
  269.  
  270.     PUBLIC    _arccosine
  271.  
  272. _arccosine     proc    far
  273.  
  274.     .386
  275.     push    ebp
  276.     mov    ebp,esp
  277.  
  278.     push    DWORD PTR x
  279.     call    far ptr _arcsine
  280.     add    esp,4
  281.  
  282.     not    dx        ; negate dx:ax
  283.     not    ax
  284.     add    ax,1
  285.     adc    dx,90        ; (90*65536L - arcsine(x));
  286.  
  287.     mov    esp,ebp
  288.     pop    ebp
  289.     ret
  290.  
  291. _arccosine    endp
  292.  
  293.  
  294. ; /********** INVERSE TANGENT **********/
  295.  
  296.              ;/* we can use a table for atan  */
  297.              ;/* as slope is fairly constant  */
  298. ;long arctan2(long y, long x)
  299.  
  300.  
  301. y    equ    [bp+8]          ; arguments
  302. x    equ    [bp+12]
  303.  
  304. sx    equ    [bp-2]        ; locals (quadrant flags)
  305. sy    equ    [bp-4]
  306. g45    equ    [bp-6]
  307.  
  308.         PUBLIC    _arctan2
  309.  
  310. _arctan2     proc    far
  311.  
  312.     .386
  313.     push    ebp
  314.     mov    ebp,esp
  315.     sub    esp,12
  316.  
  317.     push    esi
  318.     push    edi
  319.     push    ecx
  320.  
  321.     mov    WORD PTR sx,0
  322.     mov    WORD PTR sy,0
  323.     mov    WORD PTR g45,0
  324.  
  325.     mov    eax,DWORD PTR x
  326.     or    eax,eax
  327.     jge    x_pos
  328.     inc    WORD PTR sx
  329.     neg    eax
  330. x_pos:
  331.     mov    edx,DWORD PTR y
  332.     or    edx,edx
  333.     jge    y_pos
  334.     inc    WORD PTR sy
  335.     neg    edx
  336. y_pos:
  337.     or    eax,eax
  338.     jne    nonzero_x
  339.     mov    eax,005A0000h    ; 90<<16
  340.     jmp    sign_eax
  341. nonzero_x:
  342.     or    edx,edx
  343.     jne    nonzero_y
  344.     mov    eax,0        ; 0
  345.     jmp    sign_eax
  346. nonzero_y:
  347.     cmp    edx,eax
  348.     jne    not45deg
  349.     mov    eax,002D0000h    ;/* 45 degrees */
  350.     jmp    sign_eax
  351. not45deg:
  352.     jl    not_ygtx
  353.     inc    WORD PTR g45    ; which quadrant?
  354.     xchg    eax,edx
  355. not_ygtx:
  356.     mov    ebx,eax     ; save abs(x)
  357.     xor    eax,eax
  358.     shrd    eax,edx,8      ;/* shift so it will divide OK */
  359.     shr    edx,8
  360.     idiv    ebx            ;/* y/x << 24 */
  361.     mov    ebx,eax
  362.     shr    ebx,14         ;/* top 8 bits are index into table */
  363.     and    ebx,03FCh
  364.  
  365.     les    di, DWORD PTR _atantable
  366.     mov    ecx,DWORD PTR es:[bx+di]
  367.     mov    esi,DWORD PTR es:[bx+di+4]
  368.  
  369.     sub    esi,ecx
  370.     and    eax,0000FFFFh  ;/* bottom 16 bits interpolate */
  371.     mul    esi
  372.     shr    eax,16
  373.     add    eax,ecx
  374.  
  375.     test    WORD PTR g45,-1
  376.     je    firstquad
  377.     neg    eax
  378.     add    eax,005A0000h      ; result = 90*65536L - result;
  379. firstquad:
  380. sign_eax:
  381.     test    WORD PTR sx,-1    ; put result in proper quadrant
  382.     je    secquad
  383.     sub    eax,00B40000h
  384.     test    WORD PTR sy,-1
  385.     jne    out_eax
  386.     neg    eax
  387.     jmp    out_eax
  388. secquad:
  389.     test    WORD PTR sy,-1
  390.     je    out_eax
  391.     neg    eax
  392. out_eax:
  393.     shld    edx,eax,16         ; eax -> dx:ax
  394.  
  395.     pop    ecx
  396.     pop    edi
  397.     pop    esi
  398.  
  399.     mov    esp,ebp
  400.     pop    ebp
  401.     ret
  402.  
  403. _arctan2    endp
  404.  
  405.  
  406.     end