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
|
/* from: FreeBSD: head/lib/msun/src/e_coshl.c XXX */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* See e_cosh.c for complete comments.
*
* Converted to long double by Bruce D. Evans.
*/
#include <float.h>
#ifdef __i386__
#include <ieeefp.h>
#endif
#include "fpmath.h"
#include "math.h"
#include "math_private.h"
#include "k_expl.h"
#if LDBL_MAX_EXP != 0x4000
/* We also require the usual expsign encoding. */
#error "Unsupported long double format"
#endif
#define BIAS (LDBL_MAX_EXP - 1)
static const volatile long double huge = 0x1p10000L, tiny = 0x1p-10000L;
#if LDBL_MANT_DIG == 64
/*
* Domain [-1, 1], range ~[-1.8211e-21, 1.8211e-21]:
* |cosh(x) - c(x)| < 2**-68.8
*/
static const union IEEEl2bits
C4u = LD80C(0xaaaaaaaaaaaaac78, -5, 4.16666666666666682297e-2L);
#define C4 C4u.e
static const double
C2 = 0.5,
C6 = 1.3888888888888616e-3, /* 0x16c16c16c16b99.0p-62 */
C8 = 2.4801587301767953e-5, /* 0x1a01a01a027061.0p-68 */
C10 = 2.7557319163300398e-7, /* 0x127e4fb6c9b55f.0p-74 */
C12 = 2.0876768371393075e-9, /* 0x11eed99406a3f4.0p-81 */
C14 = 1.1469537039374480e-11, /* 0x1938c67cd18c48.0p-89 */
C16 = 4.8473490896852041e-14; /* 0x1b49c429701e45.0p-97 */
#elif LDBL_MANT_DIG == 113
/*
* Domain [-1, 1], range ~[-2.3194e-37, 2.3194e-37]:
* |cosh(x) - c(x)| < 2**-121.69
*/
static const long double
C4 = 4.16666666666666666666666666666666225e-2L, /* 0x1555555555555555555555555554e.0p-117L */
C6 = 1.38888888888888888888888888889434831e-3L, /* 0x16c16c16c16c16c16c16c16c1dd7a.0p-122L */
C8 = 2.48015873015873015873015871870962089e-5L, /* 0x1a01a01a01a01a01a01a017af2756.0p-128L */
C10 = 2.75573192239858906525574318600800201e-7L, /* 0x127e4fb7789f5c72ef01c8a040640.0p-134L */
C12 = 2.08767569878680989791444691755468269e-9L, /* 0x11eed8eff8d897b543d0679607399.0p-141L */
C14= 1.14707455977297247387801189650495351e-11L, /* 0x193974a8c07c9d24ae169a7fa9b54.0p-149L */
C16 = 4.77947733238737883626416876486279985e-14L; /* 0x1ae7f3e733b814d4e1b90f5727fe4.0p-157L */
static const double
C2 = 0.5,
C18 = 1.5619206968597871e-16, /* 0x16827863b9900b.0p-105 */
C20 = 4.1103176218528049e-19, /* 0x1e542ba3d3c269.0p-114 */
C22 = 8.8967926401641701e-22, /* 0x10ce399542a014.0p-122 */
C24 = 1.6116681626523904e-24, /* 0x1f2c981d1f0cb7.0p-132 */
C26 = 2.5022374732804632e-27; /* 0x18c7ecf8b2c4a0.0p-141 */
#else
#error "Unsupported long double format"
#endif /* LDBL_MANT_DIG == 64 */
/* log(2**16385 - 0.5) rounded up: */
static const float
o_threshold = 1.13572168e4; /* 0xb174de.0p-10 */
long double
coshl(long double x)
{
long double hi,lo,x2,x4;
#if LDBL_MANT_DIG == 113
double dx2;
#endif
uint16_t ix;
GET_LDBL_EXPSIGN(ix,x);
ix &= 0x7fff;
/* x is INF or NaN */
if(ix>=0x7fff) return x*x;
ENTERI();
/* |x| < 1, return 1 or c(x) */
if(ix<0x3fff) {
if (ix<BIAS-(LDBL_MANT_DIG+1)/2) /* |x| < TINY */
RETURNI(1+tiny); /* cosh(tiny) = 1(+) with inexact */
x2 = x*x;
#if LDBL_MANT_DIG == 64
x4 = x2*x2;
RETURNI(((C16*x2 + C14)*x4 + (C12*x2 + C10))*(x4*x4*x2) +
((C8*x2 + C6)*x2 + C4)*x4 + C2*x2 + 1);
#elif LDBL_MANT_DIG == 113
dx2 = x2;
RETURNI((((((((((((C26*dx2 + C24)*dx2 + C22)*dx2 +
C20)*x2 + C18)*x2 +
C16)*x2 + C14)*x2 + C12)*x2 + C10)*x2 + C8)*x2 + C6)*x2 +
C4)*(x2*x2) + C2*x2 + 1);
#endif
}
/* |x| in [1, 64), return accurate exp(|x|)/2+1/exp(|x|)/2 */
if (ix < 0x4005) {
k_hexpl(fabsl(x), &hi, &lo);
RETURNI(lo + 0.25/(hi + lo) + hi);
}
/* |x| in [64, o_threshold], return correctly-overflowing exp(|x|)/2 */
if (fabsl(x) <= o_threshold)
RETURNI(hexpl(fabsl(x)));
/* |x| > o_threshold, cosh(x) overflow */
RETURNI(huge*huge);
}
|
