sasin.S 2.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103
  1. |
  2. | sasin.sa 3.3 12/19/90
  3. |
  4. | Description: The entry point sAsin computes the inverse sine of
  5. | an input argument; sAsind does the same except for denormalized
  6. | input.
  7. |
  8. | Input: Double-extended number X in location pointed to
  9. | by address register a0.
  10. |
  11. | Output: The value arcsin(X) returned in floating-point register Fp0.
  12. |
  13. | Accuracy and Monotonicity: The returned result is within 3 ulps in
  14. | 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  15. | result is subsequently rounded to double precision. The
  16. | result is provably monotonic in double precision.
  17. |
  18. | Speed: The program sASIN takes approximately 310 cycles.
  19. |
  20. | Algorithm:
  21. |
  22. | ASIN
  23. | 1. If |X| >= 1, go to 3.
  24. |
  25. | 2. (|X| < 1) Calculate asin(X) by
  26. | z := sqrt( [1-X][1+X] )
  27. | asin(X) = atan( x / z ).
  28. | Exit.
  29. |
  30. | 3. If |X| > 1, go to 5.
  31. |
  32. | 4. (|X| = 1) sgn := sign(X), return asin(X) := sgn * Pi/2. Exit.
  33. |
  34. | 5. (|X| > 1) Generate an invalid operation by 0 * infinity.
  35. | Exit.
  36. |
  37. | Copyright (C) Motorola, Inc. 1990
  38. | All Rights Reserved
  39. |
  40. | For details on the license for this file, please see the
  41. | file, README, in this same directory.
  42. |SASIN idnt 2,1 | Motorola 040 Floating Point Software Package
  43. |section 8
  44. PIBY2: .long 0x3FFF0000,0xC90FDAA2,0x2168C235,0x00000000
  45. |xref t_operr
  46. |xref t_frcinx
  47. |xref t_extdnrm
  48. |xref satan
  49. .global sasind
  50. sasind:
  51. |--ASIN(X) = X FOR DENORMALIZED X
  52. bra t_extdnrm
  53. .global sasin
  54. sasin:
  55. fmovex (%a0),%fp0 | ...LOAD INPUT
  56. movel (%a0),%d0
  57. movew 4(%a0),%d0
  58. andil #0x7FFFFFFF,%d0
  59. cmpil #0x3FFF8000,%d0
  60. bges asinbig
  61. |--THIS IS THE USUAL CASE, |X| < 1
  62. |--ASIN(X) = ATAN( X / SQRT( (1-X)(1+X) ) )
  63. fmoves #0x3F800000,%fp1
  64. fsubx %fp0,%fp1 | ...1-X
  65. fmovemx %fp2-%fp2,-(%a7)
  66. fmoves #0x3F800000,%fp2
  67. faddx %fp0,%fp2 | ...1+X
  68. fmulx %fp2,%fp1 | ...(1+X)(1-X)
  69. fmovemx (%a7)+,%fp2-%fp2
  70. fsqrtx %fp1 | ...SQRT([1-X][1+X])
  71. fdivx %fp1,%fp0 | ...X/SQRT([1-X][1+X])
  72. fmovemx %fp0-%fp0,(%a0)
  73. bsr satan
  74. bra t_frcinx
  75. asinbig:
  76. fabsx %fp0 | ...|X|
  77. fcmps #0x3F800000,%fp0
  78. fbgt t_operr |cause an operr exception
  79. |--|X| = 1, ASIN(X) = +- PI/2.
  80. fmovex PIBY2,%fp0
  81. movel (%a0),%d0
  82. andil #0x80000000,%d0 | ...SIGN BIT OF X
  83. oril #0x3F800000,%d0 | ...+-1 IN SGL FORMAT
  84. movel %d0,-(%sp) | ...push SIGN(X) IN SGL-FMT
  85. fmovel %d1,%FPCR
  86. fmuls (%sp)+,%fp0
  87. bra t_frcinx
  88. |end