QRScanner/mobile/jni/fitz/strtod.c

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