qtrans.c

Calc Software Package for Number Calc

	zfree(X);
	qtmp = qalloc();
	k = zlowbit(sum);
	if (k) {
		zshift(sum, -k, &qtmp->num);
		zfree(sum);
	} else {
		qtmp->num = sum;
	}
	zbitvalue(m - 4 - k, &qtmp->den);
	res = qmappr(qtmp, epsilon, 24L);
	qfree(qtmp);
	return res;
}

/*
 * Inverse secant function
 */
NUMBER *
qasec(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp, *res;

	tmp = qinv(q);
	res = qacos(tmp, epsilon);
	qfree(tmp);
	return res;
}


/*
 * Inverse cosecant function
 */
NUMBER *
qacsc(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp, *res;

	tmp = qinv(q);
	res = qasin(tmp, epsilon);
	qfree(tmp);
	return res;
}


/*
 * Inverse cotangent function
 */
NUMBER *
qacot(NUMBER *q, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon for acot");
		/*NOTREACHED*/
	}
	if (qiszero(q)) {
		epsilon1 = qscale(epsilon, 1L);
		tmp1 = qpi(epsilon1);
		qfree(epsilon1);
		tmp2 = qscale(tmp1, -1L);
		qfree(tmp1);
		return tmp2;
	}
	tmp1 = qinv(q);
	if (!qisneg(q)) {
		tmp2 = qatan(tmp1, epsilon);
		qfree(tmp1);
		return tmp2;
	}
	epsilon1 = qscale(epsilon, -2L);
	tmp2 = qatan(tmp1, epsilon1);
	qfree(tmp1);
	tmp1 = qpi(epsilon1);
	qfree(epsilon1);
	tmp3 = qqadd(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qmappr(tmp3, epsilon, 24L);
	qfree(tmp3);
	return tmp1;
}


/*
 * Calculate the angle which is determined by the point (x,y).
 * This is the same as atan(y/x) for positive x, but is continuous
 * except for y = 0, x <= 0.  By convention, y is the first argument.
 * For all x, y, -pi < atan2 <= pi.  For example, qatan2(1, -1) = 3/4 * pi.
 */
NUMBER *
qatan2(NUMBER *qy, NUMBER *qx, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for atan2");
		/*NOTREACHED*/
	}
	if (qiszero(qy) && qiszero(qx)) {
		/* conform to 4.3BSD ANSI/IEEE 754-1985 math lib */
		return qlink(&_qzero_);
	}
	/*
	 * If the point is on the negative real axis, then the answer is pi.
	 */
	if (qiszero(qy) && qisneg(qx))
		return qpi(epsilon);
	/*
	 * If the point is in the right half plane, then use the normal atan.
	 */
	if (!qisneg(qx) && !qiszero(qx)) {
		if (qiszero(qy))
			return qlink(&_qzero_);
		tmp1 = qqdiv(qy, qx);
		tmp2 = qatan(tmp1, epsilon);
		qfree(tmp1);
		return tmp2;
	}
	/*
	 * The point is in the left half plane (x <= 0) with nonzero y.
	 * Calculate the angle by using the formula:
	 *	atan2(y,x) = 2 * atan(sgn(y) * sqrt((x/y)^2 + 1) - x/y).
	 */
	epsilon2 = qscale(epsilon, -4L);
	tmp1 = qqdiv(qx, qy);
	tmp2 = qsquare(tmp1);
	tmp3 = qqadd(tmp2, &_qone_);
	qfree(tmp2);
	tmp2 = qsqrt(tmp3, epsilon2, 24L | (qy->num.sign * 64));
	qfree(tmp3);
	tmp3 = qsub(tmp2, tmp1);
	qfree(tmp2);
	qfree(tmp1);
	qfree(epsilon2);
	epsilon2 = qscale(epsilon, -1L);
	tmp1 = qatan(tmp3, epsilon2);
	qfree(epsilon2);
	qfree(tmp3);
	tmp2 = qscale(tmp1, 1L);
	qfree(tmp1);
	return tmp2;
}


