grisu3_print.h (9143B)
1 /* 2 * Copyright (c) 2016 Mikkel F. Jørgensen, dvide.com 3 * Copyright author of MathGeoLib (https://github.com/juj) 4 * 5 * Licensed under the Apache License, Version 2.0 (the "License"); 6 * you may not use this file except in compliance with the License. 7 * You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. http://www.apache.org/licenses/LICENSE-2.0 16 */ 17 18 /* 19 * Extracted from MathGeoLib. 20 * 21 * mikkelfj: 22 * - Fixed final output when printing single digit negative exponent to 23 * have leading zero (important for JSON). 24 * - Changed formatting to prefer 0.012 over 1.2-e-2. 25 * 26 * Large portions of the original grisu3.c file has been moved to 27 * grisu3_math.h, the rest is placed here. 28 * 29 * See also comments in grisu3_math.h. 30 * 31 * MatGeoLib grisu3.c comment: 32 * 33 * This file is part of an implementation of the "grisu3" double to string 34 * conversion algorithm described in the research paper 35 * 36 * "Printing Floating-Point Numbers Quickly And Accurately with Integers" 37 * by Florian Loitsch, available at 38 * http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf 39 */ 40 41 #ifndef GRISU3_PRINT_H 42 #define GRISU3_PRINT_H 43 44 #ifdef __cplusplus 45 extern "C" { 46 #endif 47 48 #include <stdio.h> /* sprintf, only needed for fallback printing */ 49 #include <assert.h> /* assert */ 50 51 #include "grisu3_math.h" 52 53 /* 54 * The lightweight "portable" C library recognizes grisu3 support if 55 * included first. 56 */ 57 #define grisu3_print_double_is_defined 1 58 59 /* 60 * Not sure we have an exact definition, but we get up to 23 61 * emperically. There is some math ensuring it does not go awol though, 62 * like 18 digits + exponent or so. 63 * This max should be safe size buffer for printing, including zero term. 64 */ 65 #define GRISU3_PRINT_MAX 30 66 67 static int grisu3_round_weed(char *buffer, int len, uint64_t wp_W, uint64_t delta, uint64_t rest, uint64_t ten_kappa, uint64_t ulp) 68 { 69 uint64_t wp_Wup = wp_W - ulp; 70 uint64_t wp_Wdown = wp_W + ulp; 71 while(rest < wp_Wup && delta - rest >= ten_kappa 72 && (rest + ten_kappa < wp_Wup || wp_Wup - rest >= rest + ten_kappa - wp_Wup)) 73 { 74 --buffer[len-1]; 75 rest += ten_kappa; 76 } 77 if (rest < wp_Wdown && delta - rest >= ten_kappa 78 && (rest + ten_kappa < wp_Wdown || wp_Wdown - rest > rest + ten_kappa - wp_Wdown)) 79 return 0; 80 81 return 2*ulp <= rest && rest <= delta - 4*ulp; 82 } 83 84 static int grisu3_digit_gen(grisu3_diy_fp_t low, grisu3_diy_fp_t w, grisu3_diy_fp_t high, char *buffer, int *length, int *kappa) 85 { 86 uint64_t unit = 1; 87 grisu3_diy_fp_t too_low = { low.f - unit, low.e }; 88 grisu3_diy_fp_t too_high = { high.f + unit, high.e }; 89 grisu3_diy_fp_t unsafe_interval = grisu3_diy_fp_minus(too_high, too_low); 90 grisu3_diy_fp_t one = { 1ULL << -w.e, w.e }; 91 uint32_t p1 = (uint32_t)(too_high.f >> -one.e); 92 uint64_t p2 = too_high.f & (one.f - 1); 93 uint32_t div; 94 *kappa = grisu3_largest_pow10(p1, GRISU3_DIY_FP_FRACT_SIZE + one.e, &div); 95 *length = 0; 96 97 while(*kappa > 0) 98 { 99 uint64_t rest; 100 char digit = (char)(p1 / div); 101 buffer[*length] = '0' + digit; 102 ++*length; 103 p1 %= div; 104 --*kappa; 105 rest = ((uint64_t)p1 << -one.e) + p2; 106 if (rest < unsafe_interval.f) return grisu3_round_weed(buffer, *length, grisu3_diy_fp_minus(too_high, w).f, unsafe_interval.f, rest, (uint64_t)div << -one.e, unit); 107 div /= 10; 108 } 109 110 for(;;) 111 { 112 char digit; 113 p2 *= 10; 114 unit *= 10; 115 unsafe_interval.f *= 10; 116 /* Integer division by one. */ 117 digit = (char)(p2 >> -one.e); 118 buffer[*length] = '0' + digit; 119 ++*length; 120 p2 &= one.f - 1; /* Modulo by one. */ 121 --*kappa; 122 if (p2 < unsafe_interval.f) return grisu3_round_weed(buffer, *length, grisu3_diy_fp_minus(too_high, w).f * unit, unsafe_interval.f, p2, one.f, unit); 123 } 124 } 125 126 static int grisu3(double v, char *buffer, int *length, int *d_exp) 127 { 128 int mk, kappa, success; 129 grisu3_diy_fp_t dfp = grisu3_cast_diy_fp_from_double(v); 130 grisu3_diy_fp_t w = grisu3_diy_fp_normalize(dfp); 131 132 /* normalize boundaries */ 133 grisu3_diy_fp_t t = { (dfp.f << 1) + 1, dfp.e - 1 }; 134 grisu3_diy_fp_t b_plus = grisu3_diy_fp_normalize(t); 135 grisu3_diy_fp_t b_minus; 136 grisu3_diy_fp_t c_mk; /* Cached power of ten: 10^-k */ 137 uint64_t u64 = grisu3_cast_uint64_from_double(v); 138 assert(v > 0 && v <= 1.7976931348623157e308); /* Grisu only handles strictly positive finite numbers. */ 139 if (!(u64 & GRISU3_D64_FRACT_MASK) && (u64 & GRISU3_D64_EXP_MASK) != 0) { b_minus.f = (dfp.f << 2) - 1; b_minus.e = dfp.e - 2;} /* lower boundary is closer? */ 140 else { b_minus.f = (dfp.f << 1) - 1; b_minus.e = dfp.e - 1; } 141 b_minus.f = b_minus.f << (b_minus.e - b_plus.e); 142 b_minus.e = b_plus.e; 143 144 mk = grisu3_diy_fp_cached_pow(GRISU3_MIN_TARGET_EXP - GRISU3_DIY_FP_FRACT_SIZE - w.e, &c_mk); 145 146 w = grisu3_diy_fp_multiply(w, c_mk); 147 b_minus = grisu3_diy_fp_multiply(b_minus, c_mk); 148 b_plus = grisu3_diy_fp_multiply(b_plus, c_mk); 149 150 success = grisu3_digit_gen(b_minus, w, b_plus, buffer, length, &kappa); 151 *d_exp = kappa - mk; 152 return success; 153 } 154 155 static int grisu3_i_to_str(int val, char *str) 156 { 157 int len, i; 158 char *s; 159 char *begin = str; 160 if (val < 0) { *str++ = '-'; val = -val; } 161 s = str; 162 163 for(;;) 164 { 165 int ni = val / 10; 166 int digit = val - ni*10; 167 *s++ = (char)('0' + digit); 168 if (ni == 0) 169 break; 170 val = ni; 171 } 172 *s = '\0'; 173 len = (int)(s - str); 174 for(i = 0; i < len/2; ++i) 175 { 176 char ch = str[i]; 177 str[i] = str[len-1-i]; 178 str[len-1-i] = ch; 179 } 180 181 return (int)(s - begin); 182 } 183 184 static int grisu3_print_nan(uint64_t v, char *dst) 185 { 186 static char hexdigits[16] = "0123456789ABCDEF"; 187 int i = 0; 188 189 dst[0] = 'N'; 190 dst[1] = 'a'; 191 dst[2] = 'N'; 192 dst[3] = '('; 193 dst[20] = ')'; 194 dst[21] = '\0'; 195 dst += 4; 196 for (i = 15; i >= 0; --i) { 197 dst[i] = hexdigits[v & 0x0F]; 198 v >>= 4; 199 } 200 return 21; 201 } 202 203 static int grisu3_print_double(double v, char *dst) 204 { 205 int d_exp, len, success, decimals, i; 206 uint64_t u64 = grisu3_cast_uint64_from_double(v); 207 char *s2 = dst; 208 assert(dst); 209 210 /* Prehandle NaNs */ 211 if ((u64 << 1) > 0xFFE0000000000000ULL) return grisu3_print_nan(u64, dst); 212 /* Prehandle negative values. */ 213 if ((u64 & GRISU3_D64_SIGN) != 0) { *s2++ = '-'; v = -v; u64 ^= GRISU3_D64_SIGN; } 214 /* Prehandle zero. */ 215 if (!u64) { *s2++ = '0'; *s2 = '\0'; return (int)(s2 - dst); } 216 /* Prehandle infinity. */ 217 if (u64 == GRISU3_D64_EXP_MASK) { *s2++ = 'i'; *s2++ = 'n'; *s2++ = 'f'; *s2 = '\0'; return (int)(s2 - dst); } 218 219 success = grisu3(v, s2, &len, &d_exp); 220 /* If grisu3 was not able to convert the number to a string, then use old sprintf (suboptimal). */ 221 if (!success) return sprintf(s2, "%.17g", v) + (int)(s2 - dst); 222 223 /* We now have an integer string of form "151324135" and a base-10 exponent for that number. */ 224 /* Next, decide the best presentation for that string by whether to use a decimal point, or the scientific exponent notation 'e'. */ 225 /* We don't pick the absolute shortest representation, but pick a balance between readability and shortness, e.g. */ 226 /* 1.545056189557677e-308 could be represented in a shorter form */ 227 /* 1545056189557677e-323 but that would be somewhat unreadable. */ 228 decimals = GRISU3_MIN(-d_exp, GRISU3_MAX(1, len-1)); 229 230 /* mikkelfj: 231 * fix zero prefix .1 => 0.1, important for JSON export. 232 * prefer unscientific notation at same length: 233 * -1.2345e-4 over -1.00012345, 234 * -1.0012345 over -1.2345e-3 235 */ 236 if (d_exp < 0 && (len + d_exp) > -3 && len <= -d_exp) 237 { 238 /* mikkelfj: fix zero prefix .1 => 0.1, and short exponents 1.3e-2 => 0.013. */ 239 memmove(s2 + 2 - d_exp - len, s2, (size_t)len); 240 s2[0] = '0'; 241 s2[1] = '.'; 242 for (i = 2; i < 2-d_exp-len; ++i) s2[i] = '0'; 243 len += i; 244 } 245 else if (d_exp < 0 && len > 1) /* Add decimal point? */ 246 { 247 for(i = 0; i < decimals; ++i) s2[len-i] = s2[len-i-1]; 248 s2[len++ - decimals] = '.'; 249 d_exp += decimals; 250 /* Need scientific notation as well? */ 251 if (d_exp != 0) { s2[len++] = 'e'; len += grisu3_i_to_str(d_exp, s2+len); } 252 } 253 /* Add scientific notation? */ 254 else if (d_exp < 0 || d_exp > 2) { s2[len++] = 'e'; len += grisu3_i_to_str(d_exp, s2+len); } 255 /* Add zeroes instead of scientific notation? */ 256 else if (d_exp > 0) { while(d_exp-- > 0) s2[len++] = '0'; } 257 s2[len] = '\0'; /* grisu3 doesn't null terminate, so ensure termination. */ 258 return (int)(s2+len-dst); 259 } 260 261 #ifdef __cplusplus 262 } 263 #endif 264 265 #endif /* GRISU3_PRINT_H */