1 #include <u.h> 2 #include <libc.h> 3 #include "fmtdef.h" 4 5 static ulong 6 umuldiv(ulong a, ulong b, ulong c) 7 { 8 double d; 9 10 d = ((double)a * (double)b) / (double)c; 11 if(d >= 4294967295.) 12 d = 4294967295.; 13 return (ulong)d; 14 } 15 16 /* 17 * This routine will convert to arbitrary precision 18 * floating point entirely in multi-precision fixed. 19 * The answer is the closest floating point number to 20 * the given decimal number. Exactly half way are 21 * rounded ala ieee rules. 22 * Method is to scale input decimal between .500 and .999... 23 * with external power of 2, then binary search for the 24 * closest mantissa to this decimal number. 25 * Nmant is is the required precision. (53 for ieee dp) 26 * Nbits is the max number of bits/word. (must be <= 28) 27 * Prec is calculated - the number of words of fixed mantissa. 28 */ 29 enum 30 { 31 Nbits = 28, /* bits safely represented in a ulong */ 32 Nmant = 53, /* bits of precision required */ 33 Prec = (Nmant+Nbits+1)/Nbits, /* words of Nbits each to represent mantissa */ 34 Sigbit = 1<<(Prec*Nbits-Nmant), /* first significant bit of Prec-th word */ 35 Ndig = 1500, 36 One = (ulong)(1<<Nbits), 37 Half = (ulong)(One>>1), 38 Maxe = 310, 39 40 Fsign = 1<<0, /* found - */ 41 Fesign = 1<<1, /* found e- */ 42 Fdpoint = 1<<2, /* found . */ 43 44 S0 = 0, /* _ _S0 +S1 #S2 .S3 */ 45 S1, /* _+ #S2 .S3 */ 46 S2, /* _+# #S2 .S4 eS5 */ 47 S3, /* _+. #S4 */ 48 S4, /* _+#.# #S4 eS5 */ 49 S5, /* _+#.#e +S6 #S7 */ 50 S6, /* _+#.#e+ #S7 */ 51 S7, /* _+#.#e+# #S7 */ 52 }; 53 54 static int xcmp(char*, char*); 55 static int fpcmp(char*, ulong*); 56 static void frnorm(ulong*); 57 static void divascii(char*, int*, int*, int*); 58 static void mulascii(char*, int*, int*, int*); 59 60 typedef struct Tab Tab; 61 struct Tab 62 { 63 int bp; 64 int siz; 65 char* cmp; 66 }; 67 68 #ifndef ERANGE 69 #define ERANGE 12345 70 #endif 71 72 double 73 fmtstrtod(const char *as, char **aas) 74 { 75 int na, ex, dp, bp, c, i, flag, state; 76 ulong low[Prec], hig[Prec], mid[Prec]; 77 double d; 78 char *s, a[Ndig]; 79 80 flag = 0; /* Fsign, Fesign, Fdpoint */ 81 na = 0; /* number of digits of a[] */ 82 dp = 0; /* na of decimal point */ 83 ex = 0; /* exonent */ 84 85 state = S0; 86 for(s=(char*)as;; s++) { 87 c = *s; 88 if(c >= '0' && c <= '9') { 89 switch(state) { 90 case S0: 91 case S1: 92 case S2: 93 state = S2; 94 break; 95 case S3: 96 case S4: 97 state = S4; 98 break; 99 100 case S5: 101 case S6: 102 case S7: 103 state = S7; 104 ex = ex*10 + (c-'0'); 105 continue; 106 } 107 if(na == 0 && c == '0') { 108 dp--; 109 continue; 110 } 111 if(na < Ndig-50) 112 a[na++] = c; 113 continue; 114 } 115 switch(c) { 116 case '\t': 117 case '\n': 118 case '\v': 119 case '\f': 120 case '\r': 121 case ' ': 122 if(state == S0) 123 continue; 124 break; 125 case '-': 126 if(state == S0) 127 flag |= Fsign; 128 else 129 flag |= Fesign; 130 case '+': 131 if(state == S0) 132 state = S1; 133 else 134 if(state == S5) 135 state = S6; 136 else 137 break; /* syntax */ 138 continue; 139 case '.': 140 flag |= Fdpoint; 141 dp = na; 142 if(state == S0 || state == S1) { 143 state = S3; 144 continue; 145 } 146 if(state == S2) { 147 state = S4; 148 continue; 149 } 150 break; 151 case 'e': 152 case 'E': 153 if(state == S2 || state == S4) { 154 state = S5; 155 continue; 156 } 157 break; 158 } 159 break; 160 } 161 162 /* 163 * clean up return char-pointer 164 */ 165 switch(state) { 166 case S0: 167 if(xcmp(s, "nan") == 0) { 168 if(aas != nil) 169 *aas = s+3; 170 goto retnan; 171 } 172 case S1: 173 if(xcmp(s, "infinity") == 0) { 174 if(aas != nil) 175 *aas = s+8; 176 goto retinf; 177 } 178 if(xcmp(s, "inf") == 0) { 179 if(aas != nil) 180 *aas = s+3; 181 goto retinf; 182 } 183 case S3: 184 if(aas != nil) 185 *aas = (char*)as; 186 goto ret0; /* no digits found */ 187 case S6: 188 s--; /* back over +- */ 189 case S5: 190 s--; /* back over e */ 191 break; 192 } 193 if(aas != nil) 194 *aas = s; 195 196 if(flag & Fdpoint) 197 while(na > 0 && a[na-1] == '0') 198 na--; 199 if(na == 0) 200 goto ret0; /* zero */ 201 a[na] = 0; 202 if(!(flag & Fdpoint)) 203 dp = na; 204 if(flag & Fesign) 205 ex = -ex; 206 dp += ex; 207 if(dp < -Maxe){ 208 errno = ERANGE; 209 goto ret0; /* underflow by exp */ 210 } else 211 if(dp > +Maxe) 212 goto retinf; /* overflow by exp */ 213 214 /* 215 * normalize the decimal ascii number 216 * to range .[5-9][0-9]* e0 217 */ 218 bp = 0; /* binary exponent */ 219 while(dp > 0) 220 divascii(a, &na, &dp, &bp); 221 while(dp < 0 || a[0] < '5') 222 mulascii(a, &na, &dp, &bp); 223 224 /* close approx by naive conversion */ 225 mid[0] = 0; 226 mid[1] = 1; 227 for(i=0; (c=a[i]); i++) { 228 mid[0] = mid[0]*10 + (c-'0'); 229 mid[1] = mid[1]*10; 230 if(i >= 8) 231 break; 232 } 233 low[0] = umuldiv(mid[0], One, mid[1]); 234 hig[0] = umuldiv(mid[0]+1, One, mid[1]); 235 for(i=1; i<Prec; i++) { 236 low[i] = 0; 237 hig[i] = One-1; 238 } 239 240 /* binary search for closest mantissa */ 241 for(;;) { 242 /* mid = (hig + low) / 2 */ 243 c = 0; 244 for(i=0; i<Prec; i++) { 245 mid[i] = hig[i] + low[i]; 246 if(c) 247 mid[i] += One; 248 c = mid[i] & 1; 249 mid[i] >>= 1; 250 } 251 frnorm(mid); 252 253 /* compare */ 254 c = fpcmp(a, mid); 255 if(c > 0) { 256 c = 1; 257 for(i=0; i<Prec; i++) 258 if(low[i] != mid[i]) { 259 c = 0; 260 low[i] = mid[i]; 261 } 262 if(c) 263 break; /* between mid and hig */ 264 continue; 265 } 266 if(c < 0) { 267 for(i=0; i<Prec; i++) 268 hig[i] = mid[i]; 269 continue; 270 } 271 272 /* only hard part is if even/odd roundings wants to go up */ 273 c = mid[Prec-1] & (Sigbit-1); 274 if(c == Sigbit/2 && (mid[Prec-1]&Sigbit) == 0) 275 mid[Prec-1] -= c; 276 break; /* exactly mid */ 277 } 278 279 /* normal rounding applies */ 280 c = mid[Prec-1] & (Sigbit-1); 281 mid[Prec-1] -= c; 282 if(c >= Sigbit/2) { 283 mid[Prec-1] += Sigbit; 284 frnorm(mid); 285 } 286 goto out; 287 288 ret0: 289 return 0; 290 291 retnan: 292 return __NaN(); 293 294 retinf: 295 /* 296 * Unix strtod requires these. Plan 9 would return Inf(0) or Inf(-1). */ 297 errno = ERANGE; 298 if(flag & Fsign) 299 return -HUGE_VAL; 300 return HUGE_VAL; 301 302 out: 303 d = 0; 304 for(i=0; i<Prec; i++) 305 d = d*One + mid[i]; 306 if(flag & Fsign) 307 d = -d; 308 d = ldexp(d, bp - Prec*Nbits); 309 if(d == 0){ /* underflow */ 310 errno = ERANGE; 311 } 312 return d; 313 } 314 315 static void 316 frnorm(ulong *f) 317 { 318 int i, c; 319 320 c = 0; 321 for(i=Prec-1; i>0; i--) { 322 f[i] += c; 323 c = f[i] >> Nbits; 324 f[i] &= One-1; 325 } 326 f[0] += c; 327 } 328 329 static int 330 fpcmp(char *a, ulong* f) 331 { 332 ulong tf[Prec]; 333 int i, d, c; 334 335 for(i=0; i<Prec; i++) 336 tf[i] = f[i]; 337 338 for(;;) { 339 /* tf *= 10 */ 340 for(i=0; i<Prec; i++) 341 tf[i] = tf[i]*10; 342 frnorm(tf); 343 d = (tf[0] >> Nbits) + '0'; 344 tf[0] &= One-1; 345 346 /* compare next digit */ 347 c = *a; 348 if(c == 0) { 349 if('0' < d) 350 return -1; 351 if(tf[0] != 0) 352 goto cont; 353 for(i=1; i<Prec; i++) 354 if(tf[i] != 0) 355 goto cont; 356 return 0; 357 } 358 if(c > d) 359 return +1; 360 if(c < d) 361 return -1; 362 a++; 363 cont:; 364 } 365 } 366 367 static void 368 divby(char *a, int *na, int b) 369 { 370 int n, c; 371 char *p; 372 373 p = a; 374 n = 0; 375 while(n>>b == 0) { 376 c = *a++; 377 if(c == 0) { 378 while(n) { 379 c = n*10; 380 if(c>>b) 381 break; 382 n = c; 383 } 384 goto xx; 385 } 386 n = n*10 + c-'0'; 387 (*na)--; 388 } 389 for(;;) { 390 c = n>>b; 391 n -= c<<b; 392 *p++ = c + '0'; 393 c = *a++; 394 if(c == 0) 395 break; 396 n = n*10 + c-'0'; 397 } 398 (*na)++; 399 xx: 400 while(n) { 401 n = n*10; 402 c = n>>b; 403 n -= c<<b; 404 *p++ = c + '0'; 405 (*na)++; 406 } 407 *p = 0; 408 } 409 410 static Tab tab1[] = 411 { 412 1, 0, "", 413 3, 1, "7", 414 6, 2, "63", 415 9, 3, "511", 416 13, 4, "8191", 417 16, 5, "65535", 418 19, 6, "524287", 419 23, 7, "8388607", 420 26, 8, "67108863", 421 27, 9, "134217727", 422 }; 423 424 static void 425 divascii(char *a, int *na, int *dp, int *bp) 426 { 427 int b, d; 428 Tab *t; 429 430 d = *dp; 431 if(d >= (int)(nelem(tab1))) 432 d = (int)(nelem(tab1))-1; 433 t = tab1 + d; 434 b = t->bp; 435 if(memcmp(a, t->cmp, t->siz) > 0) 436 d--; 437 *dp -= d; 438 *bp += b; 439 divby(a, na, b); 440 } 441 442 static void 443 mulby(char *a, char *p, char *q, int b) 444 { 445 int n, c; 446 447 n = 0; 448 *p = 0; 449 for(;;) { 450 q--; 451 if(q < a) 452 break; 453 c = *q - '0'; 454 c = (c<<b) + n; 455 n = c/10; 456 c -= n*10; 457 p--; 458 *p = c + '0'; 459 } 460 while(n) { 461 c = n; 462 n = c/10; 463 c -= n*10; 464 p--; 465 *p = c + '0'; 466 } 467 } 468 469 static Tab tab2[] = 470 { 471 1, 1, "", /* dp = 0-0 */ 472 3, 3, "125", 473 6, 5, "15625", 474 9, 7, "1953125", 475 13, 10, "1220703125", 476 16, 12, "152587890625", 477 19, 14, "19073486328125", 478 23, 17, "11920928955078125", 479 26, 19, "1490116119384765625", 480 27, 19, "7450580596923828125", /* dp 8-9 */ 481 }; 482 483 static void 484 mulascii(char *a, int *na, int *dp, int *bp) 485 { 486 char *p; 487 int d, b; 488 Tab *t; 489 490 d = -*dp; 491 if(d >= (int)(nelem(tab2))) 492 d = (int)(nelem(tab2))-1; 493 t = tab2 + d; 494 b = t->bp; 495 if(memcmp(a, t->cmp, t->siz) < 0) 496 d--; 497 p = a + *na; 498 *bp -= b; 499 *dp += d; 500 *na += d; 501 mulby(a, p+d, p, b); 502 } 503 504 static int 505 xcmp(char *a, char *b) 506 { 507 int c1, c2; 508 509 while((c1 = *b++)) { 510 c2 = *a++; 511 if(isupper(c2)) 512 c2 = tolower(c2); 513 if(c1 != c2) 514 return 1; 515 } 516 return 0; 517 } 518