/*
 * Calculate the value of pi to within the required epsilon.
 * This uses the following formula which only needs integer calculations
 * except for the final operation:
 *	pi = 1 / SUMOF(comb(2 * N, N) ^ 3 * (42 * N + 5) / 2 ^ (12 * N + 4)),
 * where the summation runs from N=0.  This formula gives about 6 bits of
 * accuracy per term.  Since the denominator for each term is a power of two,
 * we can simply use shifts to sum the terms.  The combinatorial numbers
 * in the formula are calculated recursively using the formula:
 *	comb(2*(N+1), N+1) = 2 * comb(2 * N, N) * (2 * N + 1) / N.
 */
NUMBER *
qpi(NUMBER *epsilon)
{
	ZVALUE comb;			/* current combinatorial value */
	ZVALUE sum;			/* current sum */
	ZVALUE tmp1, tmp2;
	NUMBER *r, *t1, qtmp;
	long shift;			/* current shift of result */
	long N;				/* current term number */
	long bits;			/* needed number of bits of precision */
	long t;

	if (qiszero(epsilon)) {
		math_error("zero epsilon value for pi");
		/*NOTREACHED*/
	}
	if (epsilon == pivalue[0])
		return qlink(pivalue[1]);
	if (pivalue[0]) {
		qfree(pivalue[0]);
		qfree(pivalue[1]);
	}
	bits = -qilog2(epsilon) + 4;
	if (bits < 4)
		bits = 4;
	comb = _one_;
	itoz(5L, &sum);
	N = 0;
	shift = 4;
	do {
		t = 1 + (++N & 0x1);
		(void) zdivi(comb, N / (3 - t), &tmp1);
		zfree(comb);
		zmuli(tmp1, t * (2 * N - 1), &comb);
		zfree(tmp1);
		zsquare(comb, &tmp1);
		zmul(comb, tmp1, &tmp2);
		zfree(tmp1);
		zmuli(tmp2, 42 * N + 5, &tmp1);
		zfree(tmp2);
		zshift(sum, 12L, &tmp2);
		zfree(sum);
		zadd(tmp1, tmp2, &sum);
		t = zhighbit(tmp1);
		zfree(tmp1);
		zfree(tmp2);
		shift += 12;
	} while ((shift - t) < bits);
	zfree(comb);
	qtmp.num = _one_;
	qtmp.den = sum;
	t1 = qscale(&qtmp, shift);
	zfree(sum);
	r = qmappr(t1, epsilon, 24L);
	qfree(t1);
	pivalue[0] = qlink(epsilon);
	pivalue[1] = qlink(r);
	return r;
}

/*
 * Calculate the exponential function to the nearest or next to nearest
 * multiple of the positive number epsilon.
 */
NUMBER *
qexp(NUMBER *q, NUMBER *epsilon)
{
	long m, n;
	NUMBER *tmp1, *tmp2;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon value for exp");
		/*NOTREACHED*/
	}
	if (qiszero(q))
		return qlink(&_qone_);
	tmp1 = qmul(q, &_qlge_);
	m = qtoi(tmp1);		/* exp(q) < 2^(m+1) or m == MAXLONG */
	qfree(tmp1);

	if (m > (1 << 30))
		return NULL;

	n = qilog2(epsilon);	/* 2^n <= epsilon < 2^(n+1) */
	if (m < n)
		return qlink(&_qzero_);
	tmp1 = qqabs(q);
	tmp2 = qexprel(tmp1, m - n + 1);
	qfree(tmp1);
	if (tmp2 == NULL)
		return NULL;
	if (qisneg(q)) {
		tmp1 = qinv(tmp2);
		qfree(tmp2);
		tmp2 = tmp1;
	}
	tmp1 = qmappr(tmp2, epsilon, 24L);
	qfree(tmp2);
	return tmp1;
}


/*
 * Calculate the exponential function with relative error corresponding
 * to a specified number of significant bits
 * Requires *q >= 0, bitnum >= 0.
 * This returns NULL if more than 2^30 working bits would be required.
 */
