1 // Copyright 2018 Ulf Adams
2 //
3 // The contents of this file may be used under the terms of the Apache License,
4 // Version 2.0.
5 //
6 // (See accompanying file LICENSE-Apache or copy at
7 // http://www.apache.org/licenses/LICENSE-2.0)
8 //
9 // Alternatively, the contents of this file may be used under the terms of
10 // the Boost Software License, Version 1.0.
11 // (See accompanying file LICENSE-Boost or copy at
12 // https://www.boost.org/LICENSE_1_0.txt)
13 //
14 // Unless required by applicable law or agreed to in writing, this software
15 // is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16 // KIND, either express or implied.
17
18 // Runtime compiler options:
19 // -DRYU_DEBUG Generate verbose debugging output to stdout.
20
21
22
23
24 #ifdef RYU_DEBUG
25 static char* s(uint128_t v) {
26 int len = decimalLength(v);
27 char* b = (char*) malloc((len + 1) * sizeof(char));
28 for (int i = 0; i < len; i++) {
29 const uint32_t c = (uint32_t) (v % 10);
30 v /= 10;
31 b[len - 1 - i] = (char) ('0' + c);
32 }
33 b[len] = 0;
34 return b;
35 }
36 #endif
37
38 #define ONE ((uint128_t) 1)
39
40 struct floating_decimal_128 generic_binary_to_decimal(
41 const uint128_t ieeeMantissa, const uint32_t ieeeExponent, const bool ieeeSign,
42 const uint32_t mantissaBits, const uint32_t exponentBits, const bool explicitLeadingBit) {
43 #ifdef RYU_DEBUG
44 printf("IN=");
45 for (int32_t bit = 127; bit >= 0; --bit) {
46 printf("%u", (uint32_t) ((bits >> bit) & 1));
47 }
48 printf("\n");
49 #endif
50
51 const uint32_t bias = (1u << (exponentBits - 1)) - 1;
52
53 if (ieeeExponent == 0 && ieeeMantissa == 0) {
54 struct floating_decimal_128 fd;
55 fd.mantissa = 0;
56 fd.exponent = 0;
57 fd.sign = ieeeSign;
58 return fd;
59 }
60 if (ieeeExponent == ((1u << exponentBits) - 1u)) {
61 struct floating_decimal_128 fd;
62 fd.mantissa = explicitLeadingBit ? ieeeMantissa & ((ONE << (mantissaBits - 1)) - 1) : ieeeMantissa;
63 fd.exponent = FD128_EXCEPTIONAL_EXPONENT;
64 fd.sign = ieeeSign;
65 return fd;
66 }
67
68 int32_t e2;
69 uint128_t m2;
70 // We subtract 2 in all cases so that the bounds computation has 2 additional bits.
71 if (explicitLeadingBit) {
72 // mantissaBits includes the explicit leading bit, so we need to correct for that here.
73 if (ieeeExponent == 0) {
74 e2 = 1 - bias - mantissaBits + 1 - 2;
75 } else {
76 e2 = ieeeExponent - bias - mantissaBits + 1 - 2;
77 }
78 m2 = ieeeMantissa;
79 } else {
80 if (ieeeExponent == 0) {
81 e2 = 1 - bias - mantissaBits - 2;
82 m2 = ieeeMantissa;
83 } else {
84 e2 = ieeeExponent - bias - mantissaBits - 2;
85 m2 = (ONE << mantissaBits) | ieeeMantissa;
86 }
87 }
88 const bool even = (m2 & 1) == 0;
89 const bool acceptBounds = even;
90
91 #ifdef RYU_DEBUG
92 printf("-> %s %s * 2^%d\n", ieeeSign ? "-" : "+", s(m2), e2 + 2);
93 #endif
94
95 // Step 2: Determine the interval of legal decimal representations.
96 const uint128_t mv = 4 * m2;
97 // Implicit bool -> int conversion. True is 1, false is 0.
98 const uint32_t mmShift =
99 (ieeeMantissa != (explicitLeadingBit ? ONE << (mantissaBits - 1) : 0))
100 || (ieeeExponent == 0);
101
102 // Step 3: Convert to a decimal power base using 128-bit arithmetic.
103 uint128_t vr, vp, vm;
104 int32_t e10;
105 bool vmIsTrailingZeros = false;
106 bool vrIsTrailingZeros = false;
107 if (e2 >= 0) {
108 // I tried special-casing q == 0, but there was no effect on performance.
109 // This expression is slightly faster than max(0, log10Pow2(e2) - 1).
110 const uint32_t q = log10Pow2(e2) - (e2 > 3);
111 e10 = q;
112 const int32_t k = FLOAT_128_POW5_INV_BITCOUNT + pow5bits(q) - 1;
113 const int32_t i = -e2 + q + k;
114 uint64_t pow5[4];
115 generic_computeInvPow5(q, pow5);
116 vr = mulShift(4 * m2, pow5, i);
117 vp = mulShift(4 * m2 + 2, pow5, i);
118 vm = mulShift(4 * m2 - 1 - mmShift, pow5, i);
119 #ifdef RYU_DEBUG
120 printf("%s * 2^%d / 10^%d\n", s(mv), e2, q);
121 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm));
122 #endif
123 // floor(log_5(2^128)) = 55, this is very conservative
124 if (q <= 55) {
125 // Only one of mp, mv, and mm can be a multiple of 5, if any.
126 if (mv % 5 == 0) {
127 vrIsTrailingZeros = multipleOfPowerOf5(mv, q - 1);
128 } else if (acceptBounds) {
129 // Same as min(e2 + (~mm & 1), pow5Factor(mm)) >= q
130 // <=> e2 + (~mm & 1) >= q && pow5Factor(mm) >= q
131 // <=> true && pow5Factor(mm) >= q, since e2 >= q.
132 vmIsTrailingZeros = multipleOfPowerOf5(mv - 1 - mmShift, q);
133 } else {
134 // Same as min(e2 + 1, pow5Factor(mp)) >= q.
135 vp -= multipleOfPowerOf5(mv + 2, q);
136 }
137 }
138 } else {
139 // This expression is slightly faster than max(0, log10Pow5(-e2) - 1).
140 const uint32_t q = log10Pow5(-e2) - (-e2 > 1);
141 e10 = q + e2;
142 const int32_t i = -e2 - q;
143 const int32_t k = pow5bits(i) - FLOAT_128_POW5_BITCOUNT;
144 const int32_t j = q - k;
145 uint64_t pow5[4];
146 generic_computePow5(i, pow5);
147 vr = mulShift(4 * m2, pow5, j);
148 vp = mulShift(4 * m2 + 2, pow5, j);
149 vm = mulShift(4 * m2 - 1 - mmShift, pow5, j);
150 #ifdef RYU_DEBUG
151 printf("%s * 5^%d / 10^%d\n", s(mv), -e2, q);
152 printf("%d %d %d %d\n", q, i, k, j);
153 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm));
154 #endif
155 if (q <= 1) {
156 // {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q trailing 0 bits.
157 // mv = 4 m2, so it always has at least two trailing 0 bits.
158 vrIsTrailingZeros = true;
159 if (acceptBounds) {
160 // mm = mv - 1 - mmShift, so it has 1 trailing 0 bit iff mmShift == 1.
161 vmIsTrailingZeros = mmShift == 1;
162 } else {
163 // mp = mv + 2, so it always has at least one trailing 0 bit.
164 --vp;
165 }
166 } else if (q < 127) { // TODO(ulfjack): Use a tighter bound here.
167 // We need to compute min(ntz(mv), pow5Factor(mv) - e2) >= q-1
168 // <=> ntz(mv) >= q-1 && pow5Factor(mv) - e2 >= q-1
169 // <=> ntz(mv) >= q-1 (e2 is negative and -e2 >= q)
170 // <=> (mv & ((1 << (q-1)) - 1)) == 0
171 // We also need to make sure that the left shift does not overflow.
172 vrIsTrailingZeros = multipleOfPowerOf2(mv, q - 1);
173 #ifdef RYU_DEBUG
174 printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false");
175 #endif
176 }
177 }
178 #ifdef RYU_DEBUG
179 printf("e10=%d\n", e10);
180 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm));
181 printf("vm is trailing zeros=%s\n", vmIsTrailingZeros ? "true" : "false");
182 printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false");
183 #endif
184
185 // Step 4: Find the shortest decimal representation in the interval of legal representations.
186 uint32_t removed = 0;
187 uint8_t lastRemovedDigit = 0;
188 uint128_t output;
189
190 while (vp / 10 > vm / 10) {
191 vmIsTrailingZeros &= vm % 10 == 0;
192 vrIsTrailingZeros &= lastRemovedDigit == 0;
193 lastRemovedDigit = (uint8_t) (vr % 10);
194 vr /= 10;
195 vp /= 10;
196 vm /= 10;
197 ++removed;
198 }
199 #ifdef RYU_DEBUG
200 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm));
201 printf("d-10=%s\n", vmIsTrailingZeros ? "true" : "false");
202 #endif
203 if (vmIsTrailingZeros) {
204 while (vm % 10 == 0) {
205 vrIsTrailingZeros &= lastRemovedDigit == 0;
206 lastRemovedDigit = (uint8_t) (vr % 10);
207 vr /= 10;
208 vp /= 10;
209 vm /= 10;
210 ++removed;
211 }
212 }
213 #ifdef RYU_DEBUG
214 printf("%s %d\n", s(vr), lastRemovedDigit);
215 printf("vr is trailing zeros=%s\n", vrIsTrailingZeros ? "true" : "false");
216 #endif
217 if (vrIsTrailingZeros && (lastRemovedDigit == 5) && (vr % 2 == 0)) {
218 // Round even if the exact numbers is .....50..0.
219 lastRemovedDigit = 4;
220 }
221 // We need to take vr+1 if vr is outside bounds or we need to round up.
222 output = vr +
223 ((vr == vm && (!acceptBounds || !vmIsTrailingZeros)) || (lastRemovedDigit >= 5));
224 const int32_t exp = e10 + removed;
225
226 #ifdef RYU_DEBUG
227 printf("V+=%s\nV =%s\nV-=%s\n", s(vp), s(vr), s(vm));
228 printf("O=%s\n", s(output));
229 printf("EXP=%d\n", exp);
230 #endif
231
232 struct floating_decimal_128 fd;
233 fd.mantissa = output;
234 fd.exponent = exp;
235 fd.sign = ieeeSign;
236 return fd;
237 }
238
239 static inline int copy_special_str(char * const result, const struct floating_decimal_128 fd) {
240 if (fd.mantissa) {
241 memcpy(result, "NaN", 3);
242 return 3;
243 }
244 if (fd.sign) {
245 result[0] = '-';
246 }
247 memcpy(result + fd.sign, "Infinity", 8);
248 return fd.sign + 8;
249 }
250
251 int generic_to_chars(const struct floating_decimal_128 v, char* const result) {
252 if (v.exponent == FD128_EXCEPTIONAL_EXPONENT) {
253 return copy_special_str(result, v);
254 }
255
256 // Step 5: Print the decimal representation.
257 int index = 0;
258 if (v.sign) {
259 result[index++] = '-';
260 }
261
262 uint128_t output = v.mantissa;
263 const uint32_t olength = decimalLength(output);
264
265 #ifdef RYU_DEBUG
266 printf("DIGITS=%s\n", s(v.mantissa));
267 printf("OLEN=%u\n", olength);
268 printf("EXP=%u\n", v.exponent + olength);
269 #endif
270
271 for (uint32_t i = 0; i < olength - 1; ++i) {
272 const uint32_t c = (uint32_t) (output % 10);
273 output /= 10;
274 result[index + olength - i] = (char) ('0' + c);
275 }
276 result[index] = '0' + (uint32_t) (output % 10); // output should be < 10 by now.
277
278 // Print decimal point if needed.
279 if (olength > 1) {
280 result[index + 1] = '.';
281 index += olength + 1;
282 } else {
283 ++index;
284 }
285
286 // Print the exponent.
287 result[index++] = 'e';
288 int32_t exp = v.exponent + olength - 1;
289 if (exp < 0) {
290 result[index++] = '-';
291 exp = -exp;
292 } else
293 result[index++] = '+';
294
295 uint32_t elength = decimalLength(exp);
296 if (elength == 1)
297 ++elength;
298 for (uint32_t i = 0; i < elength; ++i) {
299 const uint32_t c = exp % 10;
300 exp /= 10;
301 result[index + elength - 1 - i] = (char) ('0' + c);
302 }
303 index += elength;
304 return index;
305 }