hc
2023-10-25 6c2073b7aa40e29d0eca7d571dd7bc590c7ecaa7
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
|    stanh.sa 3.1 12/10/90
|
|    The entry point sTanh computes the hyperbolic tangent of
|    an input argument; sTanhd does the same except for denormalized
|    input.
|
|    Input: Double-extended number X in location pointed to
|        by address register a0.
|
|    Output: The value tanh(X) returned in floating-point register Fp0.
|
|    Accuracy and Monotonicity: The returned result is within 3 ulps in
|        64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
|        result is subsequently rounded to double precision. The
|        result is provably monotonic in double precision.
|
|    Speed: The program stanh takes approximately 270 cycles.
|
|    Algorithm:
|
|    TANH
|    1. If |X| >= (5/2) log2 or |X| <= 2**(-40), go to 3.
|
|    2. (2**(-40) < |X| < (5/2) log2) Calculate tanh(X) by
|        sgn := sign(X), y := 2|X|, z := expm1(Y), and
|        tanh(X) = sgn*( z/(2+z) ).
|        Exit.
|
|    3. (|X| <= 2**(-40) or |X| >= (5/2) log2). If |X| < 1,
|        go to 7.
|
|    4. (|X| >= (5/2) log2) If |X| >= 50 log2, go to 6.
|
|    5. ((5/2) log2 <= |X| < 50 log2) Calculate tanh(X) by
|        sgn := sign(X), y := 2|X|, z := exp(Y),
|        tanh(X) = sgn - [ sgn*2/(1+z) ].
|        Exit.
|
|    6. (|X| >= 50 log2) Tanh(X) = +-1 (round to nearest). Thus, we
|        calculate Tanh(X) by
|        sgn := sign(X), Tiny := 2**(-126),
|        tanh(X) := sgn - sgn*Tiny.
|        Exit.
|
|    7. (|X| < 2**(-40)). Tanh(X) = X.    Exit.
|
 
|        Copyright (C) Motorola, Inc. 1990
|            All Rights Reserved
|
|       For details on the license for this file, please see the
|       file, README, in this same directory.
 
|STANH    idnt    2,1 | Motorola 040 Floating Point Software Package
 
   |section    8
 
#include "fpsp.h"
 
   .set    X,FP_SCR5
   .set    XDCARE,X+2
   .set    XFRAC,X+4
 
   .set    SGN,L_SCR3
 
   .set    V,FP_SCR6
 
BOUNDS1:    .long 0x3FD78000,0x3FFFDDCE | ... 2^(-40), (5/2)LOG2
 
   |xref    t_frcinx
   |xref    t_extdnrm
   |xref    setox
   |xref    setoxm1
 
   .global    stanhd
stanhd:
|--TANH(X) = X FOR DENORMALIZED X
 
   bra        t_extdnrm
 
   .global    stanh
stanh:
   fmovex        (%a0),%fp0    | ...LOAD INPUT
 
   fmovex        %fp0,X(%a6)
   movel        (%a0),%d0
   movew        4(%a0),%d0
   movel        %d0,X(%a6)
   andl        #0x7FFFFFFF,%d0
   cmp2l        BOUNDS1(%pc),%d0    | ...2**(-40) < |X| < (5/2)LOG2 ?
   bcss        TANHBORS
 
|--THIS IS THE USUAL CASE
|--Y = 2|X|, Z = EXPM1(Y), TANH(X) = SIGN(X) * Z / (Z+2).
 
   movel        X(%a6),%d0
   movel        %d0,SGN(%a6)
   andl        #0x7FFF0000,%d0
   addl        #0x00010000,%d0    | ...EXPONENT OF 2|X|
   movel        %d0,X(%a6)
   andl        #0x80000000,SGN(%a6)
   fmovex        X(%a6),%fp0        | ...FP0 IS Y = 2|X|
 
   movel        %d1,-(%a7)
   clrl        %d1
   fmovemx    %fp0-%fp0,(%a0)
   bsr        setoxm1        | ...FP0 IS Z = EXPM1(Y)
   movel        (%a7)+,%d1
 
   fmovex        %fp0,%fp1
   fadds        #0x40000000,%fp1    | ...Z+2
   movel        SGN(%a6),%d0
   fmovex        %fp1,V(%a6)
   eorl        %d0,V(%a6)
 
   fmovel        %d1,%FPCR        |restore users exceptions
   fdivx        V(%a6),%fp0
   bra        t_frcinx
 
TANHBORS:
   cmpl        #0x3FFF8000,%d0
   blt        TANHSM
 
   cmpl        #0x40048AA1,%d0
   bgt        TANHHUGE
 
|-- (5/2) LOG2 < |X| < 50 LOG2,
|--TANH(X) = 1 - (2/[EXP(2X)+1]). LET Y = 2|X|, SGN = SIGN(X),
|--TANH(X) = SGN -    SGN*2/[EXP(Y)+1].
 
   movel        X(%a6),%d0
   movel        %d0,SGN(%a6)
   andl        #0x7FFF0000,%d0
   addl        #0x00010000,%d0    | ...EXPO OF 2|X|
   movel        %d0,X(%a6)        | ...Y = 2|X|
   andl        #0x80000000,SGN(%a6)
   movel        SGN(%a6),%d0
   fmovex        X(%a6),%fp0        | ...Y = 2|X|
 
   movel        %d1,-(%a7)
   clrl        %d1
   fmovemx    %fp0-%fp0,(%a0)
   bsr        setox        | ...FP0 IS EXP(Y)
   movel        (%a7)+,%d1
   movel        SGN(%a6),%d0
   fadds        #0x3F800000,%fp0    | ...EXP(Y)+1
 
   eorl        #0xC0000000,%d0    | ...-SIGN(X)*2
   fmoves        %d0,%fp1        | ...-SIGN(X)*2 IN SGL FMT
   fdivx        %fp0,%fp1        | ...-SIGN(X)2 / [EXP(Y)+1 ]
 
   movel        SGN(%a6),%d0
   orl        #0x3F800000,%d0    | ...SGN
   fmoves        %d0,%fp0        | ...SGN IN SGL FMT
 
   fmovel        %d1,%FPCR        |restore users exceptions
   faddx        %fp1,%fp0
 
   bra        t_frcinx
 
TANHSM:
   movew        #0x0000,XDCARE(%a6)
 
   fmovel        %d1,%FPCR        |restore users exceptions
   fmovex        X(%a6),%fp0        |last inst - possible exception set
 
   bra        t_frcinx
 
TANHHUGE:
|---RETURN SGN(X) - SGN(X)EPS
   movel        X(%a6),%d0
   andl        #0x80000000,%d0
   orl        #0x3F800000,%d0
   fmoves        %d0,%fp0
   andl        #0x80000000,%d0
   eorl        #0x80800000,%d0    | ...-SIGN(X)*EPS
 
   fmovel        %d1,%FPCR        |restore users exceptions
   fadds        %d0,%fp0
 
   bra        t_frcinx
 
   |end