URI:
       dc.c - sbase - suckless unix tools
  HTML git clone git://git.suckless.org/sbase
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       dc.c (31359B)
       ---
            1 #include <assert.h>
            2 #include <errno.h>
            3 #include <ctype.h>
            4 #include <setjmp.h>
            5 #include <stdarg.h>
            6 #include <stdio.h>
            7 #include <stdlib.h>
            8 #include <string.h>
            9 #include <limits.h>
           10 
           11 #include "arg.h"
           12 #include "util.h"
           13 
           14 #define NDIGITS 10
           15 #define REGSIZ 16
           16 #define HASHSIZ 128
           17 
           18 enum {
           19         NVAL,
           20         STR,
           21         NUM,
           22 };
           23 
           24 enum {
           25         NE,
           26         LE,
           27         GE,
           28 };
           29 
           30 typedef struct num Num;
           31 typedef struct val Val;
           32 
           33 struct num {
           34         int begin;
           35         int scale;
           36         int room;
           37         signed char *buf, *wp, *rp;
           38 };
           39 
           40 struct digit {
           41         int val;
           42         struct digit *next;
           43 };
           44 
           45 struct val {
           46         int type;
           47         union {
           48                 Num *n;
           49                 char *s;
           50         } u;
           51 
           52         struct val *next;
           53 };
           54 
           55 struct ary {
           56         int n;
           57         Val *buf;
           58         struct ary *next;
           59 };
           60 
           61 struct reg {
           62         char *name;
           63         Val val;
           64         struct ary ary;
           65         struct reg *next;
           66 };
           67 
           68 struct input {
           69         FILE *fp;
           70         char *buf;
           71         size_t n;
           72         char *s;
           73         struct input *next;
           74 };
           75 
           76 static Val *stack;
           77 static jmp_buf env;
           78 static struct input *input;
           79 static struct reg *htbl[HASHSIZ];
           80 
           81 static signed char onestr[] = {1, 0};
           82 static Num zero, one = {.buf = onestr, .wp = onestr + 1};
           83 static char digits[] = "0123456789ABCDEF.";
           84 
           85 static int scale, ibase = 10, obase = 10;
           86 static int iflag;
           87 static int col;
           88 
           89 /*
           90  * This dc implementation follows the decription of dc found in the paper
           91  * DC - An Interactive Desk Calculator, by Robert Morris and Lorinda Cherry
           92  */
           93 
           94 static void
           95 dumpnum(char *msg, Num *num)
           96 {
           97         signed char *p;
           98 
           99         printf("Num (%s)", msg ? msg : "");
          100 
          101         if (!num) {
          102                 puts("NULL\n");
          103                 return;
          104         }
          105         printf("\n"
          106                "\tscale=%d\n"
          107                "\troom=%d\n"
          108                "\tdigits=[\n",
          109                num->scale,
          110                num->room);
          111         for (p = num->buf; p < num->wp; ++p)
          112                 printf("\t\t%02d\n", *p);
          113         printf("\t]\n");
          114 }
          115 
          116 /*
          117  * error is not called from the implementation of the
          118  * arithmetic functions because that can drive to memory
          119  * leaks very easily.
          120  */
          121 static void
          122 error(char *fmt, ...)
          123 {
          124         va_list va;
          125 
          126         va_start(va, fmt);
          127         xvprintf(fmt, va);
          128         putc('\n', stderr);
          129         va_end(va);
          130 
          131         longjmp(env, 1);
          132 }
          133 
          134 static void
          135 freenum(Num *num)
          136 {
          137         if (!num)
          138                 return;
          139         free(num->buf);
          140         free(num);
          141 }
          142 
          143 static Num *
          144 moreroom(Num *num, int more)
          145 {
          146         int ro, wo, room;
          147         signed char *p;
          148 
          149         ro = num->rp - num->buf;
          150         wo = num->wp - num->buf;
          151         room = num->room;
          152 
          153         if (room > INT_MAX - more)
          154                 eprintf("out of memory\n");
          155 
          156         room = room + more;
          157         if (room < NDIGITS)
          158                 room = NDIGITS;
          159         p = erealloc(num->buf, room);
          160         num->buf = p;
          161         num->rp = p + ro;
          162         num->wp = p + wo;
          163         num->room = room;
          164 
          165         return num;
          166 }
          167 
          168 static Num *
          169 grow(Num *num)
          170 {
          171         return moreroom(num, NDIGITS);
          172 }
          173 
          174 static Num *
          175 expand(Num *num, int min)
          176 {
          177         if (min < num->room)
          178                 return num;
          179         return moreroom(num, min - num->room);
          180 }
          181 
          182 static Num *
          183 new(int room)
          184 {
          185         Num *num = emalloc(sizeof(*num));
          186 
          187         num->rp = num->wp = num->buf = NULL;
          188         num->begin = num->room = num->scale = 0;
          189 
          190         return moreroom(num, room);
          191 }
          192 
          193 static Num *
          194 zeronum(int ndigits)
          195 {
          196         Num *num = new(ndigits);
          197 
          198         num->wp = num->buf + ndigits;
          199         memset(num->buf, 0, ndigits);
          200 
          201         return num;
          202 }
          203 
          204 static Num *
          205 wrdigit(Num *num, int d)
          206 {
          207         if (num->wp == &num->buf[num->room])
          208                 grow(num);
          209         *num->wp++ = d;
          210 
          211         return num;
          212 }
          213 
          214 static int
          215 rddigit(Num *num)
          216 {
          217         if (num->rp == num->wp)
          218                 return -1;
          219         return *num->rp++;
          220 }
          221 
          222 static int
          223 peek(Num *num)
          224 {
          225         if (num->rp == num->wp)
          226                 return -1;
          227         return *num->rp;
          228 }
          229 
          230 static Num *
          231 poke(Num *num, unsigned d)
          232 {
          233         if (num->rp == &num->buf[num->room])
          234                 grow(num);
          235         if (num->rp == num->wp)
          236                 num->wp++;
          237         *num->rp = d;
          238 
          239         return num;
          240 }
          241 
          242 static int
          243 begin(Num *num)
          244 {
          245         return num->begin != 0;
          246 }
          247 
          248 static int
          249 first(Num *num)
          250 {
          251         num->begin = 0;
          252         num->rp = num->buf;
          253         return num->rp != num->wp;
          254 }
          255 
          256 static int
          257 last(Num *num)
          258 {
          259         if (num->wp != num->buf) {
          260                 num->begin = 0;
          261                 num->rp = num->wp - 1;
          262                 return 1;
          263         }
          264         num->begin = 1;
          265         return 0;
          266 }
          267 
          268 static int
          269 prev(Num *num)
          270 {
          271         if (num->rp > num->buf) {
          272                 --num->rp;
          273                 return 1;
          274         }
          275         num->begin = 1;
          276         return 0;
          277 }
          278 
          279 static int
          280 next(Num *num)
          281 {
          282         if (num->rp != num->wp + 1) {
          283                 ++num->rp;
          284                 return 1;
          285         }
          286         return 0;
          287 }
          288 
          289 static void
          290 numtrunc(Num *num)
          291 {
          292         num->wp = num->rp;
          293         if (num->rp != num->buf)
          294                 num->rp--;
          295 }
          296 
          297 static int
          298 more(Num *num)
          299 {
          300         return (num->rp != num->wp);
          301 }
          302 
          303 static int
          304 length(Num *num)
          305 {
          306         return num->wp - num->buf;
          307 }
          308 
          309 static int
          310 tell(Num *num)
          311 {
          312         return num->rp - num->buf;
          313 }
          314 
          315 static void
          316 seek(Num *num, int pos)
          317 {
          318         num->rp = num->buf + pos;
          319 }
          320 
          321 static void
          322 rshift(Num *num, int n)
          323 {
          324         int diff;
          325 
          326         diff = length(num) - n;
          327         if (diff < 0) {
          328                 first(num);
          329                 numtrunc(num);
          330                 return;
          331         }
          332 
          333         memmove(num->buf, num->buf + n, diff);
          334         num->rp = num->buf + diff;
          335         numtrunc(num);
          336 }
          337 
          338 static void
          339 lshift(Num *num, int n)
          340 {
          341         int len;
          342 
          343         len = length(num);
          344         expand(num, len + n);
          345         memmove(num->buf + n, num->buf, len);
          346         memset(num->buf, 0, n);
          347         num->wp += n;
          348 }
          349 
          350 static Num *
          351 chsign(Num *num)
          352 {
          353         int val, d, carry;
          354 
          355         carry = 0;
          356         for (first(num); more(num); next(num)) {
          357                 d = peek(num);
          358                 val = 100 - d - carry;
          359 
          360                 carry = 1;
          361                 if (val >= 100) {
          362                         val -= 100;
          363                         carry = 0;
          364                 }
          365                 poke(num, val);
          366         }
          367 
          368         prev(num);
          369         if (carry != 0) {
          370                 if (peek(num) == 99)
          371                         poke(num, -1);
          372                 else
          373                         wrdigit(num, -1);
          374         } else {
          375                 if (peek(num) == 0)
          376                         numtrunc(num);
          377         }
          378 
          379         return num;
          380 }
          381 
          382 static Num *
          383 copy(Num *num)
          384 {
          385         Num *p;
          386         int len = length(num);
          387 
          388         p = new(len);
          389         memcpy(p->buf, num->buf, len);
          390         p->wp = p->buf + len;
          391         p->rp = p->buf;
          392         p->scale = num->scale;
          393 
          394         return p;
          395 }
          396 
          397 static int
          398 negative(Num *num)
          399 {
          400         return last(num) && peek(num) == -1;
          401 }
          402 
          403 static Num *
          404 norm(Num *n)
          405 {
          406         /* trailing 0 */
          407         for (last(n); peek(n) == 0; numtrunc(n))
          408                 ;
          409 
          410         if (negative(n)) {
          411                 for (prev(n); peek(n) == 99; prev(n)) {
          412                         poke(n, -1);
          413                         next(n);
          414                         numtrunc(n);
          415                 }
          416         }
          417 
          418         /* adjust scale for 0 case*/
          419         if (length(n) == 0)
          420                 n->scale = 0;
          421         return n;
          422 }
          423 
          424 static Num *
          425 mulnto(Num *src, Num *dst, int n)
          426 {
          427         div_t dd;
          428         int d, carry;
          429 
          430         first(dst);
          431         numtrunc(dst);
          432 
          433         carry = 0;
          434         for (first(src); more(src); next(src)) {
          435                 d = peek(src) * n + carry;
          436                 dd = div(d, 100);
          437                 carry = dd.quot;
          438                 wrdigit(dst, dd.rem);
          439         }
          440 
          441         if (carry)
          442                 wrdigit(dst, carry);
          443         return dst;
          444 }
          445 
          446 static Num *
          447 muln(Num *num, int n)
          448 {
          449         div_t dd;
          450         int d, carry;
          451 
          452         carry = 0;
          453         for (first(num); more(num); next(num)) {
          454                 d = peek(num) * n + carry;
          455                 dd = div(d, 100);
          456                 carry = dd.quot;
          457                 poke(num, dd.rem);
          458         }
          459 
          460         if (carry)
          461                 wrdigit(num, carry);
          462         return num;
          463 }
          464 
          465 static int
          466 divn(Num *num, int n)
          467 {
          468         div_t dd;
          469         int val, carry;
          470 
          471         carry = 0;
          472         for (last(num); !begin(num); prev(num)) {
          473                 val = carry * 100 + peek(num);
          474                 dd = div(val, n);
          475                 poke(num, dd.quot);
          476                 carry = dd.rem;
          477         }
          478         norm(num);
          479 
          480         return carry;
          481 }
          482 
          483 static void
          484 div10(Num *num, int n)
          485 {
          486         div_t dd = div(n, 2);
          487 
          488         if (dd.rem == 1)
          489                 divn(num, 10);
          490 
          491         rshift(num, dd.quot);
          492 }
          493 
          494 static void
          495 mul10(Num *num, int n)
          496 {
          497         div_t dd = div(n, 2);
          498 
          499         if (dd.rem == 1)
          500                 muln(num, 10);
          501 
          502         lshift(num, dd.quot);
          503 }
          504 
          505 static void
          506 align(Num *a, Num *b)
          507 {
          508         int d;
          509         Num *max, *min;
          510 
          511         d = a->scale - b->scale;
          512         if (d == 0) {
          513                 return;
          514         } else if (d > 0) {
          515                 min = b;
          516                 max = a;
          517         } else {
          518                 min = a;
          519                 max = b;
          520         }
          521 
          522         d = abs(d);
          523         mul10(min, d);
          524         min->scale += d;
          525 
          526         assert(min->scale == max->scale);
          527 }
          528 
          529 static Num *
          530 addn(Num *num, int val)
          531 {
          532         int d, carry = val;
          533 
          534         for (first(num); carry; next(num)) {
          535                 d = more(num) ? peek(num) : 0;
          536                 d += carry;
          537                 carry = 0;
          538 
          539                 if (d >= 100) {
          540                         carry = 1;
          541                         d -= 100;
          542                 }
          543                 poke(num, d);
          544         }
          545 
          546         return num;
          547 }
          548 
          549 static Num *
          550 reverse(Num *num)
          551 {
          552         int d;
          553         signed char *p, *q;
          554 
          555         for (p = num->buf, q = num->wp - 1; p < q; ++p, --q) {
          556                 d = *p;
          557                 *p = *q;
          558                 *q = d;
          559         }
          560 
          561         return num;
          562 }
          563 
          564 static Num *
          565 addnum(Num *a, Num *b)
          566 {
          567         Num *c;
          568         int len, alen, blen, carry, da, db, sum;
          569 
          570         align(a, b);
          571         alen = length(a);
          572         blen = length(b);
          573         len = (alen > blen) ? alen : blen;
          574         c = new(len);
          575 
          576         first(a);
          577         first(b);
          578         carry = 0;
          579         while (len-- > 0) {
          580                 da = (more(a)) ? rddigit(a) : 0;
          581                 db = (more(b)) ? rddigit(b) : 0;
          582 
          583                 sum = da + db + carry;
          584                 if (sum >= 100) {
          585                         carry = 1;
          586                         sum -= 100;
          587                 } else if (sum < 0) {
          588                         carry = -1;
          589                         sum += 100;
          590                 } else {
          591                         carry = 0;
          592                 }
          593 
          594                 wrdigit(c, sum);
          595         }
          596 
          597         if (carry)
          598                 wrdigit(c, carry);
          599         c->scale = a->scale;
          600 
          601         return norm(c);
          602 }
          603 
          604 static Num *
          605 subnum(Num *a, Num *b)
          606 {
          607         Num *tmp, *sum;
          608 
          609         tmp = chsign(copy(b));
          610         sum = addnum(a, tmp);
          611         freenum(tmp);
          612 
          613         return sum;
          614 }
          615 
          616 static Num *
          617 mulnum(Num *a, Num *b)
          618 {
          619         Num shadow, *c, *ca, *cb;
          620         int pos, prod, carry, dc, db, da, sc;
          621         int asign = negative(a), bsign = negative(b);
          622 
          623         c = zeronum(length(a) + length(b) + 1);
          624         c->scale = a->scale + b->scale;
          625         sc = (a->scale > b->scale) ? a->scale : b->scale;
          626 
          627         ca = a;
          628         if (asign)
          629                 ca = chsign(copy(ca));
          630         cb = b;
          631         if (bsign)
          632                 cb = chsign(copy(cb));
          633 
          634         /* avoid aliasing problems when called from expnum*/
          635         if (ca == cb) {
          636                 shadow = *cb;
          637                 b = cb = &shadow;
          638         }
          639 
          640         for (first(cb); more(cb); next(cb)) {
          641                 div_t d;
          642 
          643                 carry = 0;
          644                 db = peek(cb);
          645 
          646                 pos = tell(c);
          647                 for (first(ca); more(ca); next(ca)) {
          648                         da = peek(ca);
          649                         dc = peek(c);
          650                         prod = da * db + dc + carry;
          651                         d = div(prod, 100);
          652                         carry = d.quot;
          653                         poke(c, d.rem);
          654                         next(c);
          655                 }
          656 
          657                 for ( ; carry > 0; carry = d.quot) {
          658                         dc = peek(c) + carry;
          659                         d = div(dc, 100);
          660                         poke(c, d.rem);
          661                         next(c);
          662                 }
          663                 seek(c, pos + 1);
          664         }
          665         norm(c);
          666 
          667         if (sc < scale)
          668                 sc = scale;
          669         sc = c->scale - sc;
          670         if (sc > 0) {
          671                 div10(c, sc);
          672                 c->scale -= sc;
          673         }
          674 
          675         if (ca != a)
          676                 freenum(ca);
          677         if (cb != b)
          678                 freenum(cb);
          679 
          680         if (asign ^ bsign)
          681                 chsign(c);
          682         return c;
          683 }
          684 
          685 /*
          686  * The divmod function is implemented following the algorithm
          687  * from the plan9 version that is not exactly like the one described
          688  * in the paper. A lot of magic here.
          689  */
          690 static Num *
          691 divmod(Num *odivd, Num *odivr, Num **remp)
          692 {
          693         Num *acc, *divd, *divr, *res;
          694         int divsign, remsign;
          695         int under, magic, ndig, diff;
          696         int d, q, carry, divcarry;
          697         long dr, dd, cc;
          698 
          699         divr = odivr;
          700         acc = copy(&zero);
          701         divd = copy(odivd);
          702         res = zeronum(length(odivd));
          703 
          704         under = divcarry = divsign = remsign = 0;
          705 
          706         if (length(divr) == 0) {
          707                 weprintf("divide by 0\n");
          708                 goto ret;
          709         }
          710 
          711         divsign = negative(divd);
          712         if (divsign)
          713                 chsign(divd);
          714 
          715         remsign = negative(divr);
          716         if (remsign)
          717                 divr = chsign(copy(divr));
          718 
          719         diff = length(divd) - length(divr);
          720 
          721         seek(res, diff + 1);
          722         last(divd);
          723         last(divr);
          724 
          725         wrdigit(divd, 0);
          726 
          727         dr = peek(divr);
          728         magic = dr < 10;
          729         dr = dr * 100 + (prev(divr) ? peek(divr) : 0);
          730         if (magic) {
          731                 dr = dr * 100 + (prev(divr) ? peek(divr) : 0);
          732                 dr *= 2;
          733                 dr /= 25;
          734         }
          735 
          736         for (ndig = 0; diff >= 0; ++ndig) {
          737                 last(divd);
          738                 dd = peek(divd);
          739                 dd = dd * 100 + (prev(divd) ? peek(divd) : 0);
          740                 dd = dd * 100 + (prev(divd) ? peek(divd) : 0);
          741                 cc = dr;
          742 
          743                 if (diff == 0)
          744                         dd++;
          745                 else
          746                         cc++;
          747 
          748                 if (magic)
          749                         dd *= 8;
          750 
          751                 q = dd / cc;
          752                 under = 0;
          753                 if (q > 0 && dd % cc < 8 && magic) {
          754                         q--;
          755                         under = 1;
          756                 }
          757 
          758                 mulnto(divr, acc, q);
          759 
          760                 /* Subtract acc from dividend at offset position */
          761                 first(acc);
          762                 carry = 0;
          763                 for (seek(divd, diff); more(divd); next(divd)) {
          764                         d = peek(divd);
          765                         d = d - (more(acc) ? rddigit(acc) : 0) - carry;
          766                         carry = 0;
          767                         if (d < 0) {
          768                                 d += 100;
          769                                 carry = 1;
          770                         }
          771                         poke(divd, d);
          772                 }
          773                 divcarry = carry;
          774 
          775                 /* Store quotient digit */
          776                 prev(res);
          777                 poke(res, q);
          778 
          779                 /* Handle borrow propagation */
          780                 last(divd);
          781                 d = peek(divd);
          782                 if ((d != 0) && (diff != 0)) {
          783                         prev(divd);
          784                         d = peek(divd) + 100;
          785                         poke(divd, d);
          786                 }
          787 
          788                 /* Shorten dividend for next iteration */
          789                 if (--diff >= 0)
          790                         divd->wp--;
          791         }
          792 
          793         /*
          794          * if we have an underflow then we have to adjust
          795          * the remaining and the result
          796          */
          797         if (under) {
          798                 Num *p = subnum(divd, divr);
          799                 if (negative(p)) {
          800                         freenum(p);
          801                 } else {
          802                         freenum(divd);
          803                         poke(res, q + 1);
          804                         divd = p;
          805                 }
          806         }
          807 
          808         if (divcarry) {
          809                 Num *p;
          810 
          811                 poke(res, q - 1);
          812                 poke(divd, -1);
          813                 p = addnum(divr, divd);
          814                 freenum(divd);
          815                 divd = p;
          816         }
          817 
          818         divcarry = 0;
          819         for (first(res); more(res); next(res)) {
          820                 d = peek(res) + divcarry;
          821                 divcarry = 0;
          822                 if (d >= 100) {
          823                         d -= 100;
          824                         divcarry = 1;
          825                 }
          826                 poke(res, d);
          827         }
          828 
          829 ret:
          830         if (divsign)
          831                 chsign(divd);
          832         if (divsign ^ remsign)
          833                 chsign(res);
          834 
          835         if (remp) {
          836                 divd->scale = odivd->scale;
          837                 *remp = norm(divd);
          838         } else {
          839                 freenum(divd);
          840         }
          841 
          842         if (divr != odivr)
          843                 freenum(divr);
          844 
          845         freenum(acc);
          846 
          847         res->scale = odivd->scale - odivr->scale;
          848         if (res->scale < 0)
          849                 res->scale = 0;
          850 
          851         return norm(res);
          852 }
          853 
          854 static int
          855 divscale(Num *divd, Num *divr)
          856 {
          857         int diff;
          858 
          859         if (length(divr) == 0) {
          860                 weprintf("divide by 0\n");
          861                 return 0;
          862         }
          863 
          864         diff = scale + divr->scale - divd->scale;
          865 
          866         if (diff > 0) {
          867                 mul10(divd, diff);
          868                 divd->scale += diff;
          869         } else if (diff < 0) {
          870                 mul10(divr, -diff);
          871                 divr->scale += -diff;
          872         }
          873 
          874         return 1;
          875 }
          876 
          877 static Num *
          878 divnum(Num *a, Num *b)
          879 {
          880         Num *r;
          881         int siga, sigb;
          882 
          883         siga = negative(a);
          884         if (siga)
          885                 chsign(a);
          886 
          887         sigb = negative(b);
          888         if (sigb)
          889                 chsign(b);
          890 
          891         if (!divscale(a, b))
          892                 return copy(&zero);
          893 
          894         r = divmod(a, b, NULL);
          895         if (siga ^ sigb)
          896                 chsign(r);
          897         return r;
          898 }
          899 
          900 static Num *
          901 modnum(Num *a, Num *b)
          902 {
          903         Num *mod, *c;
          904         int siga, sigb;
          905 
          906         siga = negative(a);
          907         if (siga)
          908                 chsign(a);
          909 
          910         sigb = negative(b);
          911         if (sigb)
          912                 chsign(b);
          913 
          914         if (!divscale(a, b))
          915                 return copy(&zero);
          916 
          917         c = divmod(a, b, &mod);
          918         freenum(c);
          919 
          920         if (siga)
          921                 chsign(mod);
          922 
          923         return mod;
          924 }
          925 
          926 static Num *
          927 expnum(Num *base, Num *exp)
          928 {
          929         int neg, d;
          930         Num *res, *fact, *e, *tmp1, *tmp2;
          931 
          932         res = copy(&one);
          933         if (length(exp) == 0)
          934                 return res;
          935 
          936         e = copy(exp);
          937         if ((neg = negative(exp)) != 0)
          938                 chsign(e);
          939 
          940         if (e->scale > 0) {
          941                 div10(e, e->scale);
          942                 e->scale = 0;
          943         }
          944 
          945         fact = copy(base);
          946         while (length(e) > 0) {
          947                 first(e);
          948                 d = peek(e);
          949                 if (d % 2 == 1) {
          950                         tmp1 = mulnum(res, fact);
          951                         freenum(res);
          952                         res = tmp1;
          953                 }
          954 
          955                 /* Square fact */
          956                 tmp1 = mulnum(fact, fact);
          957                 freenum(fact);
          958                 fact = tmp1;
          959 
          960                 divn(e, 2);
          961         }
          962         freenum(fact);
          963         freenum(e);
          964 
          965         /* Handle negative exponent: 1 / res */
          966         if (neg) {
          967                 tmp2 = divnum(tmp1 = copy(&one), res);
          968                 freenum(tmp1);
          969                 freenum(res);
          970                 res = tmp2;
          971         }
          972 
          973         return res;
          974 }
          975 
          976 /*
          977  * Compare two numbers: returns <0 if a<b, 0 if a==b, >0 if a>b
          978  */
          979 static int
          980 cmpnum(Num *a, Num *b)
          981 {
          982         Num *diff;
          983         int result;
          984 
          985         diff = subnum(a, b);
          986         if (length(diff) == 0)
          987                 result = 0;
          988         else if (negative(diff))
          989                 result = -1;
          990         else
          991                 result = 1;
          992         freenum(diff);
          993 
          994         return result;
          995 }
          996 
          997 /*
          998  * Integer square root of a small integer (0-9999)
          999  * Used for initial guess in Newton's method
         1000  */
         1001 static int
         1002 isqrt(int n)
         1003 {
         1004         int x, x1;
         1005 
         1006         if (n <= 0)
         1007                 return 0;
         1008         if (n == 1)
         1009                 return 1;
         1010 
         1011         x = n;
         1012         x1 = (x + 1) / 2;
         1013         while (x1 < x) {
         1014                 x = x1;
         1015                 x1 = (x + n / x) / 2;
         1016         }
         1017         return x;
         1018 }
         1019 
         1020 /*
         1021  * Square root using Newton's method: x_{n+1} = (x_n + y/x_n) / 2
         1022  *
         1023  * Key insight: sqrt(a * 10^(2n)) = sqrt(a) * 10^n
         1024  * So we scale up the input to get the desired output precision.
         1025  *
         1026  * To compute sqrt with scale decimal places of precision:
         1027  * 1. Scale up y by 10^(2*scale + 2) (extra 2 for guard digits)
         1028  * 2. Compute integer sqrt
         1029  * 3. Result has (scale + 1) decimal places, numtrunc to scale
         1030  */
         1031 static Num *
         1032 sqrtnum(Num *oy)
         1033 {
         1034         Num *y, *x, *xprev, *q, *sum;
         1035         int top, ysc, iter;
         1036 
         1037         if (length(oy) == 0)
         1038                 return copy(&zero);
         1039 
         1040         if (negative(oy)) {
         1041                 weprintf("square root of negative number\n");
         1042                 return copy(&zero);
         1043         }
         1044 
         1045         y = copy(oy);
         1046         ysc = 2 * scale + 2 - y->scale;
         1047         if (ysc > 0)
         1048                 mul10(y, ysc);
         1049         ysc = 2 * scale + 2;
         1050 
         1051         /* Make scale even (so sqrt gives integer result) */
         1052         if (ysc % 2 == 1) {
         1053                 muln(y, 10);
         1054                 ysc++;
         1055         }
         1056         y->scale = 0;
         1057 
         1058         last(y);
         1059         top = peek(y);
         1060         if (prev(y) && length(y) > 1)
         1061                 top = top * 100 + peek(y);
         1062 
         1063         x = new(0);
         1064         wrdigit(x, isqrt(top));
         1065         x->scale = 0;
         1066 
         1067         /* Scale up the initial guess to match the magnitude of y */
         1068         lshift(x, (length(y) - 1) / 2);
         1069 
         1070 
         1071         /* Newton iteration: x = (x + y/x) / 2 */
         1072         xprev = NULL;
         1073         for (iter = 0; iter < 1000; iter++) {
         1074                 q = divmod(y, x, NULL);
         1075                 sum = addnum(x, q);
         1076                 freenum(q);
         1077                 divn(sum, 2);
         1078 
         1079                 /* Check for convergence: sum == x or sum == prev */
         1080                 if (cmpnum(sum, x) == 0) {
         1081                         freenum(sum);
         1082                         break;
         1083                 }
         1084                 if (xprev != NULL && cmpnum(sum, xprev) == 0) {
         1085                         /* Oscillating - pick smaller */
         1086                         if (cmpnum(x, sum) < 0) {
         1087                                 freenum(sum);
         1088                         } else {
         1089                                 freenum(x);
         1090                                 x = sum;
         1091                         }
         1092                         break;
         1093                 }
         1094 
         1095                 freenum(xprev);
         1096                 xprev = x;
         1097                 x = sum;
         1098         }
         1099         freenum(xprev);
         1100         freenum(y);
         1101 
         1102         /* Truncate to desired scale */
         1103         x->scale = ysc / 2;
         1104         if (x->scale > scale) {
         1105                 int diff = x->scale - scale;
         1106                 div10(x, diff);
         1107                 x->scale = scale;
         1108         }
         1109 
         1110         return norm(x);
         1111 }
         1112 
         1113 static Num *
         1114 tonum(void)
         1115 {
         1116         char *s, *t, *end, *dot;
         1117         Num *num, *denom, *numer, *frac, *q, *rem;
         1118         int sign, d, ch, nfrac;
         1119 
         1120         s = input->s;
         1121         num = new(0);
         1122         sign = 0;
         1123         if (*s == '_') {
         1124                 sign = 1;
         1125                 ++s;
         1126         }
         1127 
         1128         dot = NULL;
         1129         for (t = s; (ch = *t) > 0 || ch <= UCHAR_MAX; ++t) {
         1130                 if (!strchr(digits, ch))
         1131                         break;
         1132                 if (ch == '.') {
         1133                         if (dot)
         1134                                 break;
         1135                         dot = t;
         1136                 }
         1137         }
         1138         input->s = end = t;
         1139 
         1140         /*
         1141          * Parse integer part: process digits left-to-right.
         1142          * For each digit: num = num * ibase + digit
         1143          */
         1144         for (t = s; t < (dot ? dot : end); ++t) {
         1145                 d = strchr(digits, *t) - digits;
         1146                 muln(num, ibase);
         1147                 addn(num, d);
         1148         }
         1149         norm(num);
         1150 
         1151         if (!dot)
         1152                 goto ret;
         1153 
         1154         /*
         1155          * convert fractional digits
         1156          * Algorithm: For digits d[0], d[1], ..., d[n-1] after '.'
         1157          * Value = d[0]/ibase + d[1]/ibase^2 + ... + d[n-1]/ibase^n
         1158          *
         1159          * numerator = d[0]*ibase^(n-1) + d[1]*ibase^(n-2) + ... + d[n-1]
         1160          * denominator = ibase^n
         1161          * Then extract decimal digits by repeated: num*100/denom
         1162          */
         1163         denom = copy(&one);
         1164         numer = copy(&zero);
         1165         for (t = dot + 1; t < end; ++t) {
         1166                 d = strchr(digits, *t) - digits;
         1167                 muln(denom, ibase);
         1168                 muln(numer, ibase);
         1169                 addn(numer, d);
         1170         }
         1171 
         1172         nfrac = end - dot - 1;
         1173         frac = new(0);
         1174         d = 0;
         1175         while (frac->scale < nfrac || length(numer) > 0) {
         1176                 muln(numer, 100);
         1177                 q = divmod(numer, denom, &rem);
         1178                 freenum(numer);
         1179 
         1180                 d = first(q) ? peek(q) : 0;
         1181                 wrdigit(frac, d);
         1182                 freenum(q);
         1183                 numer = rem;
         1184                 frac->scale += 2;
         1185         }
         1186         reverse(frac);
         1187 
         1188         /* Trim to exact input scale for odd nfrac */
         1189         if (frac->scale > nfrac && d % 10 == 0) {
         1190                 divn(frac, 10);
         1191                 frac->scale--;
         1192         }
         1193 
         1194         freenum(numer);
         1195         freenum(denom);
         1196 
         1197         q = addnum(num, frac);
         1198         freenum(num);
         1199         freenum(frac);
         1200         num = q;
         1201 
         1202 ret:
         1203         if (sign)
         1204                 chsign(num);
         1205         return num;
         1206 }
         1207 
         1208 static void
         1209 prchr(int ch)
         1210 {
         1211         if (col >= 69) {
         1212                 putchar('\\');
         1213                 putchar('\n');
         1214                 col = 0;
         1215         }
         1216         putchar(ch);
         1217         col++;
         1218 }
         1219 
         1220 static void
         1221 printd(int d, int base, int space)
         1222 {
         1223         int w, n;
         1224 
         1225         if (base <= 16) {
         1226                 prchr(digits[d]);
         1227         } else {
         1228                 if (space)
         1229                         prchr(' ');
         1230 
         1231                 for (w = 1, n = base - 1; n >= 10; n /= 10)
         1232                         w++;
         1233 
         1234                 if (col + w > 69) {
         1235                         putchar('\\');
         1236                         putchar('\n');
         1237                         col = 0;
         1238                 }
         1239                 col += printf("%0*d", w, d);
         1240         }
         1241 }
         1242 
         1243 static void
         1244 pushdigit(struct digit **l, int val)
         1245 {
         1246         struct digit *it = emalloc(sizeof(*it));
         1247 
         1248         it->next = *l;
         1249         it->val = val;
         1250         *l = it;
         1251 }
         1252 
         1253 static int
         1254 popdigit(struct digit **l)
         1255 {
         1256         int val;
         1257         struct digit *next, *it = *l;
         1258 
         1259         if (it == NULL)
         1260                 return -1;
         1261 
         1262         val = it->val;
         1263         next = it->next;
         1264         free(it);
         1265         *l = next;
         1266         return val;
         1267 }
         1268 
         1269 static void
         1270 printnum(Num *onum, int base)
         1271 {
         1272         struct digit *sp;
         1273         int sc, i, sign, n;
         1274         Num *num, *inte, *frac, *opow;
         1275 
         1276         col = 0;
         1277         if (length(onum) == 0) {
         1278                 prchr('0');
         1279                 return;
         1280         }
         1281 
         1282         num = copy(onum);
         1283         if ((sign = negative(num)) != 0)
         1284                 chsign(num);
         1285 
         1286         sc = num->scale;
         1287         if (num->scale % 2 == 1) {
         1288                 muln(num, 10);
         1289                 num->scale++;
         1290         }
         1291         inte = copy(num);
         1292         rshift(inte, num->scale / 2);
         1293         inte->scale = 0;
         1294         frac = subnum(num, inte);
         1295 
         1296         sp = NULL;
         1297         for (i = 0; length(inte) > 0; ++i)
         1298                 pushdigit(&sp, divn(inte, base));
         1299         if (sign)
         1300                 prchr('-');
         1301         while (i-- > 0)
         1302                 printd(popdigit(&sp), base, 1);
         1303         assert(sp == NULL);
         1304 
         1305         if (num->scale == 0)
         1306                 goto ret;
         1307 
         1308         /*
         1309          * Print fractional part by repeated multiplication by base.
         1310          * We maintain the fraction as: frac / 10^scale
         1311          *
         1312          * Algorithm:
         1313          * 1. Multiply frac by base
         1314          * 2. Output integer part (frac / 10^scale)
         1315          * 3. Keep fractional part (frac % 10^scale)
         1316          */
         1317         prchr('.');
         1318 
         1319         opow = copy(&one);
         1320         mul10(opow, num->scale);
         1321 
         1322         for (n = 0; n < sc; ++n) {
         1323                 int d;
         1324                 Num *q, *rem;
         1325 
         1326                 muln(frac, base);
         1327                 q = divmod(frac, opow, &rem);
         1328                 d = first(q) ? peek(q) : 0;
         1329                 freenum(frac);
         1330                 freenum(q);
         1331                 frac = rem;
         1332                 printd(d, base, n > 0);
         1333         }
         1334         freenum(opow);
         1335 
         1336 ret:
         1337         freenum(num);
         1338         freenum(inte);
         1339         freenum(frac);
         1340 }
         1341 
         1342 static int
         1343 moreinput(void)
         1344 {
         1345         struct input *ip;
         1346 
         1347 repeat:
         1348         if (!input)
         1349                 return 0;
         1350 
         1351         if (input->buf != NULL && *input->s != '\0')
         1352                 return 1;
         1353 
         1354         if (input->fp) {
         1355                 if (getline(&input->buf, &input->n, input->fp) >= 0) {
         1356                         input->s = input->buf;
         1357                         return 1;
         1358                 }
         1359                 if (ferror(input->fp)) {
         1360                         eprintf("reading from file:");
         1361                         exit(1);
         1362                 }
         1363                 fclose(input->fp);
         1364         }
         1365 
         1366         ip = input;
         1367         input = ip->next;
         1368         free(ip->buf);
         1369         free(ip);
         1370         goto repeat;
         1371 }
         1372 
         1373 static void
         1374 addinput(FILE *fp, char *s)
         1375 {
         1376         struct input *ip;
         1377 
         1378         assert((!fp && !s) == 0);
         1379 
         1380         ip = emalloc(sizeof(*ip));
         1381         ip->next = input;
         1382         ip->fp = fp;
         1383         ip->n = 0;
         1384         ip->s = ip->buf = s;
         1385         input = ip;
         1386 }
         1387 
         1388 static void
         1389 delinput(int cmd, int n)
         1390 {
         1391         if (n < 0)
         1392                 error("Q command requires a number >= 0");
         1393         while (n-- > 0) {
         1394                 if (cmd == 'Q' && !input->next)
         1395                         error("Q command argument exceeded string execution depth");
         1396                 if (input->fp)
         1397                         fclose(input->fp);
         1398                 free(input->buf);
         1399                 input = input->next;
         1400                 if (!input)
         1401                         exit(0);
         1402         }
         1403 }
         1404 
         1405 static void
         1406 push(Val v)
         1407 {
         1408         Val *p = emalloc(sizeof(Val));
         1409 
         1410         *p = v;
         1411         p->next = stack;
         1412         stack = p;
         1413 }
         1414 
         1415 static void
         1416 needstack(int n)
         1417 {
         1418         Val *vp;
         1419 
         1420         for (vp = stack; n > 0 && vp; vp = vp->next)
         1421                 --n;
         1422         if (n > 0)
         1423                 error("stack empty");
         1424 }
         1425 
         1426 static Val
         1427 pop(void)
         1428 {
         1429         Val v;
         1430 
         1431         if (!stack)
         1432                 error("stack empty");
         1433         v = *stack;
         1434         free(stack);
         1435         stack = v.next;
         1436         v.next = NULL;
         1437 
         1438         return v;
         1439 }
         1440 
         1441 static Num *
         1442 popnum(void)
         1443 {
         1444         Val v = pop();
         1445 
         1446         if (v.type != NUM) {
         1447                 free(v.u.s);
         1448                 error("non-numeric value");
         1449         }
         1450         return v.u.n;
         1451 }
         1452 
         1453 static void
         1454 pushnum(Num *num)
         1455 {
         1456         push((Val) {.type = NUM, .u.n = num});
         1457 }
         1458 
         1459 static void
         1460 pushstr(char *s)
         1461 {
         1462         push((Val) {.type = STR, .u.s = s});
         1463 }
         1464 
         1465 static void
         1466 arith(Num *(*fn)(Num *, Num *))
         1467 {
         1468         Num *a, *b, *c;
         1469 
         1470         needstack(2);
         1471         b = popnum();
         1472         a = popnum();
         1473         c = (*fn)(a, b);
         1474         freenum(a);
         1475         freenum(b);
         1476         pushnum(c);
         1477 }
         1478 
         1479 static void
         1480 pushdivmod(void)
         1481 {
         1482         Num *a, *b, *q, *rem;
         1483 
         1484         needstack(2);
         1485         b = popnum();
         1486         a = popnum();
         1487 
         1488         if (!divscale(a, b)) {
         1489                 q = copy(&zero);
         1490                 rem = copy(&zero);
         1491         } else {
         1492                 q = divmod(a, b, &rem);
         1493         }
         1494 
         1495         pushnum(q);
         1496         pushnum(rem);
         1497         freenum(a);
         1498         freenum(b);
         1499 }
         1500 
         1501 static int
         1502 popint(void)
         1503 {
         1504         Num *num;
         1505         int r = -1, n, d;
         1506 
         1507         num = popnum();
         1508         if (negative(num))
         1509                 goto ret;
         1510 
         1511         /* discard fraction part */
         1512         div10(num, num->scale);
         1513 
         1514         n = 0;
         1515         for (last(num); !begin(num); prev(num)) {
         1516                 if (n > INT_MAX / 100)
         1517                         goto ret;
         1518                 n *= 100;
         1519                 d = peek(num);
         1520                 if (n > INT_MAX - d)
         1521                         goto ret;
         1522                 n += d;
         1523         }
         1524         r = n;
         1525 
         1526 ret:
         1527         freenum(num);
         1528         return r;
         1529 }
         1530 
         1531 static void
         1532 pushint(int n)
         1533 {
         1534         div_t dd;
         1535         Num *num;
         1536 
         1537         num = new(0);
         1538         for ( ; n > 0; n = dd.quot) {
         1539                 dd = div(n, 100);
         1540                 wrdigit(num, dd.rem);
         1541         }
         1542         pushnum(num);
         1543 }
         1544 
         1545 static void
         1546 printval(Val v)
         1547 {
         1548         if (v.type == STR)
         1549                 fputs(v.u.s, stdout);
         1550         else
         1551                 printnum(v.u.n, obase);
         1552 }
         1553 
         1554 static Val
         1555 dupval(Val v)
         1556 {
         1557         Val nv;
         1558 
         1559         switch (nv.type = v.type) {
         1560         case STR:
         1561                 nv.u.s = estrdup(v.u.s);
         1562                 break;
         1563         case NUM:
         1564                 nv.u.n = copy(v.u.n);
         1565                 break;
         1566         case NVAL:
         1567                 nv.type = NUM;
         1568                 nv.u.n = copy(&zero);
         1569                 break;
         1570         }
         1571         nv.next = NULL;
         1572 
         1573         return nv;
         1574 }
         1575 
         1576 static void
         1577 freeval(Val v)
         1578 {
         1579         if (v.type == STR)
         1580                 free(v.u.s);
         1581         else if (v.type == NUM)
         1582                 freenum(v.u.n);
         1583 }
         1584 
         1585 static void
         1586 dumpstack(void)
         1587 {
         1588         Val *vp;
         1589 
         1590         for (vp = stack; vp; vp = vp->next) {
         1591                 printval(*vp);
         1592                 putchar('\n');
         1593         }
         1594 }
         1595 
         1596 static void
         1597 clearstack(void)
         1598 {
         1599         Val *vp, *next;
         1600 
         1601         for (vp = stack; vp; vp = next) {
         1602                 next = vp->next;
         1603                 freeval(*vp);
         1604                 free(vp);
         1605         }
         1606         stack = NULL;
         1607 }
         1608 
         1609 static void
         1610 dupstack(void)
         1611 {
         1612         Val v;
         1613 
         1614         push(v = pop());
         1615         push(dupval(v));
         1616 }
         1617 
         1618 static void
         1619 deepstack(void)
         1620 {
         1621         int n;
         1622         Val *vp;
         1623 
         1624         n = 0;
         1625         for (vp = stack; vp; vp = vp->next) {
         1626                 if (n == INT_MAX)
         1627                         error("stack depth does not fit in a integer");
         1628                 ++n;
         1629         }
         1630         pushint(n);
         1631 }
         1632 
         1633 static void
         1634 pushfrac(void)
         1635 {
         1636         Val v = pop();
         1637 
         1638         if (v.type == STR)
         1639                 pushint(0);
         1640         else
         1641                 pushint(v.u.n->scale);
         1642         freeval(v);
         1643 }
         1644 
         1645 static void
         1646 pushlen(void)
         1647 {
         1648         int n;
         1649         Num *num;
         1650         Val v = pop();
         1651 
         1652         if (v.type == STR) {
         1653                 n = strlen(v.u.s);
         1654         } else {
         1655                 num = v.u.n;
         1656                 if (length(num) == 0) {
         1657                         n = 1;
         1658                 } else {
         1659                         n = length(num) * 2;
         1660                         n -= last(num) ? peek(num) < 10 : 0;
         1661                 }
         1662         }
         1663         pushint(n);
         1664         freeval(v);
         1665 }
         1666 
         1667 static void
         1668 setibase(void)
         1669 {
         1670         int n = popint();
         1671 
         1672         if (n < 2 || n > 16)
         1673                 error("input base must be an integer between 2 and 16");
         1674         ibase = n;
         1675 }
         1676 
         1677 static void
         1678 setobase(void)
         1679 {
         1680         int n = popint();
         1681 
         1682         if (n < 2)
         1683                 error("output base must be an integer greater than 1");
         1684         obase = n;
         1685 }
         1686 
         1687 static char *
         1688 string(char *dst, int *np)
         1689 {
         1690         int n, ch;
         1691 
         1692         n = np ? *np : 0;
         1693         for (;;) {
         1694                 ch = *input->s++;
         1695 
         1696                 switch (ch) {
         1697                 case '\0':
         1698                         if (!moreinput())
         1699                                 exit(0);
         1700                         break;
         1701                 case '\\':
         1702                         if (*input->s == '[') {
         1703                                 dst = erealloc(dst, n + 1);
         1704                                 dst[n++] = *input->s++;
         1705                                 break;
         1706                         }
         1707                         goto copy;
         1708                 case ']':
         1709                         if (!np) {
         1710                                 dst = erealloc(dst, n + 1);
         1711                                 dst[n] = '\0';
         1712                                 return dst;
         1713                         }
         1714                 case '[':
         1715                 default:
         1716                 copy:
         1717                         dst = erealloc(dst, n + 1);
         1718                         dst[n++] = ch;
         1719                         if (ch == '[')
         1720                                 dst = string(dst, &n);
         1721                         if (ch == ']') {
         1722                                 *np = n;
         1723                                 return dst;
         1724                         }
         1725                 }
         1726         }
         1727 }
         1728 
         1729 static void
         1730 setscale(void)
         1731 {
         1732         int n = popint();
         1733 
         1734         if (n < 0)
         1735                 error("scale must be a nonnegative integer");
         1736         scale = n;
         1737 }
         1738 
         1739 static unsigned
         1740 hash(char *name)
         1741 {
         1742         int c;
         1743         unsigned h = 5381;
         1744 
         1745         while (c = *name++)
         1746                 h = h*33 ^ c;
         1747 
         1748         return h;
         1749 }
         1750 
         1751 static struct reg *
         1752 lookup(char *name)
         1753 {
         1754         struct reg *rp;
         1755         int h = hash(name) & HASHSIZ-1;
         1756 
         1757         for (rp = htbl[h]; rp; rp = rp->next) {
         1758                 if (strcmp(name, rp->name) == 0)
         1759                         return rp;
         1760         }
         1761 
         1762         rp = emalloc(sizeof(*rp));
         1763         rp->next = htbl[h];
         1764         htbl[h] = rp;
         1765         rp->name = estrdup(name);
         1766 
         1767         rp->val.type = NVAL;
         1768         rp->val.next = NULL;
         1769 
         1770         rp->ary.n = 0;
         1771         rp->ary.buf = NULL;
         1772         rp->ary.next = NULL;
         1773 
         1774         return rp;
         1775 }
         1776 
         1777 static char *
         1778 regname(void)
         1779 {
         1780         int delim, ch;
         1781         char *s;
         1782         static char name[REGSIZ];
         1783 
         1784         ch = *input->s++;
         1785         if (!iflag || ch != '<' && ch != '"') {
         1786                 name[0] = ch;
         1787                 name[1] = '\0';
         1788                 return name;
         1789         }
         1790 
         1791         if ((delim = ch) == '<')
         1792                 delim = '>';
         1793 
         1794         for (s = name; s < &name[REGSIZ]; ++s) {
         1795                 ch = *input->s++;
         1796                 if (ch == '\0' ||  ch == delim) {
         1797                         *s = '\0';
         1798                         if (ch == '>') {
         1799                                 name[0] = atoi(name);
         1800                                 name[1] = '\0';
         1801                         }
         1802                         return name;
         1803                 }
         1804                 *s = ch;
         1805         }
         1806 
         1807         error("identifier too long");
         1808 }
         1809 
         1810 static void
         1811 popreg(void)
         1812 {
         1813         int i;
         1814         Val *vnext;
         1815         struct ary *anext;
         1816         char *s = regname();
         1817         struct reg *rp = lookup(s);
         1818 
         1819         if (rp->val.type == NVAL)
         1820                 error("stack register '%s' (%o) is empty", s, s[0]);
         1821 
         1822         push(rp->val);
         1823         vnext = rp->val.next;
         1824         if (!vnext) {
         1825                 rp->val.type = NVAL;
         1826         } else {
         1827                 rp->val = *vnext;
         1828                 free(vnext);
         1829         }
         1830 
         1831         for (i = 0; i < rp->ary.n; ++i)
         1832                 freeval(rp->ary.buf[i]);
         1833         free(rp->ary.buf);
         1834 
         1835         anext = rp->ary.next;
         1836         if (!anext) {
         1837                 rp->ary.n = 0;
         1838                 rp->ary.buf = NULL;
         1839         } else {
         1840                 rp->ary = *anext;
         1841                 free(anext);
         1842         }
         1843 }
         1844 
         1845 static void
         1846 pushreg(void)
         1847 {
         1848         Val v;
         1849         Val *vp;
         1850         struct ary *ap;
         1851         struct reg *rp = lookup(regname());
         1852 
         1853         v = pop();
         1854 
         1855         vp = emalloc(sizeof(Val));
         1856         *vp = rp->val;
         1857         rp->val = v;
         1858         rp->val.next = vp;
         1859 
         1860         ap = emalloc(sizeof(struct ary));
         1861         *ap = rp->ary;
         1862         rp->ary.n = 0;
         1863         rp->ary.buf = NULL;
         1864         rp->ary.next = ap;
         1865 }
         1866 
         1867 static Val *
         1868 aryidx(void)
         1869 {
         1870         int n;
         1871         int i;
         1872         Val *vp;
         1873         struct reg *rp = lookup(regname());
         1874         struct ary *ap = &rp->ary;
         1875 
         1876         n = popint();
         1877         if (n < 0)
         1878                 error("array index must fit in a positive integer");
         1879 
         1880         if (n >= ap->n) {
         1881                 ap->buf = ereallocarray(ap->buf, n+1, sizeof(Val));
         1882                 for (i = ap->n; i <= n; ++i)
         1883                         ap->buf[i].type = NVAL;
         1884                 ap->n = n + 1;
         1885         }
         1886         return &ap->buf[n];
         1887 }
         1888 
         1889 static void
         1890 aryget(void)
         1891 {
         1892         Val *vp = aryidx();
         1893 
         1894         push(dupval(*vp));
         1895 }
         1896 
         1897 static void
         1898 aryset(void)
         1899 {
         1900         Val val, *vp = aryidx();
         1901 
         1902         val = pop();
         1903         freeval(*vp);
         1904         *vp = val;
         1905 }
         1906 
         1907 static void
         1908 execmacro(void)
         1909 {
         1910         int ch;
         1911         Val v = pop();
         1912 
         1913         assert(v.type != NVAL);
         1914         if (v.type == NUM) {
         1915                 push(v);
         1916                 return;
         1917         }
         1918 
         1919         if (input->fp) {
         1920                 addinput(NULL, v.u.s);
         1921                 return;
         1922         }
         1923 
         1924         for (ch = *input->s; ch > 0 && ch <= UCHAR_MAX; ch = *input->s) {
         1925                 if (!isspace(ch))
         1926                         break;
         1927                 ++input->s;
         1928         }
         1929 
         1930         /* check for tail recursion */
         1931         if (ch == '\0' && strcmp(input->buf, v.u.s) == 0) {
         1932                 free(input->buf);
         1933                 input->buf = input->s = v.u.s;
         1934                 return;
         1935         }
         1936 
         1937         addinput(NULL, v.u.s);
         1938 }
         1939 
         1940 static void
         1941 relational(int ch)
         1942 {
         1943         int r;
         1944         char *s;
         1945         Num *a, *b;
         1946         struct reg *rp = lookup(regname());
         1947 
         1948         needstack(2);
         1949         a = popnum();
         1950         b = popnum();
         1951         r = cmpnum(a, b);
         1952         freenum(a);
         1953         freenum(b);
         1954 
         1955         switch (ch) {
         1956         case '>':
         1957                 r = r > 0;
         1958                 break;
         1959         case '<':
         1960                 r = r < 0;
         1961                 break;
         1962         case '=':
         1963                 r = r == 0;
         1964                 break;
         1965         case LE:
         1966                 r = r <= 0;
         1967                 break;
         1968         case GE:
         1969                 r = r >= 0;
         1970                 break;
         1971         case NE:
         1972                 r = r != 0;
         1973                 break;
         1974         default:
         1975                 abort();
         1976         }
         1977 
         1978         if (!r)
         1979                 return;
         1980 
         1981         push(dupval(rp->val));
         1982         execmacro();
         1983 }
         1984 
         1985 static void
         1986 printbytes(void)
         1987 {
         1988         Num *num;
         1989         Val v = pop();
         1990 
         1991         if (v.type == STR) {
         1992                 fputs(v.u.s, stdout);
         1993         } else {
         1994                 num = v.u.n;
         1995                 div10(num, num->scale);
         1996                 num->scale = 0;
         1997                 printnum(num, 100);
         1998         }
         1999         freeval(v);
         2000 }
         2001 
         2002 static void
         2003 eval(void)
         2004 {
         2005         int ch;
         2006         char *s;
         2007         Num *num;
         2008         size_t siz;
         2009         Val v1, v2;
         2010         struct reg *rp;
         2011 
         2012         if (setjmp(env))
         2013                 return;
         2014 
         2015         for (s = input->s; (ch = *s) != '\0'; ++s) {
         2016                 if (ch < 0 || ch > UCHAR_MAX || !isspace(ch))
         2017                         break;
         2018         }
         2019         input->s = s + (ch != '\0');
         2020 
         2021         switch (ch) {
         2022         case '\0':
         2023                 break;
         2024         case 'n':
         2025                 v1 = pop();
         2026                 printval(v1);
         2027                 freeval(v1);
         2028                 break;
         2029         case 'p':
         2030                 v1 = pop();
         2031                 printval(v1);
         2032                 putchar('\n');
         2033                 push(v1);
         2034                 break;
         2035         case 'P':
         2036                 printbytes();
         2037                 break;
         2038         case 'f':
         2039                 dumpstack();
         2040                 break;
         2041         case '+':
         2042                 arith(addnum);
         2043                 break;
         2044         case '-':
         2045                 arith(subnum);
         2046                 break;
         2047         case '*':
         2048                 arith(mulnum);
         2049                 break;
         2050         case '/':
         2051                 arith(divnum);
         2052                 break;
         2053         case '%':
         2054                 arith(modnum);
         2055                 break;
         2056         case '^':
         2057                 arith(expnum);
         2058                 break;
         2059         case '~':
         2060                 pushdivmod();
         2061                 break;
         2062         case 'v':
         2063                 num = popnum();
         2064                 pushnum(sqrtnum(num));
         2065                 freenum(num);
         2066                 break;
         2067         case 'c':
         2068                 clearstack();
         2069                 break;
         2070         case 'd':
         2071                 dupstack();
         2072                 break;
         2073         case 'r':
         2074                 needstack(2);
         2075                 v1 = pop();
         2076                 v2 = pop();
         2077                 push(v1);
         2078                 push(v2);
         2079                 break;
         2080         case 'S':
         2081                 pushreg();
         2082                 break;
         2083         case 'L':
         2084                 popreg();
         2085                 break;
         2086         case 's':
         2087                 rp = lookup(regname());
         2088                 v1 = pop();
         2089                 freeval(rp->val);
         2090                 rp->val.u = v1.u;
         2091                 rp->val.type = v1.type;
         2092                 break;
         2093         case 'l':
         2094                 rp = lookup(regname());
         2095                 push(dupval(rp->val));
         2096                 break;
         2097         case 'i':
         2098                 setibase();
         2099                 break;
         2100         case 'o':
         2101                 setobase();
         2102                 break;
         2103         case 'k':
         2104                 setscale();
         2105                 break;
         2106         case 'I':
         2107                 pushint(ibase);
         2108                 break;
         2109         case 'O':
         2110                 pushint(obase);
         2111                 break;
         2112         case 'K':
         2113                 pushint(scale);
         2114                 break;
         2115         case '[':
         2116                 pushstr(string(NULL, NULL));
         2117                 break;
         2118         case 'x':
         2119                 execmacro();
         2120                 break;
         2121         case '!':
         2122                 switch (*input->s++) {
         2123                 case '<':
         2124                         ch = GE;
         2125                         break;
         2126                 case '>':
         2127                         ch = LE;
         2128                         break;
         2129                 case '=':
         2130                         ch = NE;
         2131                         break;
         2132                 default:
         2133                         system(input->s-1);
         2134                         goto discard;
         2135                 }
         2136         case '>':
         2137         case '<':
         2138         case '=':
         2139                 relational(ch);
         2140                 break;
         2141         case '?':
         2142                 s = NULL;
         2143                 if (getline(&s, &siz, stdin) > 0) {
         2144                         pushstr(s);
         2145                 } else {
         2146                         free(s);
         2147                         if (ferror(stdin))
         2148                                 eprintf("reading from file\n");
         2149                 }
         2150                 break;
         2151         case 'q':
         2152                 delinput('q', 2);
         2153                 break;
         2154         case 'Q':
         2155                 delinput('Q', popint());
         2156                 break;
         2157         case 'Z':
         2158                 pushlen();
         2159                 break;
         2160         case 'X':
         2161                 pushfrac();
         2162                 break;
         2163         case 'z':
         2164                 deepstack();
         2165                 break;
         2166         case '#':
         2167         discard:
         2168                 while (*input->s)
         2169                         ++input->s;
         2170                 break;
         2171         case ':':
         2172                 aryset();
         2173                 break;
         2174         case ';':
         2175                 aryget();
         2176                 break;
         2177         default:
         2178                 if (!strchr(digits, ch))
         2179                         error("'%c' (%#o) unimplemented", ch, ch);
         2180         case '_':
         2181                 --input->s;
         2182                 pushnum(tonum());
         2183                 break;
         2184         }
         2185 }
         2186 
         2187 static void
         2188 dc(char *fname)
         2189 {
         2190         FILE *fp;
         2191 
         2192         if (strcmp(fname, "-") == 0) {
         2193                 fp = stdin;
         2194         } else {
         2195                 if ((fp = fopen(fname, "r")) == NULL)
         2196                         eprintf("opening '%s':", fname);
         2197         }
         2198         addinput(fp, NULL);
         2199 
         2200         while (moreinput())
         2201                 eval();
         2202 
         2203         free(input);
         2204         input = NULL;
         2205 }
         2206 
         2207 static void
         2208 usage(void)
         2209 {
         2210         eprintf("usage: dc [-i] [file ...]\n");
         2211 }
         2212 
         2213 int
         2214 main(int argc, char *argv[])
         2215 {
         2216         ARGBEGIN {
         2217         case 'i':
         2218                 iflag = 1;
         2219                 break;
         2220         default:
         2221                 usage();
         2222         } ARGEND
         2223 
         2224         while (*argv)
         2225                 dc(*argv++);
         2226         dc("-");
         2227 
         2228         return 0;
         2229 }