1+ /* ***
2+ DIAMOND protein sequence aligner
3+ Copyright (C) 2012-2026 Benjamin J. Buchfink
4+
5+ This program is free software: you can redistribute it and/or modify
6+ it under the terms of the GNU General Public License as published by
7+ the Free Software Foundation, either version 3 of the License, or
8+ (at your option) any later version.
9+
10+ This program is distributed in the hope that it will be useful,
11+ but WITHOUT ANY WARRANTY; without even the implied warranty of
12+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+ GNU General Public License for more details.
14+
15+ You should have received a copy of the GNU General Public License
16+ along with this program. If not, see <http://www.gnu.org/licenses/>.
17+ ****/
18+ // SPDX-License-Identifier: GPL-3.0-or-later
19+
20+ #pragma once
21+ #include < cstdint>
22+ #include < cstring>
23+
24+ namespace Math {
25+
26+ inline double from_bits (std::uint64_t b) {
27+ double d;
28+ std::memcpy (&d, &b, sizeof d);
29+ return d;
30+ }
31+
32+ inline std::uint64_t to_bits (double d) {
33+ std::uint64_t b;
34+ std::memcpy (&b, &d, sizeof d);
35+ return b;
36+ }
37+
38+ inline double quiet_nan () { return from_bits (0x7ff8000000000000ULL ); }
39+
40+ inline double exp (double x) {
41+ constexpr double LN2_HI = 6.93147180369123816490e-01 ;
42+ constexpr double LN2_LO = 1.90821492927058770002e-10 ;
43+ constexpr double INV_LN2 = 1.44269504088896338700e+00 ;
44+ constexpr double RND_SHIFT = 6755399441055744.0 ;
45+
46+ const std::uint64_t ux = to_bits (x);
47+ const std::uint64_t ax = ux & 0x7fffffffffffffffULL ;
48+
49+ if (ax > 0x7ff0000000000000ULL ) return quiet_nan ();
50+ if (ax == 0x7ff0000000000000ULL )
51+ return (ux >> 63 ) ? 0.0 : x;
52+ if (x >= 709.78271289338409 ) return from_bits (0x7ff0000000000000ULL );
53+ if (x <= -745.13321910194122 ) return 0.0 ;
54+
55+ double t = x * INV_LN2 + RND_SHIFT ;
56+ double kd = t - RND_SHIFT ;
57+ int k = static_cast <int >(kd);
58+ double r = x - kd * LN2_HI ;
59+ r = r - kd * LN2_LO ;
60+
61+ const double c2 = 5.00000000000000000000e-01 ;
62+ const double c3 = 1.66666666666666666667e-01 ;
63+ const double c4 = 4.16666666666666666667e-02 ;
64+ const double c5 = 8.33333333333333333333e-03 ;
65+ const double c6 = 1.38888888888888888889e-03 ;
66+ const double c7 = 1.98412698412698412698e-04 ;
67+ const double c8 = 2.48015873015873015873e-05 ;
68+ const double c9 = 2.75573192239858906526e-06 ;
69+ const double c10 = 2.75573192239858906526e-07 ;
70+ const double c11 = 2.50521083854417187751e-08 ;
71+ const double c12 = 2.08767569878680989792e-09 ;
72+ const double c13 = 1.60590438368216145994e-10 ;
73+
74+ double p = c13;
75+ p = c12 + r * p;
76+ p = c11 + r * p;
77+ p = c10 + r * p;
78+ p = c9 + r * p;
79+ p = c8 + r * p;
80+ p = c7 + r * p;
81+ p = c6 + r * p;
82+ p = c5 + r * p;
83+ p = c4 + r * p;
84+ p = c3 + r * p;
85+ p = c2 + r * p;
86+ double er = 1.0 + r + (r * r) * p;
87+
88+ if (k >= -1021 && k <= 1023 ) {
89+ return er * from_bits (static_cast <std::uint64_t >(1023 + k) << 52 );
90+ }
91+ if (k > 1023 ) {
92+ er *= from_bits (static_cast <std::uint64_t >(1023 + 1023 ) << 52 );
93+ k -= 1023 ;
94+ return er * from_bits (static_cast <std::uint64_t >(1023 + k) << 52 );
95+ }
96+ er *= from_bits (static_cast <std::uint64_t >(1023 - 1000 ) << 52 );
97+ k += 1000 ;
98+ return er * from_bits (static_cast <std::uint64_t >(1023 + k) << 52 );
99+ }
100+
101+ inline double log (double x) {
102+ constexpr double LN2_HI = 6.93147180369123816490e-01 ;
103+ constexpr double LN2_LO = 1.90821492927058770002e-10 ;
104+ constexpr double SQRT2 = 1.41421356237309514547e+00 ;
105+
106+ const double Lg1 = 6.666666666666735130e-01 ;
107+ const double Lg2 = 3.999999999940941908e-01 ;
108+ const double Lg3 = 2.857142874366239149e-01 ;
109+ const double Lg4 = 2.222219843214978396e-01 ;
110+ const double Lg5 = 1.818357216161805012e-01 ;
111+ const double Lg6 = 1.531383769920937332e-01 ;
112+ const double Lg7 = 1.479819860511658591e-01 ;
113+
114+ std::uint64_t u = to_bits (x);
115+
116+ if (u >> 63 ) {
117+ if ((u << 1 ) == 0 ) return from_bits (0xfff0000000000000ULL );
118+ return quiet_nan ();
119+ }
120+ if (u == 0 ) return from_bits (0xfff0000000000000ULL );
121+ if (u >= 0x7ff0000000000000ULL ) {
122+ if (u == 0x7ff0000000000000ULL ) return x;
123+ return quiet_nan ();
124+ }
125+
126+ int e = 0 ;
127+ if (u < 0x0010000000000000ULL ) {
128+ x *= 18014398509481984.0 ;
129+ u = to_bits (x);
130+ e -= 54 ;
131+ }
132+
133+ e += static_cast <int >(u >> 52 ) - 1023 ;
134+ u = (u & 0x000fffffffffffffULL ) | 0x3ff0000000000000ULL ;
135+ double m = from_bits (u);
136+ if (m > SQRT2 ) { m *= 0.5 ; e += 1 ; }
137+
138+ double f = m - 1.0 ;
139+ double s = f / (2.0 + f);
140+ double z = s * s;
141+ double w = z * z;
142+ double t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
143+ double t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
144+ double R = t2 + t1;
145+ double hfsq = 0.5 * f * f;
146+
147+ double dk = static_cast <double >(e);
148+ return dk * LN2_HI - ((hfsq - (s * (hfsq + R) + dk * LN2_LO )) - f);
149+ }
150+
151+ }
0 commit comments