File indexing completed on 2025-05-11 08:23:49
0001 #include "fpsp-namespace.h"
0002 //
0003 //
0004 // ssinh.sa 3.1 12/10/90
0005 //
0006 // The entry point sSinh computes the hyperbolic sine of
0007 // an input argument; sSinhd does the same except for denormalized
0008 // input.
0009 //
0010 // Input: Double-extended number X in location pointed to
0011 // by address register a0.
0012 //
0013 // Output: The value sinh(X) returned in floating-point register Fp0.
0014 //
0015 // Accuracy and Monotonicity: The returned result is within 3 ulps in
0016 // 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
0017 // result is subsequently rounded to double precision. The
0018 // result is provably monotonic in double precision.
0019 //
0020 // Speed: The program sSINH takes approximately 280 cycles.
0021 //
0022 // Algorithm:
0023 //
0024 // SINH
0025 // 1. If |X| > 16380 log2, go to 3.
0026 //
0027 // 2. (|X| <= 16380 log2) Sinh(X) is obtained by the formulae
0028 // y = |X|, sgn = sign(X), and z = expm1(Y),
0029 // sinh(X) = sgn*(1/2)*( z + z/(1+z) ).
0030 // Exit.
0031 //
0032 // 3. If |X| > 16480 log2, go to 5.
0033 //
0034 // 4. (16380 log2 < |X| <= 16480 log2)
0035 // sinh(X) = sign(X) * exp(|X|)/2.
0036 // However, invoking exp(|X|) may cause premature overflow.
0037 // Thus, we calculate sinh(X) as follows:
0038 // Y := |X|
0039 // sgn := sign(X)
0040 // sgnFact := sgn * 2**(16380)
0041 // Y' := Y - 16381 log2
0042 // sinh(X) := sgnFact * exp(Y').
0043 // Exit.
0044 //
0045 // 5. (|X| > 16480 log2) sinh(X) must overflow. Return
0046 // sign(X)*Huge*Huge to generate overflow and an infinity with
0047 // the appropriate sign. Huge is the largest finite number in
0048 // extended format. Exit.
0049 //
0050
0051 // Copyright (C) Motorola, Inc. 1990
0052 // All Rights Reserved
0053 //
0054 // THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA
0055 // The copyright notice above does not evidence any
0056 // actual or intended publication of such source code.
0057
0058 //SSINH idnt 2,1 | Motorola 040 Floating Point Software Package
0059
0060 |section 8
0061
0062 T1: .long 0x40C62D38,0xD3D64634 // ... 16381 LOG2 LEAD
0063 T2: .long 0x3D6F90AE,0xB1E75CC7 // ... 16381 LOG2 TRAIL
0064
0065 |xref t_frcinx
0066 |xref t_ovfl
0067 |xref t_extdnrm
0068 |xref setox
0069 |xref setoxm1
0070
0071 .global ssinhd
0072 ssinhd:
0073 //--SINH(X) = X FOR DENORMALIZED X
0074
0075 bra t_extdnrm
0076
0077 .global ssinh
0078 ssinh:
0079 fmovex (%a0),%fp0 // ...LOAD INPUT
0080
0081 movel (%a0),%d0
0082 movew 4(%a0),%d0
0083 movel %d0,%a1 // save a copy of original (compacted) operand
0084 andl #0x7FFFFFFF,%d0
0085 cmpl #0x400CB167,%d0
0086 bgts SINHBIG
0087
0088 //--THIS IS THE USUAL CASE, |X| < 16380 LOG2
0089 //--Y = |X|, Z = EXPM1(Y), SINH(X) = SIGN(X)*(1/2)*( Z + Z/(1+Z) )
0090
0091 fabsx %fp0 // ...Y = |X|
0092
0093 moveml %a1/%d1,-(%sp)
0094 fmovemx %fp0-%fp0,(%a0)
0095 clrl %d1
0096 bsr setoxm1 // ...FP0 IS Z = EXPM1(Y)
0097 fmovel #0,%fpcr
0098 moveml (%sp)+,%a1/%d1
0099
0100 fmovex %fp0,%fp1
0101 fadds #0x3F800000,%fp1 // ...1+Z
0102 fmovex %fp0,-(%sp)
0103 fdivx %fp1,%fp0 // ...Z/(1+Z)
0104 movel %a1,%d0
0105 andl #0x80000000,%d0
0106 orl #0x3F000000,%d0
0107 faddx (%sp)+,%fp0
0108 movel %d0,-(%sp)
0109
0110 fmovel %d1,%fpcr
0111 fmuls (%sp)+,%fp0 //last fp inst - possible exceptions set
0112
0113 bra t_frcinx
0114
0115 SINHBIG:
0116 cmpl #0x400CB2B3,%d0
0117 bgt t_ovfl
0118 fabsx %fp0
0119 fsubd T1(%pc),%fp0 // ...(|X|-16381LOG2_LEAD)
0120 movel #0,-(%sp)
0121 movel #0x80000000,-(%sp)
0122 movel %a1,%d0
0123 andl #0x80000000,%d0
0124 orl #0x7FFB0000,%d0
0125 movel %d0,-(%sp) // ...EXTENDED FMT
0126 fsubd T2(%pc),%fp0 // ...|X| - 16381 LOG2, ACCURATE
0127
0128 movel %d1,-(%sp)
0129 clrl %d1
0130 fmovemx %fp0-%fp0,(%a0)
0131 bsr setox
0132 fmovel (%sp)+,%fpcr
0133
0134 fmulx (%sp)+,%fp0 //possible exception
0135 bra t_frcinx
0136
0137 |end