static NUMBER *
qexprel(NUMBER *q, long bitnum)
{
	long n, m, k, h, s, t, d;
	NUMBER *qtmp1;
	ZVALUE X, B, sum, term, ztmp1, ztmp2;

	h = qilog2(q);
	k = bitnum + h + 1;
	if (k < 0)
		return qlink(&_qone_);
	s = k;
	if (k) {
		do {
			t = s;
			s = (s + k/s)/2;
		}
		while (t > s);
	}		/* s is int(sqrt(k)) */
	s++;
	if (s < -h)
		s = -h;
	n = h + s;	/* n is number of squarings that will be required */
	m = bitnum + n;
	if (m > (1 << 30))
		return NULL;
	while (s > 0) {		/* increasing m by ilog2(s) */
		s >>= 1;
		m++;
	}			/* m is working number of bits */
	qtmp1 = qscale(q, m - n);
	zquo(qtmp1->num, qtmp1->den, &X, 24);
	qfree(qtmp1);
	if (ziszero(X)) {
		zfree(X);
		return qlink(&_qone_);
	}
	zbitvalue(m, &sum);
	zcopy(X, &term);
	d = 1;
	do {
		zadd(sum, term, &ztmp1);
		zfree(sum);
		sum = ztmp1;
		zmul(term, X, &ztmp1);
		zfree(term);
		zshift(ztmp1, -m, &ztmp2);
		zfree(ztmp1);
		zdivi(ztmp2, ++d, &term);
		zfree(ztmp2);
	}
	while (!ziszero(term));
	zfree(term);
	zfree(X);
	k = 0;
	zbitvalue(2 * m + 1, &B);
	while (n-- > 0) {
		k *= 2;
		zsquare(sum, &ztmp1);
		zfree(sum);
		if (zrel(ztmp1, B) >= 0) {
			zshift(ztmp1, -m - 1, &sum);
			k++;
		} else {
			zshift(ztmp1, -m, &sum);
		}
		zfree(ztmp1);
	}
	zfree(B);
	h = zlowbit(sum);
	qtmp1 = qalloc();
	if (m > h + k) {
		zshift(sum, -h, &qtmp1->num);
		zbitvalue(m - h - k, &qtmp1->den);
	}
	else
		zshift(sum, k - m, &qtmp1->num);
	zfree(sum);
	return qtmp1;
}


/*
 * Calculate the natural logarithm of a number accurate to the specified
 * positive epsilon.
 */
NUMBER *
qln(NUMBER *q, NUMBER *epsilon)
{
	long m, n, k, h, d;
	ZVALUE term, sum, mul, pow, X, D, B, ztmp;
	NUMBER *qtmp, *res;
	BOOL neg;

	if (qiszero(q) || qiszero(epsilon)) {
		math_error("logarithm of 0");
		/*NOTREACHED*/
	}
	if (qisunit(q))
		return qlink(&_qzero_);
	q = qqabs(q);			/* Ignore sign of q */
	neg = (zrel(q->num, q->den) < 0);
	if (neg) {
		qtmp = qinv(q);
		qfree(q);
		q = qtmp;
	}
	k = qilog2(q);
	m = -qilog2(epsilon);		/* m will be number of working bits */
	if (m < 0)
		m = 0;
	h = k;
	while (h > 0) {
		h /= 2;
		m++;			/* Add 1 for each sqrt until X < 2 */
	}
	m += 18;	/* 8 more sqrts, 8 for rounding, 2 for epsilon/4 */
	qtmp = qscale(q, m - k);
	zquo(qtmp->num, qtmp->den, &X, 24L);
	qfree(q);
	qfree(qtmp);

	zbitvalue(m, &D);		/* Now "q" = X/D */
	zbitvalue(m - 8, &ztmp);
	zadd(D, ztmp, &B);		/* Will take sqrts until X <= B */
	zfree(ztmp);

	n = 1;			/* n is to count 1 + number of sqrts */

	while (k > 0 || zrel(X, B) > 0) {
		n++;
		zshift(X, m + (k & 1), &ztmp);
		zfree(X);
		zsqrt(ztmp, &X, 24);
		zfree(ztmp)
		k /= 2;
	}
	zfree(B);
	zsub(X, D, &pow);	/* pow, mul used as tmps */
	zadd(X, D, &mul);
	zfree(X);
	zfree(D);
	zshift(pow, m, &ztmp);
	zfree(pow);
	zquo(ztmp, mul, &pow, 24);	/* pow now (X - D)/(X + D) */
	zfree(ztmp);
	zfree(mul);

	zcopy(pow, &sum);	/* pow is first term of sum */
	zsquare(pow, &ztmp);
	zshift(ztmp, -m, &mul); /* mul is now multiplier for powers */
	zfree(ztmp);

	d = 1;
	for (;;) {
		zmul(pow, mul, &ztmp);
		zfree(pow);
		zshift(ztmp, -m, &pow);
		zfree(ztmp);
		d += 2;
		zdivi(pow, d, &term);	/* Round down div should be round off */
		if (ziszero(term)) {
			zfree(term);
			break;
		}
		zadd(sum, term, &ztmp);
		zfree(term);
		zfree(sum);
		sum = ztmp;
	}
	zfree(pow);
	zfree(mul);
	qtmp = qalloc();	/* qtmp is to be 2^n * sum / 2^m */
	k = zlowbit(sum);
	sum.sign = neg;
	if (k + n >= m) {
		zshift(sum, n - m, &qtmp->num);
	} else {
		if (k) {
			zshift(sum, -k, &qtmp->num);
			zfree(sum);
		} else {
			qtmp->num = sum;
		}
		zbitvalue(m - k - n, &qtmp->den);
	}
	res = qmappr(qtmp, epsilon, 24L);
	qfree(qtmp);
	return res;
}


