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