sbase/dc.c

#include <assert.h>
#include <errno.h>
#include <ctype.h>
#include <setjmp.h>
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>

#include "arg.h"
#include "util.h"

#define NDIGITS 10
#define REGSIZ 16
#define HASHSIZ 128

enum {
	NVAL,
	STR,
	NUM,
};

enum {
	NE,
	LE,
	GE,
};

typedef struct num Num;
typedef struct val Val;

struct num {
	int begin;
	int scale;
	int room;
	signed char *buf, *wp, *rp;
};

struct digit {
	int val;
	struct digit *next;
};

struct val {
	int type;
	union {
		Num *n;
		char *s;
	} u;

	struct val *next;
};

struct ary {
	int n;
	Val *buf;
	struct ary *next;
};

struct reg {
	char *name;
	Val val;
	struct ary ary;
	struct reg *next;
};

struct input {
	FILE *fp;
	char *buf;
	size_t n;
	char *s;
	struct input *next;
};

static Val *stack;
static jmp_buf env;
static struct input *input;
static struct reg *htbl[HASHSIZ];

static signed char onestr[] = {1, 0};
static Num zero, one = {.buf = onestr, .wp = onestr + 1};
static char digits[] = "0123456789ABCDEF.";

static int scale, ibase = 10, obase = 10;
static int iflag;
static int col;

/*
 * This dc implementation follows the decription of dc found in the paper
 * DC - An Interactive Desk Calculator, by Robert Morris and Lorinda Cherry
 */

static void
dumpnum(char *msg, Num *num)
{
	signed char *p;

	printf("Num (%s)", msg ? msg : "");

	if (!num) {
		puts("NULL\n");
		return;
	}
	printf("\n"
	       "\tscale=%d\n"
	       "\troom=%d\n"
	       "\tdigits=[\n",
	       num->scale,
	       num->room);
	for (p = num->buf; p < num->wp; ++p)
		printf("\t\t%02d\n", *p);
	printf("\t]\n");
}

/*
 * error is not called from the implementation of the
 * arithmetic functions because that can drive to memory
 * leaks very easily.
 */
static void
error(char *fmt, ...)
{
	va_list va;

	va_start(va, fmt);
	xvprintf(fmt, va);
	putc('\n', stderr);
	va_end(va);

	longjmp(env, 1);
}

static void
freenum(Num *num)
{
	if (!num)
		return;
	free(num->buf);
	free(num);
}

static Num *
moreroom(Num *num, int more)
{
	int ro, wo, room;
	signed char *p;

	ro = num->rp - num->buf;
	wo = num->wp - num->buf;
	room = num->room;

	if (room > INT_MAX - more)
		eprintf("out of memory\n");

	room = room + more;
	if (room < NDIGITS)
		room = NDIGITS;
	p = erealloc(num->buf, room);
	num->buf = p;
	num->rp = p + ro;
	num->wp = p + wo;
	num->room = room;

	return num;
}

static Num *
grow(Num *num)
{
	return moreroom(num, NDIGITS);
}

static Num *
expand(Num *num, int min)
{
	if (min < num->room)
		return num;
	return moreroom(num, min - num->room);
}

static Num *
new(int room)
{
	Num *num = emalloc(sizeof(*num));

	num->rp = num->wp = num->buf = NULL;
	num->begin = num->room = num->scale = 0;

	return moreroom(num, room);
}

static Num *
zeronum(int ndigits)
{
	Num *num = new(ndigits);

	num->wp = num->buf + ndigits;
	memset(num->buf, 0, ndigits);

	return num;
}

static Num *
wrdigit(Num *num, int d)
{
	if (num->wp == &num->buf[num->room])
		grow(num);
	*num->wp++ = d;

	return num;
}

static int
rddigit(Num *num)
{
	if (num->rp == num->wp)
		return -1;
	return *num->rp++;
}

static int
peek(Num *num)
{
	if (num->rp == num->wp)
		return -1;
	return *num->rp;
}

static Num *
poke(Num *num, unsigned d)
{
	if (num->rp == &num->buf[num->room])
		grow(num);
	if (num->rp == num->wp)
		num->wp++;
	*num->rp = d;

	return num;
}

static int
begin(Num *num)
{
	return num->begin != 0;
}

static int
first(Num *num)
{
	num->begin = 0;
	num->rp = num->buf;
	return num->rp != num->wp;
}

static int
last(Num *num)
{
	if (num->wp != num->buf) {
		num->begin = 0;
		num->rp = num->wp - 1;
		return 1;
	}
	num->begin = 1;
	return 0;
}

static int
prev(Num *num)
{
	if (num->rp > num->buf) {
		--num->rp;
		return 1;
	}
	num->begin = 1;
	return 0;
}

static int
next(Num *num)
{
	if (num->rp != num->wp + 1) {
		++num->rp;
		return 1;
	}
	return 0;
}

static void
numtrunc(Num *num)
{
	num->wp = num->rp;
	if (num->rp != num->buf)
		num->rp--;
}

static int
more(Num *num)
{
	return (num->rp != num->wp);
}

static int
length(Num *num)
{
	return num->wp - num->buf;
}

static int
tell(Num *num)
{
	return num->rp - num->buf;
}

