dec.c (10901B)
1 /* dec.c - aritmetica decimal de 18 digitos (ver dec.h). 2 * 3 * Todo con uint64: las mantisas tienen 18 digitos (< 10^18 < 2^63), el producto 4 * se arma en "limbs" de 10^9 y la division es la larga de la escuela. El 5 * redondeo es a la mitad par usando la informacion de lo descartado (ri). */ 6 #include "dec.h" 7 #include "str.h" 8 9 #define E9 1000000000ull 10 #define E17 100000000000000000ull 11 #define E18 1000000000000000000ull 12 13 const Dec DEC_ZERO = { 0, 0, 0, DE_OK }; 14 const Dec DEC_ONE = { E17, 0, 0, DE_OK }; 15 16 static const uint64_t P10[20] = { 17 1ull, 10ull, 100ull, 1000ull, 10000ull, 100000ull, 1000000ull, 10000000ull, 18 100000000ull, 1000000000ull, 10000000000ull, 100000000000ull, 1000000000000ull, 19 10000000000000ull, 100000000000000ull, 1000000000000000ull, 10000000000000000ull, 20 100000000000000000ull, 1000000000000000000ull, 10000000000000000000ull 21 }; 22 23 Dec dec_err(void) { Dec d = { 0, 0, 0, DE_MATH }; return d; } 24 25 /* Informacion de redondeo de un resto r de una unidad u: 26 * 0 exacto, 1 menos de media, 2 justo media, 3 mas de media. */ 27 static int rinfo(uint64_t r, uint64_t u) 28 { 29 if (!r) return 0; 30 uint64_t h = u / 2; /* u siempre es par (potencia de 10 >= 10) */ 31 return r < h ? 1 : r == h ? 2 : 3; 32 } 33 34 /* Arma un Dec con valor ~ (m + lo descartado) * 10^p y lo redondea a 18 digitos. */ 35 static Dec pack(int neg, uint64_t m, int p, int ri) 36 { 37 Dec d = { 0, 0, (uint8_t)neg, DE_OK }; 38 if (!m) { d.neg = 0; return d; } 39 while (m >= E18) { 40 int dg = (int)(m % 10); 41 m /= 10; p++; 42 if (dg > 5) ri = 3; 43 else if (dg == 5) ri = ri ? 3 : 2; 44 else if (dg > 0) ri = 1; 45 else ri = ri ? 1 : 0; 46 } 47 while (m < E17) { m *= 10; p--; } /* quien llama garantiza 18 digitos si ri != 0 */ 48 if (ri == 3 || (ri == 2 && (m & 1))) { 49 m++; 50 if (m == E18) { m = E17; p++; } 51 } 52 int e = p + 17; 53 if (e > DEC_EMAX) return dec_err(); 54 if (e < -DEC_EMAX - 1) return DEC_ZERO; /* subdesborde: cero */ 55 d.m = m; d.e = (int16_t)e; 56 return d; 57 } 58 59 Dec dec_int(int64_t v) 60 { 61 uint64_t m = v < 0 ? (uint64_t)0 - (uint64_t)v : (uint64_t)v; 62 return pack(v < 0, m, 0, 0); 63 } 64 65 Dec dec_neg(Dec a) { if (a.m) a.neg ^= 1; return a; } 66 Dec dec_abs(Dec a) { a.neg = 0; return a; } 67 68 int dec_cmp_abs(Dec a, Dec b) 69 { 70 if (!a.m || !b.m) return a.m ? 1 : b.m ? -1 : 0; 71 if (a.e != b.e) return a.e > b.e ? 1 : -1; 72 return a.m > b.m ? 1 : a.m < b.m ? -1 : 0; 73 } 74 75 int dec_cmp(Dec a, Dec b) 76 { 77 int sa = !a.m ? 0 : a.neg ? -1 : 1, sb = !b.m ? 0 : b.neg ? -1 : 1; 78 if (sa != sb) return sa > sb ? 1 : -1; 79 int c = dec_cmp_abs(a, b); 80 return sa < 0 ? -c : c; 81 } 82 83 Dec dec_add(Dec a, Dec b) 84 { 85 if (a.err || b.err) return dec_err(); 86 if (!b.m) return a; 87 if (!a.m) return b; 88 if (dec_cmp_abs(a, b) < 0) { Dec t = a; a = b; b = t; } 89 int d = a.e - b.e; 90 if (a.neg == b.neg) { 91 if (d > 19) return a; 92 uint64_t B = d ? b.m / P10[d] : b.m; 93 int rb = d ? rinfo(b.m % P10[d], P10[d]) : 0; 94 return pack(a.neg, a.m + B, a.e - 17, rb); 95 } 96 /* signos distintos, |a| >= |b|: un digito extra de guarda */ 97 uint64_t A = a.m * 10, B; 98 int rb = 0; 99 if (d == 0) B = b.m * 10; 100 else if (d == 1) B = b.m; 101 else if (d <= 19) { B = b.m / P10[d - 1]; rb = rinfo(b.m % P10[d - 1], P10[d - 1]); } 102 else { B = 0; rb = 1; } 103 uint64_t D = A - B; 104 if (rb) { D--; rb = 4 - rb; } /* restar una fraccion: queda D-1 + (1-frac) */ 105 if (rb == 4) rb = 0; 106 return pack(a.neg, D, a.e - 18, rb); 107 } 108 109 Dec dec_sub(Dec a, Dec b) { return dec_add(a, dec_neg(b)); } 110 111 Dec dec_mul(Dec a, Dec b) 112 { 113 if (a.err || b.err) return dec_err(); 114 if (!a.m || !b.m) return DEC_ZERO; 115 uint64_t a1 = a.m / E9, a0 = a.m % E9, b1 = b.m / E9, b0 = b.m % E9; 116 uint64_t p0 = a0 * b0, p1 = a1 * b0 + a0 * b1, p2 = a1 * b1; 117 uint64_t L0 = p0 % E9, t = p1 + p0 / E9; 118 uint64_t L1 = t % E9; 119 t = p2 + t / E9; /* t = L3*10^9 + L2: los 18 digitos altos */ 120 uint64_t low = L1 * E9 + L0; /* producto = t*10^18 + low */ 121 int neg = a.neg ^ b.neg, p = a.e + b.e - 34; 122 if (t >= E17) return pack(neg, t, p + 18, rinfo(low, E18)); 123 return pack(neg, t * 10 + low / E17, p + 17, rinfo(low % E17, E17)); 124 } 125 126 Dec dec_div(Dec a, Dec b) 127 { 128 if (a.err || b.err || !b.m) return dec_err(); 129 if (!a.m) return DEC_ZERO; 130 uint64_t q = a.m / b.m, r = a.m % b.m; 131 for (int i = 0; i < 18; i++) { 132 r *= 10; 133 q = q * 10 + r / b.m; 134 r %= b.m; 135 } 136 int ri = !r ? 0 : 2 * r < b.m ? 1 : 2 * r == b.m ? 2 : 3; 137 return pack(a.neg ^ b.neg, q, a.e - b.e - 18, ri); 138 } 139 140 Dec dec_mul_int(Dec a, int64_t k) { return dec_mul(a, dec_int(k)); } 141 Dec dec_div_int(Dec a, int64_t k) { return dec_div(a, dec_int(k)); } 142 143 Dec dec_scale10(Dec a, int k) 144 { 145 if (a.err || !a.m) return a; 146 int e = a.e + k; 147 if (e > DEC_EMAX) return dec_err(); 148 if (e < -DEC_EMAX - 1) return DEC_ZERO; 149 a.e = (int16_t)e; 150 return a; 151 } 152 153 static uint64_t isqrt64(uint64_t v) 154 { 155 uint64_t r = 0, bit = 1ull << 62; 156 while (bit > v) bit >>= 2; 157 while (bit) { 158 if (v >= r + bit) { v -= r + bit; r = (r >> 1) + bit; } 159 else r >>= 1; 160 bit >>= 2; 161 } 162 return r; 163 } 164 165 Dec dec_sqrt(Dec a) 166 { 167 if (a.err || a.neg) return a.m ? dec_err() : DEC_ZERO; 168 if (!a.m) return DEC_ZERO; 169 /* a = M * 10^(2k) con M de 17 o 18 digitos: sqrt(a) ~ isqrt(M) * 10^k (9 digitos) */ 170 int p = a.e - 17; 171 uint64_t M = a.m; 172 if (p & 1) { M /= 10; p++; } 173 Dec x = pack(0, isqrt64(M), p / 2, 0); 174 /* Newton: cada paso duplica los digitos (9 -> 18 -> 36) */ 175 for (int i = 0; i < 3; i++) { 176 Dec y = dec_div(a, x); 177 x = dec_div_int(dec_add(x, y), 2); 178 } 179 return x; 180 } 181 182 Dec dec_pow_int(Dec a, int64_t n) 183 { 184 if (a.err) return a; 185 bool inv = n < 0; 186 uint64_t k = inv ? (uint64_t)0 - (uint64_t)n : (uint64_t)n; 187 Dec r = DEC_ONE, b = a; 188 while (k) { 189 if (k & 1) r = dec_mul(r, b); 190 k >>= 1; 191 if (k) b = dec_mul(b, b); 192 if (r.err || b.err) return dec_err(); 193 } 194 return inv ? dec_div(DEC_ONE, r) : r; 195 } 196 197 bool dec_is_int(Dec a) 198 { 199 if (a.err) return false; 200 if (!a.m || a.e >= 17) return true; 201 if (a.e < 0) return false; 202 return a.m % P10[17 - a.e] == 0; 203 } 204 205 bool dec_to_int(Dec a, int64_t *out) 206 { 207 if (!dec_is_int(a) || a.e > 18) return false; 208 uint64_t v = !a.m ? 0 : a.e >= 17 ? a.m * P10[a.e - 17] : a.m / P10[17 - a.e]; 209 if (v > 9223372036854775807ull) return false; 210 *out = a.neg ? -(int64_t)v : (int64_t)v; 211 return true; 212 } 213 214 Dec dec_trunc(Dec a) 215 { 216 if (a.err || !a.m || a.e >= 17) return a; 217 if (a.e < 0) return DEC_ZERO; 218 uint64_t p = P10[17 - a.e]; 219 a.m = a.m / p * p; 220 return a; 221 } 222 223 Dec dec_floor(Dec a) 224 { 225 Dec t = dec_trunc(a); 226 if (a.neg && dec_cmp(t, a) != 0) t = dec_sub(t, DEC_ONE); 227 return t; 228 } 229 230 /* Redondea dejando 'keep' digitos de la mantisa (0..18), mitad lejos de cero. */ 231 static Dec round_keep(Dec a, int keep) 232 { 233 if (a.err || !a.m || keep >= 18) return a; 234 if (keep < 0) return DEC_ZERO; 235 uint64_t p = P10[18 - keep]; 236 uint64_t q = a.m / p, r = a.m % p; 237 if (r >= p / 2) q++; 238 if (keep == 0) { /* todo era redondeo: queda 0 o 1 en el digito de arriba */ 239 if (!q) return DEC_ZERO; 240 a.m = E17; a.e++; 241 return a.e > DEC_EMAX ? dec_err() : a; 242 } 243 q *= p; 244 if (q >= E18) { q /= 10; a.e++; if (a.e > DEC_EMAX) return dec_err(); } 245 a.m = q; 246 return a; 247 } 248 249 Dec dec_round_int(Dec a) { return a.e >= 17 ? a : round_keep(a, a.e + 1); } 250 Dec dec_round_sig(Dec a, int n) { return round_keep(a, n); } 251 Dec dec_round_dec(Dec a, int n) { return a.e + 1 + n >= 18 ? a : round_keep(a, a.e + 1 + n); } 252 253 Dec dec_fmod(Dec a, Dec b) 254 { 255 if (a.err || b.err || !b.m) return dec_err(); 256 Dec q = dec_trunc(dec_div(a, b)); 257 Dec r = dec_sub(a, dec_mul(q, b)); 258 /* el cociente redondeado puede pasarse en uno */ 259 if (r.m && r.neg != a.neg) r = dec_add(r, a.neg == b.neg ? b : dec_neg(b)); 260 if (dec_cmp_abs(r, b) >= 0) r = dec_sub(r, a.neg == b.neg ? b : dec_neg(b)); 261 return r; 262 } 263 264 Dec dec_snap(Dec a, int tol) 265 { 266 if (a.err || !a.m) return a; 267 if (a.e < -tol) return DEC_ZERO; 268 if (a.e == -1 || a.e == 0) { 269 Dec d = dec_sub(dec_abs(a), DEC_ONE); 270 if (!d.m || d.e < -tol) { Dec one = DEC_ONE; one.neg = a.neg; return one; } 271 } 272 return a; 273 } 274 275 void dec_digits(Dec a, char d[DEC_DIGITS + 1]) 276 { 277 uint64_t m = a.m; 278 for (int i = DEC_DIGITS - 1; i >= 0; i--) { d[i] = (char)('0' + m % 10); m /= 10; } 279 d[DEC_DIGITS] = 0; 280 } 281 282 Dec dec_parse(const char *s, const char **end) 283 { 284 const char *p = s; 285 int neg = 0; 286 if (*p == '-' || *p == '+') { neg = *p == '-'; p++; } 287 uint64_t m = 0; 288 int nd = 0, exp = 0, ri = 0, any = 0; 289 bool dot = false; 290 for (;; p++) { 291 if (*p == '.' && !dot) { dot = true; continue; } 292 if (*p < '0' || *p > '9') break; 293 any = 1; 294 int dg = *p - '0'; 295 if (!m && !dg) { if (dot) exp--; continue; } /* ceros a la izquierda */ 296 if (nd < 19) { m = m * 10 + (uint64_t)dg; nd++; if (dot) exp--; } 297 else { /* digitos de mas: solo redondeo */ 298 if (!dot) exp++; 299 if (nd == 19) ri = dg > 5 ? 3 : dg == 5 ? 2 : dg ? 1 : 0; 300 else if (dg && ri != 3) ri = ri == 2 ? 3 : ri ? ri : 1; 301 nd++; 302 } 303 } 304 if (!any) { if (end) *end = 0; return DEC_ZERO; } 305 if ((*p == 'e' || *p == 'E') && ((p[1] >= '0' && p[1] <= '9') || 306 ((p[1] == '-' || p[1] == '+') && p[2] >= '0' && p[2] <= '9'))) { 307 long x; 308 const char *q = str_scan_int(p + 1, &x); 309 if (q) { p = q; if (x > 9999) x = 9999; if (x < -9999) x = -9999; exp += (int)x; } 310 } 311 if (end) *end = p; 312 if (!m) return DEC_ZERO; 313 return pack(neg, m, exp, ri); 314 } 315 316 int dec_to_text(Dec a, char *buf, int n) 317 { 318 if (a.err) return str_put(buf, n, 0, "Error"); 319 if (!a.m) return str_put(buf, n, 0, "0"); 320 char d[DEC_DIGITS + 1]; 321 dec_digits(a, d); 322 int nd = DEC_DIGITS; 323 while (nd > 1 && d[nd - 1] == '0') nd--; 324 int p = 0; 325 if (a.neg) p = str_put(buf, n, p, "-"); 326 char c[2] = { d[0], 0 }; 327 p = str_put(buf, n, p, c); 328 if (nd > 1) { 329 p = str_put(buf, n, p, "."); 330 d[nd] = 0; 331 p = str_put(buf, n, p, d + 1); 332 } 333 if (a.e) { 334 p = str_put(buf, n, p, "E"); 335 p = str_int(buf, n, p, a.e); 336 } 337 return p; 338 } 339 340 int i64_text(char *buf, int n, int pos, int64_t v) 341 { 342 char t[24]; 343 int k = 0; 344 uint64_t u = v < 0 ? (uint64_t)0 - (uint64_t)v : (uint64_t)v; 345 do { t[k++] = (char)('0' + u % 10); u /= 10; } while (u); 346 if (v < 0) t[k++] = '-'; 347 char r[24]; 348 for (int i = 0; i < k; i++) r[i] = t[k - 1 - i]; 349 r[k] = 0; 350 return str_put(buf, n, pos, r); 351 }