1
2 /*******************************************************************************
3 MIT License
4 -----------
5
6 Copyright (c) 2002-2019 Advanced Micro Devices, Inc.
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this Software and associated documentaon files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in
16 all copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24 THE SOFTWARE.
25 *******************************************************************************/
26
27 #include "libm.h"
28 #include "libm_util.h"
29
30 #define USE_SPLITEXP
31 #define USE_SCALEDOUBLE_1
32 #define USE_SCALEDOUBLE_2
33 #define USE_INFINITY_WITH_FLAGS
34 #define USE_VAL_WITH_FLAGS
35 #define USE_HANDLE_ERROR
36 #include "libm_inlines.h"
37 #undef USE_SPLITEXP
38 #undef USE_SCALEDOUBLE_1
39 #undef USE_SCALEDOUBLE_2
40 #undef USE_INFINITY_WITH_FLAGS
41 #undef USE_VAL_WITH_FLAGS
42 #undef USE_HANDLE_ERROR
43
44 #include "libm_errno.h"
45
46 #ifdef _MSC_VER
47 #pragma function(sinh)
48 #endif
49
sinh(double x)50 double sinh(double x)
51 {
52 /*
53 After dealing with special cases the computation is split into
54 regions as follows:
55
56 abs(x) >= max_sinh_arg:
57 sinh(x) = sign(x)*Inf
58
59 abs(x) >= small_threshold:
60 sinh(x) = sign(x)*exp(abs(x))/2 computed using the
61 splitexp and scaleDouble functions as for exp_amd().
62
63 abs(x) < small_threshold:
64 compute p = exp(y) - 1 and then z = 0.5*(p+(p/(p+1.0)))
65 sinh(x) is then sign(x)*z. */
66
67 static const double
68 max_sinh_arg = 7.10475860073943977113e+02, /* 0x408633ce8fb9f87e */
69 thirtytwo_by_log2 = 4.61662413084468283841e+01, /* 0x40471547652b82fe */
70 log2_by_32_lead = 2.16608493356034159660e-02, /* 0x3f962e42fe000000 */
71 log2_by_32_tail = 5.68948749532545630390e-11, /* 0x3dcf473de6af278e */
72 small_threshold = 8*BASEDIGITS_DP64*0.30102999566398119521373889;
73 /* (8*BASEDIGITS_DP64*log10of2) ' exp(-x) insignificant c.f. exp(x) */
74
75 /* Lead and tail tabulated values of sinh(i) and cosh(i)
76 for i = 0,...,36. The lead part has 26 leading bits. */
77
78 static const double sinh_lead[37] = {
79 0.00000000000000000000e+00, /* 0x0000000000000000 */
80 1.17520117759704589844e+00, /* 0x3ff2cd9fc0000000 */
81 3.62686038017272949219e+00, /* 0x400d03cf60000000 */
82 1.00178747177124023438e+01, /* 0x40240926e0000000 */
83 2.72899169921875000000e+01, /* 0x403b4a3800000000 */
84 7.42032089233398437500e+01, /* 0x40528d0160000000 */
85 2.01713153839111328125e+02, /* 0x406936d228000000 */
86 5.48316116333007812500e+02, /* 0x4081228768000000 */
87 1.49047882080078125000e+03, /* 0x409749ea50000000 */
88 4.05154187011718750000e+03, /* 0x40afa71570000000 */
89 1.10132326660156250000e+04, /* 0x40c5829dc8000000 */
90 2.99370708007812500000e+04, /* 0x40dd3c4488000000 */
91 8.13773945312500000000e+04, /* 0x40f3de1650000000 */
92 2.21206695312500000000e+05, /* 0x410b00b590000000 */
93 6.01302140625000000000e+05, /* 0x412259ac48000000 */
94 1.63450865625000000000e+06, /* 0x4138f0cca8000000 */
95 4.44305525000000000000e+06, /* 0x4150f2ebd0000000 */
96 1.20774762500000000000e+07, /* 0x4167093488000000 */
97 3.28299845000000000000e+07, /* 0x417f4f2208000000 */
98 8.92411500000000000000e+07, /* 0x419546d8f8000000 */
99 2.42582596000000000000e+08, /* 0x41aceb0888000000 */
100 6.59407856000000000000e+08, /* 0x41c3a6e1f8000000 */
101 1.79245641600000000000e+09, /* 0x41dab5adb8000000 */
102 4.87240166400000000000e+09, /* 0x41f226af30000000 */
103 1.32445608960000000000e+10, /* 0x4208ab7fb0000000 */
104 3.60024494080000000000e+10, /* 0x4220c3d390000000 */
105 9.78648043520000000000e+10, /* 0x4236c93268000000 */
106 2.66024116224000000000e+11, /* 0x424ef822f0000000 */
107 7.23128516608000000000e+11, /* 0x42650bba30000000 */
108 1.96566712320000000000e+12, /* 0x427c9aae40000000 */
109 5.34323724288000000000e+12, /* 0x4293704708000000 */
110 1.45244246507520000000e+13, /* 0x42aa6b7658000000 */
111 3.94814795284480000000e+13, /* 0x42c1f43fc8000000 */
112 1.07321789251584000000e+14, /* 0x42d866f348000000 */
113 2.91730863685632000000e+14, /* 0x42f0953e28000000 */
114 7.93006722514944000000e+14, /* 0x430689e220000000 */
115 2.15561576592179200000e+15}; /* 0x431ea215a0000000 */
116
117 static const double sinh_tail[37] = {
118 0.00000000000000000000e+00, /* 0x0000000000000000 */
119 1.60467555584448807892e-08, /* 0x3e513ae6096a0092 */
120 2.76742892754807136947e-08, /* 0x3e5db70cfb79a640 */
121 2.09697499555224576530e-07, /* 0x3e8c2526b66dc067 */
122 2.04940252448908240062e-07, /* 0x3e8b81b18647f380 */
123 1.65444891522700935932e-06, /* 0x3ebbc1cdd1e1eb08 */
124 3.53116789999998198721e-06, /* 0x3ecd9f201534fb09 */
125 6.94023870987375490695e-06, /* 0x3edd1c064a4e9954 */
126 4.98876893611587449271e-06, /* 0x3ed4eca65d06ea74 */
127 3.19656024605152215752e-05, /* 0x3f00c259bcc0ecc5 */
128 2.08687768377236501204e-04, /* 0x3f2b5a6647cf9016 */
129 4.84668088325403796299e-05, /* 0x3f09691adefb0870 */
130 1.17517985422733832468e-03, /* 0x3f53410fc29cde38 */
131 6.90830086959560562415e-04, /* 0x3f46a31a50b6fb3c */
132 1.45697262451506548420e-03, /* 0x3f57defc71805c40 */
133 2.99859023684906737806e-02, /* 0x3f9eb49fd80e0bab */
134 1.02538800507941396667e-02, /* 0x3f84fffc7bcd5920 */
135 1.26787628407699110022e-01, /* 0x3fc03a93b6c63435 */
136 6.86652479544033744752e-02, /* 0x3fb1940bb255fd1c */
137 4.81593627621056619148e-01, /* 0x3fded26e14260b50 */
138 1.70489513795397629181e+00, /* 0x3ffb47401fc9f2a2 */
139 1.12416073482258713767e+01, /* 0x40267bb3f55634f1 */
140 7.06579578070110514432e+00, /* 0x401c435ff8194ddc */
141 5.91244512999659974639e+01, /* 0x404d8fee052ba63a */
142 1.68921736147050694399e+02, /* 0x40651d7edccde3f6 */
143 2.60692936262073658327e+02, /* 0x40704b1644557d1a */
144 3.62419382134885609048e+02, /* 0x4076a6b5ca0a9dc4 */
145 4.07689930834187271103e+03, /* 0x40afd9cc72249aba */
146 1.55377375868385224749e+04, /* 0x40ce58de693edab5 */
147 2.53720210371943067003e+04, /* 0x40d8c70158ac6363 */
148 4.78822310734952334315e+04, /* 0x40e7614764f43e20 */
149 1.81871712615542812273e+05, /* 0x4106337db36fc718 */
150 5.62892347580489004031e+05, /* 0x41212d98b1f611e2 */
151 6.41374032312148716301e+05, /* 0x412392bc108b37cc */
152 7.57809544070145115256e+06, /* 0x415ce87bdc3473dc */
153 3.64177136406482197344e+06, /* 0x414bc8d5ae99ad14 */
154 7.63580561355670914054e+06}; /* 0x415d20d76744835c */
155
156 static const double cosh_lead[37] = {
157 1.00000000000000000000e+00, /* 0x3ff0000000000000 */
158 1.54308062791824340820e+00, /* 0x3ff8b07550000000 */
159 3.76219564676284790039e+00, /* 0x400e18fa08000000 */
160 1.00676617622375488281e+01, /* 0x402422a490000000 */
161 2.73082327842712402344e+01, /* 0x403b4ee858000000 */
162 7.42099475860595703125e+01, /* 0x40528d6fc8000000 */
163 2.01715633392333984375e+02, /* 0x406936e678000000 */
164 5.48317031860351562500e+02, /* 0x4081228948000000 */
165 1.49047915649414062500e+03, /* 0x409749eaa8000000 */
166 4.05154199218750000000e+03, /* 0x40afa71580000000 */
167 1.10132329101562500000e+04, /* 0x40c5829dd0000000 */
168 2.99370708007812500000e+04, /* 0x40dd3c4488000000 */
169 8.13773945312500000000e+04, /* 0x40f3de1650000000 */
170 2.21206695312500000000e+05, /* 0x410b00b590000000 */
171 6.01302140625000000000e+05, /* 0x412259ac48000000 */
172 1.63450865625000000000e+06, /* 0x4138f0cca8000000 */
173 4.44305525000000000000e+06, /* 0x4150f2ebd0000000 */
174 1.20774762500000000000e+07, /* 0x4167093488000000 */
175 3.28299845000000000000e+07, /* 0x417f4f2208000000 */
176 8.92411500000000000000e+07, /* 0x419546d8f8000000 */
177 2.42582596000000000000e+08, /* 0x41aceb0888000000 */
178 6.59407856000000000000e+08, /* 0x41c3a6e1f8000000 */
179 1.79245641600000000000e+09, /* 0x41dab5adb8000000 */
180 4.87240166400000000000e+09, /* 0x41f226af30000000 */
181 1.32445608960000000000e+10, /* 0x4208ab7fb0000000 */
182 3.60024494080000000000e+10, /* 0x4220c3d390000000 */
183 9.78648043520000000000e+10, /* 0x4236c93268000000 */
184 2.66024116224000000000e+11, /* 0x424ef822f0000000 */
185 7.23128516608000000000e+11, /* 0x42650bba30000000 */
186 1.96566712320000000000e+12, /* 0x427c9aae40000000 */
187 5.34323724288000000000e+12, /* 0x4293704708000000 */
188 1.45244246507520000000e+13, /* 0x42aa6b7658000000 */
189 3.94814795284480000000e+13, /* 0x42c1f43fc8000000 */
190 1.07321789251584000000e+14, /* 0x42d866f348000000 */
191 2.91730863685632000000e+14, /* 0x42f0953e28000000 */
192 7.93006722514944000000e+14, /* 0x430689e220000000 */
193 2.15561576592179200000e+15}; /* 0x431ea215a0000000 */
194
195 static const double cosh_tail[37] = {
196 0.00000000000000000000e+00, /* 0x0000000000000000 */
197 6.89700037027478056904e-09, /* 0x3e3d9f5504c2bd28 */
198 4.43207835591715833630e-08, /* 0x3e67cb66f0a4c9fd */
199 2.33540217013828929694e-07, /* 0x3e8f58617928e588 */
200 5.17452463948269748331e-08, /* 0x3e6bc7d000c38d48 */
201 9.38728274131605919153e-07, /* 0x3eaf7f9d4e329998 */
202 2.73012191010840495544e-06, /* 0x3ec6e6e464885269 */
203 3.29486051438996307950e-06, /* 0x3ecba3a8b946c154 */
204 4.75803746362771416375e-06, /* 0x3ed3f4e76110d5a4 */
205 3.33050940471947692369e-05, /* 0x3f017622515a3e2b */
206 9.94707313972136215365e-06, /* 0x3ee4dc4b528af3d0 */
207 6.51685096227860253398e-05, /* 0x3f11156278615e10 */
208 1.18132406658066663359e-03, /* 0x3f535ad50ed821f5 */
209 6.93090416366541877541e-04, /* 0x3f46b61055f2935c */
210 1.45780415323416845386e-03, /* 0x3f57e2794a601240 */
211 2.99862082708111758744e-02, /* 0x3f9eb4b45f6aadd3 */
212 1.02539925859688602072e-02, /* 0x3f85000b967b3698 */
213 1.26787669807076286421e-01, /* 0x3fc03a940fadc092 */
214 6.86652631843830962843e-02, /* 0x3fb1940bf3bf874c */
215 4.81593633223853068159e-01, /* 0x3fded26e1a2a2110 */
216 1.70489514001513020602e+00, /* 0x3ffb4740205796d6 */
217 1.12416073489841270572e+01, /* 0x40267bb3f55cb85d */
218 7.06579578098005001152e+00, /* 0x401c435ff81e18ac */
219 5.91244513000686140458e+01, /* 0x404d8fee052bdea4 */
220 1.68921736147088438429e+02, /* 0x40651d7edccde926 */
221 2.60692936262087528121e+02, /* 0x40704b1644557e0e */
222 3.62419382134890611269e+02, /* 0x4076a6b5ca0a9e1c */
223 4.07689930834187453002e+03, /* 0x40afd9cc72249abe */
224 1.55377375868385224749e+04, /* 0x40ce58de693edab5 */
225 2.53720210371943103382e+04, /* 0x40d8c70158ac6364 */
226 4.78822310734952334315e+04, /* 0x40e7614764f43e20 */
227 1.81871712615542812273e+05, /* 0x4106337db36fc718 */
228 5.62892347580489004031e+05, /* 0x41212d98b1f611e2 */
229 6.41374032312148716301e+05, /* 0x412392bc108b37cc */
230 7.57809544070145115256e+06, /* 0x415ce87bdc3473dc */
231 3.64177136406482197344e+06, /* 0x414bc8d5ae99ad14 */
232 7.63580561355670914054e+06}; /* 0x415d20d76744835c */
233
234 unsigned long long ux, aux, xneg;
235 double y, z, z1, z2;
236 int m;
237
238 /* Special cases */
239
240 GET_BITS_DP64(x, ux);
241 aux = ux & ~SIGNBIT_DP64;
242 if (aux < 0x3e30000000000000) /* |x| small enough that sinh(x) = x */
243 {
244 if (aux == 0)
245 /* with no inexact */
246 return x;
247 else
248 return val_with_flags(x, AMD_F_INEXACT);
249 }
250 else if (aux >= 0x7ff0000000000000) /* |x| is NaN or Inf */
251 {
252 if (aux > 0x7ff0000000000000)
253 /* x is NaN */
254 return _handle_error("sinh", OP_SINH, ux|0x0008000000000000, _DOMAIN,
255 0, EDOM, x, 0.0, 1);
256 else
257 return x + x;
258 }
259
260
261 xneg = (aux != ux);
262
263 y = x;
264 if (xneg) y = -x;
265
266 if (y >= max_sinh_arg)
267 {
268 if (xneg)
269 return _handle_error("sinh", OP_SINH, NINFBITPATT_DP64, _OVERFLOW,
270 AMD_F_OVERFLOW, ERANGE, x, 0.0, 1);
271 else
272 return _handle_error("sinh", OP_SINH, PINFBITPATT_DP64, _OVERFLOW,
273 AMD_F_OVERFLOW, ERANGE, x, 0.0, 1);
274 }
275 else if (y >= small_threshold)
276 {
277 /* In this range y is large enough so that
278 the negative exponential is negligible,
279 so sinh(y) is approximated by sign(x)*exp(y)/2. The
280 code below is an inlined version of that from
281 exp() with two changes (it operates on
282 y instead of x, and the division by 2 is
283 done by reducing m by 1). */
284
285 splitexp(y, 1.0, thirtytwo_by_log2, log2_by_32_lead,
286 log2_by_32_tail, &m, &z1, &z2);
287 m -= 1;
288
289 if (m >= EMIN_DP64 && m <= EMAX_DP64)
290 z = scaleDouble_1((z1+z2),m);
291 else
292 z = scaleDouble_2((z1+z2),m);
293 }
294 else
295 {
296 /* In this range we find the integer part y0 of y
297 and the increment dy = y - y0. We then compute
298
299 z = sinh(y) = sinh(y0)cosh(dy) + cosh(y0)sinh(dy)
300
301 where sinh(y0) and cosh(y0) are tabulated above. */
302
303 int ind;
304 double dy, dy2, sdy, cdy, sdy1, sdy2;
305
306 ind = (int)y;
307 dy = y - ind;
308
309 dy2 = dy*dy;
310 sdy = dy*dy2*(0.166666666666666667013899e0 +
311 (0.833333333333329931873097e-2 +
312 (0.198412698413242405162014e-3 +
313 (0.275573191913636406057211e-5 +
314 (0.250521176994133472333666e-7 +
315 (0.160576793121939886190847e-9 +
316 0.7746188980094184251527126e-12*dy2)*dy2)*dy2)*dy2)*dy2)*dy2);
317
318 cdy = dy2*(0.500000000000000005911074e0 +
319 (0.416666666666660876512776e-1 +
320 (0.138888888889814854814536e-2 +
321 (0.248015872460622433115785e-4 +
322 (0.275573350756016588011357e-6 +
323 (0.208744349831471353536305e-8 +
324 0.1163921388172173692062032e-10*dy2)*dy2)*dy2)*dy2)*dy2)*dy2);
325
326 /* At this point sinh(dy) is approximated by dy + sdy.
327 Shift some significant bits from dy to sdy. */
328
329 GET_BITS_DP64(dy, ux);
330 ux &= 0xfffffffff8000000;
331 PUT_BITS_DP64(ux, sdy1);
332 sdy2 = sdy + (dy - sdy1);
333
334 z = ((((((cosh_tail[ind]*sdy2 + sinh_tail[ind]*cdy)
335 + cosh_tail[ind]*sdy1) + sinh_tail[ind])
336 + cosh_lead[ind]*sdy2) + sinh_lead[ind]*cdy)
337 + cosh_lead[ind]*sdy1) + sinh_lead[ind];
338 }
339
340 if (xneg) z = - z;
341 return z;
342 }
343