Newsgroups: comp.sources.unix
From: dbell@canb.auug.org.au (David I. Bell)
Subject: v27i135: calc-2.9.0 - arbitrary precision C-like programmable calculator, Part08/19
References: <1.755316719.21314@gw.home.vix.com>
Sender: unix-sources-moderator@gw.home.vix.com
Approved: vixie@gw.home.vix.com

Submitted-By: dbell@canb.auug.org.au (David I. Bell)
Posting-Number: Volume 27, Issue 135
Archive-Name: calc-2.9.0/part08

#!/bin/sh
# this is part 8 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc2.9.0/opcodes.c continued
#
CurArch=8
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
     exit 1; fi
( read Scheck
  if test "$Scheck" != $CurArch
  then echo "Please unpack part $Scheck next!"
       exit 1;
  else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file calc2.9.0/opcodes.c"
sed 's/^X//' << 'SHAR_EOF' >> calc2.9.0/opcodes.c
X#endif
X	vsprintf(buf, fmt, ap);
X	va_end(ap);
X	fprintf(stderr, "%s\n", buf);
X	funcname = NULL;
X	longjmp(jmpbuf, 1);
X}
X
X/* END CODE */
SHAR_EOF
echo "File calc2.9.0/opcodes.c is complete"
chmod 0644 calc2.9.0/opcodes.c || echo "restore of calc2.9.0/opcodes.c fails"
set `wc -c calc2.9.0/opcodes.c`;Sum=$1
if test "$Sum" != "52101"
then echo original size 52101, current size $Sum;fi
echo "x - extracting calc2.9.0/opcodes.h (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/opcodes.h &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X */
X
X#ifndef	OPCODES_H
X#define	OPCODES_H
X
X
X/*
X * Opcodes
X */
X#define OP_NOP		0L	/* no operation */
X#define OP_LOCALADDR	1L	/* load address of local variable */
X#define OP_GLOBALADDR	2L	/* load address of global variable */
X#define OP_PARAMADDR	3L	/* load address of paramater variable */
X#define OP_LOCALVALUE	4L	/* load value of local variable */
X#define OP_GLOBALVALUE	5L	/* load value of global variable */
X#define OP_PARAMVALUE	6L	/* load value of paramater variable */
X#define OP_NUMBER	7L	/* load constant real numeric value */
X#define OP_INDEXADDR	8L	/* load array index address */
X#define	OP_PRINTRESULT	9L	/* print result of top-level expression */
X#define OP_ASSIGN	10L	/* assign value to variable */
X#define OP_ADD		11L	/* add top two values */
X#define OP_SUB		12L	/* subtract top two values */
X#define OP_MUL		13L	/* multiply top two values */
X#define OP_DIV		14L	/* divide top two values */
X#define OP_MOD		15L	/* take mod of top two values */
X#define OP_SAVE		16L	/* save value for later use */
X#define OP_NEGATE	17L	/* negate top value */
X#define OP_INVERT	18L	/* invert top value */
X#define OP_INT		19L	/* take integer part of top value */
X#define OP_FRAC		20L	/* take fraction part of top value */
X#define OP_NUMERATOR	21L	/* take numerator of top value */
X#define OP_DENOMINATOR	22L	/* take denominator of top value */
X#define OP_DUPLICATE	23L	/* duplicate top value on stack */
X#define OP_POP		24L	/* pop top value from stack */
X#define OP_RETURN	25L	/* return value of function */
X#define OP_JUMPEQ	26L	/* jump if top value is zero */
X#define OP_JUMPNE	27L	/* jump if top value is nonzero */
X#define OP_JUMP		28L	/* jump unconditionally */
X#define OP_USERCALL	29L	/* call a user-defined function */
X#define OP_GETVALUE	30L	/* convert address to value */
X#define OP_EQ		31L	/* test top two elements for equality */
X#define OP_NE		32L	/* test top two elements for inequality */
X#define OP_LE		33L	/* test top two elements for <= */
X#define OP_GE		34L	/* test top two elements for >= */
X#define OP_LT		35L	/* test top two elements for < */
X#define OP_GT		36L	/* test top two elements for > */
X#define OP_PREINC	37L	/* add one to variable (++x) */
X#define OP_PREDEC	38L	/* subtract one from variable (--x) */
X#define OP_POSTINC	39L	/* add one to variable (x++) */
X#define OP_POSTDEC	40L	/* subtract one from variable (x--) */
X#define OP_DEBUG	41L	/* debugging point */
X#define OP_PRINT	42L	/* print value */
X#define OP_ASSIGNPOP	43L	/* assign to variable and remove it */
X#define OP_ZERO		44L	/* put zero on the stack */
X#define OP_ONE		45L	/* put one on the stack */
X#define OP_PRINTEOL	46L	/* print end of line */
X#define OP_PRINTSPACE	47L	/* print a space */
X#define OP_PRINTSTRING	48L	/* print constant string */
X#define OP_DUPVALUE	49L	/* duplicate value of top value */
X#define OP_OLDVALUE	50L	/* old calculation value */
X#define OP_QUO		51L	/* integer quotient of top two values */
X#define OP_POWER	52L	/* number raised to a power */
X#define OP_QUIT		53L	/* quit program */
X#define OP_CALL		54L	/* call built-in routine */
X#define OP_GETEPSILON	55L	/* get allowed error for calculations */
X#define OP_AND		56L	/* arithmetic and */
X#define OP_OR		57L	/* arithmetic or */
X#define OP_NOT		58L	/* logical not */
X#define OP_ABS		59L	/* absolute value */
X#define OP_SGN		60L	/* sign of number */
X#define OP_ISINT	61L	/* whether top value is integer */
X#define OP_CONDORJUMP	62L	/* conditional or jump */
X#define OP_CONDANDJUMP	63L	/* conditional and jump */
X#define OP_SQUARE	64L	/* square top value */
X#define OP_STRING	65L	/* load constant string value */
X#define OP_ISNUM	66L	/* whether top value is a number */
X#define OP_UNDEF	67L	/* load undefined value on stack */
X#define OP_ISNULL	68L	/* whether variable is the null value */
X#define OP_ARGVALUE	69L	/* load value of argument (parameter) n */
X#define OP_MATCREATE	70L	/* create matrix */
X#define OP_ISMAT	71L	/* whether variable is a matrix */
X#define OP_ISSTR	72L	/* whether variable is a string */
X#define OP_GETCONFIG	73L	/* get value of configuration parameter */
X#define OP_LEFTSHIFT	74L	/* left shift of integer */
X#define OP_RIGHTSHIFT	75L	/* right shift of integer */
X#define OP_CASEJUMP	76L	/* test case and jump if not matched */
X#define OP_ISODD	77L	/* whether value is an odd integer */
X#define OP_ISEVEN	78L	/* whether value is even integer */
X#define OP_FIADDR	79L	/* 'fast index' matrix value address */
X#define OP_FIVALUE	80L	/* 'fast index' matrix value */
X#define OP_ISREAL	81L	/* test value for real number */
X#define OP_IMAGINARY	82L	/* load imaginary numeric constant */
X#define OP_RE		83L	/* real part of complex number */
X#define OP_IM		84L	/* imaginary part of complex number */
X#define OP_CONJUGATE	85L	/* complex conjugate of complex number */
X#define OP_OBJCREATE	86L	/* create object */
X#define OP_ISOBJ	87L	/* whether value is an object */
X#define OP_NORM		88L	/* norm of value (square of abs) */
X#define OP_ELEMADDR	89L	/* address of element of object */
X#define OP_ELEMVALUE	90L	/* value of element of object */
X#define OP_ISTYPE	91L	/* whether two values are the same type */
X#define OP_SCALE	92L	/* scale value by a power of two */
X#define	OP_ISLIST	93L	/* whether value is a list */
X#define	OP_SWAP		94L	/* swap values of two variables */
X#define	OP_ISSIMPLE	95L	/* whether value is a simple type */
X#define	OP_CMP		96L	/* compare values returning -1, 0, or 1 */
X#define	OP_QUOMOD	97L	/* calculate quotient and remainder */
X#define	OP_SETCONFIG	98L	/* set configuration parameter */
X#define	OP_SETEPSILON	99L	/* set allowed error for calculations */
X#define	OP_ISFILE	100L	/* whether value is a file */
X#define	OP_ISASSOC	101L	/* whether value is an association */
X#define	OP_INITSTATIC	102L	/* once only code for static initialization */
X#define	OP_ELEMINIT	103L	/* assign element of matrix or object */
X#define MAX_OPCODE	103L	/* highest legal opcode */
X
X#endif
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/opcodes.h || echo "restore of calc2.9.0/opcodes.h fails"
set `wc -c calc2.9.0/opcodes.h`;Sum=$1
if test "$Sum" != "6052"
then echo original size 6052, current size $Sum;fi
echo "x - extracting calc2.9.0/qfunc.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/qfunc.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Extended precision rational arithmetic non-primitive functions
X */
X
X#include "qmath.h"
X
X
XNUMBER *_epsilon_;	/* default precision for real functions */
Xlong _epsilonprec_;	/* bits of precision for epsilon */
X
X
X/*
X * Set the default precision for real calculations.
X * The precision must be between zero and one.
X */
Xvoid
Xsetepsilon(q)
X	NUMBER *q;		/* number to be set as the new epsilon */
X{
X	NUMBER *old;
X
X	if (qisneg(q) || qiszero(q) || (qreli(q, 1L) >= 0))
X		math_error("Epsilon value must be between zero and one");
X	old = _epsilon_;
X	_epsilonprec_ = qprecision(q);
X	_epsilon_ = qlink(q);
X	if (old)
X		qfree(old);
X}
X
X
X/*
X * Return the inverse of one number modulo another.
X * That is, find x such that:
X *	Ax = 1 (mod B)
X * Returns zero if the numbers are not relatively prime (temporary hack).
X */
XNUMBER *
Xqminv(q1, q2)
X	NUMBER *q1, *q2;
X{
X	NUMBER *r;
X
X	if (qisfrac(q1) || qisfrac(q2))
X		math_error("Non-integers for minv");
X	r = qalloc();
X	if (zmodinv(q1->num, q2->num, &r->num)) {
X		qfree(r);
X		return qlink(&_qzero_);
X	}
X	return r;
X}
X
X
X/*
X * Return the modulo of one number raised to another.
X * Here q1 is the number to be raised, q2 is the power to raise it to,
X * and q3 is the number to take the modulo with the result.
X * The second and third numbers are assumed nonnegative.
X * Returned answer is non-negative.
X *	q4 = qpowermod(q1, q2, q3);
X */
XNUMBER *
Xqpowermod(q1, q2, q3)
X	NUMBER *q1, *q2, *q3;
X{
X	NUMBER *r;
X
X	if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
X		math_error("Non-integers for powermod");
X	r = qalloc();
X	zpowermod(q1->num, q2->num, q3->num, &r->num);
X	return r;
X}
X
X
X/*
X * Return the power function of one number with another.
X * The power must be integral.
X *	q3 = qpowi(q1, q2);
X */
XNUMBER *
Xqpowi(q1, q2)
X	NUMBER *q1, *q2;
X{
X	register NUMBER *r;
X	BOOL invert, sign;
X	ZVALUE num, den, z2;
X
X	if (qisfrac(q2))
X		math_error("Raising number to fractional power");
X	num = q1->num;
X	den = q1->den;
X	z2 = q2->num;
X	sign = (num.sign && zisodd(z2));
X	invert = z2.sign;
X	num.sign = 0;
X	z2.sign = 0;
X	/*
X	* Check for trivial cases first.
X	*/
X	if (ziszero(num)) {	/* zero raised to a power */
X		if (invert || ziszero(z2))
X			math_error("Zero raised to non-positive power");
X		return qlink(&_qzero_);
X	}
X	if (zisunit(num) && zisunit(den)) {	/* 1 or -1 raised to a power */
X		r = (sign ? q1 : &_qone_);
X		r->links++;
X		return r;
X	}
X	if (ziszero(z2))	/* raising to zeroth power */
X		return qlink(&_qone_);
X	if (zisunit(z2)) {	/* raising to power 1 or -1 */
X		if (invert)
X			return qinv(q1);
X		return qlink(q1);
X	}
X	/*
X	 * Not a trivial case.  Do the real work.
X	 */
X	r = qalloc();
X	if (!zisunit(num))
X		zpowi(num, z2, &r->num);
X	if (!zisunit(den))
X		zpowi(den, z2, &r->den);
X	if (invert) {
X		z2 = r->num;
X		r->num = r->den;
X		r->den = z2;
X	}
X	r->num.sign = sign;
X	return r;
X}
X
X
X/*
X * Given the legs of a right triangle, compute its hypothenuse within
X * the specified error.  This is sqrt(a^2 + b^2).
X */
XNUMBER *
Xqhypot(q1, q2, epsilon)
X	NUMBER *q1, *q2, *epsilon;
X{
X	NUMBER *tmp1, *tmp2, *tmp3;
X
X	if (qisneg(epsilon) || qiszero(epsilon))
X		math_error("Bad epsilon value for hypot");
X	if (qiszero(q1))
X		return qabs(q2);
X	if (qiszero(q2))
X		return qabs(q1);
X	tmp1 = qsquare(q1);
X	tmp2 = qsquare(q2);
X	tmp3 = qadd(tmp1, tmp2);
X	qfree(tmp1);
X	qfree(tmp2);
X	tmp1 = qsqrt(tmp3, epsilon);
X	qfree(tmp3);
X	return tmp1;
X}
X
X
X/*
X * Given one leg of a right triangle with unit hypothenuse, calculate
X * the other leg within the specified error.  This is sqrt(1 - a^2).
X * If the wantneg flag is nonzero, then negative square root is returned.
X */
XNUMBER *
Xqlegtoleg(q, epsilon, wantneg)
X	NUMBER *q, *epsilon;
X	BOOL wantneg;
X{
X	NUMBER *qt, *res, qtmp;
X	ZVALUE num, ztmp1, ztmp2;
X
X	if (qisneg(epsilon) || qiszero(epsilon))
X		math_error("Bad epsilon value for legtoleg");
X	if (qisunit(q))
X		return qlink(&_qzero_);
X	if (qiszero(q)) {
X		if (wantneg)
X			return qlink(&_qnegone_);
X		return qlink(&_qone_);
X	}
X	num = q->num;
X	num.sign = 0;
X	if (zrel(num, q->den) >= 0)
X		math_error("Leg too large in legtoleg");
X	zsquare(q->den, &ztmp1);
X	zsquare(num, &ztmp2);
X	zsub(ztmp1, ztmp2, &qtmp.num);
X	zfree(ztmp1);
X	zfree(ztmp2);
X	qtmp.den = _one_;
X	qt = qsqrt(&qtmp, epsilon);
X	zfree(qtmp.num);
X	qtmp.num = q->den;
X	res = qdiv(qt, &qtmp);
X	qfree(qt);
X	qt = qbappr(res, epsilon);
X	qfree(res);
X	if (wantneg) {
X		res = qneg(qt);
X		qfree(qt);
X		qt = res;
X	}
X	return qt;
X}
X
X
X/*
X * Compute the square root of a number to within the specified error.
X * If the number is an exact square, the exact result is returned.
X *	q3 = qsqrt(q1, q2);
X */
XNUMBER *
Xqsqrt(q1, epsilon)
X	register NUMBER *q1, *epsilon;
X{
X	long bits, bits2;	/* number of bits of precision */
X	int exact;
X	NUMBER *r;
X	ZVALUE t1, t2;
X
X	if (qisneg(q1))
X		math_error("Square root of negative number");
X	if (qisneg(epsilon) || qiszero(epsilon))
X		math_error("Bad epsilon value for sqrt");
X	if (qiszero(q1))
X		return qlink(&_qzero_);
X	if (qisunit(q1))
X		return qlink(&_qone_);
X	if (qiszero(epsilon) && qisint(q1) && zistiny(q1->num) && (*q1->num.v <= 3))
X		return qlink(&_qone_);
X	bits = zhighbit(epsilon->den) - zhighbit(epsilon->num) + 1;
X	if (bits < 0)
X		bits = 0;
X	bits2 = zhighbit(q1->num) - zhighbit(q1->den) + 1;
X	if (bits2 > 0)
X		bits += bits2;
X	r = qalloc();
X	zshift(q1->num, bits * 2, &t2);
X	zmul(q1->den, t2, &t1);
X	zfree(t2);
X	exact = zsqrt(t1, &t2);
X	zfree(t1);
X	if (exact) {
X		zshift(q1->den, bits, &t1);
X		zreduce(t2, t1, &r->num, &r->den);
X	} else {
X		zquo(t2, q1->den, &t1);
X		zfree(t2);
X		zbitvalue(bits, &t2);
X		zreduce(t1, t2, &r->num, &r->den);
X	}
X	zfree(t1);
X	zfree(t2);
X	if (qiszero(r)) {
X		qfree(r);
X		r = qlink(&_qzero_);
X	}
X	return r;
X}
X
X
X/*
X * Calculate the integral part of the square root of a number.
X * Example:  qisqrt(13) = 3.
X */
XNUMBER *
Xqisqrt(q)
X	NUMBER *q;
X{
X	NUMBER *r;
X	ZVALUE tmp;
X
X	if (qisneg(q))
X		math_error("Square root of negative number");
X	if (qiszero(q))
X		return qlink(&_qzero_);
X	if (qisint(q) && zistiny(q->num) && (z1tol(q->num) < 4))
X		return qlink(&_qone_);
X	r = qalloc();
X	if (qisint(q)) {
X		(void) zsqrt(q->num, &r->num);
X		return r;
X	}
X	zquo(q->num, q->den, &tmp);
X	(void) zsqrt(tmp, &r->num);
X	zfree(tmp);
X	return r;
X}
X
X
X/*
X * Return whether or not a number is an exact square.
X */
XBOOL
Xqissquare(q)
X	NUMBER *q;
X{
X	BOOL flag;
X
X	flag = zissquare(q->num);
X	if (qisint(q) || !flag)
X		return flag;
X	return zissquare(q->den);
X}
X
X
X/*
X * Compute the greatest integer of the Kth root of a number.
X * Example:  qiroot(85, 3) = 4.
X */
XNUMBER *
Xqiroot(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	NUMBER *r;
X	ZVALUE tmp;
X
X	if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
X		math_error("Taking number to bad root value");
X	if (qiszero(q1))
X		return qlink(&_qzero_);
X	if (qisone(q1) || qisone(q2))
X		return qlink(q1);
X	if (qistwo(q2))
X		return qisqrt(q1);
X	r = qalloc();
X	if (qisint(q1)) {
X		zroot(q1->num, q2->num, &r->num);
X		return r;
X	}
X	zquo(q1->num, q1->den, &tmp);
X	zroot(tmp, q2->num, &r->num);
X	zfree(tmp);
X	return r;
X}
X
X
X/*
X * Return the greatest integer of the base 2 log of a number.
X * This is the number such that  1 <= x / log2(x) < 2.
X * Examples:  qilog2(8) = 3, qilog2(1.3) = 1, qilog2(1/7) = -3.
X */
Xlong
Xqilog2(q)
X	NUMBER *q;		/* number to take log of */
X{
X	long n;			/* power of two */
X	int c;			/* result of comparison */
X	ZVALUE tmp;		/* temporary value */
X
X	if (qisneg(q) || qiszero(q))
X		math_error("Non-positive number for log2");
X	if (qisint(q))
X		return zhighbit(q->num);
X	n = zhighbit(q->num) - zhighbit(q->den);
X	if (n == 0)
X		c = zrel(q->num, q->den);
X	else if (n > 0) {
X		zshift(q->den, n, &tmp);
X		c = zrel(q->num, tmp);
X	} else {
X		zshift(q->num, n, &tmp);
X		c = zrel(tmp, q->den);
X	}
X	if (n)
X		zfree(tmp);
X	if (c < 0)
X		n--;
X	return n;
X}
X
X
X/*
X * Return the greatest integer of the base 10 log of a number.
X * This is the number such that  1 <= x / log10(x) < 10.
X * Examples:  qilog10(100) = 2, qilog10(12.3) = 1, qilog10(.023) = -2.
X */
Xlong
Xqilog10(q)
X	NUMBER *q;		/* number to take log of */
X{
X	long n;			/* log value */
X	ZVALUE temp;		/* temporary value */
X
X	if (qisneg(q) || qiszero(q))
X		math_error("Non-positive number for log10");
X	if (qisint(q))
X		return zlog10(q->num);
X	/*
X	 * The number is not an integer.
X	 * Compute the result if the number is greater than one.
X	 */
X	if ((q->num.len > q->den.len) ||
X		((q->num.len == q->den.len) && (zrel(q->num, q->den) > 0))) {
X			zquo(q->num, q->den, &temp);
X			n = zlog10(temp);
X			zfree(temp);
X			return n;
X	}
X	/*
X	 * Here if the number is less than one.
X	 * If the number is the inverse of a power of ten, then the obvious answer
X	 * will be off by one.  Subtracting one if the number is the inverse of an
X	 * integer will fix it.
X	 */
X	if (zisunit(q->num))
X		zsub(q->den, _one_, &temp);
X	else
X		zquo(q->den, q->num, &temp);
X	n = -zlog10(temp) - 1;
X	zfree(temp);
X	return n;
X}
X
X
X/*
X * Return the number of digits in a number, ignoring the sign.
X * For fractions, this is the number of digits of its greatest integer.
X * Examples: qdigits(3456) = 4, qdigits(-23.45) = 2, qdigits(.0120) = 1.
X */
Xlong
Xqdigits(q)
X	NUMBER *q;		/* number to count digits of */
X{
X	long n;			/* number of digits */
X	ZVALUE temp;		/* temporary value */
X
X	if (qisint(q))
X		return zdigits(q->num);
X	zquo(q->num, q->den, &temp);
X	n = zdigits(temp);
X	zfree(temp);
X	return n;
X}
X
X
X/*
X * Return the digit at the specified decimal place of a number represented
X * in floating point.  The lowest digit of the integral part of a number
X * is the zeroth decimal place.  Negative decimal places indicate digits
X * to the right of the decimal point.  Examples: qdigit(1234.5678, 1) = 3,
X * qdigit(1234.5678, -3) = 7.
X */
XFLAG
Xqdigit(q, n)
X	NUMBER *q;
X	long n;
X{
X	ZVALUE tenpow, tmp1, tmp2;
X	FLAG res;
X
X	/*
X	 * Zero number or negative decimal place of integer is trivial.
X	 */
X	if (qiszero(q) || (qisint(q) && (n < 0)))
X		return 0;
X	/*
X	 * For non-negative decimal places, answer is easy.
X	 */
X	if (n >= 0) {
X		if (qisint(q))
X			return zdigit(q->num, n);
X		zquo(q->num, q->den, &tmp1);
X		res = zdigit(tmp1, n);
X		zfree(tmp1);
X		return res;
X	}
X	/*
X	 * Fractional value and want negative digit, must work harder.
X	 */
X	ztenpow(-n, &tenpow);
X	zmul(q->num, tenpow, &tmp1);
X	zfree(tenpow);
X	zquo(tmp1, q->den, &tmp2);
X	res = zmodi(tmp2, 10L);
X	zfree(tmp1);
X	zfree(tmp2);
X	return res;
X}
X
X
X/*
X * Return whether or not a bit is set at the specified bit position in a
X * number.  The lowest bit of the integral part of a number is the zeroth
X * bit position.  Negative bit positions indicate bits to the right of the
X * binary decimal point.  Examples: qdigit(17.1, 0) = 1, qdigit(17.1, -1) = 0.
X */
XBOOL
Xqisset(q, n)
X	NUMBER *q;
X	long n;
X{
X	NUMBER *qtmp1, *qtmp2;
X	ZVALUE ztmp;
X	BOOL res;
X
X	/*
X	 * Zero number or negative bit position place of integer is trivial.
X	 */
X	if (qiszero(q) || (qisint(q) && (n < 0)))
X		return FALSE;
X	/*
X	 * For non-negative bit positions, answer is easy.
X	 */
X	if (n >= 0) {
X		if (qisint(q))
X			return zisset(q->num, n);
X		zquo(q->num, q->den, &ztmp);
X		res = zisset(ztmp, n);
X		zfree(ztmp);
X		return res;
X	}
X	/*
X	 * Fractional value and want negative bit position, must work harder.
X	 */
X	qtmp1 = qscale(q, -n);
X	qtmp2 = qint(qtmp1);
X	qfree(qtmp1);
X	res = ((qtmp2->num.v[0] & 0x01) != 0);
X	qfree(qtmp2);
X	return res;
X}
X
X
X/*
X * Compute the factorial of an integer.
X *	q2 = qfact(q1);
X */
XNUMBER *
Xqfact(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (qisfrac(q))
X		math_error("Non-integral factorial");
X	if (qiszero(q) || zisone(q->num))
X		return qlink(&_qone_);
X	r = qalloc();
X	zfact(q->num, &r->num);
X	return r;
X}
X
X
X/*
X * Compute the product of the primes less than or equal to a number.
X *	q2 = qpfact(q1);
X */
XNUMBER *
Xqpfact(q)
X	register NUMBER *q;
X{
X	NUMBER *r;
X
X	if (qisfrac(q))
X		math_error("Non-integral factorial");
X	r = qalloc();
X	zpfact(q->num, &r->num);
X	return r;
X}
X
X
X/*
X * Compute the lcm of all the numbers less than or equal to a number.
X *	q2 = qlcmfact(q1);
X */
XNUMBER *
Xqlcmfact(q)
X	register NUMBER *q;
X{
X	NUMBER *r;
X
X	if (qisfrac(q))
X		math_error("Non-integral lcmfact");
X	r = qalloc();
X	zlcmfact(q->num, &r->num);
X	return r;
X}
X
X
X/*
X * Compute the permutation function  M! / (M - N)!.
X */
XNUMBER *
Xqperm(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	register NUMBER *r;
X
X	if (qisfrac(q1) || qisfrac(q2))
X		math_error("Non-integral arguments for permutation");
X	r = qalloc();
X	zperm(q1->num, q2->num, &r->num);
X	return r;
X}
X
X
X/*
X * Compute the combinatorial function  M! / (N! * (M - N)!).
X */
XNUMBER *
Xqcomb(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	register NUMBER *r;
X
X	if (qisfrac(q1) || qisfrac(q2))
X		math_error("Non-integral arguments for combinatorial");
X	r = qalloc();
X	zcomb(q1->num, q2->num, &r->num);
X	return r;
X}
X
X
X/*
X * Compute the Jacobi function (a / b).
X * -1 => a is not quadratic residue mod b
X * 1 => b is composite, or a is quad residue of b
X * 0 => b is even or b < 0
X */
XNUMBER *
Xqjacobi(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	if (qisfrac(q1) || qisfrac(q2))
X		math_error("Non-integral arguments for jacobi");
X	return itoq((long) zjacobi(q1->num, q2->num));
X}
X
X
X/*
X * Compute the Fibonacci number F(n).
X */
XNUMBER *
Xqfib(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (qisfrac(q))
X		math_error("Non-integral Fibonacci number");
X	r = qalloc();
X	zfib(q->num, &r->num);
X	return r;
X}
X
X
X/*
X * Truncate a number to the specified number of decimal places.
X * Specifying zero places makes the result identical to qint.
X * Example: qtrunc(2/3, 3) = .666
X */
XNUMBER *
Xqtrunc(q1, q2)
X	NUMBER *q1, *q2;
X{
X	long places;
X	NUMBER *r;
X	ZVALUE tenpow, tmp1, tmp2;
X
X	if (qisfrac(q2) || qisneg(q2) || !zistiny(q2->num))
X		math_error("Bad number of places for qtrunc");
X	if (qisint(q1))
X		return qlink(q1);
X	r = qalloc();
X	places = z1tol(q2->num);
X	/*
X	 * Ok, produce the result.
X	 * First see if we want no places, in which case just take integer part.
X	 */
X	if (places == 0) {
X		zquo(q1->num, q1->den, &r->num);
X		return r;
X	}
X	ztenpow(places, &tenpow);
X	zmul(q1->num, tenpow, &tmp1);
X	zquo(tmp1, q1->den, &tmp2);
X	zfree(tmp1);
X	if (ziszero(tmp2)) {
X		zfree(tmp2);
X		return qlink(&_qzero_);
X	}
X	/*
X	 * Now reduce the result to the lowest common denominator.
X	 */
X	zgcd(tmp2, tenpow, &tmp1);
X	if (zisunit(tmp1)) {
X		zfree(tmp1);
X		r->num = tmp2;
X		r->den = tenpow;
X		return r;
X	}
X	zquo(tmp2, tmp1, &r->num);
X	zquo(tenpow, tmp1, &r->den);
X	zfree(tmp1);
X	zfree(tmp2);
X	zfree(tenpow);
X	return r;
X}
X
X
X/*
X * Round a number to the specified number of decimal places.
X * Zero decimal places means round to the nearest integer.
X * Example: qround(2/3, 3) = .667
X */
XNUMBER *
Xqround(q, places)
X	NUMBER *q;		/* number to be rounded */
X	long places;		/* number of decimal places to round to */
X{
X	NUMBER *r;
X	ZVALUE tenpow, roundval, tmp1, tmp2;
X
X	if (places < 0)
X		math_error("Negative places for qround");
X	if (qisint(q))
X		return qlink(q);
X	/*
X	 * Calculate one half of the denominator, ignoring fractional results.
X	 * This is the value we will add in order to cause rounding.
X	 */
X	zshift(q->den, -1L, &roundval);
X	roundval.sign = q->num.sign;
X	/*
X	 * Ok, now do the actual work to produce the result.
X	 */
X	r = qalloc();
X	ztenpow(places, &tenpow);
X	zmul(q->num, tenpow, &tmp2);
X	zadd(tmp2, roundval, &tmp1);
X	zfree(tmp2);
X	zfree(roundval);
X	zquo(tmp1, q->den, &tmp2);
X	zfree(tmp1);
X	if (ziszero(tmp2)) {
X		zfree(tmp2);
X		return qlink(&_qzero_);
X	}
X	/*
X	 * Now reduce the result to the lowest common denominator.
X	 */
X	zgcd(tmp2, tenpow, &tmp1);
X	if (zisunit(tmp1)) {
X		zfree(tmp1);
X		r->num = tmp2;
X		r->den = tenpow;
X		return r;
X	}
X	zquo(tmp2, tmp1, &r->num);
X	zquo(tenpow, tmp1, &r->den);
X	zfree(tmp1);
X	zfree(tmp2);
X	zfree(tenpow);
X	return r;
X}
X
X
X/*
X * Truncate a number to the specified number of binary places.
X * Specifying zero places makes the result identical to qint.
X */
XNUMBER *
Xqbtrunc(q1, q2)
X	NUMBER *q1, *q2;
X{
X	long places, twopow;
X	NUMBER *r;
X	ZVALUE tmp1, tmp2;
X
X	if (qisfrac(q2) || qisneg(q2) || !zistiny(q2->num))
X		math_error("Bad number of places for qtrunc");
X	if (qisint(q1))
X		return qlink(q1);
X	r = qalloc();
X	places = z1tol(q2->num);
X	/*
X	 * Ok, produce the result.
X	 * First see if we want no places, in which case just take integer part.
X	 */
X	if (places == 0) {
X		zquo(q1->num, q1->den, &r->num);
X		return r;
X	}
X	zshift(q1->num, places, &tmp1);
X	zquo(tmp1, q1->den, &tmp2);
X	zfree(tmp1);
X	if (ziszero(tmp2)) {
X		zfree(tmp2);
X		return qlink(&_qzero_);
X	}
X	/*
X	 * Now reduce the result to the lowest common denominator.
X	 */
X	if (zisodd(tmp2)) {
X		r->num = tmp2;
X		zbitvalue(places, &r->den);
X		return r;
X	}
X	twopow = zlowbit(tmp2);
X	if (twopow > places)
X		twopow = places;
X	places -= twopow;
X	zshift(tmp2, -twopow, &r->num);
X	zfree(tmp2);
X	zbitvalue(places, &r->den);
X	return r;
X}
X
X
X/*
X * Round a number to the specified number of binary places.
X * Zero binary places means round to the nearest integer.
X */
XNUMBER *
Xqbround(q, places)
X	NUMBER *q;		/* number to be rounded */
X	long places;		/* number of binary places to round to */
X{
X	long twopow;
X	NUMBER *r;
X	ZVALUE roundval, tmp1, tmp2;
X
X	if (places < 0)
X		math_error("Negative places for qbround");
X	if (qisint(q))
X		return qlink(q);
X	r = qalloc();
X	/*
X	 * Calculate one half of the denominator, ignoring fractional results.
X	 * This is the value we will add in order to cause rounding.
X	 */
X	zshift(q->den, -1L, &roundval);
X	roundval.sign = q->num.sign;
X	/*
X	 * Ok, produce the result.
X	 */
X	zshift(q->num, places, &tmp1);
X	zadd(tmp1, roundval, &tmp2);
X	zfree(roundval);
X	zfree(tmp1);
X	zquo(tmp2, q->den, &tmp1);
X	zfree(tmp2);
X	if (ziszero(tmp1)) {
X		zfree(tmp1);
X		return qlink(&_qzero_);
X	}
X	/*
X	 * Now reduce the result to the lowest common denominator.
X	 */
X	if (zisodd(tmp1)) {
X		r->num = tmp1;
X		zbitvalue(places, &r->den);
X		return r;
X	}
X	twopow = zlowbit(tmp1);
X	if (twopow > places)
X		twopow = places;
X	places -= twopow;
X	zshift(tmp1, -twopow, &r->num);
X	zfree(tmp1);
X	zbitvalue(places, &r->den);
X	return r;
X}
X
X
X/*
X * Approximate a number by using binary rounding with the minimum number
X * of binary places so that the resulting number is within the specified
X * epsilon of the original number.
X */
XNUMBER *
Xqbappr(q, e)
X	NUMBER *q, *e;
X{
X	long bits;
X
X	if (qisneg(e) || qiszero(e))
X		math_error("Bad epsilon value for qbappr");
X	if (e == _epsilon_)
X		bits = _epsilonprec_ + 1;
X	else
X		bits = qprecision(e) + 1;
X	return qbround(q, bits);
X}
X
X
X/*
X * Approximate a number using continued fractions until the approximation
X * error is less than the specified value.  If a NULL pointer is given
X * for the error value, then the closest simpler fraction is returned.
X */
XNUMBER *
Xqcfappr(q, e)
X	NUMBER *q, *e;
X{
X	NUMBER qtest, *qtmp;
X	ZVALUE u1, u2, u3, v1, v2, v3, t1, t2, t3, qq, tt;
X	int i;
X	BOOL haveeps;
X
X	haveeps = TRUE;
X	if (e == NULL) {
X		haveeps = FALSE;
X		e = &_qzero_;
X	}
X	if (qisneg(e))
X		math_error("Negative epsilon for cfappr");
X	if (qisint(q) || zisunit(q->num) || (haveeps && qiszero(e)))
X		return qlink(q);
X	u1 = _one_;
X	u2 = _zero_;
X	u3 = q->num;
X	u3.sign = 0;
X	v1 = _zero_;
X	v2 = _one_;
X	v3 = q->den;
X	while (!ziszero(v3)) {
X		if (!ziszero(u2) && !ziszero(u1)) {
X			qtest.num = u2;
X			qtest.den = u1;
X			qtest.den.sign = 0;
X			qtest.num.sign = q->num.sign;
X			qtmp = qsub(q, &qtest);
X			qtest = *qtmp;
X			qtest.num.sign = 0;
X			i = qrel(&qtest, e);
X			qfree(qtmp);
X			if (i <= 0)
X				break;
X		}
X		zquo(u3, v3, &qq);
X		zmul(qq, v1, &tt); zsub(u1, tt, &t1); zfree(tt);
X		zmul(qq, v2, &tt); zsub(u2, tt, &t2); zfree(tt);
X		zmul(qq, v3, &tt); zsub(u3, tt, &t3); zfree(tt);
X		zfree(qq); zfree(u1); zfree(u2);
X		if ((u3.v != q->num.v) && (u3.v != q->den.v))
X			zfree(u3);
X		u1 = v1; u2 = v2; u3 = v3;
X		v1 = t1; v2 = t2; v3 = t3;
X	}
X	if (u3.v != q->den.v)
X		zfree(u3);
X	zfree(v1);
X	zfree(v2);
X	i = ziszero(v3);
X	zfree(v3);
X	if (i && haveeps) {
X		zfree(u1);
X		zfree(u2);
X		return qlink(q);
X	}
X	qtest.num = u2;
X	qtest.den = u1;
X	qtest.den.sign = 0;
X	qtest.num.sign = q->num.sign;
X	qtmp = qcopy(&qtest);
X	zfree(u1);
X	zfree(u2);
X	return qtmp;
X}
X
X
X/*
X * Return an indication on whether or not two fractions are approximately
X * equal within the specified epsilon. Returns negative if the absolute value
X * of the difference between the two numbers is less than epsilon, zero if
X * the difference is equal to epsilon, and positive if the difference is
X * greater than epsilon.
X */
XFLAG
Xqnear(q1, q2, epsilon)
X	NUMBER *q1, *q2, *epsilon;
X{
X	int res;
X	NUMBER qtemp, *qq;
X
X	if (qisneg(epsilon))
X		math_error("Negative epsilon for near");
X	if (q1 == q2) {
X		if (qiszero(epsilon))
X			return 0;
X		return -1;
X	}
X	if (qiszero(epsilon))
X		return qcmp(q1, q2);
X	if (qiszero(q2)) {
X		qtemp = *q1;
X		qtemp.num.sign = 0;
X		return qrel(&qtemp, epsilon);
X	}
X	if (qiszero(q1)) {
X		qtemp = *q2;
X		qtemp.num.sign = 0;
X		return qrel(&qtemp, epsilon);
X	}
X	qq = qsub(q1, q2);
X	qtemp = *qq;
X	qtemp.num.sign = 0;
X	res = qrel(&qtemp, epsilon);
X	qfree(qq);
X	return res;
X}
X
X
X/*
X * Compute the gcd (greatest common divisor) of two numbers.
X *	q3 = qgcd(q1, q2);
X */
XNUMBER *
Xqgcd(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	ZVALUE z;
X	NUMBER *q;
X
X	if (q1 == q2)
X		return qabs(q1);
X	if (qisfrac(q1) || qisfrac(q2)) {
X		q = qalloc();
X		zgcd(q1->num, q2->num, &q->num);
X		zlcm(q1->den, q2->den, &q->den);
X		return q;
X	}
X	if (qiszero(q1))
X		return qabs(q2);
X	if (qiszero(q2))
X		return qabs(q1);
X	if (qisunit(q1) || qisunit(q2))
X		return qlink(&_qone_);
X	zgcd(q1->num, q2->num, &z);
X	if (zisunit(z)) {
X		zfree(z);
X		return qlink(&_qone_);
X	}
X	q = qalloc();
X	q->num = z;
X	return q;
X}
X
X
X/*
X * Compute the lcm (least common multiple) of two numbers.
X *	q3 = qlcm(q1, q2);
X */
XNUMBER *
Xqlcm(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	NUMBER *q;
X
X	if (qiszero(q1) || qiszero(q2))
X		return qlink(&_qzero_);
X	if (q1 == q2)
X		return qabs(q1);
X	if (qisunit(q1))
X		return qabs(q2);
X	if (qisunit(q2))
X		return qabs(q1);
X	q = qalloc();
X	zlcm(q1->num, q2->num, &q->num);
X	if (qisfrac(q1) || qisfrac(q2))
X		zgcd(q1->den, q2->den, &q->den);
X	return q;
X}
X
X
X/*
X * Remove all occurances of the specified factor from a number.
X * Returned number is always positive.
X */
XNUMBER *
Xqfacrem(q1, q2)
X	NUMBER *q1, *q2;
X{
X	long count;
X	ZVALUE tmp;
X	NUMBER *r;
X
X	if (qisfrac(q1) || qisfrac(q2))
X		math_error("Non-integers for factor removal");
X	count = zfacrem(q1->num, q2->num, &tmp);
X	if (zisunit(tmp)) {
X		zfree(tmp);
X		return qlink(&_qone_);
X	}
X	if (count == 0) {
X		zfree(tmp);
X		return qlink(q1);
X	}
X	r = qalloc();
X	r->num = tmp;
X	return r;
X}
X
X
X/*
X * Divide one number by the gcd of it with another number repeatedly until
X * the number is relatively prime.
X */
XNUMBER *
Xqgcdrem(q1, q2)
X	NUMBER *q1, *q2;
X{
X	ZVALUE tmp;
X	NUMBER *r;
X
X	if (qisfrac(q1) || qisfrac(q2))
X		math_error("Non-integers for gcdrem");
X	zgcdrem(q1->num, q2->num, &tmp);
X	if (zisunit(tmp)) {
X		zfree(tmp);
X		return qlink(&_qone_);
X	}
X	if (zcmp(q1->num, tmp) == 0) {
X		zfree(tmp);
X		return qlink(q1);
X	}
X	r = qalloc();
X	r->num = tmp;
X	return r;
X}
X
X
X/*
X * Return the lowest prime factor of a number.
X * Search is conducted for the specified number of primes.
X * Returns one if no factor was found.
X */
XNUMBER *
Xqlowfactor(q1, q2)
X	NUMBER *q1, *q2;
X{
X	if (qisfrac(q1) || qisfrac(q2))
X		math_error("Non-integers for lowfactor");
X	return itoq(zlowfactor(q1->num, ztoi(q2->num)));
X}
X
X
X/*
X * Return the number of places after the decimal point needed to exactly
X * represent the specified number as a real number.  Integers return zero,
X * and non-terminating decimals return minus one.  Examples:
X *	qplaces(1/7)=-1, qplaces(3/10)= 1, qplaces(1/8)=3, qplaces(4)=0.
X */
Xlong
Xqplaces(q)
X	NUMBER *q;
X{
X	long twopow, fivepow;
X	HALF fiveval[2];
X	ZVALUE five;
X	ZVALUE tmp;
X
X	if (qisint(q))	/* no decimal places if number is integer */
X		return 0;
X	/*
X	 * The number of decimal places of a fraction in lowest terms is finite
X	 * if an only if the denominator is of the form 2^A * 5^B, and then the
X	 * number of decimal places is equal to MAX(A, B).
X	 */
X	five.sign = 0;
X	five.len = 1;
X	five.v = fiveval;
X	fiveval[0] = 5;
X	fivepow = zfacrem(q->den, five, &tmp);
X	if (!zisonebit(tmp)) {
X		zfree(tmp);
X		return -1;
X	}
X	twopow = zlowbit(tmp);
X	zfree(tmp);
X	if (twopow < fivepow)
X		twopow = fivepow;
X	return twopow;
X}
X
X
X/*
X * Perform a probabilistic primality test (algorithm P in Knuth).
X * Returns FALSE if definitely not prime, or TRUE if probably prime.
X * Second arg determines how many times to check for primality.
X */
XBOOL
Xqprimetest(q1, q2)
X	NUMBER *q1, *q2;
X{
X	if (qisfrac(q1) || qisfrac(q2) || qisneg(q2))
X		math_error("Bad arguments for qprimetest");
X	return zprimetest(q1->num, qtoi(q2));
X}
X
X
X/*
X * Return a trivial hash value for a number.
X */
XHASH
Xqhash(q)
X	NUMBER *q;
X{
X	HASH hash;
X
X	hash = zhash(q->num);
X	if (qisfrac(q))
X		hash += zhash(q->den) * 2000003;
X	return hash;
X}
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/qfunc.c || echo "restore of calc2.9.0/qfunc.c fails"
set `wc -c calc2.9.0/qfunc.c`;Sum=$1
if test "$Sum" != "24730"
then echo original size 24730, current size $Sum;fi
echo "x - extracting calc2.9.0/qio.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/qio.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Scanf and printf routines for arbitrary precision rational numbers
X */
X
X#include "stdarg.h"
X#include "qmath.h"
X
X
X#define	PUTCHAR(ch)		math_chr(ch)
X#define	PUTSTR(str)		math_str(str)
X#define	PRINTF1(fmt, a1)	math_fmt(fmt, a1)
X#define	PRINTF2(fmt, a1, a2)	math_fmt(fmt, a1, a2)
X
X
X#if 0
Xstatic long	etoalen;
Xstatic char	*etoabuf = NULL;
X#endif
X
Xstatic long	scalefactor;
Xstatic ZVALUE	scalenumber = { 0, 0, 0 };
X
X
X/*
X * Print a formatted string containing arbitrary numbers, similar to printf.
X * ALL numeric arguments to this routine are rational NUMBERs.
X * Various forms of printing such numbers are supplied, in addition
X * to strings and characters.  Output can actually be to any FILE
X * stream or a string.
X */
X#ifdef VARARGS
X# define VA_ALIST1 fmt, va_alist
X# define VA_DCL1 char *fmt; va_dcl
X#else
X# ifdef __STDC__
X#  define VA_ALIST1 char *fmt, ...
X#  define VA_DCL1
X# else
X#  define VA_ALIST1 fmt
X#  define VA_DCL1 char *fmt;
X# endif
X#endif
X/*VARARGS*/
Xvoid
Xqprintf(VA_ALIST1)
X	VA_DCL1
X{
X	va_list ap;
X	NUMBER *q;
X	int ch, sign;
X	long width, precision;
X
X#ifdef VARARGS
X	va_start(ap);
X#else
X	va_start(ap, fmt);
X#endif
X	while ((ch = *fmt++) != '\0') {
X		if (ch == '\\') {
X			ch = *fmt++;
X			switch (ch) {
X				case 'n': ch = '\n'; break;
X				case 'r': ch = '\r'; break;
X				case 't': ch = '\t'; break;
X				case 'f': ch = '\f'; break;
X				case 'v': ch = '\v'; break;
X				case 'b': ch = '\b'; break;
X				case 0:
X					va_end(ap);
X					return;
X			}
X			PUTCHAR(ch);
X			continue;
X		}
X		if (ch != '%') {
X			PUTCHAR(ch);
X			continue;
X		}
X		ch = *fmt++;
X		width = 0; precision = 8; sign = 1;
Xpercent:	;
X		switch (ch) {
X			case 'd':
X				q = va_arg(ap, NUMBER *);
X				qprintfd(q, width);
X				break;
X			case 'f':
X				q = va_arg(ap, NUMBER *);
X				qprintff(q, width, precision);
X				break;
X			case 'e':
X				q = va_arg(ap, NUMBER *);
X				qprintfe(q, width, precision);
X				break;
X			case 'r':
X			case 'R':
X				q = va_arg(ap, NUMBER *);
X				qprintfr(q, width, (BOOL) (ch == 'R'));
X				break;
X			case 'N':
X				q = va_arg(ap, NUMBER *);
X				zprintval(q->num, 0L, width);
X				break;
X			case 'D':
X				q = va_arg(ap, NUMBER *);
X				zprintval(q->den, 0L, width);
X				break;
X			case 'o':
X				q = va_arg(ap, NUMBER *);
X				qprintfo(q, width);
X				break;
X			case 'x':
X				q = va_arg(ap, NUMBER *);
X				qprintfx(q, width);
X				break;
X			case 'b':
X				q = va_arg(ap, NUMBER *);
X				qprintfb(q, width);
X				break;
X			case 's':
X				PUTSTR(va_arg(ap, char *));
X				break;
X			case 'c':
X				PUTCHAR(va_arg(ap, int));
X				break;
X			case 0:
X				va_end(ap);
X				return;
X			case '-':
X				sign = -1;
X				ch = *fmt++;
X			default:
X		if (('0' <= ch && ch <= '9') || ch == '.' || ch == '*') {
X			if (ch == '*') {
X				q = va_arg(ap, NUMBER *);
X				width = sign * qtoi(q);
X				ch = *fmt++;
X			} else if (ch != '.') {
X				width = ch - '0';
X				while ('0' <= (ch = *fmt++) && ch <= '9')
X					width = width * 10 + ch - '0';
X				width *= sign;
X			}
X			if (ch == '.') {
X				if ((ch = *fmt++) == '*') {
X					q = va_arg(ap, NUMBER *);
X					precision = qtoi(q);
X					ch = *fmt++;
X				} else {
X					precision = 0;
X					while ('0' <= (ch = *fmt++) && ch <= '9')
X						precision = precision * 10 + ch - '0';
X				}
X			}
X			goto percent;
X		}
X		}
X	}
X	va_end(ap);
X}
X
X
X#if 0
X/*
X * Read a number from the specified FILE stream (NULL means stdin).
X * The number can be an integer, a fraction, a real number, an
X * exponential number, or a hex, octal or binary number.  Leading blanks
X * are skipped.  Illegal numbers return NULL.  Unrecognized characters
X * remain to be read on the line.
X *	q = qreadval(fp);
X */
XNUMBER *
Xqreadval(fp)
X	FILE *fp;		/* file stream to read from (or NULL) */
X{
X	NUMBER *r;		/* returned number */
X	char *cp; 		/* current buffer location */
X	long savecc;		/* characters saved in buffer */
X	long scancc; 		/* characters parsed correctly */
X	int ch;			/* current character */
X
X	if (fp == NULL)
X		fp = stdin;
X	if (etoabuf == NULL) {
X		etoabuf = (char *)malloc(OUTBUFSIZE + 2);
X		if (etoabuf == NULL)
X			return NULL;
X		etoalen = OUTBUFSIZE;
X	}
X	cp = etoabuf;
X	ch = fgetc(fp);
X	while ((ch == ' ') || (ch == '\t'))
X		ch = fgetc(fp);
X	savecc = 0;
X	for (;;) {
X		if (ch == EOF)
X			return NULL;
X		if (savecc >= etoalen)
X		{
X			cp = (char *)realloc(etoabuf, etoalen + OUTBUFSIZE + 2);
X			if (cp == NULL)
X				return NULL;
X			etoabuf = cp;
X			etoalen += OUTBUFSIZE;
X			cp += savecc;
X		}
X		*cp++ = (char)ch;
X		*cp = '\0';
X		scancc = qparse(etoabuf, QPF_SLASH);
X		if (scancc != ++savecc)
X			break;
X		ch = fgetc(fp);
X	}
X	ungetc(ch, fp);
X	if (scancc < 0)
X		return NULL;
X	r = atoq(etoabuf);
X	if (ziszero(r->den)) {
X		qfree(r);
X		r = NULL;
X	}
X	return r;
X}
X#endif
X
X
X/*
X * Print a number in the specified output mode.
X * If MODE_DEFAULT is given, then the default output mode is used.
X * Any approximate output is flagged with a leading tilde.
X * Integers are always printed as themselves.
X */
Xvoid
Xqprintnum(q, outmode)
X	NUMBER *q;
X{
X	NUMBER tmpval;
X	long prec, exp;
X
X	if (outmode == MODE_DEFAULT)
X		outmode = _outmode_;
X	if ((outmode == MODE_FRAC) || ((outmode == MODE_REAL) && qisint(q))) {
X		qprintfr(q, 0L, FALSE);
X		return;
X	}
X	switch (outmode) {
X		case MODE_INT:
X			if (qisfrac(q))
X				PUTCHAR('~');
X			qprintfd(q, 0L);
X			break;
X
X		case MODE_REAL:
X			prec = qplaces(q);
X			if ((prec < 0) || (prec > _outdigits_)) {
X				prec = _outdigits_;
X				PUTCHAR('~');
X			}
X			qprintff(q, 0L, prec);
X			break;
X
X		case MODE_EXP:
X			if (qiszero(q)) {
X				PUTCHAR('0');
X				return;
X			}
X			tmpval = *q;
X			tmpval.num.sign = 0;
X			exp = qilog10(&tmpval);
X			if (exp == 0) {		/* in range to output as real */
X				qprintnum(q, MODE_REAL);
X				return;
X			}
X			tmpval.num = _one_;
X			tmpval.den = _one_;
X			if (exp > 0)
X				ztenpow(exp, &tmpval.den);
X			else
X				ztenpow(-exp, &tmpval.num);
X			q = qmul(q, &tmpval);
X			zfree(tmpval.num);
X			zfree(tmpval.den);
X			qprintnum(q, MODE_REAL);
X			qfree(q);
X			PRINTF1("e%ld", exp);
X			break;
X
X		case MODE_HEX:
X			qprintfx(q, 0L);
X			break;
X
X		case MODE_OCTAL:
X			qprintfo(q, 0L);
X			break;
X
X		case MODE_BINARY:
X			qprintfb(q, 0L);
X			break;
X
X		default:
X			math_error("Bad mode for print");
X	}
X}
X
X
X/*
X * Print a number in floating point representation.
X * Example:  193.784
X */
Xvoid
Xqprintff(q, width, precision)
X	NUMBER *q;
X	long width;
X	long precision;
X{
X	ZVALUE z, z1;
X
X	if (precision != scalefactor) {
X		if (scalenumber.v)
X			zfree(scalenumber);
X		ztenpow(precision, &scalenumber);
X		scalefactor = precision;
X	}
X	if (scalenumber.v)
X		zmul(q->num, scalenumber, &z);
X	else
X		z = q->num;
X	if (qisfrac(q)) {
X		zquo(z, q->den, &z1);
X		if (z.v != q->num.v)
X			zfree(z);
X		z = z1;
X	}
X	if (qisneg(q) && ziszero(z))
X		PUTCHAR('-');
X	zprintval(z, precision, width);
X	if (z.v != q->num.v)
X		zfree(z);
X}
X
X
X/*
X * Print a number in exponential notation.
X * Example: 4.1856e34
X */
X/*ARGSUSED*/
Xvoid
Xqprintfe(q, width, precision)
X	register NUMBER *q;
X	long width;
X	long precision;
X{
X	long exponent;
X	NUMBER q2;
X	ZVALUE num, den, tenpow, tmp;
X
X	if (qiszero(q)) {
X		PUTSTR("0.0");
X		return;
X	}
X	num = q->num;
X	den = q->den;
X	num.sign = 0;
X	exponent = zdigits(num) - zdigits(den);
X	if (exponent > 0) {
X		ztenpow(exponent, &tenpow);
X		zmul(den, tenpow, &tmp);
X		zfree(tenpow);
X		den = tmp;
X	}
X	if (exponent < 0) {
X		ztenpow(-exponent, &tenpow);
X		zmul(num, tenpow, &tmp);
X		zfree(tenpow);
X		num = tmp;
X	}
X	if (zrel(num, den) < 0) {
X		zmuli(num, 10L, &tmp);
X		if (num.v != q->num.v)
X			zfree(num);
X		num = tmp;
X		exponent--;
X	}
X	q2.num = num;
X	q2.den = den;
X	q2.num.sign = q->num.sign;
X	qprintff(&q2, 0L, precision);
X	if (exponent)
X		PRINTF1("e%ld", exponent);
X	if (num.v != q->num.v)
X		zfree(num);
X	if (den.v != q->den.v)
X		zfree(den);
X}
X
X
X/*
X * Print a number in rational representation.
X * Example: 397/37
X */
Xvoid
Xqprintfr(q, width, force)
X	NUMBER *q;
X	long width;
X	BOOL force;
X{
X	zprintval(q->num, 0L, width);
X	if (force || qisfrac(q)) {
X		PUTCHAR('/');
X		zprintval(q->den, 0L, width);
X	}
X}
X
X
X/*
X * Print a number as an integer (truncating fractional part).
X * Example: 958421
X */
Xvoid
Xqprintfd(q, width)
X	NUMBER *q;
X	long width;
X{
X	ZVALUE z;
X
X	if (qisfrac(q)) {
X		zquo(q->num, q->den, &z);
X		zprintval(z, 0L, width);
X		zfree(z);
X	} else
X		zprintval(q->num, 0L, width);
X}
X
X
X/*
X * Print a number in hex.
X * This prints the numerator and denominator in hex.
X */
Xvoid
Xqprintfx(q, width)
X	NUMBER *q;
X	long width;
X{
X	zprintx(q->num, width);
X	if (qisfrac(q)) {
X		PUTCHAR('/');
X		zprintx(q->den, 0L);
X	}
X}
X
X
X/*
X * Print a number in binary.
X * This prints the numerator and denominator in binary.
X */
Xvoid
Xqprintfb(q, width)
X	NUMBER *q;
X	long width;
X{
X	zprintb(q->num, width);
X	if (qisfrac(q)) {
X		PUTCHAR('/');
X		zprintb(q->den, 0L);
X	}
X}
X
X
X/*
X * Print a number in octal.
X * This prints the numerator and denominator in octal.
X */
Xvoid
Xqprintfo(q, width)
X	NUMBER *q;
X	long width;
X{
X	zprinto(q->num, width);
X	if (qisfrac(q)) {
X		PUTCHAR('/');
X		zprinto(q->den, 0L);
X	}
X}
X
X
X/*
X * Convert a string to a number in rational, floating point,
X * exponential notation, hex, or octal.
X *	q = atoq(string);
X */
XNUMBER *
Xatoq(s)
X	register char *s;
X{
X	register NUMBER *q;
X	register char *t;
X	ZVALUE div, newnum, newden, tmp;
X	long decimals, exp;
X	BOOL hex, negexp;
X
X	q = qalloc();
X	decimals = 0;
X	exp = 0;
X	negexp = FALSE;
X	hex = FALSE;
X	t = s;
X	if ((*t == '+') || (*t == '-'))
X		t++;
X	if ((*t == '0') && ((t[1] == 'x') || (t[1] == 'X'))) {
X		hex = TRUE;
X		t += 2;
X	}
X	while (((*t >= '0') && (*t <= '9')) || (hex &&
X		(((*t >= 'a') && (*t <= 'f')) || ((*t >= 'A') && (*t <= 'F')))))
X			t++;
X	if (*t == '/') {
X		t++;
X		atoz(t, &q->den);
X	} else if ((*t == '.') || (*t == 'e') || (*t == 'E')) {
X		if (*t == '.') {
X			t++;
X			while ((*t >= '0') && (*t <= '9')) {
X				t++;
X				decimals++;
X			}
X		}
X		/*
X		 * Parse exponent if any
X		 */
X		if ((*t == 'e') || (*t == 'E')) {
X			t++;
X			if (*t == '+')
X				t++;
X			else if (*t == '-') {
X				negexp = TRUE;
X				t++;
X			}
X			while ((*t >= '0') && (*t <= '9')) {
X				exp = (exp * 10) + *t++ - '0';
X				if (exp > 1000000)
X					math_error("Exponent too large");
X			}
X		}
X		ztenpow(decimals, &q->den);
X	}
X	atoz(s, &q->num);
X	if (qiszero(q)) {
X		qfree(q);
X		return qlink(&_qzero_);
X	}
X	/*
X	 * Apply the exponential if any
X	 */
X	if (exp) {
X		ztenpow(exp, &tmp);
X		if (negexp) {
X			zmul(q->den, tmp, &newden);
X			zfree(q->den);
X			q->den = newden;
X		} else {
X			zmul(q->num, tmp, &newnum);
X			zfree(q->num);
X			q->num = newnum;
X		}
X		zfree(tmp);
X	}
X	/*
X	 * Reduce the fraction to lowest terms
X	 */
X	if (zisunit(q->num) || zisunit(q->den))
X		return q;
X	zgcd(q->num, q->den, &div);
X	if (zisunit(div))
X		return q;
X	zquo(q->num, div, &newnum);
X	zfree(q->num);
X	zquo(q->den, div, &newden);
X	zfree(q->den);
X	q->num = newnum;
X	q->den = newden;
X	return q;
X}
X
X
X/*
X * Parse a number in any of the various legal forms, and return the count
X * of characters that are part of a legal number.  Numbers can be either a
X * decimal integer, possibly two decimal integers separated with a slash, a
X * floating point or exponential number, a hex number beginning with "0x",
X * a binary number beginning with "0b", or an octal number beginning with "0".
X * The flags argument modifies the end of number testing for ease in handling
X * fractions or complex numbers.  Minus one is returned if the number format
X * is definitely illegal.
X */
Xlong
Xqparse(cp, flags)
X	register char *cp;
X{
X	char *oldcp;
X
X	oldcp = cp;
X	if ((*cp == '+') || (*cp == '-'))
X		cp++;
X	if ((*cp == '+') || (*cp == '-'))
X		return -1;
X	if ((*cp == '0') && ((cp[1] == 'x') || (cp[1] == 'X'))) {	/* hex */
X		cp += 2;
X		while (((*cp >= '0') && (*cp <= '9')) ||
X			((*cp >= 'a') && (*cp <= 'f')) ||
X			((*cp >= 'A') && (*cp <= 'F')))
X				cp++;
X		goto done;
X	}
X	if ((*cp == '0') && ((cp[1] == 'b') || (cp[1] == 'B'))) {	/* binary */
X		cp += 2;
X		while ((*cp == '0') || (*cp == '1'))
X			cp++;
X		goto done;
X	}
X	if ((*cp == '0') && (cp[1] >= '0') && (cp[1] <= '9')) { /* octal */
X		while ((*cp >= '0') && (*cp <= '7'))
X			cp++;
X		goto done;
X	}
X	/*
X	 * Number is decimal, but can still be a fraction or real or exponential.
X	 */
X	while ((*cp >= '0') && (*cp <= '9'))
X		cp++;
X	if (*cp == '/' && flags & QPF_SLASH) {	/* fraction */
X		cp++;
X		while ((*cp >= '0') && (*cp <= '9'))
X			cp++;
X		goto done;
X	}
X	if (*cp == '.') {	/* floating point */
X		cp++;
X		while ((*cp >= '0') && (*cp <= '9'))
X			cp++;
X	}
X	if ((*cp == 'e') || (*cp == 'E')) {	/* exponential */
X		cp++;
X		if ((*cp == '+') || (*cp == '-'))
X			cp++;
X		if ((*cp == '+') || (*cp == '-'))
X			return -1;
X		while ((*cp >= '0') && (*cp <= '9'))
X			cp++;
X	}
X
Xdone:
X	if (((*cp == 'i') || (*cp == 'I')) && (flags & QPF_IMAG))
X		cp++;
X	if ((*cp == '.') || ((*cp == '/') && (flags & QPF_SLASH)) ||
X		((*cp >= '0') && (*cp <= '9')) ||
X		((*cp >= 'a') && (*cp <= 'z')) ||
X		((*cp >= 'A') && (*cp <= 'Z')))
X			return -1;
X	return (cp - oldcp);
X}
X
X/* END CODE */
SHAR_EOF
chmod 0644 calc2.9.0/qio.c || echo "restore of calc2.9.0/qio.c fails"
set `wc -c calc2.9.0/qio.c`;Sum=$1
if test "$Sum" != "12895"
then echo original size 12895, current size $Sum;fi
echo "x - extracting calc2.9.0/qmath.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > calc2.9.0/qmath.c &&
X/*
X * Copyright (c) 1993 David I. Bell
X * Permission is granted to use, distribute, or modify this source,
X * provided that this copyright notice remains intact.
X *
X * Extended precision rational arithmetic primitive routines
X */
X
X#include "qmath.h"
X
X
XNUMBER _qzero_ =	{ { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
XNUMBER _qone_ =		{ { _oneval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
Xstatic NUMBER _qtwo_ =	{ { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
Xstatic NUMBER _qten_ =	{ { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1 };
XNUMBER _qnegone_ =	{ { _oneval_, 1, 1 }, { _oneval_, 1, 0 }, 1 };
XNUMBER _qonehalf_ =	{ { _oneval_, 1, 0 }, { _twoval_, 1, 0 }, 1 };
X
X
X/*
X * Create another copy of a number.
X *	q2 = qcopy(q1);
X */
XNUMBER *
Xqcopy(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	r = qalloc();
X	r->num.sign = q->num.sign;
X	if (!zisunit(q->num)) {
X		r->num.len = q->num.len;
X		r->num.v = alloc(r->num.len);
X		zcopyval(q->num, r->num);
X	}
X	if (!zisunit(q->den)) {
X		r->den.len = q->den.len;
X		r->den.v = alloc(r->den.len);
X		zcopyval(q->den, r->den);
X	}
X	return r;
X}
X
X
X/*
X * Convert a number to a normal integer.
X *	i = qtoi(q);
X */
Xlong
Xqtoi(q)
X	register NUMBER *q;
X{
X	long i;
X	ZVALUE res;
X
X	if (qisint(q))
X		return ztoi(q->num);
X	zquo(q->num, q->den, &res);
X	i = ztoi(res);
X	zfree(res);
X	return i;
X}
X
X
X/*
X * Convert a normal integer into a number.
X *	q = itoq(i);
X */
XNUMBER *
Xitoq(i)
X	long i;
X{
X	register NUMBER *q;
X
X	if ((i >= -1) && (i <= 10)) {
X		switch ((int) i) {
X			case 0: q = &_qzero_; break;
X			case 1: q = &_qone_; break;
X			case 2: q = &_qtwo_; break;
X			case 10: q = &_qten_; break;
X			case -1: q = &_qnegone_; break;
X			default: q = NULL;
X		}
X		if (q)
X			return qlink(q);
X	}
X	q = qalloc();
X	itoz(i, &q->num);
X	return q;
X}
X
X
X/*
X * Create a number from the given integral numerator and denominator.
X *	q = iitoq(inum, iden);
X */
XNUMBER *
Xiitoq(inum, iden)
X	long inum, iden;
X{
X	register NUMBER *q;
X	long d;
X	BOOL sign;
X
X	if (iden == 0)
X		math_error("Division by zero");
X	if (inum == 0)
X		return qlink(&_qzero_);
X	sign = 0;
X	if (inum < 0) {
X		sign = 1;
X		inum = -inum;
X	}
X	if (iden < 0) {
X		sign = 1 - sign;
X		iden = -iden;
X	}
X	d = iigcd(inum, iden);
X	inum /= d;
X	iden /= d;
X	if (iden == 1)
X		return itoq(sign ? -inum : inum);
X	q = qalloc();
X	if (inum != 1)
X		itoz(inum, &q->num);
X	itoz(iden, &q->den);
X	q->num.sign = sign;
X	return q;
X}
X
X
X/*
X * Add two numbers to each other.
X *	q3 = qadd(q1, q2);
X */
XNUMBER *
Xqadd(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	NUMBER *r;
X	ZVALUE t1, t2, temp, d1, d2, vpd1, upd1;
X
X	if (qiszero(q1))
X		return qlink(q2);
X	if (qiszero(q2))
X		return qlink(q1);
X	r = qalloc();
X	/*
X	 * If either number is an integer, then the result is easy.
X	 */
X	if (qisint(q1) && qisint(q2)) {
X		zadd(q1->num, q2->num, &r->num);
X		return r;
X	}
X	if (qisint(q2)) {
X		zmul(q1->den, q2->num, &temp);
X		zadd(q1->num, temp, &r->num);
X		zfree(temp);
X		zcopy(q1->den, &r->den);
X		return r;
X	}
X	if (qisint(q1)) {
X		zmul(q2->den, q1->num, &temp);
X		zadd(q2->num, temp, &r->num);
X		zfree(temp);
X		zcopy(q2->den, &r->den);
X		return r;
X	}
X	/*
X	 * Both arguments are true fractions, so we need more work.
X	 * If the denominators are relatively prime, then the answer is the
X	 * straightforward cross product result with no need for reduction.
X	 */
X	zgcd(q1->den, q2->den, &d1);
X	if (zisunit(d1)) {
X		zfree(d1);
X		zmul(q1->num, q2->den, &t1);
X		zmul(q1->den, q2->num, &t2);
X		zadd(t1, t2, &r->num);
X		zfree(t1);
X		zfree(t2);
X		zmul(q1->den, q2->den, &r->den);
X		return r;
X	}
X	/*
X	 * The calculation is now more complicated.
X	 * See Knuth Vol 2 for details.
X	 */
X	zquo(q2->den, d1, &vpd1);
X	zquo(q1->den, d1, &upd1);
X	zmul(q1->num, vpd1, &t1);
X	zmul(q2->num, upd1, &t2);
X	zadd(t1, t2, &temp);
X	zfree(t1);
X	zfree(t2);
X	zfree(vpd1);
X	zgcd(temp, d1, &d2);
X	zfree(d1);
X	if (zisunit(d2)) {
X		zfree(d2);
X		r->num = temp;
X		zmul(upd1, q2->den, &r->den);
X		zfree(upd1);
X		return r;
X	}
X	zquo(temp, d2, &r->num);
X	zfree(temp);
X	zquo(q2->den, d2, &temp);
X	zfree(d2);
X	zmul(temp, upd1, &r->den);
X	zfree(temp);
X	zfree(upd1);
X	return r;
X}
X
X
X/*
X * Subtract one number from another.
X *	q3 = qsub(q1, q2);
X */
XNUMBER *
Xqsub(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	NUMBER t, *r;
X
X	if (q1 == q2)
X		return qlink(&_qzero_);
X	if (qiszero(q2))
X		return qlink(q1);
X	if (qisint(q1) && qisint(q2)) {
X		r = qalloc();
X		zsub(q1->num, q2->num, &r->num);
X		return r;
X	}
X	t = *q2;
X	t.num.sign = !t.num.sign;
X	return qadd(q1, &t);
X}
X
X
X/*
X * Increment a number by one.
X */
XNUMBER *
Xqinc(q)
X	NUMBER *q;
X{
X	NUMBER *r;
X
X	r = qalloc();
X	if (qisint(q)) {
X		zadd(q->num, _one_, &r->num);
X		return r;
X	}
X	zadd(q->num, q->den, &r->num);
X	zcopy(q->den, &r->den);
X	return r;
X}
X
X
X/*
X * Decrement a number by one.
X */
XNUMBER *
Xqdec(q)
X	NUMBER *q;
X{
X	NUMBER *r;
X
X	r = qalloc();
X	if (qisint(q)) {
X		zsub(q->num, _one_, &r->num);
X		return r;
X	}
X	zsub(q->num, q->den, &r->num);
X	zcopy(q->den, &r->den);
X	return r;
X}
X
X
X/*
X * Add a normal small integer value to an arbitrary number.
X */
XNUMBER *
Xqaddi(q1, n)
X	NUMBER *q1;
X	long n;
X{
X	NUMBER addnum;		/* temporary number */
X	HALF addval[2];		/* value of small number */
X	BOOL neg;		/* TRUE if number is neg */
X
X	if (n == 0)
X		return qlink(q1);
X	if (n == 1)
X		return qinc(q1);
X	if (n == -1)
X		return qdec(q1);
X	if (qiszero(q1))
X		return itoq(n);
X	addnum.num.sign = 0;
X	addnum.num.len = 1;
X	addnum.num.v = addval;
X	addnum.den = _one_;
X	neg = (n < 0);
X	if (neg)
X		n = -n;
X	addval[0] = (HALF) n;
X	n = (((FULL) n) >> BASEB);
X	if (n) {
X		addval[1] = (HALF) n;
X		addnum.num.len = 2;
X	}
X	if (neg)
X		return qsub(q1, &addnum);
X	else
X		return qadd(q1, &addnum);
X}
X
X
X/*
X * Multiply two numbers.
X *	q3 = qmul(q1, q2);
X */
XNUMBER *
Xqmul(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	NUMBER *r;			/* returned value */
X	ZVALUE n1, n2, d1, d2;		/* numerators and denominators */
X	ZVALUE tmp;
X
X	if (qiszero(q1) || qiszero(q2))
X		return qlink(&_qzero_);
X	if (qisone(q1))
X		return qlink(q2);
X	if (qisone(q2))
X		return qlink(q1);
X	if (qisint(q1) && qisint(q2)) {	/* easy results if integers */
X		r = qalloc();
X		zmul(q1->num, q2->num, &r->num);
X		return r;
X	}
X	n1 = q1->num;
X	n2 = q2->num;
X	d1 = q1->den;
X	d2 = q2->den;
X	if (ziszero(d1) || ziszero(d2))
X		math_error("Division by zero");
X	if (ziszero(n1) || ziszero(n2))
X		return qlink(&_qzero_);
X	if (!zisunit(n1) && !zisunit(d2)) {	/* possibly reduce */
X		zgcd(n1, d2, &tmp);
X		if (!zisunit(tmp)) {
X			zquo(q1->num, tmp, &n1);
X			zquo(q2->den, tmp, &d2);
X		}
X		zfree(tmp);
X	}
X	if (!zisunit(n2) && !zisunit(d1)) {	/* again possibly reduce */
X		zgcd(n2, d1, &tmp);
X		if (!zisunit(tmp)) {
X			zquo(q2->num, tmp, &n2);
X			zquo(q1->den, tmp, &d1);
X		}
X		zfree(tmp);
X	}
X	r = qalloc();
X	zmul(n1, n2, &r->num);
X	zmul(d1, d2, &r->den);
X	if (q1->num.v != n1.v)
X		zfree(n1);
X	if (q1->den.v != d1.v)
X		zfree(d1);
X	if (q2->num.v != n2.v)
X		zfree(n2);
X	if (q2->den.v != d2.v)
X		zfree(d2);
X	return r;
X}
X
X
X/*
X * Multiply a number by a small integer.
X *	q2 = qmuli(q1, n);
X */
XNUMBER *
Xqmuli(q, n)
X	NUMBER *q;
X	long n;
X{
X	NUMBER *r;
X	long d;			/* gcd of multiplier and denominator */
X	int sign;
X
X	if ((n == 0) || qiszero(q))
X		return qlink(&_qzero_);
X	if (n == 1)
X		return qlink(q);
X	r = qalloc();
X	if (qisint(q)) {
X		zmuli(q->num, n, &r->num);
X		return r;
X	}
X	sign = 1;
X	if (n < 0) {
X		n = -n;
X		sign = -1;
X	}
X	d = zmodi(q->den, n);
X	d = iigcd(d, n);
X	zmuli(q->num, (n * sign) / d, &r->num);
X	(void) zdivi(q->den, d, &r->den);
X	return r;
X}
X
X
X/*
X * Divide two numbers (as fractions).
X *	q3 = qdiv(q1, q2);
X */
XNUMBER *
Xqdiv(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	NUMBER temp;
X
X	if (qiszero(q2))
X		math_error("Division by zero");
X	if ((q1 == q2) || !qcmp(q1, q2))
X		return qlink(&_qone_);
X	if (qisone(q1))
X		return qinv(q2);
X	temp.num = q2->den;
X	temp.den = q2->num;
X	temp.num.sign = temp.den.sign;
X	temp.den.sign = 0;
X	temp.links = 1;
X	return qmul(q1, &temp);
X}
X
X
X/*
X * Divide a number by a small integer.
X *	q2 = qdivi(q1, n);
X */
XNUMBER *
Xqdivi(q, n)
X	NUMBER *q;
X	long n;
X{
X	NUMBER *r;
X	long d;			/* gcd of divisor and numerator */
X	int sign;
X
X	if (n == 0)
X		math_error("Division by zero");
X	if ((n == 1) || qiszero(q))
X		return qlink(q);
X	sign = 1;
X	if (n < 0) {
X		n = -n;
X		sign = -1;
X	}
X	r = qalloc();
X	d = zmodi(q->num, n);
X	d = iigcd(d, n);
X	(void) zdivi(q->num, d * sign, &r->num);
X	zmuli(q->den, n / d, &r->den);
X	return r;
X}
X
X
X/*
X * Return the quotient when one number is divided by another.
X * This works for fractional values also, and in all cases:
X *	qquo(q1, q2) = int(q1 / q2).
X */
XNUMBER *
Xqquo(q1, q2)
X	register NUMBER *q1, *q2;
X{
X	ZVALUE num, den, res;
X	NUMBER *q;
X
X	if (zisunit(q1->num))
X		num = q2->den;
X	else if (zisunit(q2->den))
X		num = q1->num;
X	else
X		zmul(q1->num, q2->den, &num);
X	if (zisunit(q1->den))
X		den = q2->num;
X	else if (zisunit(q2->num))
X		den = q1->den;
X	else
X		zmul(q1->den, q2->num, &den);
X	zquo(num, den, &res);
X	if ((num.v != q2->den.v) && (num.v != q1->num.v))
X		zfree(num);
X	if ((den.v != q2->num.v) && (den.v != q1->den.v))
X		zfree(den);
X	if (ziszero(res)) {
X		zfree(res);
X		return qlink(&_qzero_);
X	}
X	res.sign = (q1->num.sign != q2->num.sign);
X	if (zisunit(res)) {
X		q = (res.sign ? &_qnegone_ : &_qone_);
X		zfree(res);
X		return qlink(q);
X	}
X	q = qalloc();
X	q->num = res;
X	return q;
X}
X
X
X/*
X * Return the absolute value of a number.
X *	q2 = qabs(q1);
X */
XNUMBER *
Xqabs(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (q->num.sign == 0)
X		return qlink(q);
X	r = qalloc();
X	if (!zisunit(q->num))
X		zcopy(q->num, &r->num);
X	if (!zisunit(q->den))
X		zcopy(q->den, &r->den);
X	r->num.sign = 0;
X	return r;
X}
X
X
X/*
X * Negate a number.
X *	q2 = qneg(q1);
X */
XNUMBER *
Xqneg(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (qiszero(q))
X		return qlink(&_qzero_);
X	r = qalloc();
X	if (!zisunit(q->num))
X		zcopy(q->num, &r->num);
X	if (!zisunit(q->den))
X		zcopy(q->den, &r->den);
X	r->num.sign = !q->num.sign;
X	return r;
X}
X
X
X/*
X * Return the sign of a number (-1, 0 or 1)
X */
XNUMBER *
Xqsign(q)
X	NUMBER *q;
X{
X	if (qiszero(q))
X		return qlink(&_qzero_);
X	if (qisneg(q))
X		return qlink(&_qnegone_);
X	return qlink(&_qone_);
X}
X
X
X/*
X * Invert a number.
X *	q2 = qinv(q1);
X */
XNUMBER *
Xqinv(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (qisunit(q)) {
X		r = (qisneg(q) ? &_qnegone_ : &_qone_);
X		return qlink(r);
X	}
X	if (qiszero(q))
X		math_error("Division by zero");
X	r = qalloc();
X	if (!zisunit(q->num))
X		zcopy(q->num, &r->den);
X	if (!zisunit(q->den))
X		zcopy(q->den, &r->num);
X	r->num.sign = q->num.sign;
X	r->den.sign = 0;
X	return r;
X}
X
X
X/*
X * Return just the numerator of a number.
X *	q2 = qnum(q1);
X */
XNUMBER *
Xqnum(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (qisint(q))
X		return qlink(q);
X	if (zisunit(q->num)) {
X		r = (qisneg(q) ? &_qnegone_ : &_qone_);
X		return qlink(r);
X	}
X	r = qalloc();
X	zcopy(q->num, &r->num);
X	return r;
X}
X
X
X/*
X * Return just the denominator of a number.
X *	q2 = qden(q1);
X */
XNUMBER *
Xqden(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (qisint(q))
X		return qlink(&_qone_);
X	r = qalloc();
X	zcopy(q->den, &r->num);
X	return r;
X}
X
X
X/*
X * Return the fractional part of a number.
X *	q2 = qfrac(q1);
X */
XNUMBER *
Xqfrac(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (qisint(q))
X		return qlink(&_qzero_);
X	if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
X		(q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
X			return qlink(q);
X	r = qalloc();
X	zmod(q->num, q->den, &r->num);
X	zcopy(q->den, &r->den);
X	r->num.sign = q->num.sign;
X	return r;
X}
X
X
X/*
X * Return the integral part of a number.
X *	q2 = qint(q1);
X */
XNUMBER *
Xqint(q)
X	register NUMBER *q;
X{
X	register NUMBER *r;
X
X	if (qisint(q))
X		return qlink(q);
X	if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
X		(q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
SHAR_EOF
echo "End of part 8"
echo "File calc2.9.0/qmath.c is continued in part 9"
echo "9" > s2_seq_.tmp
exit 0