ssinh.S 3.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134
  1. |
  2. | ssinh.sa 3.1 12/10/90
  3. |
  4. | The entry point sSinh computes the hyperbolic sine of
  5. | an input argument; sSinhd 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 sinh(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 sSINH takes approximately 280 cycles.
  19. |
  20. | Algorithm:
  21. |
  22. | SINH
  23. | 1. If |X| > 16380 log2, go to 3.
  24. |
  25. | 2. (|X| <= 16380 log2) Sinh(X) is obtained by the formulae
  26. | y = |X|, sgn = sign(X), and z = expm1(Y),
  27. | sinh(X) = sgn*(1/2)*( z + z/(1+z) ).
  28. | Exit.
  29. |
  30. | 3. If |X| > 16480 log2, go to 5.
  31. |
  32. | 4. (16380 log2 < |X| <= 16480 log2)
  33. | sinh(X) = sign(X) * exp(|X|)/2.
  34. | However, invoking exp(|X|) may cause premature overflow.
  35. | Thus, we calculate sinh(X) as follows:
  36. | Y := |X|
  37. | sgn := sign(X)
  38. | sgnFact := sgn * 2**(16380)
  39. | Y' := Y - 16381 log2
  40. | sinh(X) := sgnFact * exp(Y').
  41. | Exit.
  42. |
  43. | 5. (|X| > 16480 log2) sinh(X) must overflow. Return
  44. | sign(X)*Huge*Huge to generate overflow and an infinity with
  45. | the appropriate sign. Huge is the largest finite number in
  46. | extended format. Exit.
  47. |
  48. | Copyright (C) Motorola, Inc. 1990
  49. | All Rights Reserved
  50. |
  51. | For details on the license for this file, please see the
  52. | file, README, in this same directory.
  53. |SSINH idnt 2,1 | Motorola 040 Floating Point Software Package
  54. |section 8
  55. T1: .long 0x40C62D38,0xD3D64634 | ... 16381 LOG2 LEAD
  56. T2: .long 0x3D6F90AE,0xB1E75CC7 | ... 16381 LOG2 TRAIL
  57. |xref t_frcinx
  58. |xref t_ovfl
  59. |xref t_extdnrm
  60. |xref setox
  61. |xref setoxm1
  62. .global ssinhd
  63. ssinhd:
  64. |--SINH(X) = X FOR DENORMALIZED X
  65. bra t_extdnrm
  66. .global ssinh
  67. ssinh:
  68. fmovex (%a0),%fp0 | ...LOAD INPUT
  69. movel (%a0),%d0
  70. movew 4(%a0),%d0
  71. movel %d0,%a1 | save a copy of original (compacted) operand
  72. andl #0x7FFFFFFF,%d0
  73. cmpl #0x400CB167,%d0
  74. bgts SINHBIG
  75. |--THIS IS THE USUAL CASE, |X| < 16380 LOG2
  76. |--Y = |X|, Z = EXPM1(Y), SINH(X) = SIGN(X)*(1/2)*( Z + Z/(1+Z) )
  77. fabsx %fp0 | ...Y = |X|
  78. moveml %a1/%d1,-(%sp)
  79. fmovemx %fp0-%fp0,(%a0)
  80. clrl %d1
  81. bsr setoxm1 | ...FP0 IS Z = EXPM1(Y)
  82. fmovel #0,%fpcr
  83. moveml (%sp)+,%a1/%d1
  84. fmovex %fp0,%fp1
  85. fadds #0x3F800000,%fp1 | ...1+Z
  86. fmovex %fp0,-(%sp)
  87. fdivx %fp1,%fp0 | ...Z/(1+Z)
  88. movel %a1,%d0
  89. andl #0x80000000,%d0
  90. orl #0x3F000000,%d0
  91. faddx (%sp)+,%fp0
  92. movel %d0,-(%sp)
  93. fmovel %d1,%fpcr
  94. fmuls (%sp)+,%fp0 |last fp inst - possible exceptions set
  95. bra t_frcinx
  96. SINHBIG:
  97. cmpl #0x400CB2B3,%d0
  98. bgt t_ovfl
  99. fabsx %fp0
  100. fsubd T1(%pc),%fp0 | ...(|X|-16381LOG2_LEAD)
  101. movel #0,-(%sp)
  102. movel #0x80000000,-(%sp)
  103. movel %a1,%d0
  104. andl #0x80000000,%d0
  105. orl #0x7FFB0000,%d0
  106. movel %d0,-(%sp) | ...EXTENDED FMT
  107. fsubd T2(%pc),%fp0 | ...|X| - 16381 LOG2, ACCURATE
  108. movel %d1,-(%sp)
  109. clrl %d1
  110. fmovemx %fp0-%fp0,(%a0)
  111. bsr setox
  112. fmovel (%sp)+,%fpcr
  113. fmulx (%sp)+,%fp0 |possible exception
  114. bra t_frcinx
  115. |end