static void
seek(Num *num, int pos)
{
	num->rp = num->buf + pos;
}

static void
rshift(Num *num, int n)
{
	int diff;

	diff = length(num) - n;
	if (diff < 0) {
		first(num);
		numtrunc(num);
		return;
	}

	memmove(num->buf, num->buf + n, diff);
	num->rp = num->buf + diff;
	numtrunc(num);
}

static void
lshift(Num *num, int n)
{
	int len;

	len = length(num);
	expand(num, len + n);
	memmove(num->buf + n, num->buf, len);
	memset(num->buf, 0, n);
	num->wp += n;
}

static Num *
chsign(Num *num)
{
	int val, d, carry;

	carry = 0;
	for (first(num); more(num); next(num)) {
		d = peek(num);
		val = 100 - d - carry;

		carry = 1;
		if (val >= 100) {
			val -= 100;
			carry = 0;
		}
		poke(num, val);
	}

	prev(num);
	if (carry != 0) {
		if (peek(num) == 99)
			poke(num, -1);
		else
			wrdigit(num, -1);
	} else {
		if (peek(num) == 0)
			numtrunc(num);
	}

	return num;
}

static Num *
copy(Num *num)
{
	Num *p;
	int len = length(num);

	p = new(len);
	memcpy(p->buf, num->buf, len);
	p->wp = p->buf + len;
	p->rp = p->buf;
	p->scale = num->scale;

	return p;
}

static int
negative(Num *num)
{
	return last(num) && peek(num) == -1;
}

static Num *
norm(Num *n)
{
	/* trailing 0 */
	for (last(n); peek(n) == 0; numtrunc(n))
		;

	if (negative(n)) {
		for (prev(n); peek(n) == 99; prev(n)) {
			poke(n, -1);
			next(n);
			numtrunc(n);
		}
	}

	/* adjust scale for 0 case*/
	if (length(n) == 0)
		n->scale = 0;
	return n;
}

static Num *
mulnto(Num *src, Num *dst, int n)
{
	div_t dd;
	int d, carry;

	first(dst);
	numtrunc(dst);

	carry = 0;
	for (first(src); more(src); next(src)) {
		d = peek(src) * n + carry;
		dd = div(d, 100);
		carry = dd.quot;
		wrdigit(dst, dd.rem);
	}

	if (carry)
		wrdigit(dst, carry);
	return dst;
}

static Num *
muln(Num *num, int n)
{
	div_t dd;
	int d, carry;

	carry = 0;
	for (first(num); more(num); next(num)) {
		d = peek(num) * n + carry;
		dd = div(d, 100);
		carry = dd.quot;
		poke(num, dd.rem);
	}

	if (carry)
		wrdigit(num, carry);
	return num;
}

static int
divn(Num *num, int n)
{
	div_t dd;
	int val, carry;

	carry = 0;
	for (last(num); !begin(num); prev(num)) {
		val = carry * 100 + peek(num);
		dd = div(val, n);
		poke(num, dd.quot);
		carry = dd.rem;
	}
	norm(num);

	return carry;
}

static void
div10(Num *num, int n)
{
	div_t dd = div(n, 2);

	if (dd.rem == 1)
		divn(num, 10);

	rshift(num, dd.quot);
}

static void
mul10(Num *num, int n)
{
	div_t dd = div(n, 2);

	if (dd.rem == 1)
		muln(num, 10);

	lshift(num, dd.quot);
}

static void
align(Num *a, Num *b)
{
	int d;
	Num *max, *min;

	d = a->scale - b->scale;
	if (d == 0) {
		return;
	} else if (d > 0) {
		min = b;
		max = a;
	} else {
		min = a;
		max = b;
	}

	d = abs(d);
	mul10(min, d);
	min->scale += d;

	assert(min->scale == max->scale);
}

static Num *
addn(Num *num, int val)
{
	int d, carry = val;

	for (first(num); carry; next(num)) {
		d = more(num) ? peek(num) : 0;
		d += carry;
		carry = 0;

		if (d >= 100) {
			carry = 1;
			d -= 100;
		}
		poke(num, d);
	}

	return num;
}

static Num *
reverse(Num *num)
{
	int d;
	signed char *p, *q;

	for (p = num->buf, q = num->wp - 1; p < q; ++p, --q) {
		d = *p;
		*p = *q;
		*q = d;
	}

	return num;
}

static Num *
addnum(Num *a, Num *b)
{
	Num *c;
	int len, alen, blen, carry, da, db, sum;

	align(a, b);
	alen = length(a);
	blen = length(b);
	len = (alen > blen) ? alen : blen;
	c = new(len);

	first(a);
	first(b);
	carry = 0;
	while (len-- > 0) {
		da = (more(a)) ? rddigit(a) : 0;
		db = (more(b)) ? rddigit(b) : 0;

		sum = da + db + carry;
		if (sum >= 100) {
			carry = 1;
			sum -= 100;
		} else if (sum < 0) {
			carry = -1;
			sum += 100;
		} else {
			carry = 0;
		}

		wrdigit(c, sum);
	}

	if (carry)
		wrdigit(c, carry);
	c->scale = a->scale;

	return norm(c);
}

static Num *
subnum(Num *a, Num *b)
{
	Num *tmp, *sum;

	tmp = chsign(copy(b));
	sum = addnum(a, tmp);
	freenum(tmp);

	return sum;
}

