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 (31011B)
       ---
            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         /* Adjust scale to max(a->scale, b->scale) while c is still positive */
          668         sc = c->scale - sc;
          669         if (sc > 0) {
          670                 div10(c, sc);
          671                 c->scale -= sc;
          672         }
          673 
          674         if (ca != a)
          675                 freenum(ca);
          676         if (cb != b)
          677                 freenum(cb);
          678 
          679         if (asign ^ bsign)
          680                 chsign(c);
          681         return c;
          682 }
          683 
          684 /*
          685  * The divmod function is implemented following the algorithm
          686  * from the plan9 version that is not exactly like the one described
          687  * in the paper. A lot of magic here.
          688  */
          689 static Num *
          690 divmod(Num *odivd, Num *odivr, Num **remp)
          691 {
          692         Num *acc, *divd, *divr, *res;
          693         int divsign, remsign;
          694         int under, magic, ndig, diff;
          695         int d, q, carry, divcarry;
          696         long dr, dd, cc;
          697 
          698         divr = odivr;
          699         acc = copy(&zero);
          700         divd = copy(odivd);
          701         res = zeronum(length(odivd));
          702 
          703         under = divcarry = divsign = remsign = 0;
          704 
          705         if (length(divr) == 0) {
          706                 weprintf("divide by 0\n");
          707                 goto ret;
          708         }
          709 
          710         divsign = negative(divd);
          711         if (divsign)
          712                 chsign(divd);
          713 
          714         remsign = negative(divr);
          715         if (remsign)
          716                 divr = chsign(copy(divr));
          717 
          718         diff = length(divd) - length(divr);
          719 
          720         seek(res, diff + 1);
          721         last(divd);
          722         last(divr);
          723 
          724         wrdigit(divd, 0);
          725 
          726         dr = peek(divr);
          727         magic = dr < 10;
          728         dr = dr * 100 + (prev(divr) ? peek(divr) : 0);
          729         if (magic) {
          730                 dr = dr * 100 + (prev(divr) ? peek(divr) : 0);
          731                 dr *= 2;
          732                 dr /= 25;
          733         }
          734 
          735         for (ndig = 0; diff >= 0; ++ndig) {
          736                 last(divd);
          737                 dd = peek(divd);
          738                 dd = dd * 100 + (prev(divd) ? peek(divd) : 0);
          739                 dd = dd * 100 + (prev(divd) ? peek(divd) : 0);
          740                 cc = dr;
          741 
          742                 if (diff == 0)
          743                         dd++;
          744                 else
          745                         cc++;
          746 
          747                 if (magic)
          748                         dd *= 8;
          749 
          750                 q = dd / cc;
          751                 under = 0;
          752                 if (q > 0 && dd % cc < 8 && magic) {
          753                         q--;
          754                         under = 1;
          755                 }
          756 
          757                 mulnto(divr, acc, q);
          758 
          759                 /* Subtract acc from dividend at offset position */
          760                 first(acc);
          761                 carry = 0;
          762                 for (seek(divd, diff); more(divd); next(divd)) {
          763                         d = peek(divd);
          764                         d = d - (more(acc) ? rddigit(acc) : 0) - carry;
          765                         carry = 0;
          766                         if (d < 0) {
          767                                 d += 100;
          768                                 carry = 1;
          769                         }
          770                         poke(divd, d);
          771                 }
          772                 divcarry = carry;
          773 
          774                 /* Store quotient digit */
          775                 prev(res);
          776                 poke(res, q);
          777 
          778                 /* Handle borrow propagation */
          779                 last(divd);
          780                 d = peek(divd);
          781                 if ((d != 0) && (diff != 0)) {
          782                         prev(divd);
          783                         d = peek(divd) + 100;
          784                         poke(divd, d);
          785                 }
          786 
          787                 /* Shorten dividend for next iteration */
          788                 if (--diff >= 0)
          789                         divd->wp--;
          790         }
          791 
          792         /*
          793          * if we have an underflow then we have to adjust
          794          * the remaining and the result
          795          */
          796         if (under) {
          797                 Num *p = subnum(divd, divr);
          798                 if (negative(p)) {
          799                         freenum(p);
          800                 } else {
          801                         freenum(divd);
          802                         poke(res, q + 1);
          803                         divd = p;
          804                 }
          805         }
          806 
          807         if (divcarry) {
          808                 Num *p;
          809 
          810                 poke(res, q - 1);
          811                 poke(divd, -1);
          812                 p = addnum(divr, divd);
          813                 freenum(divd);
          814                 divd = p;
          815         }
          816 
          817         divcarry = 0;
          818         for (first(res); more(res); next(res)) {
          819                 d = peek(res) + divcarry;
          820                 divcarry = 0;
          821                 if (d >= 100) {
          822                         d -= 100;
          823                         divcarry = 1;
          824                 }
          825                 poke(res, d);
          826         }
          827 
          828 ret:
          829         if (divsign)
          830                 chsign(divd);
          831         if (divsign ^ remsign)
          832                 chsign(res);
          833 
          834         if (remp) {
          835                 divd->scale = odivd->scale;
          836                 *remp = norm(divd);
          837         } else {
          838                 freenum(divd);
          839         }
          840 
          841         if (divr != odivr)
          842                 freenum(divr);
          843 
          844         freenum(acc);
          845 
          846         res->scale = odivd->scale - odivr->scale;
          847         if (res->scale < 0)
          848                 res->scale = 0;
          849 
          850         return norm(res);
          851 }
          852 
          853 static int
          854 divscale(Num *divd, Num *divr)
          855 {
          856         int diff;
          857 
          858         if (length(divr) == 0) {
          859                 weprintf("divide by 0\n");
          860                 return 0;
          861         }
          862 
          863         diff = scale + divr->scale - divd->scale;
          864 
          865         if (diff > 0) {
          866                 mul10(divd, diff);
          867                 divd->scale += diff;
          868         } else if (diff < 0) {
          869                 mul10(divr, -diff);
          870                 divr->scale += -diff;
          871         }
          872 
          873         return 1;
          874 }
          875 
          876 static Num *
          877 divnum(Num *a, Num *b)
          878 {
          879         if (!divscale(a, b))
          880                 return copy(&zero);
          881 
          882         return divmod(a, b, NULL);
          883 }
          884 
          885 static Num *
          886 modnum(Num *a, Num *b)
          887 {
          888         Num *mod, *c;
          889 
          890         if (!divscale(a, b))
          891                 return copy(&zero);
          892 
          893         c = divmod(a, b, &mod);
          894         freenum(c);
          895 
          896         return mod;
          897 }
          898 
          899 static Num *
          900 expnum(Num *base, Num *exp)
          901 {
          902         int neg, d;
          903         Num *res, *fact, *e, *tmp1, *tmp2;
          904 
          905         res = copy(&one);
          906         if (length(exp) == 0)
          907                 return res;
          908 
          909         e = copy(exp);
          910         if ((neg = negative(exp)) != 0)
          911                 chsign(e);
          912 
          913         if (e->scale > 0) {
          914                 div10(e, e->scale);
          915                 e->scale = 0;
          916         }
          917 
          918         fact = copy(base);
          919         while (length(e) > 0) {
          920                 first(e);
          921                 d = peek(e);
          922                 if (d % 2 == 1) {
          923                         tmp1 = mulnum(res, fact);
          924                         freenum(res);
          925                         res = tmp1;
          926                 }
          927 
          928                 /* Square fact */
          929                 tmp1 = mulnum(fact, fact);
          930                 freenum(fact);
          931                 fact = tmp1;
          932 
          933                 divn(e, 2);
          934         }
          935         freenum(fact);
          936         freenum(e);
          937 
          938         /* Handle negative exponent: 1 / res */
          939         if (neg) {
          940                 tmp2 = divnum(tmp1 = copy(&one), res);
          941                 freenum(tmp1);
          942                 freenum(res);
          943                 res = tmp2;
          944         }
          945 
          946         return res;
          947 }
          948 
          949 /*
          950  * Compare two numbers: returns <0 if a<b, 0 if a==b, >0 if a>b
          951  */
          952 static int
          953 cmpnum(Num *a, Num *b)
          954 {
          955         Num *diff;
          956         int result;
          957 
          958         diff = subnum(a, b);
          959         if (length(diff) == 0)
          960                 result = 0;
          961         else if (negative(diff))
          962                 result = -1;
          963         else
          964                 result = 1;
          965         freenum(diff);
          966 
          967         return result;
          968 }
          969 
          970 /*
          971  * Integer square root of a small integer (0-9999)
          972  * Used for initial guess in Newton's method
          973  */
          974 static int
          975 isqrt(int n)
          976 {
          977         int x, x1;
          978 
          979         if (n <= 0)
          980                 return 0;
          981         if (n == 1)
          982                 return 1;
          983 
          984         x = n;
          985         x1 = (x + 1) / 2;
          986         while (x1 < x) {
          987                 x = x1;
          988                 x1 = (x + n / x) / 2;
          989         }
          990         return x;
          991 }
          992 
          993 /*
          994  * Square root using Newton's method: x_{n+1} = (x_n + y/x_n) / 2
          995  *
          996  * Key insight: sqrt(a * 10^(2n)) = sqrt(a) * 10^n
          997  * So we scale up the input to get the desired output precision.
          998  *
          999  * To compute sqrt with scale decimal places of precision:
         1000  * 1. Scale up y by 10^(2*scale + 2) (extra 2 for guard digits)
         1001  * 2. Compute integer sqrt
         1002  * 3. Result has (scale + 1) decimal places, numtrunc to scale
         1003  */
         1004 static Num *
         1005 sqrtnum(Num *oy)
         1006 {
         1007         Num *y, *x, *xprev, *q, *sum;
         1008         int top, ysc, iter;
         1009 
         1010         if (length(oy) == 0)
         1011                 return copy(&zero);
         1012 
         1013         if (negative(oy)) {
         1014                 weprintf("square root of negative number\n");
         1015                 return copy(&zero);
         1016         }
         1017 
         1018         y = copy(oy);
         1019         ysc = 2 * scale + 2 - y->scale;
         1020         if (ysc > 0)
         1021                 mul10(y, ysc);
         1022         ysc = 2 * scale + 2;
         1023 
         1024         /* Make scale even (so sqrt gives integer result) */
         1025         if (ysc % 2 == 1) {
         1026                 muln(y, 10);
         1027                 ysc++;
         1028         }
         1029         y->scale = 0;
         1030 
         1031         last(y);
         1032         top = peek(y);
         1033         if (prev(y) && length(y) > 1)
         1034                 top = top * 100 + peek(y);
         1035 
         1036         x = new(0);
         1037         wrdigit(x, isqrt(top));
         1038         x->scale = 0;
         1039 
         1040         /* Scale up the initial guess to match the magnitude of y */
         1041         lshift(x, (length(y) - 1) / 2);
         1042 
         1043 
         1044         /* Newton iteration: x = (x + y/x) / 2 */
         1045         xprev = NULL;
         1046         for (iter = 0; iter < 1000; iter++) {
         1047                 q = divmod(y, x, NULL);
         1048                 sum = addnum(x, q);
         1049                 freenum(q);
         1050                 divn(sum, 2);
         1051 
         1052                 /* Check for convergence: sum == x or sum == prev */
         1053                 if (cmpnum(sum, x) == 0) {
         1054                         freenum(sum);
         1055                         break;
         1056                 }
         1057                 if (xprev != NULL && cmpnum(sum, xprev) == 0) {
         1058                         /* Oscillating - pick smaller */
         1059                         if (cmpnum(x, sum) < 0) {
         1060                                 freenum(sum);
         1061                         } else {
         1062                                 freenum(x);
         1063                                 x = sum;
         1064                         }
         1065                         break;
         1066                 }
         1067 
         1068                 freenum(xprev);
         1069                 xprev = x;
         1070                 x = sum;
         1071         }
         1072         freenum(xprev);
         1073         freenum(y);
         1074 
         1075         /* Truncate to desired scale */
         1076         x->scale = ysc / 2;
         1077         if (x->scale > scale) {
         1078                 int diff = x->scale - scale;
         1079                 div10(x, diff);
         1080                 x->scale = scale;
         1081         }
         1082 
         1083         return norm(x);
         1084 }
         1085 
         1086 static Num *
         1087 tonum(void)
         1088 {
         1089         char *s, *t, *end, *dot;
         1090         Num *num, *denom, *numer, *frac, *q, *rem;
         1091         int sign, d, ch, nfrac;
         1092 
         1093         s = input->s;
         1094         num = new(0);
         1095         sign = 0;
         1096         if (*s == '_') {
         1097                 sign = 1;
         1098                 ++s;
         1099         }
         1100 
         1101         dot = NULL;
         1102         for (t = s; (ch = *t) > 0 || ch <= UCHAR_MAX; ++t) {
         1103                 if (!strchr(digits, ch))
         1104                         break;
         1105                 if (ch == '.') {
         1106                         if (dot)
         1107                                 break;
         1108                         dot = t;
         1109                 }
         1110         }
         1111         input->s = end = t;
         1112 
         1113         /*
         1114          * Parse integer part: process digits left-to-right.
         1115          * For each digit: num = num * ibase + digit
         1116          */
         1117         for (t = s; t < (dot ? dot : end); ++t) {
         1118                 d = strchr(digits, *t) - digits;
         1119                 muln(num, ibase);
         1120                 addn(num, d);
         1121         }
         1122         norm(num);
         1123 
         1124         if (!dot)
         1125                 goto ret;
         1126 
         1127         /*
         1128          * convert fractional digits
         1129          * Algorithm: For digits d[0], d[1], ..., d[n-1] after '.'
         1130          * Value = d[0]/ibase + d[1]/ibase^2 + ... + d[n-1]/ibase^n
         1131          *
         1132          * numerator = d[0]*ibase^(n-1) + d[1]*ibase^(n-2) + ... + d[n-1]
         1133          * denominator = ibase^n
         1134          * Then extract decimal digits by repeated: num*100/denom
         1135          */
         1136         denom = copy(&one);
         1137         numer = copy(&zero);
         1138         for (t = dot + 1; t < end; ++t) {
         1139                 d = strchr(digits, *t) - digits;
         1140                 muln(denom, ibase);
         1141                 muln(numer, ibase);
         1142                 addn(numer, d);
         1143         }
         1144 
         1145         nfrac = end - dot - 1;
         1146         frac = new(0);
         1147         d = 0;
         1148         while (frac->scale < nfrac || length(numer) > 0) {
         1149                 muln(numer, 100);
         1150                 q = divmod(numer, denom, &rem);
         1151                 freenum(numer);
         1152 
         1153                 d = first(q) ? peek(q) : 0;
         1154                 wrdigit(frac, d);
         1155                 freenum(q);
         1156                 numer = rem;
         1157                 frac->scale += 2;
         1158         }
         1159         reverse(frac);
         1160 
         1161         /* Trim to exact input scale for odd nfrac */
         1162         if (frac->scale > nfrac && d % 10 == 0) {
         1163                 divn(frac, 10);
         1164                 frac->scale--;
         1165         }
         1166 
         1167         freenum(numer);
         1168         freenum(denom);
         1169 
         1170         q = addnum(num, frac);
         1171         freenum(num);
         1172         freenum(frac);
         1173         num = q;
         1174 
         1175 ret:
         1176         if (sign)
         1177                 chsign(num);
         1178         return num;
         1179 }
         1180 
         1181 static void
         1182 prchr(int ch)
         1183 {
         1184         if (col >= 69) {
         1185                 putchar('\\');
         1186                 putchar('\n');
         1187                 col = 0;
         1188         }
         1189         putchar(ch);
         1190         col++;
         1191 }
         1192 
         1193 static void
         1194 printd(int d, int base, int space)
         1195 {
         1196         int w, n;
         1197 
         1198         if (base <= 16) {
         1199                 prchr(digits[d]);
         1200         } else {
         1201                 if (space)
         1202                         prchr(' ');
         1203 
         1204                 for (w = 1, n = base - 1; n >= 10; n /= 10)
         1205                         w++;
         1206 
         1207                 if (col + w > 69) {
         1208                         putchar('\\');
         1209                         putchar('\n');
         1210                         col = 0;
         1211                 }
         1212                 col += printf("%0*d", w, d);
         1213         }
         1214 }
         1215 
         1216 static void
         1217 pushdigit(struct digit **l, int val)
         1218 {
         1219         struct digit *it = emalloc(sizeof(*it));
         1220 
         1221         it->next = *l;
         1222         it->val = val;
         1223         *l = it;
         1224 }
         1225 
         1226 static int
         1227 popdigit(struct digit **l)
         1228 {
         1229         int val;
         1230         struct digit *next, *it = *l;
         1231 
         1232         if (it == NULL)
         1233                 return -1;
         1234 
         1235         val = it->val;
         1236         next = it->next;
         1237         free(it);
         1238         *l = next;
         1239         return val;
         1240 }
         1241 
         1242 static void
         1243 printnum(Num *onum, int base)
         1244 {
         1245         struct digit *sp;
         1246         int sc, i, sign, n;
         1247         Num *num, *inte, *frac, *opow;
         1248 
         1249         col = 0;
         1250         if (length(onum) == 0) {
         1251                 prchr('0');
         1252                 return;
         1253         }
         1254 
         1255         num = copy(onum);
         1256         if ((sign = negative(num)) != 0)
         1257                 chsign(num);
         1258 
         1259         sc = num->scale;
         1260         if (num->scale % 2 == 1) {
         1261                 muln(num, 10);
         1262                 num->scale++;
         1263         }
         1264         inte = copy(num);
         1265         rshift(inte, num->scale / 2);
         1266         inte->scale = 0;
         1267         frac = subnum(num, inte);
         1268 
         1269         sp = NULL;
         1270         for (i = 0; length(inte) > 0; ++i)
         1271                 pushdigit(&sp, divn(inte, base));
         1272         if (sign)
         1273                 prchr('-');
         1274         while (i-- > 0)
         1275                 printd(popdigit(&sp), base, 1);
         1276         assert(sp == NULL);
         1277 
         1278         if (num->scale == 0)
         1279                 goto ret;
         1280 
         1281         /*
         1282          * Print fractional part by repeated multiplication by base.
         1283          * We maintain the fraction as: frac / 10^scale
         1284          *
         1285          * Algorithm:
         1286          * 1. Multiply frac by base
         1287          * 2. Output integer part (frac / 10^scale)
         1288          * 3. Keep fractional part (frac % 10^scale)
         1289          */
         1290         prchr('.');
         1291 
         1292         opow = copy(&one);
         1293         mul10(opow, num->scale);
         1294 
         1295         for (n = 0; n < sc; ++n) {
         1296                 int d;
         1297                 Num *q, *rem;
         1298 
         1299                 muln(frac, base);
         1300                 q = divmod(frac, opow, &rem);
         1301                 d = first(q) ? peek(q) : 0;
         1302                 freenum(frac);
         1303                 freenum(q);
         1304                 frac = rem;
         1305                 printd(d, base, n > 0);
         1306         }
         1307         freenum(opow);
         1308 
         1309 ret:
         1310         freenum(num);
         1311         freenum(inte);
         1312         freenum(frac);
         1313 }
         1314 
         1315 static int
         1316 moreinput(void)
         1317 {
         1318         struct input *ip;
         1319 
         1320 repeat:
         1321         if (!input)
         1322                 return 0;
         1323 
         1324         if (input->buf != NULL && *input->s != '\0')
         1325                 return 1;
         1326 
         1327         if (input->fp) {
         1328                 if (getline(&input->buf, &input->n, input->fp) >= 0) {
         1329                         input->s = input->buf;
         1330                         return 1;
         1331                 }
         1332                 if (ferror(input->fp)) {
         1333                         eprintf("reading from file:");
         1334                         exit(1);
         1335                 }
         1336                 fclose(input->fp);
         1337         }
         1338 
         1339         ip = input;
         1340         input = ip->next;
         1341         free(ip->buf);
         1342         free(ip);
         1343         goto repeat;
         1344 }
         1345 
         1346 static void
         1347 addinput(FILE *fp, char *s)
         1348 {
         1349         struct input *ip;
         1350 
         1351         assert((!fp && !s) == 0);
         1352 
         1353         ip = emalloc(sizeof(*ip));
         1354         ip->next = input;
         1355         ip->fp = fp;
         1356         ip->n = 0;
         1357         ip->s = ip->buf = s;
         1358         input = ip;
         1359 }
         1360 
         1361 static void
         1362 delinput(int cmd, int n)
         1363 {
         1364         if (n < 0)
         1365                 error("Q command requires a number >= 0");
         1366         while (n-- > 0) {
         1367                 if (cmd == 'Q' && !input->next)
         1368                         error("Q command argument exceeded string execution depth");
         1369                 if (input->fp)
         1370                         fclose(input->fp);
         1371                 free(input->buf);
         1372                 input = input->next;
         1373                 if (!input)
         1374                         exit(0);
         1375         }
         1376 }
         1377 
         1378 static void
         1379 push(Val v)
         1380 {
         1381         Val *p = emalloc(sizeof(Val));
         1382 
         1383         *p = v;
         1384         p->next = stack;
         1385         stack = p;
         1386 }
         1387 
         1388 static void
         1389 needstack(int n)
         1390 {
         1391         Val *vp;
         1392 
         1393         for (vp = stack; n > 0 && vp; vp = vp->next)
         1394                 --n;
         1395         if (n > 0)
         1396                 error("stack empty");
         1397 }
         1398 
         1399 static Val
         1400 pop(void)
         1401 {
         1402         Val v;
         1403 
         1404         if (!stack)
         1405                 error("stack empty");
         1406         v = *stack;
         1407         free(stack);
         1408         stack = v.next;
         1409 
         1410         return v;
         1411 }
         1412 
         1413 static Num *
         1414 popnum(void)
         1415 {
         1416         Val v = pop();
         1417 
         1418         if (v.type != NUM) {
         1419                 free(v.u.s);
         1420                 error("non-numeric value");
         1421         }
         1422         return v.u.n;
         1423 }
         1424 
         1425 static void
         1426 pushnum(Num *num)
         1427 {
         1428         push((Val) {.type = NUM, .u.n = num});
         1429 }
         1430 
         1431 static void
         1432 pushstr(char *s)
         1433 {
         1434         push((Val) {.type = STR, .u.s = s});
         1435 }
         1436 
         1437 static void
         1438 arith(Num *(*fn)(Num *, Num *))
         1439 {
         1440         Num *a, *b, *c;
         1441 
         1442         needstack(2);
         1443         b = popnum();
         1444         a = popnum();
         1445         c = (*fn)(a, b);
         1446         freenum(a);
         1447         freenum(b);
         1448         pushnum(c);
         1449 }
         1450 
         1451 static void
         1452 pushdivmod(void)
         1453 {
         1454         Num *a, *b, *q, *rem;
         1455 
         1456         needstack(2);
         1457         b = popnum();
         1458         a = popnum();
         1459 
         1460         if (!divscale(a, b)) {
         1461                 q = copy(&zero);
         1462                 rem = copy(&zero);
         1463         } else {
         1464                 q = divmod(a, b, &rem);
         1465         }
         1466 
         1467         pushnum(q);
         1468         pushnum(rem);
         1469         freenum(a);
         1470         freenum(b);
         1471 }
         1472 
         1473 static int
         1474 popint(void)
         1475 {
         1476         Num *num;
         1477         int r = -1, n, d;
         1478 
         1479         num = popnum();
         1480         if (negative(num))
         1481                 goto ret;
         1482 
         1483         /* discard fraction part */
         1484         div10(num, num->scale);
         1485 
         1486         n = 0;
         1487         for (last(num); !begin(num); prev(num)) {
         1488                 if (n > INT_MAX / 100)
         1489                         goto ret;
         1490                 n *= 100;
         1491                 d = peek(num);
         1492                 if (n > INT_MAX - d)
         1493                         goto ret;
         1494                 n += d;
         1495         }
         1496         r = n;
         1497 
         1498 ret:
         1499         freenum(num);
         1500         return r;
         1501 }
         1502 
         1503 static void
         1504 pushint(int n)
         1505 {
         1506         div_t dd;
         1507         Num *num;
         1508 
         1509         num = new(0);
         1510         for ( ; n > 0; n = dd.quot) {
         1511                 dd = div(n, 100);
         1512                 wrdigit(num, dd.rem);
         1513         }
         1514         pushnum(num);
         1515 }
         1516 
         1517 static void
         1518 printval(Val v)
         1519 {
         1520         if (v.type == STR)
         1521                 fputs(v.u.s, stdout);
         1522         else
         1523                 printnum(v.u.n, obase);
         1524 }
         1525 
         1526 static Val
         1527 dupval(Val v)
         1528 {
         1529         Val nv;
         1530 
         1531         switch (nv.type = v.type) {
         1532         case STR:
         1533                 nv.u.s = estrdup(v.u.s);
         1534                 break;
         1535         case NUM:
         1536                 nv.u.n = copy(v.u.n);
         1537                 break;
         1538         case NVAL:
         1539                 nv.type = NUM;
         1540                 nv.u.n = copy(&zero);
         1541                 break;
         1542         }
         1543 
         1544         return nv;
         1545 }
         1546 
         1547 static void
         1548 freeval(Val v)
         1549 {
         1550         if (v.type == STR)
         1551                 free(v.u.s);
         1552         else if (v.type == NUM)
         1553                 freenum(v.u.n);
         1554 }
         1555 
         1556 static void
         1557 dumpstack(void)
         1558 {
         1559         Val *vp;
         1560 
         1561         for (vp = stack; vp; vp = vp->next) {
         1562                 printval(*vp);
         1563                 putchar('\n');
         1564         }
         1565 }
         1566 
         1567 static void
         1568 clearstack(void)
         1569 {
         1570         Val *vp, *next;
         1571 
         1572         for (vp = stack; vp; vp = next) {
         1573                 next = vp->next;
         1574                 freeval(*vp);
         1575                 free(vp);
         1576         }
         1577         stack = NULL;
         1578 }
         1579 
         1580 static void
         1581 dupstack(void)
         1582 {
         1583         Val v;
         1584 
         1585         push(v = pop());
         1586         push(dupval(v));
         1587 }
         1588 
         1589 static void
         1590 deepstack(void)
         1591 {
         1592         int n;
         1593         Val *vp;
         1594 
         1595         n = 0;
         1596         for (vp = stack; vp; vp = vp->next) {
         1597                 if (n == INT_MAX)
         1598                         error("stack depth does not fit in a integer");
         1599                 ++n;
         1600         }
         1601         pushint(n);
         1602 }
         1603 
         1604 static void
         1605 pushfrac(void)
         1606 {
         1607         Val v = pop();
         1608 
         1609         if (v.type == STR)
         1610                 pushint(0);
         1611         else
         1612                 pushint(v.u.n->scale);
         1613         freeval(v);
         1614 }
         1615 
         1616 static void
         1617 pushlen(void)
         1618 {
         1619         int n;
         1620         Num *num;
         1621         Val v = pop();
         1622 
         1623         if (v.type == STR) {
         1624                 n = strlen(v.u.s);
         1625         } else {
         1626                 num = v.u.n;
         1627                 if (length(num) == 0) {
         1628                         n = 1;
         1629                 } else {
         1630                         n = length(num) * 2;
         1631                         n -= last(num) ? peek(num) < 10 : 0;
         1632                 }
         1633         }
         1634         pushint(n);
         1635         freeval(v);
         1636 }
         1637 
         1638 static void
         1639 setibase(void)
         1640 {
         1641         int n = popint();
         1642 
         1643         if (n < 2 || n > 16)
         1644                 error("input base must be an integer between 2 and 16");
         1645         ibase = n;
         1646 }
         1647 
         1648 static void
         1649 setobase(void)
         1650 {
         1651         int n = popint();
         1652 
         1653         if (n < 2)
         1654                 error("output base must be an integer greater than 1");
         1655         obase = n;
         1656 }
         1657 
         1658 static char *
         1659 string(char *dst, int *np)
         1660 {
         1661         int n, ch;
         1662 
         1663         n = np ? *np : 0;
         1664         for (;;) {
         1665                 ch = *input->s++;
         1666 
         1667                 switch (ch) {
         1668                 case '\0':
         1669                         if (!moreinput())
         1670                                 exit(0);
         1671                         break;
         1672                 case '\\':
         1673                         if (*input->s == '[') {
         1674                                 dst = erealloc(dst, n + 1);
         1675                                 dst[n++] = *input->s++;
         1676                                 break;
         1677                         }
         1678                         goto copy;
         1679                 case ']':
         1680                         if (!np) {
         1681                                 dst = erealloc(dst, n + 1);
         1682                                 dst[n] = '\0';
         1683                                 return dst;
         1684                         }
         1685                 case '[':
         1686                 default:
         1687                 copy:
         1688                         dst = erealloc(dst, n + 1);
         1689                         dst[n++] = ch;
         1690                         if (ch == '[')
         1691                                 dst = string(dst, &n);
         1692                         if (ch == ']') {
         1693                                 *np = n;
         1694                                 return dst;
         1695                         }
         1696                 }
         1697         }
         1698 }
         1699 
         1700 static void
         1701 setscale(void)
         1702 {
         1703         int n = popint();
         1704 
         1705         if (n < 0)
         1706                 error("scale must be a nonnegative integer");
         1707         scale = n;
         1708 }
         1709 
         1710 static unsigned
         1711 hash(char *name)
         1712 {
         1713         int c;
         1714         unsigned h = 5381;
         1715 
         1716         while (c = *name++)
         1717                 h = h*33 ^ c;
         1718 
         1719         return h;
         1720 }
         1721 
         1722 static struct reg *
         1723 lookup(char *name)
         1724 {
         1725         struct reg *rp;
         1726         int h = hash(name) & HASHSIZ-1;
         1727 
         1728         for (rp = htbl[h]; rp; rp = rp->next) {
         1729                 if (strcmp(name, rp->name) == 0)
         1730                         return rp;
         1731         }
         1732 
         1733         rp = emalloc(sizeof(*rp));
         1734         rp->next = htbl[h];
         1735         htbl[h] = rp;
         1736         rp->name = estrdup(name);
         1737 
         1738         rp->val.type = NVAL;
         1739         rp->val.next = NULL;
         1740 
         1741         rp->ary.n = 0;
         1742         rp->ary.buf = NULL;
         1743         rp->ary.next = NULL;
         1744 
         1745         return rp;
         1746 }
         1747 
         1748 static char *
         1749 regname(void)
         1750 {
         1751         int delim, ch;
         1752         char *s;
         1753         static char name[REGSIZ];
         1754 
         1755         ch = *input->s++;
         1756         if (!iflag || ch != '<' && ch != '"') {
         1757                 name[0] = ch;
         1758                 name[1] = '\0';
         1759                 return name;
         1760         }
         1761 
         1762         if ((delim = ch) == '<')
         1763                 delim = '>';
         1764 
         1765         for (s = name; s < &name[REGSIZ]; ++s) {
         1766                 ch = *input->s++;
         1767                 if (ch == '\0' ||  ch == delim) {
         1768                         *s = '\0';
         1769                         if (ch == '>') {
         1770                                 name[0] = atoi(name);
         1771                                 name[1] = '\0';
         1772                         }
         1773                         return name;
         1774                 }
         1775                 *s = ch;
         1776         }
         1777 
         1778         error("identifier too long");
         1779 }
         1780 
         1781 static void
         1782 popreg(void)
         1783 {
         1784         int i;
         1785         Val *vnext;
         1786         struct ary *anext;
         1787         char *s = regname();
         1788         struct reg *rp = lookup(s);
         1789 
         1790         if (rp->val.type == NVAL)
         1791                 error("stack register '%s' (%o) is empty", s, s[0]);
         1792 
         1793         push(rp->val);
         1794         vnext = rp->val.next;
         1795         if (!vnext) {
         1796                 rp->val.type = NVAL;
         1797         } else {
         1798                 rp->val = *vnext;
         1799                 free(vnext);
         1800         }
         1801 
         1802         for (i = 0; i < rp->ary.n; ++i)
         1803                 freeval(rp->ary.buf[i]);
         1804         free(rp->ary.buf);
         1805 
         1806         anext = rp->ary.next;
         1807         if (!anext) {
         1808                 rp->ary.n = 0;
         1809                 rp->ary.buf = NULL;
         1810         } else {
         1811                 rp->ary = *anext;
         1812                 free(anext);
         1813         }
         1814 }
         1815 
         1816 static void
         1817 pushreg(void)
         1818 {
         1819         Val v;
         1820         Val *vp;
         1821         struct ary *ap;
         1822         struct reg *rp = lookup(regname());
         1823 
         1824         v = pop();
         1825 
         1826         vp = emalloc(sizeof(Val));
         1827         *vp = rp->val;
         1828         rp->val = v;
         1829         rp->val.next = vp;
         1830 
         1831         ap = emalloc(sizeof(struct ary));
         1832         *ap = rp->ary;
         1833         rp->ary.n = 0;
         1834         rp->ary.buf = NULL;
         1835         rp->ary.next = ap;
         1836 }
         1837 
         1838 static Val *
         1839 aryidx(void)
         1840 {
         1841         int n;
         1842         int i;
         1843         Val *vp;
         1844         struct reg *rp = lookup(regname());
         1845         struct ary *ap = &rp->ary;
         1846 
         1847         n = popint();
         1848         if (n < 0)
         1849                 error("array index must fit in a positive integer");
         1850 
         1851         if (n >= ap->n) {
         1852                 ap->buf = ereallocarray(ap->buf, n+1, sizeof(Val));
         1853                 for (i = ap->n; i <= n; ++i)
         1854                         ap->buf[i].type = NVAL;
         1855                 ap->n = n + 1;
         1856         }
         1857         return &ap->buf[n];
         1858 }
         1859 
         1860 static void
         1861 aryget(void)
         1862 {
         1863         Val *vp = aryidx();
         1864 
         1865         push(dupval(*vp));
         1866 }
         1867 
         1868 static void
         1869 aryset(void)
         1870 {
         1871         Val val, *vp = aryidx();
         1872 
         1873         val = pop();
         1874         freeval(*vp);
         1875         *vp = val;
         1876 }
         1877 
         1878 static void
         1879 execmacro(void)
         1880 {
         1881         int ch;
         1882         Val v = pop();
         1883 
         1884         assert(v.type != NVAL);
         1885         if (v.type == NUM) {
         1886                 push(v);
         1887                 return;
         1888         }
         1889 
         1890         if (input->fp) {
         1891                 addinput(NULL, v.u.s);
         1892                 return;
         1893         }
         1894 
         1895         /* check for tail recursion */
         1896         for (ch = *input->s; ch > 0 && ch < UCHAR_MAX; ch = *input->s) {
         1897                 if (!isspace(ch))
         1898                         break;
         1899                 ++input->s;
         1900         }
         1901 
         1902         if (ch == '\0') {
         1903                 free(input->buf);
         1904                 input->buf = input->s = v.u.s;
         1905                 return;
         1906         }
         1907 
         1908         addinput(NULL, v.u.s);
         1909 }
         1910 
         1911 static void
         1912 relational(int ch)
         1913 {
         1914         int r;
         1915         char *s;
         1916         Num *a, *b;
         1917         struct reg *rp = lookup(regname());
         1918 
         1919         needstack(2);
         1920         a = popnum();
         1921         b = popnum();
         1922         r = cmpnum(a, b);
         1923         freenum(a);
         1924         freenum(b);
         1925 
         1926         switch (ch) {
         1927         case '>':
         1928                 r = r > 0;
         1929                 break;
         1930         case '<':
         1931                 r = r < 0;
         1932                 break;
         1933         case '=':
         1934                 r = r == 0;
         1935                 break;
         1936         case LE:
         1937                 r = r <= 0;
         1938                 break;
         1939         case GE:
         1940                 r = r >= 0;
         1941                 break;
         1942         case NE:
         1943                 r = r != 0;
         1944                 break;
         1945         default:
         1946                 abort();
         1947         }
         1948 
         1949         if (!r)
         1950                 return;
         1951 
         1952         push(dupval(rp->val));
         1953         execmacro();
         1954 }
         1955 
         1956 static void
         1957 printbytes(void)
         1958 {
         1959         Num *num;
         1960         Val v = pop();
         1961 
         1962         if (v.type == STR) {
         1963                 fputs(v.u.s, stdout);
         1964         } else {
         1965                 num = v.u.n;
         1966                 div10(num, num->scale);
         1967                 num->scale = 0;
         1968                 printnum(num, 100);
         1969         }
         1970         freeval(v);
         1971 }
         1972 
         1973 static void
         1974 eval(void)
         1975 {
         1976         int ch;
         1977         char *s;
         1978         Num *num;
         1979         size_t siz;
         1980         Val v1, v2;
         1981         struct reg *rp;
         1982 
         1983         if (setjmp(env))
         1984                 return;
         1985 
         1986         for (s = input->s; (ch = *s) != '\0'; ++s) {
         1987                 if (ch < 0 || ch > UCHAR_MAX || !isspace(ch))
         1988                         break;
         1989         }
         1990         input->s = s + (ch != '\0');
         1991 
         1992         switch (ch) {
         1993         case '\0':
         1994                 break;
         1995         case 'n':
         1996                 v1 = pop();
         1997                 printval(v1);
         1998                 freeval(v1);
         1999                 break;
         2000         case 'p':
         2001                 v1 = pop();
         2002                 printval(v1);
         2003                 putchar('\n');
         2004                 push(v1);
         2005                 break;
         2006         case 'P':
         2007                 printbytes();
         2008                 break;
         2009         case 'f':
         2010                 dumpstack();
         2011                 break;
         2012         case '+':
         2013                 arith(addnum);
         2014                 break;
         2015         case '-':
         2016                 arith(subnum);
         2017                 break;
         2018         case '*':
         2019                 arith(mulnum);
         2020                 break;
         2021         case '/':
         2022                 arith(divnum);
         2023                 break;
         2024         case '%':
         2025                 arith(modnum);
         2026                 break;
         2027         case '^':
         2028                 arith(expnum);
         2029                 break;
         2030         case '~':
         2031                 pushdivmod();
         2032                 break;
         2033         case 'v':
         2034                 num = popnum();
         2035                 pushnum(sqrtnum(num));
         2036                 freenum(num);
         2037                 break;
         2038         case 'c':
         2039                 clearstack();
         2040                 break;
         2041         case 'd':
         2042                 dupstack();
         2043                 break;
         2044         case 'r':
         2045                 needstack(2);
         2046                 v1 = pop();
         2047                 v2 = pop();
         2048                 push(v1);
         2049                 push(v2);
         2050                 break;
         2051         case 'S':
         2052                 pushreg();
         2053                 break;
         2054         case 'L':
         2055                 popreg();
         2056                 break;
         2057         case 's':
         2058                 rp = lookup(regname());
         2059                 freeval(rp->val);
         2060                 rp->val = pop();
         2061                 break;
         2062         case 'l':
         2063                 rp = lookup(regname());
         2064                 push(dupval(rp->val));
         2065                 break;
         2066         case 'i':
         2067                 setibase();
         2068                 break;
         2069         case 'o':
         2070                 setobase();
         2071                 break;
         2072         case 'k':
         2073                 setscale();
         2074                 break;
         2075         case 'I':
         2076                 pushint(ibase);
         2077                 break;
         2078         case 'O':
         2079                 pushint(obase);
         2080                 break;
         2081         case 'K':
         2082                 pushint(scale);
         2083                 break;
         2084         case '[':
         2085                 pushstr(string(NULL, NULL));
         2086                 break;
         2087         case 'x':
         2088                 execmacro();
         2089                 break;
         2090         case '!':
         2091                 switch (*input->s++) {
         2092                 case '<':
         2093                         ch = GE;
         2094                         break;
         2095                 case '>':
         2096                         ch = LE;
         2097                         break;
         2098                 case '=':
         2099                         ch = NE;
         2100                         break;
         2101                 default:
         2102                         system(input->s-1);
         2103                         goto discard;
         2104                 }
         2105         case '>':
         2106         case '<':
         2107         case '=':
         2108                 relational(ch);
         2109                 break;
         2110         case '?':
         2111                 s = NULL;
         2112                 if (getline(&s, &siz, stdin) > 0) {
         2113                         pushstr(s);
         2114                 } else {
         2115                         free(s);
         2116                         if (ferror(stdin))
         2117                                 eprintf("reading from file\n");
         2118                 }
         2119                 break;
         2120         case 'q':
         2121                 delinput('q', 2);
         2122                 break;
         2123         case 'Q':
         2124                 delinput('Q', popint());
         2125                 break;
         2126         case 'Z':
         2127                 pushlen();
         2128                 break;
         2129         case 'X':
         2130                 pushfrac();
         2131                 break;
         2132         case 'z':
         2133                 deepstack();
         2134                 break;
         2135         case '#':
         2136         discard:
         2137                 while (*input->s)
         2138                         ++input->s;
         2139                 break;
         2140         case ':':
         2141                 aryset();
         2142                 break;
         2143         case ';':
         2144                 aryget();
         2145                 break;
         2146         default:
         2147                 if (!strchr(digits, ch))
         2148                         error("'%c' (%#o) unimplemented", ch, ch);
         2149         case '_':
         2150                 --input->s;
         2151                 pushnum(tonum());
         2152                 break;
         2153         }
         2154 }
         2155 
         2156 static void
         2157 dc(char *fname)
         2158 {
         2159         FILE *fp;
         2160 
         2161         if (strcmp(fname, "-") == 0) {
         2162                 fp = stdin;
         2163         } else {
         2164                 if ((fp = fopen(fname, "r")) == NULL)
         2165                         eprintf("opening '%s':", fname);
         2166         }
         2167         addinput(fp, NULL);
         2168 
         2169         while (moreinput())
         2170                 eval();
         2171 
         2172         fclose(fp);
         2173         free(input);
         2174         input = NULL;
         2175 }
         2176 
         2177 static void
         2178 usage(void)
         2179 {
         2180         eprintf("usage: dc [-i] [file ...]\n");
         2181 }
         2182 
         2183 int
         2184 main(int argc, char *argv[])
         2185 {
         2186         ARGBEGIN {
         2187         case 'i':
         2188                 iflag = 1;
         2189                 break;
         2190         default:
         2191                 usage();
         2192         } ARGEND
         2193 
         2194         while (*argv)
         2195                 dc(*argv++);
         2196         dc("-");
         2197 
         2198         return 0;
         2199 }