QRScanner/mobile/jni/thirdparty/mujs/jsdtoa.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 <stdio.h>
  15 #include <math.h>
  16 #include <float.h>
  17 #include <string.h>
  18 #include <stdlib.h>
  19 #include <errno.h>
  20
  21 #include "jsi.h"
  22
  23 typedef unsigned long ulong;
  24
  25 enum { NSIGNIF  = 17 };
  26
  27 /*
  28  * first few powers of 10, enough for about 1/2 of the
  29  * total space for doubles.
  30  */
  31 static double pows10[] =
  32 {
  33         1e0,    1e1,    1e2,    1e3,    1e4,    1e5,    1e6,    1e7,    1e8,    1e9,
  34         1e10,   1e11,   1e12,   1e13,   1e14,   1e15,   1e16,   1e17,   1e18,   1e19,
  35         1e20,   1e21,   1e22,   1e23,   1e24,   1e25,   1e26,   1e27,   1e28,   1e29,
  36         1e30,   1e31,   1e32,   1e33,   1e34,   1e35,   1e36,   1e37,   1e38,   1e39,
  37         1e40,   1e41,   1e42,   1e43,   1e44,   1e45,   1e46,   1e47,   1e48,   1e49,
  38         1e50,   1e51,   1e52,   1e53,   1e54,   1e55,   1e56,   1e57,   1e58,   1e59,
  39         1e60,   1e61,   1e62,   1e63,   1e64,   1e65,   1e66,   1e67,   1e68,   1e69,
  40         1e70,   1e71,   1e72,   1e73,   1e74,   1e75,   1e76,   1e77,   1e78,   1e79,
  41         1e80,   1e81,   1e82,   1e83,   1e84,   1e85,   1e86,   1e87,   1e88,   1e89,
  42         1e90,   1e91,   1e92,   1e93,   1e94,   1e95,   1e96,   1e97,   1e98,   1e99,
  43         1e100,  1e101,  1e102,  1e103,  1e104,  1e105,  1e106,  1e107,  1e108,  1e109,
  44         1e110,  1e111,  1e112,  1e113,  1e114,  1e115,  1e116,  1e117,  1e118,  1e119,
  45         1e120,  1e121,  1e122,  1e123,  1e124,  1e125,  1e126,  1e127,  1e128,  1e129,
  46         1e130,  1e131,  1e132,  1e133,  1e134,  1e135,  1e136,  1e137,  1e138,  1e139,
  47         1e140,  1e141,  1e142,  1e143,  1e144,  1e145,  1e146,  1e147,  1e148,  1e149,
  48         1e150,  1e151,  1e152,  1e153,  1e154,  1e155,  1e156,  1e157,  1e158,  1e159,
  49 };
  50 #define npows10 ((int)(sizeof(pows10)/sizeof(pows10[0])))
  51 #define pow10(x) fmtpow10(x)
  52
  53 static double
  54 pow10(int n)
  55 {
  56         double d;
  57         int neg;
  58
  59         neg = 0;
  60         if(n < 0){
  61                 neg = 1;
  62                 n = -n;
  63         }
  64
  65         if(n < npows10)
  66                 d = pows10[n];
  67         else{
  68                 d = pows10[npows10-1];
  69                 for(;;){
  70                         n -= npows10 - 1;
  71                         if(n < npows10){
  72                                 d *= pows10[n];
  73                                 break;
  74                         }
  75                         d *= pows10[npows10 - 1];
  76                 }
  77         }
  78         if(neg)
  79                 return 1./d;
  80         return d;
  81 }
  82
  83 /*
  84  * add 1 to the decimal integer string a of length n.
  85  * if 99999 overflows into 10000, return 1 to tell caller
  86  * to move the virtual decimal point.
  87  */
  88 static int
  89 xadd1(char *a, int n)
  90 {
  91         char *b;
  92         int c;
  93
  94         if(n < 0 || n > NSIGNIF)
  95                 return 0;
  96         for(b = a+n-1; b >= a; b--) {
  97                 c = *b + 1;
  98                 if(c <= '9') {
  99                         *b = c;
 100                         return 0;
 101                 }
 102                 *b = '0';
 103         }
 104         /*
 105          * need to overflow adding digit.
 106          * shift number down and insert 1 at beginning.
 107          * decimal is known to be 0s or we wouldn't
 108          * have gotten this far. (e.g., 99999+1 => 00000)
 109          */
 110         a[0] = '1';
 111         return 1;
 112 }
 113
 114 /*
 115  * subtract 1 from the decimal integer string a.
 116  * if 10000 underflows into 09999, make it 99999
 117  * and return 1 to tell caller to move the virtual
 118  * decimal point. this way, xsub1 is inverse of xadd1.
 119  */
 120 static int
 121 xsub1(char *a, int n)
 122 {
 123         char *b;
 124         int c;
 125
 126         if(n < 0 || n > NSIGNIF)
 127                 return 0;
 128         for(b = a+n-1; b >= a; b--) {
 129                 c = *b - 1;
 130                 if(c >= '0') {
 131                         if(c == '0' && b == a) {
 132                                 /*
 133                                  * just zeroed the top digit; shift everyone up.
 134                                  * decimal is known to be 9s or we wouldn't
 135                                  * have gotten this far. (e.g., 10000-1 => 09999)
 136                                  */
 137                                 *b = '9';
 138                                 return 1;
 139                         }
 140                         *b = c;
 141                         return 0;
 142                 }
 143                 *b = '9';
 144         }
 145         /*
 146          * can't get here. the number a is always normalized
 147          * so that it has a nonzero first digit.
 148          */
 149         return 0;
 150 }
 151
 152 /*
 153  * format exponent like sprintf(p, "e%+d", e)
 154  */
 155 void
 156 js_fmtexp(char *p, int e)
 157 {
 158         char se[9];
 159         int i;
 160
 161         *p++ = 'e';
 162         if(e < 0) {
 163                 *p++ = '-';
 164                 e = -e;
 165         } else
 166                 *p++ = '+';
 167         i = 0;
 168         while(e) {
 169                 se[i++] = e % 10 + '0';
 170                 e /= 10;
 171         }
 172         while(i < 1)
 173                 se[i++] = '0';
 174         while(i > 0)
 175                 *p++ = se[--i];
 176         *p++ = '\0';
 177 }
 178
 179 /*
 180  * compute decimal integer m, exp such that:
 181  *      f = m*10^exp
 182  *      m is as short as possible with losing exactness
 183  * assumes special cases (NaN, +Inf, -Inf) have been handled.
 184  */
 185 void
 186 js_dtoa(double f, char *s, int *exp, int *neg, int *ns)
 187 {
 188         int c, d, e2, e, ee, i, ndigit, oerrno;
 189         char tmp[NSIGNIF+10];
 190         double g;
 191
 192         oerrno = errno; /* in case strtod smashes errno */
 193
 194         /*
 195          * make f non-negative.
 196          */
 197         *neg = 0;
 198         if(f < 0) {
 199                 f = -f;
 200                 *neg = 1;
 201         }
 202
 203         /*
 204          * must handle zero specially.
 205          */
 206         if(f == 0){
 207                 *exp = 0;
 208                 s[0] = '0';
 209                 s[1] = '\0';
 210                 *ns = 1;
 211                 return;
 212         }
 213
 214         /*
 215          * find g,e such that f = g*10^e.
 216          * guess 10-exponent using 2-exponent, then fine tune.
 217          */
 218         frexp(f, &e2);
 219         e = (int)(e2 * .301029995664);
 220         g = f * pow10(-e);
 221         while(g < 1) {
 222                 e--;
 223                 g = f * pow10(-e);
 224         }
 225         while(g >= 10) {
 226                 e++;
 227                 g = f * pow10(-e);
 228         }
 229
 230         /*
 231          * convert NSIGNIF digits as a first approximation.
 232          */
 233         for(i=0; i<NSIGNIF; i++) {
 234                 d = (int)g;
 235                 s[i] = d+'0';
 236                 g = (g-d) * 10;
 237         }
 238         s[i] = 0;
 239
 240         /*
 241          * adjust e because s is 314159... not 3.14159...
 242          */
 243         e -= NSIGNIF-1;
 244         js_fmtexp(s+NSIGNIF, e);
 245
 246         /*
 247          * adjust conversion until strtod(s) == f exactly.
 248          */
 249         for(i=0; i<10; i++) {
 250                 g = js_strtod(s, NULL);
 251                 if(f > g) {
 252                         if(xadd1(s, NSIGNIF)) {
 253                                 /* gained a digit */
 254                                 e--;
 255                                 js_fmtexp(s+NSIGNIF, e);
 256                         }
 257                         continue;
 258                 }
 259                 if(f < g) {
 260                         if(xsub1(s, NSIGNIF)) {
 261                                 /* lost a digit */
 262                                 e++;
 263                                 js_fmtexp(s+NSIGNIF, e);
 264                         }
 265                         continue;
 266                 }
 267                 break;
 268         }
 269
 270         /*
 271          * play with the decimal to try to simplify.
 272          */
 273
 274         /*
 275          * bump last few digits up to 9 if we can
 276          */
 277         for(i=NSIGNIF-1; i>=NSIGNIF-3; i--) {
 278                 c = s[i];
 279                 if(c != '9') {
 280                         s[i] = '9';
 281                         g = js_strtod(s, NULL);
 282                         if(g != f) {
 283                                 s[i] = c;
 284                                 break;
 285                         }
 286                 }
 287         }
 288
 289         /*
 290          * add 1 in hopes of turning 9s to 0s
 291          */
 292         if(s[NSIGNIF-1] == '9') {
 293                 strcpy(tmp, s);
 294                 ee = e;
 295                 if(xadd1(tmp, NSIGNIF)) {
 296                         ee--;
 297                         js_fmtexp(tmp+NSIGNIF, ee);
 298                 }
 299                 g = js_strtod(tmp, NULL);
 300                 if(g == f) {
 301                         strcpy(s, tmp);
 302                         e = ee;
 303                 }
 304         }
 305
 306         /*
 307          * bump last few digits down to 0 as we can.
 308          */
 309         for(i=NSIGNIF-1; i>=NSIGNIF-3; i--) {
 310                 c = s[i];
 311                 if(c != '0') {
 312                         s[i] = '0';
 313                         g = js_strtod(s, NULL);
 314                         if(g != f) {
 315                                 s[i] = c;
 316                                 break;
 317                         }
 318                 }
 319         }
 320
 321         /*
 322          * remove trailing zeros.
 323          */
 324         ndigit = NSIGNIF;
 325         while(ndigit > 1 && s[ndigit-1] == '0'){
 326                 e++;
 327                 --ndigit;
 328         }
 329         s[ndigit] = 0;
 330         *exp = e;
 331         *ns = ndigit;
 332         errno = oerrno;
 333 }
 334
 335 static ulong
 336 umuldiv(ulong a, ulong b, ulong c)
 337 {
 338         double d;
 339
 340         d = ((double)a * (double)b) / (double)c;
 341         if(d >= 4294967295.)
 342                 d = 4294967295.;
 343         return (ulong)d;
 344 }
 345
 346 /*
 347  * This routine will convert to arbitrary precision
 348  * floating point entirely in multi-precision fixed.
 349  * The answer is the closest floating point number to
 350  * the given decimal number. Exactly half way are
 351  * rounded ala ieee rules.
 352  * Method is to scale input decimal between .500 and .999...
 353  * with external power of 2, then binary search for the
 354  * closest mantissa to this decimal number.
 355  * Nmant is is the required precision. (53 for ieee dp)
 356  * Nbits is the max number of bits/word. (must be <= 28)
 357  * Prec is calculated - the number of words of fixed mantissa.
 358  */
 359 enum
 360 {
 361         Nbits   = 28,                           /* bits safely represented in a ulong */
 362         Nmant   = 53,                           /* bits of precision required */
 363         Prec    = (Nmant+Nbits+1)/Nbits,        /* words of Nbits each to represent mantissa */
 364         Sigbit  = 1<<(Prec*Nbits-Nmant),        /* first significant bit of Prec-th word */
 365         Ndig    = 1500,
 366         One     = (ulong)(1<<Nbits),
 367         Half    = (ulong)(One>>1),
 368         Maxe    = 310,
 369
 370         Fsign   = 1<<0,         /* found - */
 371         Fesign  = 1<<1,         /* found e- */
 372         Fdpoint = 1<<2,         /* found . */
 373
 374         S0      = 0,            /* _            _S0     +S1     #S2     .S3 */
 375         S1,                     /* _+           #S2     .S3 */
 376         S2,                     /* _+#          #S2     .S4     eS5 */
 377         S3,                     /* _+.          #S4 */
 378         S4,                     /* _+#.#        #S4     eS5 */
 379         S5,                     /* _+#.#e       +S6     #S7 */
 380         S6,                     /* _+#.#e+      #S7 */
 381         S7                      /* _+#.#e+#     #S7 */
 382 };
 383
 384 static  int     xcmp(char*, char*);
 385 static  int     fpcmp(char*, ulong*);
 386 static  void    frnorm(ulong*);
 387 static  void    divascii(char*, int*, int*, int*);
 388 static  void    mulascii(char*, int*, int*, int*);
 389
 390 typedef struct  Tab     Tab;
 391 struct  Tab
 392 {
 393         int     bp;
 394         int     siz;
 395         char*   cmp;
 396 };
 397
 398 double
 399 js_strtod(const char *as, char **aas)
 400 {
 401         int na, ex, dp, bp, c, i, flag, state;
 402         ulong low[Prec], hig[Prec], mid[Prec];
 403         double d;
 404         char *s, a[Ndig];
 405
 406         flag = 0;       /* Fsign, Fesign, Fdpoint */
 407         na = 0;         /* number of digits of a[] */
 408         dp = 0;         /* na of decimal point */
 409         ex = 0;         /* exonent */
 410
 411         state = S0;
 412         for(s=(char*)as;; s++) {
 413                 c = *s;
 414                 if(c >= '0' && c <= '9') {
 415                         switch(state) {
 416                         case S0:
 417                         case S1:
 418                         case S2:
 419                                 state = S2;
 420                                 break;
 421                         case S3:
 422                         case S4:
 423                                 state = S4;
 424                                 break;
 425
 426                         case S5:
 427                         case S6:
 428                         case S7:
 429                                 state = S7;
 430                                 ex = ex*10 + (c-'0');
 431                                 continue;
 432                         }
 433                         if(na == 0 && c == '0') {
 434                                 dp--;
 435                                 continue;
 436                         }
 437                         if(na < Ndig-50)
 438                                 a[na++] = c;
 439                         continue;
 440                 }
 441                 switch(c) {
 442                 case '\t':
 443                 case '\n':
 444                 case '\v':
 445                 case '\f':
 446                 case '\r':
 447                 case ' ':
 448                         if(state == S0)
 449                                 continue;
 450                         break;
 451                 case '-':
 452                         if(state == S0)
 453                                 flag |= Fsign;
 454                         else
 455                                 flag |= Fesign;
 456                         /* fall through */
 457                 case '+':
 458                         if(state == S0)
 459                                 state = S1;
 460                         else
 461                         if(state == S5)
 462                                 state = S6;
 463                         else
 464                                 break;  /* syntax */
 465                         continue;
 466                 case '.':
 467                         flag |= Fdpoint;
 468                         dp = na;
 469                         if(state == S0 || state == S1) {
 470                                 state = S3;
 471                                 continue;
 472                         }
 473                         if(state == S2) {
 474                                 state = S4;
 475                                 continue;
 476                         }
 477                         break;
 478                 case 'e':
 479                 case 'E':
 480                         if(state == S2 || state == S4) {
 481                                 state = S5;
 482                                 continue;
 483                         }
 484                         break;
 485                 }
 486                 break;
 487         }
 488
 489         /*
 490          * clean up return char-pointer
 491          */
 492         switch(state) {
 493         case S0:
 494                 if(xcmp(s, "nan") == 0) {
 495                         if(aas != NULL)
 496                                 *aas = s+3;
 497                         goto retnan;
 498                 }
 499                 /* fall through */
 500         case S1:
 501                 if(xcmp(s, "infinity") == 0) {
 502                         if(aas != NULL)
 503                                 *aas = s+8;
 504                         goto retinf;
 505                 }
 506                 if(xcmp(s, "inf") == 0) {
 507                         if(aas != NULL)
 508                                 *aas = s+3;
 509                         goto retinf;
 510                 }
 511                 /* fall through */
 512         case S3:
 513                 if(aas != NULL)
 514                         *aas = (char*)as;
 515                 goto ret0;      /* no digits found */
 516         case S6:
 517                 s--;            /* back over +- */
 518                 /* fall through */
 519         case S5:
 520                 s--;            /* back over e */
 521                 break;
 522         }
 523         if(aas != NULL)
 524                 *aas = s;
 525
 526         if(flag & Fdpoint)
 527         while(na > 0 && a[na-1] == '0')
 528                 na--;
 529         if(na == 0)
 530                 goto ret0;      /* zero */
 531         a[na] = 0;
 532         if(!(flag & Fdpoint))
 533                 dp = na;
 534         if(flag & Fesign)
 535                 ex = -ex;
 536         dp += ex;
 537         if(dp < -Maxe){
 538                 errno = ERANGE;
 539                 goto ret0;      /* underflow by exp */
 540         } else
 541         if(dp > +Maxe)
 542                 goto retinf;    /* overflow by exp */
 543
 544         /*
 545          * normalize the decimal ascii number
 546          * to range .[5-9][0-9]* e0
 547          */
 548         bp = 0;         /* binary exponent */
 549         while(dp > 0)
 550                 divascii(a, &na, &dp, &bp);
 551         while(dp < 0 || a[0] < '5')
 552                 mulascii(a, &na, &dp, &bp);
 553
 554         /* close approx by naive conversion */
 555         mid[0] = 0;
 556         mid[1] = 1;
 557         for(i=0; (c=a[i]) != '\0'; i++) {
 558                 mid[0] = mid[0]*10 + (c-'0');
 559                 mid[1] = mid[1]*10;
 560                 if(i >= 8)
 561                         break;
 562         }
 563         low[0] = umuldiv(mid[0], One, mid[1]);
 564         hig[0] = umuldiv(mid[0]+1, One, mid[1]);
 565         for(i=1; i<Prec; i++) {
 566                 low[i] = 0;
 567                 hig[i] = One-1;
 568         }
 569
 570         /* binary search for closest mantissa */
 571         for(;;) {
 572                 /* mid = (hig + low) / 2 */
 573                 c = 0;
 574                 for(i=0; i<Prec; i++) {
 575                         mid[i] = hig[i] + low[i];
 576                         if(c)
 577                                 mid[i] += One;
 578                         c = mid[i] & 1;
 579                         mid[i] >>= 1;
 580                 }
 581                 frnorm(mid);
 582
 583                 /* compare */
 584                 c = fpcmp(a, mid);
 585                 if(c > 0) {
 586                         c = 1;
 587                         for(i=0; i<Prec; i++)
 588                                 if(low[i] != mid[i]) {
 589                                         c = 0;
 590                                         low[i] = mid[i];
 591                                 }
 592                         if(c)
 593                                 break;  /* between mid and hig */
 594                         continue;
 595                 }
 596                 if(c < 0) {
 597                         for(i=0; i<Prec; i++)
 598                                 hig[i] = mid[i];
 599                         continue;
 600                 }
 601
 602                 /* only hard part is if even/odd roundings wants to go up */
 603                 c = mid[Prec-1] & (Sigbit-1);
 604                 if(c == Sigbit/2 && (mid[Prec-1]&Sigbit) == 0)
 605                         mid[Prec-1] -= c;
 606                 break;  /* exactly mid */
 607         }
 608
 609         /* normal rounding applies */
 610         c = mid[Prec-1] & (Sigbit-1);
 611         mid[Prec-1] -= c;
 612         if(c >= Sigbit/2) {
 613                 mid[Prec-1] += Sigbit;
 614                 frnorm(mid);
 615         }
 616         goto out;
 617
 618 ret0:
 619         if(flag & Fsign)
 620                 return -0.0;
 621         return 0;
 622
 623 retnan:
 624         return NAN;
 625
 626 retinf:
 627         /*
 628          * Unix strtod requires these. Plan 9 would return Inf(0) or Inf(-1). */
 629         errno = ERANGE;
 630         if(flag & Fsign)
 631                 return -HUGE_VAL;
 632         return HUGE_VAL;
 633
 634 out:
 635         d = 0;
 636         for(i=0; i<Prec; i++)
 637                 d = d*One + mid[i];
 638         if(flag & Fsign)
 639                 d = -d;
 640         d = ldexp(d, bp - Prec*Nbits);
 641         if(d == 0){     /* underflow */
 642                 errno = ERANGE;
 643         }
 644         return d;
 645 }
 646
 647 static void
 648 frnorm(ulong *f)
 649 {
 650         int i, c;
 651
 652         c = 0;
 653         for(i=Prec-1; i>0; i--) {
 654                 f[i] += c;
 655                 c = f[i] >> Nbits;
 656                 f[i] &= One-1;
 657         }
 658         f[0] += c;
 659 }
 660
 661 static int
 662 fpcmp(char *a, ulong* f)
 663 {
 664         ulong tf[Prec];
 665         int i, d, c;
 666
 667         for(i=0; i<Prec; i++)
 668                 tf[i] = f[i];
 669
 670         for(;;) {
 671                 /* tf *= 10 */
 672                 for(i=0; i<Prec; i++)
 673                         tf[i] = tf[i]*10;
 674                 frnorm(tf);
 675                 d = (tf[0] >> Nbits) + '0';
 676                 tf[0] &= One-1;
 677
 678                 /* compare next digit */
 679                 c = *a;
 680                 if(c == 0) {
 681                         if('0' < d)
 682                                 return -1;
 683                         if(tf[0] != 0)
 684                                 goto cont;
 685                         for(i=1; i<Prec; i++)
 686                                 if(tf[i] != 0)
 687                                         goto cont;
 688                         return 0;
 689                 }
 690                 if(c > d)
 691                         return +1;
 692                 if(c < d)
 693                         return -1;
 694                 a++;
 695         cont:;
 696         }
 697 }
 698
 699 static void
 700 divby(char *a, int *na, int b)
 701 {
 702         int n, c;
 703         char *p;
 704
 705         p = a;
 706         n = 0;
 707         while(n>>b == 0) {
 708                 c = *a++;
 709                 if(c == 0) {
 710                         while(n) {
 711                                 c = n*10;
 712                                 if(c>>b)
 713                                         break;
 714                                 n = c;
 715                         }
 716                         goto xx;
 717                 }
 718                 n = n*10 + c-'0';
 719                 (*na)--;
 720         }
 721         for(;;) {
 722                 c = n>>b;
 723                 n -= c<<b;
 724                 *p++ = c + '0';
 725                 c = *a++;
 726                 if(c == 0)
 727                         break;
 728                 n = n*10 + c-'0';
 729         }
 730         (*na)++;
 731 xx:
 732         while(n) {
 733                 n = n*10;
 734                 c = n>>b;
 735                 n -= c<<b;
 736                 *p++ = c + '0';
 737                 (*na)++;
 738         }
 739         *p = 0;
 740 }
 741
 742 static  Tab     tab1[] =
 743 {
 744         { 1, 0, "" },
 745         { 3, 1, "7" },
 746         { 6, 2, "63" },
 747         { 9, 3, "511" },
 748         { 13, 4, "8191" },
 749         { 16, 5, "65535" },
 750         { 19, 6, "524287" },
 751         { 23, 7, "8388607" },
 752         { 26, 8, "67108863" },
 753         { 27, 9, "134217727" },
 754 };
 755
 756 static void
 757 divascii(char *a, int *na, int *dp, int *bp)
 758 {
 759         int b, d;
 760         Tab *t;
 761
 762         d = *dp;
 763         if(d >= (int)(nelem(tab1)))
 764                 d = (int)(nelem(tab1))-1;
 765         t = tab1 + d;
 766         b = t->bp;
 767         if(memcmp(a, t->cmp, t->siz) > 0)
 768                 d--;
 769         *dp -= d;
 770         *bp += b;
 771         divby(a, na, b);
 772 }
 773
 774 static void
 775 mulby(char *a, char *p, char *q, int b)
 776 {
 777         int n, c;
 778
 779         n = 0;
 780         *p = 0;
 781         for(;;) {
 782                 q--;
 783                 if(q < a)
 784                         break;
 785                 c = *q - '0';
 786                 c = (c<<b) + n;
 787                 n = c/10;
 788                 c -= n*10;
 789                 p--;
 790                 *p = c + '0';
 791         }
 792         while(n) {
 793                 c = n;
 794                 n = c/10;
 795                 c -= n*10;
 796                 p--;
 797                 *p = c + '0';
 798         }
 799 }
 800
 801 static  Tab     tab2[] =
 802 {
 803         { 1, 1, "" },                           /* dp = 0-0 */
 804         { 3, 3, "125" },
 805         { 6, 5, "15625" },
 806         { 9, 7, "1953125" },
 807         { 13, 10, "1220703125" },
 808         { 16, 12, "152587890625" },
 809         { 19, 14, "19073486328125" },
 810         { 23, 17, "11920928955078125" },
 811         { 26, 19, "1490116119384765625" },
 812         { 27, 19, "7450580596923828125" },              /* dp 8-9 */
 813 };
 814
 815 static void
 816 mulascii(char *a, int *na, int *dp, int *bp)
 817 {
 818         char *p;
 819         int d, b;
 820         Tab *t;
 821
 822         d = -*dp;
 823         if(d >= (int)(nelem(tab2)))
 824                 d = (int)(nelem(tab2))-1;
 825         t = tab2 + d;
 826         b = t->bp;
 827         if(memcmp(a, t->cmp, t->siz) < 0)
 828                 d--;
 829         p = a + *na;
 830         *bp -= b;
 831         *dp += d;
 832         *na += d;
 833         mulby(a, p+d, p, b);
 834 }
 835
 836 static int
 837 xcmp(char *a, char *b)
 838 {
 839         int c1, c2;
 840
 841         while((c1 = *b++) != '\0') {
 842                 c2 = *a++;
 843                 if(c2 >= 'A' && c2 <= 'Z')
 844                         c2 = c2 - 'A' + 'a';
 845                 if(c1 != c2)
 846                         return 1;
 847         }
 848         return 0;
 849 }