static Num *
mulnum(Num *a, Num *b)
{
	Num shadow, *c, *ca, *cb;
	int pos, prod, carry, dc, db, da, sc;
	int asign = negative(a), bsign = negative(b);

	c = zeronum(length(a) + length(b) + 1);
	c->scale = a->scale + b->scale;
	sc = (a->scale > b->scale) ? a->scale : b->scale;

	ca = a;
	if (asign)
		ca = chsign(copy(ca));
	cb = b;
	if (bsign)
		cb = chsign(copy(cb));

	/* avoid aliasing problems when called from expnum*/
	if (ca == cb) {
		shadow = *cb;
		b = cb = &shadow;
	}

	for (first(cb); more(cb); next(cb)) {
		div_t d;

		carry = 0;
		db = peek(cb);

		pos = tell(c);
		for (first(ca); more(ca); next(ca)) {
			da = peek(ca);
			dc = peek(c);
			prod = da * db + dc + carry;
			d = div(prod, 100);
			carry = d.quot;
			poke(c, d.rem);
			next(c);
		}

		for ( ; carry > 0; carry = d.quot) {
			dc = peek(c) + carry;
			d = div(dc, 100);
			poke(c, d.rem);
			next(c);
		}
		seek(c, pos + 1);
	}
	norm(c);

	/* Adjust scale to max(a->scale, b->scale) while c is still positive */
	sc = c->scale - sc;
	if (sc > 0) {
		div10(c, sc);
		c->scale -= sc;
	}

	if (ca != a)
		freenum(ca);
	if (cb != b)
		freenum(cb);

	if (asign ^ bsign)
		chsign(c);
	return c;
}

/*
 * The divmod function is implemented following the algorithm
 * from the plan9 version that is not exactly like the one described
 * in the paper. A lot of magic here.
 */
static Num *
divmod(Num *odivd, Num *odivr, Num **remp)
{
	Num *acc, *divd, *divr, *res;
	int divsign, remsign;
	int under, magic, ndig, diff;
	int d, q, carry, divcarry;
	long dr, dd, cc;

	divr = odivr;
	acc = copy(&zero);
	divd = copy(odivd);
	res = zeronum(length(odivd));

	under = divcarry = divsign = remsign = 0;

	if (length(divr) == 0) {
		weprintf("divide by 0\n");
		goto ret;
	}

	divsign = negative(divd);
	if (divsign)
		chsign(divd);

	remsign = negative(divr);
	if (remsign)
		divr = chsign(copy(divr));

	diff = length(divd) - length(divr);

	seek(res, diff + 1);
	last(divd);
	last(divr);

	wrdigit(divd, 0);

	dr = peek(divr);
	magic = dr < 10;
	dr = dr * 100 + (prev(divr) ? peek(divr) : 0);
	if (magic) {
		dr = dr * 100 + (prev(divr) ? peek(divr) : 0);
		dr *= 2;
		dr /= 25;
	}

	for (ndig = 0; diff >= 0; ++ndig) {
		last(divd);
		dd = peek(divd);
		dd = dd * 100 + (prev(divd) ? peek(divd) : 0);
		dd = dd * 100 + (prev(divd) ? peek(divd) : 0);
		cc = dr;

		if (diff == 0)
			dd++;
		else
			cc++;

		if (magic)
			dd *= 8;

		q = dd / cc;
		under = 0;
		if (q > 0 && dd % cc < 8 && magic) {
			q--;
			under = 1;
		}

		mulnto(divr, acc, q);

		/* Subtract acc from dividend at offset position */
		first(acc);
		carry = 0;
		for (seek(divd, diff); more(divd); next(divd)) {
			d = peek(divd);
			d = d - (more(acc) ? rddigit(acc) : 0) - carry;
			carry = 0;
			if (d < 0) {
				d += 100;
				carry = 1;
			}
			poke(divd, d);
		}
		divcarry = carry;

		/* Store quotient digit */
		prev(res);
		poke(res, q);

		/* Handle borrow propagation */
		last(divd);
		d = peek(divd);
		if ((d != 0) && (diff != 0)) {
			prev(divd);
			d = peek(divd) + 100;
			poke(divd, d);
		}

		/* Shorten dividend for next iteration */
		if (--diff >= 0)
			divd->wp--;
	}

	/*
	 * if we have an underflow then we have to adjust
	 * the remaining and the result
	 */
	if (under) {
		Num *p = subnum(divd, divr);
		if (negative(p)) {
			freenum(p);
		} else {
			freenum(divd);
			poke(res, q + 1);
			divd = p;
		}
	}

	if (divcarry) {
		Num *p;

		poke(res, q - 1);
		poke(divd, -1);
		p = addnum(divr, divd);
		freenum(divd);
		divd = p;
	}

	divcarry = 0;
	for (first(res); more(res); next(res)) {
		d = peek(res) + divcarry;
		divcarry = 0;
		if (d >= 100) {
			d -= 100;
			divcarry = 1;
		}
		poke(res, d);
	}

ret:
	if (divsign)
		chsign(divd);
	if (divsign ^ remsign)
		chsign(res);

	if (remp) {
		divd->scale = odivd->scale;
		*remp = norm(divd);
	} else {
		freenum(divd);
	}

	if (divr != odivr)
		freenum(divr);

	freenum(acc);

	res->scale = odivd->scale - odivr->scale;
	if (res->scale < 0)
		res->scale = 0;

	return norm(res);
}

