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