home *** CD-ROM | disk | FTP | other *** search
/ EuroCD 3 / EuroCD 3.iso / Programming / vbcc / machines / amiga68k / libsrc / math / math_040 / acos.s next >
Encoding:
Text File  |  1998-06-24  |  4.1 KB  |  101 lines
  1. *
  2. *   $VER: acos.s 33.1 (22.1.97)
  3. *
  4. *   Calculates the arccosine of the source
  5. *
  6. *   Version history:
  7. *
  8. *   33.1    22.1.97 (c) Motorola
  9. *
  10. *           - snipped from M68060SP sources
  11. *
  12.  
  13.     machine 68040
  14.     fpu     1
  15.  
  16.     XDEF    _acos
  17.     XDEF    @acos
  18.  
  19.     XREF    @atan
  20.  
  21. *************************************************************************
  22. * acos():  computes the inverse cosine of a normalized input            *
  23. *                                                                       *
  24. * INPUT *************************************************************** *
  25. *       fp0 = extended precision input                                  *
  26. *                                                                       *
  27. * OUTPUT ************************************************************** *
  28. *       fp0 = arccos(X)                                                 *
  29. *                                                                       *
  30. * ACCURACY and MONOTONICITY ******************************************* *
  31. *       The returned result is within 3 ulps in 64 significant bit,     *
  32. *       i.e. within 0.5001 ulp to 53 bits if the result is subsequently *
  33. *       rounded to double precision. The result is provably monotonic   *
  34. *       in double precision.                                            *
  35. *                                                                       *
  36. * ALGORITHM *********************************************************** *
  37. *                                                                       *
  38. *       ACOS                                                            *
  39. *       1. If |X| >= 1, go to 3.                                        *
  40. *                                                                       *
  41. *       2. (|X| < 1) Calculate acos(X) by                               *
  42. *               z := (1-X) / (1+X)                                      *
  43. *               acos(X) = 2 * atan( sqrt(z) ).                          *
  44. *               Exit.                                                   *
  45. *                                                                       *
  46. *       3. If |X| > 1, go to 5.                                         *
  47. *                                                                       *
  48. *       4. (|X| = 1) If X > 0, return 0. Otherwise, return Pi. Exit.    *
  49. *                                                                       *
  50. *       5. (|X| > 1) Generate an invalid operation by 0 * infinity.     *
  51. *               Exit.                                                   *
  52. *                                                                       *
  53. *************************************************************************
  54.  
  55. PI      dc.l            $40000000,$C90FDAA2,$2168C235,$00000000
  56.  
  57. _acos
  58.         fmove.d         (4,sp),fp0
  59. @acos
  60.         fmove.x         fp0,-(sp)
  61.         move.l          (sp),d1                 ; pack exp w/ upper 16 fraction
  62.         move.w          (4,sp),d1
  63.         and.l           #$7FFFFFFF,d1
  64.         cmp.l           #$3FFF8000,d1
  65.         bge.b           .ACOSBIG
  66.  
  67. ;--THIS IS THE USUAL CASE, |X| < 1
  68. ;--ACOS(X) = 2 * ATAN(  SQRT( (1-X)/(1+X) ) )
  69.  
  70. .ACOSMAIN
  71.         lea             (12,sp),sp
  72.         fmove.s         #$3F800000,fp1
  73.         fadd.x          fp0,fp1                 ; 1+X
  74.         fneg.x          fp0                     ; -X
  75.         fadd.s          #$3F800000,fp0          ; 1-X
  76.         fdiv.x          fp1,fp0                 ; (1-X)/(1+X)
  77.         fsqrt.x         fp0                     ; SQRT((1-X)/(1+X))
  78.         jsr             @atan                   ; ATAN(SQRT([1-X]/[1+X]))
  79.  
  80.         fadd.x          fp0,fp0                 ; 2 * ATAN( STUFF )
  81.         rts
  82.  
  83. .ACOSBIG
  84.         fabs.x          fp0
  85.         fcmp.s          #$3f800000,fp0
  86.         fbgt            .operr
  87.  
  88. ;--|X| = 1, ACOS(X) = 0 OR PI
  89.         tst.b           (sp)                    ; is X positive or negative?
  90.         bpl.b           .ACOSP1
  91.  
  92. ;--X = -1
  93. ;Returns PI and inexact exception
  94. .ACOSM1
  95.         fmove.x         (PI,pc),fp0             ; load PI
  96. .operr  lea             (12,sp),sp
  97.         rts
  98. .ACOSP1 fmove.s         #0,fp0
  99.         lea             (12,sp),sp
  100.         rts
  101.