4.3BSD/usr/contrib/B/src/bsmall/B1num.c

/* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1984. */
/* $Header: B1num.c,v 1.1 84/06/28 00:48:56 timo Exp $ */

/* B numbers, small version */

	/*
	 * THIS VERSION SHOULD ONLY BE USED IF
	 * THE SYSTEM IS TOO LARGE OTHERWISE.
	 * IT USES FLOATING POINT ARITHMETIC FOR EXACT NUMBERS
	 * INSTEAD OF ARBITRARY LENGTH RATIONAL ARITHMETIC.
	 */

#include "b.h"
#include "b0con.h"
#include "b1obj.h"
#include "b2syn.h" /* for def of Keymark() only */
#include "B1num.h"

value numerator(v) value v; {
	Checknum(v);
	if (!Exact(v)) error("*/ on approximate number");
	return mk_int(Numerator(v));
}

value denominator(v) value v; {
	Checknum(v);
	if (!Exact(v)) error("/* on approximate number"); /* */
	return mk_int(Denominator(v));
}

double numval(v) value v; {
	Checknum(v);
	return Numval(v);
}

checkint(v) value v; {
	Checknum(v);
	if (Denominator(v) != One) error("number not an integer");
}

bool large(v) value v; {
	checkint(v);
	if (Numerator(v) < -Maxint || Numerator(v) > Maxint) return Yes;
	return No;
}
	
int intval(v) value v; {
	checkint(v);
	return (int)Numerator(v);
}

intlet propintlet(i) int i; {
	if (i < -Maxintlet || i > Maxintlet)
		error("exceedingly large integer");
	return i;
}

integer gcd(i, j) integer i, j; {
	integer k;
	if (i == Zero && j == Zero) syserr("gcd(0, 0)");
	if (i != floor(i) || j != floor(j))
		syserr("gcd called with non-integer");
	if (i < Zero) i= -i; if (j < Zero) j= -j;
	if (i < j) {
		k= i; i= j; j= k;
	}
	while (j >= One) {
		k= i-j*floor(i/j);
		i= j; j= k;
	}
	if (j != Zero) error(
		"arithmetic overflow while simplifying exact number");
	if (i != floor(i)) syserr("gcd returns non-integer");
	return i;
}

value b_zero, b_one, b_minus_one, zero, one;

value mk_exact(p, q, len) register integer p, q; intlet len; {
	value v; integer d;
	if (q == One && len ==0) {
		if (p == Zero) return copy(b_zero);
		if (p == One)  return copy(b_one);
		if (p == -One) return copy(b_minus_one);
	}
	v= grab_num(len);
	if (q == One) {
		Numerator(v)= p; Denominator(v)= q;
		return v;
	}
	if (q == Zero) error("attempt to make exact number with denominator 0");
	if (q < Zero) {p= -p; q= -q;}
	d= (q == One ? One : p == One ? One : gcd(p, q));
	Numerator(v)= p/d; Denominator(v)= q/d;
	return v;
}

bool integral(v) value v; {
	return Integral(v);
}

value mk_integer(p) int p; {
	return mk_exact((integer)p, One, 0);
}

value mk_int(p) integer p; {
	return mk_exact(p, One, 0);
}

value mk_approx(x) register double x; {
	value v= grab_num(0);
	Approxval(v)= x; Denominator(v)= Zero;
	return v;
}

initnum() {
	b_zero= grab_num(0);
		Numerator(b_zero)= Zero; Denominator(b_zero)= One;
	b_one= grab_num(0);
		Numerator(b_one)= One; Denominator(b_one)= One;
	b_minus_one= grab_num(0);
		Numerator(b_minus_one)= -One; Denominator(b_minus_one)= One;
	zero= mk_integer(0);
	one= mk_integer(1);
}

value approximate(v) value v; {
	if (!Exact(v)) return copy(v);
	return mk_approx(Numerator(v)/Denominator(v));
}

numcomp(v, w) value v, w; {
	double vv= Numval(v), ww= Numval(w);
	if (vv < ww) return -1;
	if (vv > ww) return  1;
	if (Exact(v) && Exact(w)) return 0;
	if (Exact(v)) return -1; /* 1 < 1E0 */
	if (Exact(w)) return  1; /* 1E0 > 1 */
	return 0;
}