static int
divscale(Num *divd, Num *divr)
{
	int diff;

	if (length(divr) == 0) {
		weprintf("divide by 0\n");
		return 0;
	}

	diff = scale + divr->scale - divd->scale;

	if (diff > 0) {
		mul10(divd, diff);
		divd->scale += diff;
	} else if (diff < 0) {
		mul10(divr, -diff);
		divr->scale += -diff;
	}

	return 1;
}

static Num *
divnum(Num *a, Num *b)
{
	if (!divscale(a, b))
		return copy(&zero);

	return divmod(a, b, NULL);
}

static Num *
modnum(Num *a, Num *b)
{
	Num *mod, *c;

	if (!divscale(a, b))
		return copy(&zero);

	c = divmod(a, b, &mod);
	freenum(c);

	return mod;
}

static Num *
expnum(Num *base, Num *exp)
{
	int neg, d;
	Num *res, *fact, *e, *tmp1, *tmp2;

	res = copy(&one);
	if (length(exp) == 0)
		return res;

	e = copy(exp);
	if ((neg = negative(exp)) != 0)
		chsign(e);

	if (e->scale > 0) {
		div10(e, e->scale);
		e->scale = 0;
	}

	fact = copy(base);
	while (length(e) > 0) {
		first(e);
		d = peek(e);
		if (d % 2 == 1) {
			tmp1 = mulnum(res, fact);
			freenum(res);
			res = tmp1;
		}

		/* Square fact */
		tmp1 = mulnum(fact, fact);
		freenum(fact);
		fact = tmp1;

		divn(e, 2);
	}
	freenum(fact);
	freenum(e);

	/* Handle negative exponent: 1 / res */
	if (neg) {
		tmp2 = divnum(tmp1 = copy(&one), res);
		freenum(tmp1);
		freenum(res);
		res = tmp2;
	}

	return res;
}

/*
 * Compare two numbers: returns <0 if a<b, 0 if a==b, >0 if a>b
 */
static int
cmpnum(Num *a, Num *b)
{
	Num *diff;
	int result;

	diff = subnum(a, b);
	if (length(diff) == 0)
		result = 0;
	else if (negative(diff))
		result = -1;
	else
		result = 1;
	freenum(diff);

	return result;
}

/*
 * Integer square root of a small integer (0-9999)
 * Used for initial guess in Newton's method
 */
static int
isqrt(int n)
{
	int x, x1;

	if (n <= 0)
		return 0;
	if (n == 1)
		return 1;

	x = n;
	x1 = (x + 1) / 2;
	while (x1 < x) {
		x = x1;
		x1 = (x + n / x) / 2;
	}
	return x;
}

/*
 * Square root using Newton's method: x_{n+1} = (x_n + y/x_n) / 2
 *
 * Key insight: sqrt(a * 10^(2n)) = sqrt(a) * 10^n
 * So we scale up the input to get the desired output precision.
 *
 * To compute sqrt with scale decimal places of precision:
 * 1. Scale up y by 10^(2*scale + 2) (extra 2 for guard digits)
 * 2. Compute integer sqrt
 * 3. Result has (scale + 1) decimal places, numtrunc to scale
 */
static Num *
sqrtnum(Num *oy)
{
	Num *y, *x, *xprev, *q, *sum;
	int top, ysc, iter;

	if (length(oy) == 0)
		return copy(&zero);

	if (negative(oy)) {
		weprintf("square root of negative number\n");
		return copy(&zero);
	}

	y = copy(oy);
	ysc = 2 * scale + 2 - y->scale;
	if (ysc > 0)
		mul10(y, ysc);
	ysc = 2 * scale + 2;

	/* Make scale even (so sqrt gives integer result) */
	if (ysc % 2 == 1) {
		muln(y, 10);
		ysc++;
	}
	y->scale = 0;

	last(y);
	top = peek(y);
	if (prev(y) && length(y) > 1)
		top = top * 100 + peek(y);

	x = new(0);
	wrdigit(x, isqrt(top));
	x->scale = 0;

	/* Scale up the initial guess to match the magnitude of y */
	lshift(x, (length(y) - 1) / 2);


	/* Newton iteration: x = (x + y/x) / 2 */
	xprev = NULL;
	for (iter = 0; iter < 1000; iter++) {
		q = divmod(y, x, NULL);
		sum = addnum(x, q);
		freenum(q);
		divn(sum, 2);

		/* Check for convergence: sum == x or sum == prev */
		if (cmpnum(sum, x) == 0) {
			freenum(sum);
			break;
		}
		if (xprev != NULL && cmpnum(sum, xprev) == 0) {
			/* Oscillating - pick smaller */
			if (cmpnum(x, sum) < 0) {
				freenum(sum);
			} else {
				freenum(x);
				x = sum;
			}
			break;
		}

		freenum(xprev);
		xprev = x;
		x = sum;
	}
	freenum(xprev);
	freenum(y);

	/* Truncate to desired scale */
	x->scale = ysc / 2;
	if (x->scale > scale) {
		int diff = x->scale - scale;
		div10(x, diff);
		x->scale = scale;
	}

	return norm(x);
}

