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