double numhash(v) value v; {
	number *n= (number *)Ats(v);
	return .123*n->p + .777*n->q;
}

#define CONVBUFSIZ 100
char convbuf[CONVBUFSIZ];

string convnum(v) value v; {
	double x; string bp; bool prec_loss= No;
	Checknum(v);
	x= Numval(v);
 conv:	if (!prec_loss && Exact(v) && fabs(x) <= LONG &&
		fabs(Numerator(v)) < BIG && fabs(Denominator(v)) < BIG) {
		intlet len= 0 < Length(v) && Length(v) <= MAXNUMDIG ? Length(v) : 0;
		intlet dcnt, sigcnt; bool sig;
		if (Denominator(v) != One) {
			intlet k; double p= 1.0, q;
			prec_loss= Yes;
			for (k= 1; k < MAXNUMDIG; k++) {
				p*= 10.0;
				q= p/Denominator(v);
				if (k >= len && q == floor(q)) {
					prec_loss= No;
					break;
				}
			}
			len= k;
		}
	convex:	sprintf(convbuf, "%.*f", len, x);
		dcnt= sigcnt= 0; sig= No;
		for (bp= convbuf; *bp != '\0'; bp++) 
			if ('0' <= *bp && *bp <= '9') {
				dcnt++;
				if (*bp != '0') sig= Yes;
				if (sig) sigcnt++;
			}
		if (sigcnt < MINNUMDIG && prec_loss) goto conv;
		if (dcnt > MAXNUMDIG) {
			if (len <= 0) syserr("conversion error 1");
			if (Denominator(v) == One) len= 0;
			else len-= dcnt-MAXNUMDIG;
			if (len < 0) syserr("conversion error 2");
			goto convex;
		}
	} else { /*approx etc*/
		sprintf(convbuf, "%.*e", MAXNUMDIG-5, x);
		for (bp= convbuf; *bp != '\0'; bp++)
		if (*bp == 'e') {
			*bp= 'E';
			break;
		}
	}
	return convbuf;
}

value numconst(tx, q) txptr tx, q; {
	bool dig= No; double ex= 0, ap= 1; intlet ndap, len= 0;
	while (tx < q && '0' <= *tx && *tx <= '9') {
		dig= Yes;
		ex= 10*ex+(*tx++ - '0');
	}
	if (tx < q && *tx == '.') {
		tx++; ndap= 0;
		while (tx < q && '0' <= *tx && *tx <= '9') {
			dig= Yes; ndap++;
			len= *tx == '0' ? ndap : 0;
			ex= 10*ex+(*tx++ - '0'); ap*= 10;
		}
		if (!dig) syserr("numconst[1]");
	}
	if (tx < q && *tx == 'E') {
		intlet sign= 1; double expo= 0;
		tx++;
		if (!('0' <= *tx && *tx <= '9') && Keymark(*tx)) {
			tx--;
			goto exact;
		}
		if (!dig) ex= 1;
		if (tx < q && (*tx == '+' || *tx == '-'))
			if (*tx++ == '-') sign= -1;
		dig= No;
		while (tx < q && '0' <= *tx && *tx <= '9') {
			dig= Yes;
			expo= 10*expo+(*tx++ - '0');
		}
		if (!dig) syserr("numconst[2]");
		return mk_approx(ex/ap*exp(sign*expo*log(10.0)));
	}
exact:	return mk_exact(ex, ap, len);
}

printnum(f1, v) FILE *f1; value v; {
	FILE *f= f1 ? f1 : stdout;
	if (!Exact(v) || Denominator(v) == One) {
		if (!Exact(v))
			fputc('~', f);
		fputs(convnum(v), f);
	}
	else {
		value w = numerator(v);
		fputs(convnum(w), f);
		release(w);
		fputc('/', f);
		w = denominator(v);
		fputs(convnum(w), f);
		release(w);
	}
	if (!f1) fputc('\n', f); /* Flush buffer for sdb */
}