static Num *
tonum(void)
{
	char *s, *t, *end, *dot;
	Num *num, *denom, *numer, *frac, *q, *rem;
	int sign, d, ch, nfrac;

	s = input->s;
	num = new(0);
	sign = 0;
	if (*s == '_') {
		sign = 1;
		++s;
	}

	dot = NULL;
	for (t = s; (ch = *t) > 0 || ch <= UCHAR_MAX; ++t) {
		if (!strchr(digits, ch))
			break;
		if (ch == '.') {
			if (dot)
				break;
			dot = t;
		}
	}
	input->s = end = t;

	/*
	 * Parse integer part: process digits left-to-right.
	 * For each digit: num = num * ibase + digit
	 */
	for (t = s; t < (dot ? dot : end); ++t) {
		d = strchr(digits, *t) - digits;
		muln(num, ibase);
		addn(num, d);
	}
	norm(num);

	if (!dot)
		goto ret;

	/*
	 * convert fractional digits
	 * Algorithm: For digits d[0], d[1], ..., d[n-1] after '.'
	 * Value = d[0]/ibase + d[1]/ibase^2 + ... + d[n-1]/ibase^n
	 *
	 * numerator = d[0]*ibase^(n-1) + d[1]*ibase^(n-2) + ... + d[n-1]
	 * denominator = ibase^n
	 * Then extract decimal digits by repeated: num*100/denom
	 */
	denom = copy(&one);
	numer = copy(&zero);
	for (t = dot + 1; t < end; ++t) {
		d = strchr(digits, *t) - digits;
		muln(denom, ibase);
		muln(numer, ibase);
		addn(numer, d);
	}

	nfrac = end - dot - 1;
	frac = new(0);
	d = 0;
	while (frac->scale < nfrac || length(numer) > 0) {
		muln(numer, 100);
		q = divmod(numer, denom, &rem);
		freenum(numer);

		d = first(q) ? peek(q) : 0;
		wrdigit(frac, d);
		freenum(q);
		numer = rem;
		frac->scale += 2;
	}
	reverse(frac);

	/* Trim to exact input scale for odd nfrac */
	if (frac->scale > nfrac && d % 10 == 0) {
		divn(frac, 10);
		frac->scale--;
	}

	freenum(numer);
	freenum(denom);

	q = addnum(num, frac);
	freenum(num);
	freenum(frac);
	num = q;

ret:
	if (sign)
		chsign(num);
	return num;
}

static void
prchr(int ch)
{
	if (col >= 69) {
		putchar('\\');
		putchar('\n');
		col = 0;
	}
	putchar(ch);
	col++;
}

static void
printd(int d, int base, int space)
{
	int w, n;

	if (base <= 16) {
		prchr(digits[d]);
	} else {
		if (space)
			prchr(' ');

		for (w = 1, n = base - 1; n >= 10; n /= 10)
			w++;

		if (col + w > 69) {
			putchar('\\');
			putchar('\n');
			col = 0;
		}
		col += printf("%0*d", w, d);
	}
}

static void
pushdigit(struct digit **l, int val)
{
	struct digit *it = emalloc(sizeof(*it));

	it->next = *l;
	it->val = val;
	*l = it;
}

static int
popdigit(struct digit **l)
{
	int val;
	struct digit *next, *it = *l;

	if (it == NULL)
		return -1;

	val = it->val;
	next = it->next;
	free(it);
	*l = next;
	return val;
}

static void
printnum(Num *onum, int base)
{
	struct digit *sp;
	int sc, i, sign, n;
	Num *num, *inte, *frac, *opow;

	col = 0;
	if (length(onum) == 0) {
		prchr('0');
		return;
	}

	num = copy(onum);
	if ((sign = negative(num)) != 0)
		chsign(num);

	sc = num->scale;
	if (num->scale % 2 == 1) {
		muln(num, 10);
		num->scale++;
	}
	inte = copy(num);
	rshift(inte, num->scale / 2);
	inte->scale = 0;
	frac = subnum(num, inte);

	sp = NULL;
	for (i = 0; length(inte) > 0; ++i)
		pushdigit(&sp, divn(inte, base));
	if (sign)
		prchr('-');
	while (i-- > 0)
		printd(popdigit(&sp), base, 1);
	assert(sp == NULL);

	if (num->scale == 0)
		goto ret;

        /*
         * Print fractional part by repeated multiplication by base.
         * We maintain the fraction as: frac / 10^scale
         *
         * Algorithm:
         * 1. Multiply frac by base
         * 2. Output integer part (frac / 10^scale)
         * 3. Keep fractional part (frac % 10^scale)
         */
	prchr('.');

	opow = copy(&one);
	mul10(opow, num->scale);

	for (n = 0; n < sc; ++n) {
		int d;
		Num *q, *rem;

		muln(frac, base);
		q = divmod(frac, opow, &rem);
		d = first(q) ? peek(q) : 0;
		freenum(frac);
		freenum(q);
		frac = rem;
		printd(d, base, n > 0);
	}
	freenum(opow);

ret:
	freenum(num);
	freenum(inte);
	freenum(frac);
}

