scosh.S 2.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131
  1. |
  2. | scosh.sa 3.1 12/10/90
  3. |
  4. | The entry point sCosh computes the hyperbolic cosine of
  5. | an input argument; sCoshd 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 cosh(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 sCOSH takes approximately 250 cycles.
  19. |
  20. | Algorithm:
  21. |
  22. | COSH
  23. | 1. If |X| > 16380 log2, go to 3.
  24. |
  25. | 2. (|X| <= 16380 log2) Cosh(X) is obtained by the formulae
  26. | y = |X|, z = exp(Y), and
  27. | cosh(X) = (1/2)*( z + 1/z ).
  28. | Exit.
  29. |
  30. | 3. (|X| > 16380 log2). If |X| > 16480 log2, go to 5.
  31. |
  32. | 4. (16380 log2 < |X| <= 16480 log2)
  33. | cosh(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. | Fact := 2**(16380)
  38. | Y' := Y - 16381 log2
  39. | cosh(X) := Fact * exp(Y').
  40. | Exit.
  41. |
  42. | 5. (|X| > 16480 log2) sinh(X) must overflow. Return
  43. | Huge*Huge to generate overflow and an infinity with
  44. | the appropriate sign. Huge is the largest finite number in
  45. | extended format. Exit.
  46. |
  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. |SCOSH idnt 2,1 | Motorola 040 Floating Point Software Package
  54. |section 8
  55. |xref t_ovfl
  56. |xref t_frcinx
  57. |xref setox
  58. T1: .long 0x40C62D38,0xD3D64634 | ... 16381 LOG2 LEAD
  59. T2: .long 0x3D6F90AE,0xB1E75CC7 | ... 16381 LOG2 TRAIL
  60. TWO16380: .long 0x7FFB0000,0x80000000,0x00000000,0x00000000
  61. .global scoshd
  62. scoshd:
  63. |--COSH(X) = 1 FOR DENORMALIZED X
  64. fmoves #0x3F800000,%fp0
  65. fmovel %d1,%FPCR
  66. fadds #0x00800000,%fp0
  67. bra t_frcinx
  68. .global scosh
  69. scosh:
  70. fmovex (%a0),%fp0 | ...LOAD INPUT
  71. movel (%a0),%d0
  72. movew 4(%a0),%d0
  73. andil #0x7FFFFFFF,%d0
  74. cmpil #0x400CB167,%d0
  75. bgts COSHBIG
  76. |--THIS IS THE USUAL CASE, |X| < 16380 LOG2
  77. |--COSH(X) = (1/2) * ( EXP(X) + 1/EXP(X) )
  78. fabsx %fp0 | ...|X|
  79. movel %d1,-(%sp)
  80. clrl %d1
  81. fmovemx %fp0-%fp0,(%a0) |pass parameter to setox
  82. bsr setox | ...FP0 IS EXP(|X|)
  83. fmuls #0x3F000000,%fp0 | ...(1/2)EXP(|X|)
  84. movel (%sp)+,%d1
  85. fmoves #0x3E800000,%fp1 | ...(1/4)
  86. fdivx %fp0,%fp1 | ...1/(2 EXP(|X|))
  87. fmovel %d1,%FPCR
  88. faddx %fp1,%fp0
  89. bra t_frcinx
  90. COSHBIG:
  91. cmpil #0x400CB2B3,%d0
  92. bgts COSHHUGE
  93. fabsx %fp0
  94. fsubd T1(%pc),%fp0 | ...(|X|-16381LOG2_LEAD)
  95. fsubd T2(%pc),%fp0 | ...|X| - 16381 LOG2, ACCURATE
  96. movel %d1,-(%sp)
  97. clrl %d1
  98. fmovemx %fp0-%fp0,(%a0)
  99. bsr setox
  100. fmovel (%sp)+,%fpcr
  101. fmulx TWO16380(%pc),%fp0
  102. bra t_frcinx
  103. COSHHUGE:
  104. fmovel #0,%fpsr |clr N bit if set by source
  105. bclrb #7,(%a0) |always return positive value
  106. fmovemx (%a0),%fp0-%fp0
  107. bra t_ovfl
  108. |end