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 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