static int
moreinput(void)
{
	struct input *ip;

repeat:
	if (!input)
		return 0;

	if (input->buf != NULL && *input->s != '\0')
		return 1;

	if (input->fp) {
		if (getline(&input->buf, &input->n, input->fp) >= 0) {
			input->s = input->buf;
			return 1;
		}
		if (ferror(input->fp)) {
			eprintf("reading from file:");
			exit(1);
		}
		fclose(input->fp);
	}

	ip = input;
	input = ip->next;
	free(ip->buf);
	free(ip);
	goto repeat;
}

static void
addinput(FILE *fp, char *s)
{
	struct input *ip;

	assert((!fp && !s) == 0);

	ip = emalloc(sizeof(*ip));
	ip->next = input;
	ip->fp = fp;
	ip->n = 0;
	ip->s = ip->buf = s;
	input = ip;
}

static void
delinput(int cmd, int n)
{
	if (n < 0)
		error("Q command requires a number >= 0");
	while (n-- > 0) {
		if (cmd == 'Q' && !input->next)
			error("Q command argument exceeded string execution depth");
		if (input->fp)
			fclose(input->fp);
		free(input->buf);
		input = input->next;
		if (!input)
			exit(0);
	}
}

static void
push(Val v)
{
	Val *p = emalloc(sizeof(Val));

	*p = v;
	p->next = stack;
	stack = p;
}

static void
needstack(int n)
{
	Val *vp;

	for (vp = stack; n > 0 && vp; vp = vp->next)
		--n;
	if (n > 0)
		error("stack empty");
}

static Val
pop(void)
{
	Val v;

	if (!stack)
		error("stack empty");
	v = *stack;
	free(stack);
	stack = v.next;

	return v;
}

static Num *
popnum(void)
{
	Val v = pop();

	if (v.type != NUM) {
		free(v.u.s);
		error("non-numeric value");
	}
	return v.u.n;
}

static void
pushnum(Num *num)
{
	push((Val) {.type = NUM, .u.n = num});
}

static void
pushstr(char *s)
{
	push((Val) {.type = STR, .u.s = s});
}

static void
arith(Num *(*fn)(Num *, Num *))
{
	Num *a, *b, *c;

	needstack(2);
	b = popnum();
	a = popnum();
	c = (*fn)(a, b);
	freenum(a);
	freenum(b);
	pushnum(c);
}

static void
pushdivmod(void)
{
	Num *a, *b, *q, *rem;

	needstack(2);
	b = popnum();
	a = popnum();

	if (!divscale(a, b)) {
		q = copy(&zero);
		rem = copy(&zero);
	} else {
		q = divmod(a, b, &rem);
	}

	pushnum(q);
	pushnum(rem);
	freenum(a);
	freenum(b);
}

static int
popint(void)
{
	Num *num;
	int r = -1, n, d;

	num = popnum();
	if (negative(num))
		goto ret;

	/* discard fraction part */
	div10(num, num->scale);

	n = 0;
	for (last(num); !begin(num); prev(num)) {
		if (n > INT_MAX / 100)
			goto ret;
		n *= 100;
		d = peek(num);
		if (n > INT_MAX - d)
			goto ret;
		n += d;
	}
	r = n;

ret:
	freenum(num);
	return r;
}

static void
pushint(int n)
{
	div_t dd;
	Num *num;

	num = new(0);
	for ( ; n > 0; n = dd.quot) {
		dd = div(n, 100);
		wrdigit(num, dd.rem);
	}
	pushnum(num);
}

static void
printval(Val v)
{
	if (v.type == STR)
		fputs(v.u.s, stdout);
	else
		printnum(v.u.n, obase);
}

static Val
dupval(Val v)
{
	Val nv;

	switch (nv.type = v.type) {
	case STR:
		nv.u.s = estrdup(v.u.s);
		break;
	case NUM:
		nv.u.n = copy(v.u.n);
		break;
	case NVAL:
		nv.type = NUM;
		nv.u.n = copy(&zero);
		break;
	}

	return nv;
}

static void
freeval(Val v)
{
	if (v.type == STR)
		free(v.u.s);
	else if (v.type == NUM)
		freenum(v.u.n);
}

static void
dumpstack(void)
{
	Val *vp;

	for (vp = stack; vp; vp = vp->next) {
		printval(*vp);
		putchar('\n');
	}
}

static void
clearstack(void)
{
	Val *vp, *next;

	for (vp = stack; vp; vp = next) {
		next = vp->next;
		freeval(*vp);
		free(vp);
	}
	stack = NULL;
}

static void
dupstack(void)
{
	Val v;

	push(v = pop());
	push(dupval(v));
}

static void
deepstack(void)
{
	int n;
	Val *vp;

	n = 0;
	for (vp = stack; vp; vp = vp->next) {
		if (n == INT_MAX)
			error("stack depth does not fit in a integer");
		++n;
	}
	pushint(n);
}

static void
pushfrac(void)
{
	Val v = pop();

	if (v.type == STR)
		pushint(0);
	else
		pushint(v.u.n->scale);
	freeval(v);
}

