home *** CD-ROM | disk | FTP | other *** search
/ EuroCD 3 / EuroCD 3.iso / Programming / vbcc / machines / amiga68k / libsrc / math / math_040 / asin.s < prev    next >
Encoding:
Text File  |  1998-06-24  |  4.7 KB  |  114 lines
  1. *
  2. *   $VER: asin.s 33.1 (22.1.97)
  3. *
  4. *   Calculates the arcsine 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    _asin
  17.     XDEF    @asin
  18.  
  19.     XREF    @atan
  20.  
  21. *************************************************************************
  22. * asin():  computes the inverse sine of a normalized input              *
  23. *                                                                       *
  24. * INPUT *************************************************************** *
  25. *       fp0 = extended precision input                                  *
  26. *                                                                       *
  27. * OUTPUT ************************************************************** *
  28. *       fp0 = arcsin(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. *       ASIN                                                            *
  39. *       1. If |X| >= 1, go to 3.                                        *
  40. *                                                                       *
  41. *       2. (|X| < 1) Calculate asin(X) by                               *
  42. *               z := sqrt( [1-X][1+X] )                                 *
  43. *               asin(X) = atan( x / z ).                                *
  44. *               Exit.                                                   *
  45. *                                                                       *
  46. *       3. If |X| > 1, go to 5.                                         *
  47. *                                                                       *
  48. *       4. (|X| = 1) sgn := sign(X), return asin(X) := sgn * Pi/2. Exit.*
  49. *                                                                       *
  50. *       5. (|X| > 1) Generate an invalid operation by 0 * infinity.     *
  51. *               Exit.                                                   *
  52. *                                                                       *
  53. *************************************************************************
  54.  
  55. PIBY2   dc.l            $3FFF0000,$C90FDAA2,$2168C235,$00000000
  56.  
  57. _asin
  58.         fmove.d         (4,sp),fp0
  59. @asin
  60.         fmove.x         fp0,-(sp)
  61.  
  62.         move.l          (sp),d1
  63.         move.w          (4,sp),d1
  64.         and.l           #$7FFFFFFF,d1
  65.         cmp.l           #$3FFF8000,d1
  66.         bge.b           .ASINBIG
  67.  
  68. ; This catch is added here for the '060 QSP. Originally, the call to
  69. ; satan() would handle this case by causing the exception which would
  70. ; not be caught until gen_except(). Now, with the exceptions being
  71. ; detected inside of satan(), the exception would have been handled there
  72. ; instead of inside sasin() as expected.
  73.         cmp.l           #$3FD78000,d1
  74.         blt.w           .ASINTINY
  75.  
  76. ;--THIS IS THE USUAL CASE, |X| < 1
  77. ;--ASIN(X) = ATAN( X / SQRT( (1-X)(1+X) ) )
  78.  
  79. .ASINMAIN
  80.         lea             (12,sp),sp
  81.         fmove.s         #$3F800000,fp1
  82.         fsub.x          fp0,fp1                 ; 1-X
  83.         fmovem.x        fp2,-(sp)               ;  {fp2}
  84.         fmove.s         #$3F800000,fp2
  85.         fadd.x          fp0,fp2                 ; 1+X
  86.         fmul.x          fp2,fp1                 ; (1+X)(1-X)
  87.         fmovem.x        (sp)+,fp2               ;  {fp2}
  88.         fsqrt.x         fp1                     ; SQRT([1-X][1+X])
  89.         fdiv.x          fp1,fp0                 ; X/SQRT([1-X][1+X])
  90.         jsr             @atan
  91.         rts
  92.  
  93. .ASINBIG
  94.         fabs.x          fp0                     ; |X|
  95.         fcmp.s          #$3F800000,fp0
  96.         fbgt            .operr                  ; cause an operr exception
  97.  
  98. ;--|X| = 1, ASIN(X) = +- PI/2.
  99. .ASINONE
  100.         fmove.x         (PIBY2,pc),fp0
  101.         move.l          (sp),d1
  102.         and.l           #$80000000,d1           ; SIGN BIT OF X
  103.         lea             (12,sp),sp
  104.         or.l            #$3F800000,d1           ; +-1 IN SGL FORMAT
  105.         move.l          d1,-(sp)                ; push SIGN(X) IN SGL-FMT
  106.         fmul.s          (sp)+,fp0
  107.         rts
  108.  
  109. ;--|X| < 2^(-40), ATAN(X) = X
  110. .ASINTINY
  111. .operr
  112.         lea             (12,sp),sp
  113.         rts
  114.