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 49| Copyright (C) Motorola, Inc. 1990 50| All Rights Reserved 51| 52| THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA 53| The copyright notice above does not evidence any 54| actual or intended publication of such source code. 55 56|SSINH idnt 2,1 | Motorola 040 Floating Point Software Package 57 58 |section 8 59 60T1: .long 0x40C62D38,0xD3D64634 | ... 16381 LOG2 LEAD 61T2: .long 0x3D6F90AE,0xB1E75CC7 | ... 16381 LOG2 TRAIL 62 63 |xref t_frcinx 64 |xref t_ovfl 65 |xref t_extdnrm 66 |xref setox 67 |xref setoxm1 68 69 .global ssinhd 70ssinhd: 71|--SINH(X) = X FOR DENORMALIZED X 72 73 bra t_extdnrm 74 75 .global ssinh 76ssinh: 77 fmovex (%a0),%fp0 | ...LOAD INPUT 78 79 movel (%a0),%d0 80 movew 4(%a0),%d0 81 movel %d0,%a1 | save a copy of original (compacted) operand 82 andl #0x7FFFFFFF,%d0 83 cmpl #0x400CB167,%d0 84 bgts SINHBIG 85 86|--THIS IS THE USUAL CASE, |X| < 16380 LOG2 87|--Y = |X|, Z = EXPM1(Y), SINH(X) = SIGN(X)*(1/2)*( Z + Z/(1+Z) ) 88 89 fabsx %fp0 | ...Y = |X| 90 91 moveml %a1/%d1,-(%sp) 92 fmovemx %fp0-%fp0,(%a0) 93 clrl %d1 94 bsr setoxm1 | ...FP0 IS Z = EXPM1(Y) 95 fmovel #0,%fpcr 96 moveml (%sp)+,%a1/%d1 97 98 fmovex %fp0,%fp1 99 fadds #0x3F800000,%fp1 | ...1+Z 100 fmovex %fp0,-(%sp) 101 fdivx %fp1,%fp0 | ...Z/(1+Z) 102 movel %a1,%d0 103 andl #0x80000000,%d0 104 orl #0x3F000000,%d0 105 faddx (%sp)+,%fp0 106 movel %d0,-(%sp) 107 108 fmovel %d1,%fpcr 109 fmuls (%sp)+,%fp0 |last fp inst - possible exceptions set 110 111 bra t_frcinx 112 113SINHBIG: 114 cmpl #0x400CB2B3,%d0 115 bgt t_ovfl 116 fabsx %fp0 117 fsubd T1(%pc),%fp0 | ...(|X|-16381LOG2_LEAD) 118 movel #0,-(%sp) 119 movel #0x80000000,-(%sp) 120 movel %a1,%d0 121 andl #0x80000000,%d0 122 orl #0x7FFB0000,%d0 123 movel %d0,-(%sp) | ...EXTENDED FMT 124 fsubd T2(%pc),%fp0 | ...|X| - 16381 LOG2, ACCURATE 125 126 movel %d1,-(%sp) 127 clrl %d1 128 fmovemx %fp0-%fp0,(%a0) 129 bsr setox 130 fmovel (%sp)+,%fpcr 131 132 fmulx (%sp)+,%fp0 |possible exception 133 bra t_frcinx 134 135 |end 136