static void
pushlen(void)
{
	int n;
	Num *num;
	Val v = pop();

	if (v.type == STR) {
		n = strlen(v.u.s);
	} else {
		num = v.u.n;
		if (length(num) == 0) {
			n = 1;
		} else {
			n = length(num) * 2;
			n -= last(num) ? peek(num) < 10 : 0;
		}
	}
	pushint(n);
	freeval(v);
}

static void
setibase(void)
{
	int n = popint();

	if (n < 2 || n > 16)
		error("input base must be an integer between 2 and 16");
	ibase = n;
}

static void
setobase(void)
{
	int n = popint();

	if (n < 2)
		error("output base must be an integer greater than 1");
	obase = n;
}

static char *
string(char *dst, int *np)
{
	int n, ch;

	n = np ? *np : 0;
	for (;;) {
		ch = *input->s++;

		switch (ch) {
		case '\0':
			if (!moreinput())
				exit(0);
			break;
		case '\\':
			if (*input->s == '[') {
				dst = erealloc(dst, n + 1);
				dst[n++] = *input->s++;
				break;
			}
			goto copy;
		case ']':
			if (!np) {
				dst = erealloc(dst, n + 1);
				dst[n] = '\0';
				return dst;
			}
		case '[':
		default:
		copy:
			dst = erealloc(dst, n + 1);
			dst[n++] = ch;
			if (ch == '[')
				dst = string(dst, &n);
			if (ch == ']') {
				*np = n;
				return dst;
			}
		}
	}
}

static void
setscale(void)
{
	int n = popint();

	if (n < 0)
		error("scale must be a nonnegative integer");
	scale = n;
}

static unsigned
hash(char *name)
{
	int c;
	unsigned h = 5381;

	while (c = *name++)
		h = h*33 ^ c;

	return h;
}

static struct reg *
lookup(char *name)
{
	struct reg *rp;
	int h = hash(name) & HASHSIZ-1;

	for (rp = htbl[h]; rp; rp = rp->next) {
		if (strcmp(name, rp->name) == 0)
			return rp;
	}

	rp = emalloc(sizeof(*rp));
	rp->next = htbl[h];
	htbl[h] = rp;
	rp->name = estrdup(name);

	rp->val.type = NVAL;
	rp->val.next = NULL;

	rp->ary.n = 0;
	rp->ary.buf = NULL;
	rp->ary.next = NULL;

	return rp;
}

static char *
regname(void)
{
	int delim, ch;
	char *s;
	static char name[REGSIZ];

	ch = *input->s++;
	if (!iflag || ch != '<' && ch != '"') {
		name[0] = ch;
		name[1] = '\0';
		return name;
	}

	if ((delim = ch) == '<')
		delim = '>';

	for (s = name; s < &name[REGSIZ]; ++s) {
		ch = *input->s++;
		if (ch == '\0' ||  ch == delim) {
			*s = '\0';
			if (ch == '>') {
				name[0] = atoi(name);
				name[1] = '\0';
			}
			return name;
		}
		*s = ch;
	}

	error("identifier too long");
}

static void
popreg(void)
{
	int i;
	Val *vnext;
	struct ary *anext;
	char *s = regname();
	struct reg *rp = lookup(s);

	if (rp->val.type == NVAL)
		error("stack register '%s' (%o) is empty", s, s[0]);

	push(rp->val);
	vnext = rp->val.next;
	if (!vnext) {
		rp->val.type = NVAL;
	} else {
		rp->val = *vnext;
		free(vnext);
	}

	for (i = 0; i < rp->ary.n; ++i)
		freeval(rp->ary.buf[i]);
	free(rp->ary.buf);

	anext = rp->ary.next;
	if (!anext) {
		rp->ary.n = 0;
		rp->ary.buf = NULL;
	} else {
		rp->ary = *anext;
		free(anext);
	}
}

static void
pushreg(void)
{
	Val v;
	Val *vp;
	struct ary *ap;
	struct reg *rp = lookup(regname());

	v = pop();

	vp = emalloc(sizeof(Val));
	*vp = rp->val;
	rp->val = v;
	rp->val.next = vp;

	ap = emalloc(sizeof(struct ary));
	*ap = rp->ary;
	rp->ary.n = 0;
	rp->ary.buf = NULL;
	rp->ary.next = ap;
}

static Val *
aryidx(void)
{
	int n;
	int i;
	Val *vp;
	struct reg *rp = lookup(regname());
	struct ary *ap = &rp->ary;

	n = popint();
	if (n < 0)
		error("array index must fit in a positive integer");

	if (n >= ap->n) {
		ap->buf = ereallocarray(ap->buf, n+1, sizeof(Val));
		for (i = ap->n; i <= n; ++i)
			ap->buf[i].type = NVAL;
		ap->n = n + 1;
	}
	return &ap->buf[n];
}

static void
aryget(void)
{
	Val *vp = aryidx();

	push(dupval(*vp));
}

static void
aryset(void)
{
	Val val, *vp = aryidx();

	val = pop();
	freeval(*vp);
	*vp = val;
}