/*
 * Calculate the base 10 logarithm
 *
 *	log(q) = ln(q) / ln(10)
 */
NUMBER *
qlog(NUMBER *q, NUMBER *epsilon)
{
	int need_new_ln_10 = TRUE;	/* FALSE => use cached ln_10 value */
	NUMBER *ln_q;			/* ln(x) */
	NUMBER *ret;			/* base 10 logarithm of x */

	/* firewall */
	if (qiszero(q) || qiszero(epsilon)) {
		math_error("logarithm of 0");
		/*NOTREACHED*/
	}

	/*
	 * shortcut for small integer powers of 10
	 */
	if (qisint(q) && qispos(q) && !zge8192b(q->num) && ziseven(q->num)) {
		BOOL is_10_power;	/* TRUE ==> q is a power of 10 */
		long ilog_10;		/* integer log base 10 */

		/* try for a quick small power of 10 log */
		ilog_10 = zlog10(q->num, &is_10_power );
		if (is_10_power == TRUE) {
			/* is small power of 10, return log */
			return itoq(ilog_10);
		}
		/* q is an even integer that is not a power of 10 */
	}

	/*
	 * compute ln(c) first
	 */
	ln_q = qln(q, epsilon);
	/* log(1) == 0 */
	if (qiszero(ln_q)) {
		return ln_q;
	}

	/*
	 * save epsilon for ln(10) if needed
	 */
	if (ln_10_epsilon == NULL) {
		/* first time call */
		ln_10_epsilon = qcopy(epsilon);
	} else if (qcmp(ln_10_epsilon, epsilon) == FALSE) {
		/* replaced cacheed value with epsilon arg */
		qfree(ln_10_epsilon);
		ln_10_epsilon = qcopy(epsilon);
	} else if (ln_10 != NULL) {
		/* the previously computed ln(2) is OK to use */
		need_new_ln_10 = FALSE;
	}

	/*
	 * compute ln(10) if needed
	 */
	if (need_new_ln_10 == TRUE) {
		if (ln_10 != NULL) {
			qfree(ln_10);
		}
		ln_10 = qln(&_qten_, ln_10_epsilon);
	}

	/*
	 * return ln(q) / ln(10)
	 */
	ret = qqdiv(ln_q, ln_10);
	qfree(ln_q);
	return ret;
}


/*
 * Calculate the result of raising one number to the power of another.
 * The result is calculated to the nearest or next to nearest multiple of
 * epsilon.
 */
NUMBER *
qpower(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
{
	NUMBER *tmp1, *tmp2, *epsilon2;
	NUMBER *q1tmp, *q2tmp;
	long m, n;

	if (qiszero(epsilon)) {
		math_error("Zero epsilon for power");
		/*NOTREACHED*/
	}
	if (qiszero(q1) && qisneg(q2)) {
		math_error("Negative power of zero");
		/*NOTREACHED*/
	}
	if (qiszero(q2) || qisone(q1))
		return qlink(&_qone_);
	if (qiszero(q1))
		return qlink(&_qzero_);