home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 364b.lha / PCQ_v1.1 / Runtime / PCQ / SinCos.asm < prev    next >
Encoding:
Assembly Source File  |  1990-04-08  |  6.3 KB  |  218 lines
  1. *
  2. *    SinCos.asm (of PCQ Pascal)
  3. *
  4. *    These routines replace earlier attempts at the sine and
  5. *    cosine functions.  They are based on the series
  6. *    representations of the functions, so they're much more
  7. *    accurate for most of the curve.
  8. *
  9. *    Martin Combs sent me a Pascal version of the implementation
  10. *    used here.  I had been looking at various series for exp()
  11. *    and ln(), but I didn't like the algorithms I was coming up
  12. *    with for sine and cosine.  Thus I appreciate Martin's
  13. *    contribution.
  14. *
  15. *    Martin's original Pascal functions duplicated by the assembly
  16. *    below are:
  17. *
  18. * FUNCTION sine(x:REAL):REAL;
  19. * 
  20. *   FUNCTION GetVal : REAL;
  21. *     CONST rfac2 = 0.5;        {rfac stands for reciprocal of factorial}
  22. *           rfac3 = 0.1666667;
  23. *           rfac4 = 0.0416667;
  24. *           rfac5 = 0.0083333;
  25. *           sincosfourthpi = 0.7071068; {sine or cosine of pi/4}
  26. *     VAR   x2 : REAL;
  27. * {This variation gives four significant figure accuracy}
  28. *   BEGIN
  29. *     IF x>fourthpi THEN BEGIN   { if x in second octant }
  30. *       x:=x-fourthpi;           { reduce x to first octant }
  31. *       x2:=sqr(x);
  32. *       x2:=x2*(rfac4*x2-rfac2)+1.0; {calculates cosine in first octant}
  33. *    { above is cos(x)=1-(x^2)/2!+(x^4)/4! }
  34. *       GETVAL:=sincosfourthpi*(sine(x)+x2)  { note the recursion }
  35. * { Note: sin(x)=sin((x-pi/4)+pi/4)=sin(x-pi/4)cos(pi/4)+cos(x-pi/4)sin(pi/4)}
  36. * { =.7071068(sin(x-pi/4)+cos(x-pi/4)) }
  37. *     END ELSE
  38. *     x2:=sqr(x);
  39. *     GetVal:=x*(x2*(rfac5*x2-rfac3)+1.0)   {calculates sine in first octant}
  40. *  {above is sin(x)=x-(x^3)/3!+(x^5)/5! }
  41. *   END;
  42. * 
  43. * {This variation gives two significant figure accuracy, but faster}
  44. * {  BEGIN
  45. *     IF x>fourthpi THEN BEGIN
  46. *       x:=x-fourthpi;
  47. *       x2:=sqr(x);
  48. *       GetVal:=sincosfourthpi*(sine(x)+1.0-rfac2*x2)  }
  49. * {cosine portion of above is 1-(x^2)/2! }
  50. * {    END ELSE
  51. *     x2:=sqr(x);
  52. *     GetVal:=x*(1.0-rfac3*x2)
  53. *   END;   }
  54. * 
  55. * 
  56. * BEGIN
  57. *   IF (x>twopi) OR (x<0.0) THEN  { faster than always doing it }
  58. *     x:= x-Floor(x/twopi)*twopi;
  59. *   IF x<halfpi THEN          { first quadrant }
  60. *     sine:=GetVal
  61. *   ELSE IF x<pi THEN BEGIN   { second quadrant }
  62. *     x:=pi-x;
  63. *     sine:=GetVal
  64. *   END ELSE IF x<threehalvespi THEN BEGIN   { third quadrant }
  65. *     x:=x-pi;
  66. *     sine:=-GetVal
  67. *   END ELSE BEGIN            { fourth quadrant }
  68. *     x:=twopi-x;
  69. *     sine:=-GetVal
  70. *   END
  71. * END;
  72. * 
  73. * FUNCTION cosine(x:REAL):REAL;
  74. * BEGIN
  75. *   cosine:=sine(halfpi-x)
  76. * END;
  77. *
  78.  
  79.     XREF    _p%MathBase
  80.  
  81.     XREF    _LVOSPAdd
  82.     XREF    _LVOSPSub
  83.     XREF    _LVOSPMul
  84.     XREF    _LVOSPFloor
  85.     XREF    _LVOSPCmp
  86.     XREF    _LVOSPNeg
  87.     XREF    _LVOSPDiv
  88.  
  89. *
  90. *    GetVal returns the sine of theta, held in d0, which is
  91. *    an angle between 0 and Pi/2
  92. *
  93.  
  94. GetVal
  95.     move.l    d0,-(sp)    ; save theta
  96.     move.l    #$C90FD940,d1    ; d1 := pi/4
  97.     jsr    _LVOSPCmp(a6)    ; compare 'em
  98.     ble    FirstOctant    ; if theta <= pi/4, skip this
  99.  
  100.     move.l    (sp),d0        ; d0 := theta
  101.     move.l    #$C90FD940,d1    ; d1 := pi/4
  102.     jsr    _LVOSPSub(a6)    ; d0 := theta - pi/4
  103.     move.l    d0,(sp)        ; save new theta
  104.     move.l    d0,d1        ; set up for square
  105.     jsr    _LVOSPMul(a6)    ; d0 := theta squared
  106.     move.l    d0,-(sp)    ; save theta2
  107.  
  108.     move.l    #$AAAAB33C,d1    ; d1 := 1/4!
  109.     jsr    _LVOSPMul(a6)    ; d0 := theta^2 * 1/4!
  110.     move.l    #$80000040,d1    ; d1 := 1/2!
  111.     jsr    _LVOSPSub(a6)    ; d0 := theta^2 * 1/4! - 1/2!
  112.     move.l    (sp)+,d1    ; d1 := theta^2
  113.     jsr    _LVOSPMul(a6)    ; d0 := t^2 * (t^2/4! - 1/2!) = t^4/4! - t^2/2!
  114.     move.l    #$80000041,d1    ; d1 := 1.0
  115.     jsr    _LVOSPAdd(a6)    ; d0 := everything + 1.0
  116.  
  117.     move.l    d0,-(sp)    ; save this result
  118.     move.l    4(sp),-(sp)    ; put theta back on stack
  119.     bsr    _p%sin        ; call sine recursively
  120.     addq.l    #4,sp        ; and pop it back off
  121.     move.l    (sp)+,d1    ; get theta^2
  122.     jsr    _LVOSPAdd(a6)    ; add them
  123.     move.l    #$B504F440,d1    ; d1 := sin(pi/4)
  124.     jsr    _LVOSPMul(a6)    ; multiply by sin(pi/4)
  125.     addq.l    #4,sp        ; pop previous theta
  126.     rts            ; and return
  127.  
  128. FirstOctant
  129.     move.l    (sp),d0        ; d0 := theta
  130.     move.l    d0,d1        ; so does d1
  131.     jsr    _LVOSPMul(a6)    ; d0 := theta^2
  132.     move.l    d0,-(sp)    ; save this
  133.  
  134.     move.l    #$8888653A,d1    ; d1 := 1/5! = 1/96
  135.     jsr    _LVOSPMul(a6)    ; d0 := theta^2 * 1/5!
  136.     move.l    #$AAAAAD3E,d1    ; d1 := 1/3!
  137.     jsr    _LVOSPSub(a6)    ; d0 := theta^2/5! - 1/3!
  138.     move.l    (sp)+,d1    ; pop theta^2
  139.     jsr    _LVOSPMul(a6)    ; d0 := theta^4/5! - theta^2/3!
  140.     move.l    #$80000041,d1    ; d1 := 1.0
  141.     jsr    _LVOSPAdd(a6)    ; d0 := above + 1.0
  142.     move.l    (sp)+,d1    ; d1 := theta
  143.     jmp    _LVOSPMul(a6)    ; d0 := theta - theta^3/3! + theta^5/5!
  144.                 ; return with that value
  145.  
  146.  
  147. * d0 contains the angle, which is also on the stack
  148.  
  149.     XDEF    _p%sin
  150. _p%sin
  151.     move.l    _p%MathBase,a6    ; we'll use this all over
  152.     move.l    #$C90FD943,d1    ; d1 := 2Pi
  153.     jsr    _LVOSPCmp(a6)    ; compare theta and 2Pi
  154.     bgt.s    ShrinkIt    ; normalize to 0 <= theta <= 2Pi
  155.     move.l    4(sp),d0    ; d0 := theta
  156.     moveq.l    #0,d1        ; d1 := 0.0
  157.     jsr    _LVOSPCmp(a6)    ; compare theta & 0.0
  158.     bge.s    AfterShrink    ; don't bother with normalization
  159. ShrinkIt
  160.     move.l    4(sp),d0    ; d0 := theta
  161.     move.l    #$C90FD943,d1    ; d1 := 2Pi
  162.     jsr    _LVOSPDiv(a6)    ; d0 := theta / 2Pi
  163.     jsr    _LVOSPFloor(a6)    ; d0 := |_theta / 2Pi_| ( floor )
  164.     move.l    #$C90FD943,d1    ; d1 := 2Pi
  165.     jsr    _LVOSPMul(a6)    ; d0 := floor * 2Pi
  166.     move.l    d0,d1        ; shift for upcoming subtract
  167.     move.l    4(sp),d0    ; d0 := theta
  168.     jsr    _LVOSPSub(a6)    ; d0 := x - the rest
  169.     move.l    d0,4(sp)    ; save as new theta, and fall through
  170.  
  171. AfterShrink
  172.     move.l    4(sp),d0    ; for not normalized case
  173.     move.l    #$C90FD941,d1    ; d1 := Pi/2
  174.     jsr    _LVOSPCmp(a6)    ; compare theta and Pi/2
  175.     bge.s    NotFirst    ; it's not in first quadrant
  176.     move.l    4(sp),d0    ; get theta
  177.     bra    GetVal        ; get value and leave
  178.  
  179. NotFirst
  180.     move.l    #$C90FD942,d1    ; d1 := Pi
  181.     move.l    4(sp),d0    ; d0 := theta
  182.     jsr    _LVOSPCmp(a6)    ; compare these
  183.     bge.s    NotSecond    ; it's not in first or second
  184.  
  185.     move.l    #$C90FD942,d0    ; d1 := Pi
  186.     move.l    4(sp),d1    ; d0 := theta
  187.     jsr    _LVOSPSub(a6)    ; d0 := Pi - theta
  188.     bra    GetVal        ; get value and leave
  189.  
  190. NotSecond
  191.     move.l    #$96CBE343,d1    ; d1 := 3Pi/2
  192.     move.l    4(sp),d0    ; d0 := theta
  193.     jsr    _LVOSPCmp(a6)    ; compare theta & 3Pi/2
  194.     bge    FourthQuad    ; definitely
  195.  
  196.     move.l    4(sp),d1    ; d1 := theta
  197.     move.l    #$C90FD942,d0    ; d0 := Pi
  198.     jsr    _LVOSPSub(a6)    ; d0 := Pi - theta
  199.     bra    GetVal        ; get value, and return
  200.  
  201. FourthQuad
  202.     move.l    #$C90FD943,d0    ; d0 := 2Pi
  203.     move.l    4(sp),d1    ; d1 := theta
  204.     jsr    _LVOSPSub(a6)    ; d0 := 2Pi - theta
  205.     bsr    GetVal        ; get sine here
  206.     jmp    _LVOSPNeg(a6)    ; return negative
  207.  
  208.     XDEF    _p%cos
  209. _p%cos
  210.     move.l    _p%MathBase,a6
  211.     move.l    4(sp),d1    ; d1 := theta
  212.     move.l    #$C90FD941,d0    ; d0 := Pi/2
  213.     jsr    _LVOSPSub(a6)    ; d0 := Pi/2 - theta
  214.     move.l    d0,4(sp)    ; save as new theta
  215.     bra    _p%sin        ; call sin, leave from there
  216.  
  217.     END
  218.