static void
execmacro(void)
{
	int ch;
	Val v = pop();

	assert(v.type != NVAL);
	if (v.type == NUM) {
		push(v);
		return;
	}

	if (input->fp) {
		addinput(NULL, v.u.s);
		return;
	}

	/* check for tail recursion */
	for (ch = *input->s; ch > 0 && ch < UCHAR_MAX; ch = *input->s) {
		if (!isspace(ch))
			break;
		++input->s;
	}

	if (ch == '\0') {
		free(input->buf);
		input->buf = input->s = v.u.s;
		return;
	}

	addinput(NULL, v.u.s);
}

static void
relational(int ch)
{
	int r;
	char *s;
	Num *a, *b;
	struct reg *rp = lookup(regname());

	needstack(2);
	a = popnum();
	b = popnum();
	r = cmpnum(a, b);
	freenum(a);
	freenum(b);

	switch (ch) {
	case '>':
		r = r > 0;
		break;
	case '<':
		r = r < 0;
		break;
	case '=':
		r = r == 0;
		break;
	case LE:
		r = r <= 0;
		break;
	case GE:
		r = r >= 0;
		break;
	case NE:
		r = r != 0;
		break;
	default:
		abort();
	}

	if (!r)
		return;

	push(dupval(rp->val));
	execmacro();
}

static void
printbytes(void)
{
	Num *num;
	Val v = pop();

	if (v.type == STR) {
		fputs(v.u.s, stdout);
	} else {
		num = v.u.n;
		div10(num, num->scale);
		num->scale = 0;
		printnum(num, 100);
	}
	freeval(v);
}

static void
eval(void)
{
	int ch;
	char *s;
	Num *num;
	size_t siz;
	Val v1, v2;
	struct reg *rp;

	if (setjmp(env))
		return;

	for (s = input->s; (ch = *s) != '\0'; ++s) {
		if (ch < 0 || ch > UCHAR_MAX || !isspace(ch))
			break;
	}
	input->s = s + (ch != '\0');

	switch (ch) {
	case '\0':
		break;
	case 'n':
		v1 = pop();
		printval(v1);
		freeval(v1);
		break;
	case 'p':
		v1 = pop();
		printval(v1);
		putchar('\n');
		push(v1);
		break;
	case 'P':
		printbytes();
		break;
	case 'f':
		dumpstack();
		break;
	case '+':
		arith(addnum);
		break;
	case '-':
		arith(subnum);
		break;
	case '*':
		arith(mulnum);
		break;
	case '/':
		arith(divnum);
		break;
	case '%':
		arith(modnum);
		break;
	case '^':
		arith(expnum);
		break;
	case '~':
		pushdivmod();
		break;
	case 'v':
		num = popnum();
		pushnum(sqrtnum(num));
		freenum(num);
		break;
	case 'c':
		clearstack();
		break;
	case 'd':
		dupstack();
		break;
	case 'r':
		needstack(2);
		v1 = pop();
		v2 = pop();
		push(v1);
		push(v2);
		break;
	case 'S':
		pushreg();
		break;
	case 'L':
		popreg();
		break;
	case 's':
		rp = lookup(regname());
		freeval(rp->val);
		rp->val = pop();
		break;
	case 'l':
		rp = lookup(regname());
		push(dupval(rp->val));
		break;
	case 'i':
		setibase();
		break;
	case 'o':
		setobase();
		break;
	case 'k':
		setscale();
		break;
	case 'I':
		pushint(ibase);
		break;
	case 'O':
		pushint(obase);
		break;
	case 'K':
		pushint(scale);
		break;
	case '[':
		pushstr(string(NULL, NULL));
		break;
	case 'x':
		execmacro();
		break;
	case '!':
		switch (*input->s++) {
		case '<':
			ch = GE;
			break;
		case '>':
			ch = LE;
			break;
		case '=':
			ch = NE;
			break;
		default:
			system(input->s-1);
			goto discard;
		}
	case '>':
	case '<':
	case '=':
		relational(ch);
		break;
	case '?':
		s = NULL;
		if (getline(&s, &siz, stdin) > 0) {
			pushstr(s);
		} else {
			free(s);
			if (ferror(stdin))
				eprintf("reading from file\n");
		}
		break;
	case 'q':
		delinput('q', 2);
		break;
	case 'Q':
		delinput('Q', popint());
		break;
	case 'Z':
		pushlen();
		break;
	case 'X':
		pushfrac();
		break;
	case 'z':
		deepstack();
		break;
	case '#':
	discard:
		while (*input->s)
			++input->s;
		break;
	case ':':
		aryset();
		break;
	case ';':
		aryget();
		break;
	default:
		if (!strchr(digits, ch))
			error("'%c' (%#o) unimplemented", ch, ch);
	case '_':
		--input->s;
		pushnum(tonum());
		break;
	}
}

static void
dc(char *fname)
{
	FILE *fp;

	if (strcmp(fname, "-") == 0) {
		fp = stdin;
	} else {
		if ((fp = fopen(fname, "r")) == NULL)
			eprintf("opening '%s':", fname);
	}
	addinput(fp, NULL);

	while (moreinput())
		eval();

	fclose(fp);
	free(input);
	input = NULL;
}

static void
usage(void)
{
	eprintf("usage: dc [-i] [file ...]\n");
}

int
main(int argc, char *argv[])
{
	ARGBEGIN {
	case 'i':
		iflag = 1;
		break;
	default:
		usage();
	} ARGEND

	while (*argv)
		dc(*argv++);
	dc("-");

	return 0;
}