1 /*
2 * Copyright (C) Alexander Borisov
3 *
4 * Based on nxt_diyfp.h from NGINX NJS project
5 *
6 * Copyright (C) Dmitry Volyntsev
7 * Copyright (C) NGINX, Inc.
8 *
9 * An internal diy_fp implementation.
10 * For details, see Loitsch, Florian. "Printing floating-point numbers quickly
11 * and accurately with integers." ACM Sigplan Notices 45.6 (2010): 233-243.
12 */
13
14 #ifndef LEXBOR_DIYFP_H
15 #define LEXBOR_DIYFP_H
16
17 #ifdef __cplusplus
18 extern "C" {
19 #endif
20
21 #include "lexbor/core/base.h"
22
23 #include <math.h>
24
25
26 #ifdef __cplusplus
27 #define lexbor_diyfp(_s, _e) { .significand = (_s), .exp = (int) (_e) }
28 #else
29 #define lexbor_diyfp(_s, _e) (lexbor_diyfp_t) \
30 { .significand = (_s), .exp = (int) (_e) }
31 #endif
32 #define lexbor_uint64_hl(h, l) (((uint64_t) (h) << 32) + (l))
33
34
35 #define LEXBOR_DBL_SIGNIFICAND_SIZE 52
36 #define LEXBOR_DBL_EXPONENT_BIAS (0x3FF + LEXBOR_DBL_SIGNIFICAND_SIZE)
37 #define LEXBOR_DBL_EXPONENT_MIN (-LEXBOR_DBL_EXPONENT_BIAS)
38 #define LEXBOR_DBL_EXPONENT_MAX (0x7FF - LEXBOR_DBL_EXPONENT_BIAS)
39 #define LEXBOR_DBL_EXPONENT_DENORMAL (-LEXBOR_DBL_EXPONENT_BIAS + 1)
40
41 #define LEXBOR_DBL_SIGNIFICAND_MASK lexbor_uint64_hl(0x000FFFFF, 0xFFFFFFFF)
42 #define LEXBOR_DBL_HIDDEN_BIT lexbor_uint64_hl(0x00100000, 0x00000000)
43 #define LEXBOR_DBL_EXPONENT_MASK lexbor_uint64_hl(0x7FF00000, 0x00000000)
44
45 #define LEXBOR_DIYFP_SIGNIFICAND_SIZE 64
46
47 #define LEXBOR_SIGNIFICAND_SIZE 53
48 #define LEXBOR_SIGNIFICAND_SHIFT (LEXBOR_DIYFP_SIGNIFICAND_SIZE \
49 - LEXBOR_DBL_SIGNIFICAND_SIZE)
50
51 #define LEXBOR_DECIMAL_EXPONENT_OFF 348
52 #define LEXBOR_DECIMAL_EXPONENT_MIN (-348)
53 #define LEXBOR_DECIMAL_EXPONENT_MAX 340
54 #define LEXBOR_DECIMAL_EXPONENT_DIST 8
55
56
57 typedef struct {
58 uint64_t significand;
59 int exp;
60 }
61 lexbor_diyfp_t;
62
63
64 LXB_API lexbor_diyfp_t
65 lexbor_cached_power_dec(int exp, int *dec_exp);
66
67 LXB_API lexbor_diyfp_t
68 lexbor_cached_power_bin(int exp, int *dec_exp);
69
70
71 /*
72 * Inline functions
73 */
74 #if (LEXBOR_HAVE_BUILTIN_CLZLL)
75 #define nxt_leading_zeros64(x) (((x) == 0) ? 64 : __builtin_clzll(x))
76
77 #else
78
79 lxb_inline uint64_t
lexbor_diyfp_leading_zeros64(uint64_t x)80 lexbor_diyfp_leading_zeros64(uint64_t x)
81 {
82 uint64_t n;
83
84 if (x == 0) {
85 return 64;
86 }
87
88 n = 0;
89
90 while ((x & 0x8000000000000000) == 0) {
91 n++;
92 x <<= 1;
93 }
94
95 return n;
96 }
97
98 #endif
99
100 lxb_inline lexbor_diyfp_t
lexbor_diyfp_from_d2(double d)101 lexbor_diyfp_from_d2(double d)
102 {
103 int biased_exp;
104 uint64_t significand;
105 lexbor_diyfp_t r;
106
107 union {
108 double d;
109 uint64_t u64;
110 }
111 u;
112
113 u.d = d;
114
115 biased_exp = (u.u64 & LEXBOR_DBL_EXPONENT_MASK)
116 >> LEXBOR_DBL_SIGNIFICAND_SIZE;
117 significand = u.u64 & LEXBOR_DBL_SIGNIFICAND_MASK;
118
119 if (biased_exp != 0) {
120 r.significand = significand + LEXBOR_DBL_HIDDEN_BIT;
121 r.exp = biased_exp - LEXBOR_DBL_EXPONENT_BIAS;
122 }
123 else {
124 r.significand = significand;
125 r.exp = LEXBOR_DBL_EXPONENT_MIN + 1;
126 }
127
128 return r;
129 }
130
131 lxb_inline double
lexbor_diyfp_2d(lexbor_diyfp_t v)132 lexbor_diyfp_2d(lexbor_diyfp_t v)
133 {
134 int exp;
135 uint64_t significand, biased_exp;
136
137 union {
138 double d;
139 uint64_t u64;
140 }
141 u;
142
143 exp = v.exp;
144 significand = v.significand;
145
146 while (significand > LEXBOR_DBL_HIDDEN_BIT + LEXBOR_DBL_SIGNIFICAND_MASK) {
147 significand >>= 1;
148 exp++;
149 }
150
151 if (exp >= LEXBOR_DBL_EXPONENT_MAX) {
152 return INFINITY;
153 }
154
155 if (exp < LEXBOR_DBL_EXPONENT_DENORMAL) {
156 return 0.0;
157 }
158
159 while (exp > LEXBOR_DBL_EXPONENT_DENORMAL
160 && (significand & LEXBOR_DBL_HIDDEN_BIT) == 0)
161 {
162 significand <<= 1;
163 exp--;
164 }
165
166 if (exp == LEXBOR_DBL_EXPONENT_DENORMAL
167 && (significand & LEXBOR_DBL_HIDDEN_BIT) == 0)
168 {
169 biased_exp = 0;
170
171 } else {
172 biased_exp = (uint64_t) (exp + LEXBOR_DBL_EXPONENT_BIAS);
173 }
174
175 u.u64 = (significand & LEXBOR_DBL_SIGNIFICAND_MASK)
176 | (biased_exp << LEXBOR_DBL_SIGNIFICAND_SIZE);
177
178 return u.d;
179 }
180
181 lxb_inline lexbor_diyfp_t
lexbor_diyfp_shift_left(lexbor_diyfp_t v,unsigned shift)182 lexbor_diyfp_shift_left(lexbor_diyfp_t v, unsigned shift)
183 {
184 return lexbor_diyfp(v.significand << shift, v.exp - shift);
185 }
186
187 lxb_inline lexbor_diyfp_t
lexbor_diyfp_shift_right(lexbor_diyfp_t v,unsigned shift)188 lexbor_diyfp_shift_right(lexbor_diyfp_t v, unsigned shift)
189 {
190 return lexbor_diyfp(v.significand >> shift, v.exp + shift);
191 }
192
193 lxb_inline lexbor_diyfp_t
lexbor_diyfp_sub(lexbor_diyfp_t lhs,lexbor_diyfp_t rhs)194 lexbor_diyfp_sub(lexbor_diyfp_t lhs, lexbor_diyfp_t rhs)
195 {
196 return lexbor_diyfp(lhs.significand - rhs.significand, lhs.exp);
197 }
198
199 lxb_inline lexbor_diyfp_t
lexbor_diyfp_mul(lexbor_diyfp_t lhs,lexbor_diyfp_t rhs)200 lexbor_diyfp_mul(lexbor_diyfp_t lhs, lexbor_diyfp_t rhs)
201 {
202 #if (LEXBOR_HAVE_UNSIGNED_INT128)
203
204 uint64_t l, h;
205 lxb_uint128_t u128;
206
207 u128 = (lxb_uint128_t) (lhs.significand)
208 * (lxb_uint128_t) (rhs.significand);
209
210 h = u128 >> 64;
211 l = (uint64_t) u128;
212
213 /* rounding. */
214
215 if (l & ((uint64_t) 1 << 63)) {
216 h++;
217 }
218
219 return lexbor_diyfp(h, lhs.exp + rhs.exp + 64);
220
221 #else
222
223 uint64_t a, b, c, d, ac, bc, ad, bd, tmp;
224
225 a = lhs.significand >> 32;
226 b = lhs.significand & 0xffffffff;
227 c = rhs.significand >> 32;
228 d = rhs.significand & 0xffffffff;
229
230 ac = a * c;
231 bc = b * c;
232 ad = a * d;
233 bd = b * d;
234
235 tmp = (bd >> 32) + (ad & 0xffffffff) + (bc & 0xffffffff);
236
237 /* mult_round. */
238
239 tmp += 1U << 31;
240
241 return lexbor_diyfp(ac + (ad >> 32) + (bc >> 32) + (tmp >> 32),
242 lhs.exp + rhs.exp + 64);
243 #endif
244 }
245
246 lxb_inline lexbor_diyfp_t
lexbor_diyfp_normalize(lexbor_diyfp_t v)247 lexbor_diyfp_normalize(lexbor_diyfp_t v)
248 {
249 return lexbor_diyfp_shift_left(v,
250 (unsigned) lexbor_diyfp_leading_zeros64(v.significand));
251 }
252
253
254 #ifdef __cplusplus
255 } /* extern "C" */
256 #endif
257
258 #endif /* LEXBOR_DIYFP_H */
259
