qtrans.c

早期freebsd实现

	 * Also make the argument non-negative.
	 */
	qs = qabs(q);
	scale = zhighbit(q->num) - zhighbit(q->den) + 1;
	if (scale < 0)
		scale = 0;
	if (scale > 0) {
		if (scale > 100000)
			math_error("Very large argument for exp");
		tmp = qscale(qs, -scale);
		qfree(qs);
		qs = tmp;
		tmp = qscale(epsilon, -scale);
		qfree(epsilon);
		epsilon = tmp;
	}
	epsilon2 = qscale(epsilon, -4L);
	bits = qprecision(epsilon) + 1;
	bits2 = bits + 10;
	qfree(epsilon);
	/*
	 * Now use the Taylor series expansion to calculate the exponential.
	 * Keep using approximations so that the fractions don't get too large.
	 */
	sum = qlink(&_qone_);
	term = qlink(&_qone_);
	n = 0;
	while (qrel(term, epsilon2) > 0) {
		n++;
		tmp = qmul(term, qs);
		qfree(term);
		term = qdivi(tmp, (long) n);
		qfree(tmp);
		tmp = qbround(term, bits2);
		qfree(term);
		term = tmp;
		tmp = qadd(sum, term);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	qfree(term);
	qfree(qs);
	qfree(epsilon2);
	/*
	 * Now repeatedly square the answer to get the final result.
	 * Then invert it if the original argument was negative.
	 */
	while (--scale >= 0) {
		tmp = qsquare(sum);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	tmp = qbround(sum, bits);
	qfree(sum);
	if (qisneg(q)) {
		sum = qinv(tmp);
		qfree(tmp);
		tmp = sum;
	}
	return tmp;
}


/*
 * Calculate the natural logarithm of a number accurate to the specified
 * epsilon.
 */
NUMBER *
qln(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *term, *term2, *sum, *epsilon2, *tmp1, *tmp2, *maxr;
	long shift, bits, bits2;
	int j, k;
	FULL n;
	BOOL neg;

	if (qisneg(q) || qiszero(q))
		math_error("log of non-positive number");
	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon for ln");
	if (qisone(q))
		return qlink(&_qzero_);
	/*
	 * If the number is less than one, invert it and remember that
	 * the result is to be negative.
	 */
	neg = FALSE;
	if (zrel(q->num, q->den) < 0) {
		neg = TRUE;
		q = qinv(q);
	} else
		q = qlink(q);
	j = 16;
	k = zhighbit(q->num) - zhighbit(q->den) + 1;
	while (k >>= 1)
		j++;
	epsilon2 = qscale(epsilon, -j);
	bits = qprecision(epsilon) + 1;
	bits2 = qprecision(epsilon2) + 5;
	/*
	 * By repeated square-roots scale number down to a value close
	 * to 1 so that Taylor series to be used will converge rapidly.
	 * The effect of scaling will be reversed by a later shift.
	 */
	maxr = iitoq(BASE + 1, BASE);
	shift = 1;
	while (qrel(q, maxr) > 0) {
		tmp1 = qsqrt(q, epsilon2);
		qfree(q);
		q = tmp1;
		shift++;
	}
	qfree(maxr);
	/*
	 * Calculate a value which will always converge using the formula:
	 *	ln((1+x)/(1-x)) = ln(1+x) - ln(1-x).
	 */
	tmp1 = qdec(q);
	tmp2 = qinc(q);
	term = qdiv(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	qfree(q);
	/*
	 * Now use the Taylor series expansion to calculate the result.
	 */
	n = 1;
	term2 = qsquare(term);
	sum = qlink(term);
	while (qrel(term, epsilon2) > 0) {
		n += 2;
		tmp1 = qmul(term, term2);
		qfree(term);
		term = qbround(tmp1, bits2);
		qfree(tmp1);
		tmp1 = qdivi(term, (long) n);
		tmp2 = qadd(sum, tmp1);
		qfree(tmp1);
		qfree(sum);
		sum = qbround(tmp2, bits2);
	}
	qfree(epsilon2);
	qfree(term);
	qfree(term2);
	/*
	 * Calculate the final result by multiplying by the proper power
	 * of two to undo the square roots done at the top, and possibly
	 * negating the result.
	 */
	tmp1 = qscale(sum, shift);
	qfree(sum);
	sum = qbround(tmp1, bits);
	qfree(tmp1);
	if (neg) {
		tmp1 = qneg(sum);
		qfree(sum);
		sum = tmp1;
	}
	return sum;
}


/*
 * Calculate the result of raising one number to the power of another.
 * The result is calculated to within the specified relative error.
 */
NUMBER *
qpower(q1, q2, epsilon)
	NUMBER *q1, *q2, *epsilon;
{
	NUMBER *tmp1, *tmp2, *epsilon2;

	if (qisint(q2))
		return qpowi(q1, q2);
	epsilon2 = qscale(epsilon, -4L);
	tmp1 = qln(q1, epsilon2);
	tmp2 = qmul(tmp1, q2);
	qfree(tmp1);
	tmp1 = qexp(tmp2, epsilon);
	qfree(tmp2);
	qfree(epsilon2);
	return tmp1;
}


/*
 * Calculate the Kth root of a number to within the specified accuracy.
 */
NUMBER *
qroot(q1, q2, epsilon)
	NUMBER *q1, *q2, *epsilon;
{
	NUMBER *tmp1, *tmp2, *epsilon2;
	int neg;

	if (qisneg(q2) || qiszero(q2) || qisfrac(q2))
		math_error("Taking bad root of number");
	if (qiszero(q1) || qisone(q1) || qisone(q2))
		return qlink(q1);
	if (qistwo(q2))
		return qsqrt(q1, epsilon);
	neg = qisneg(q1);
	if (neg) {
		if (ziseven(q2->num))
			math_error("Taking even root of negative number");
		q1 = qabs(q1);
	}
	epsilon2 = qscale(epsilon, -4L);
	tmp1 = qln(q1, epsilon2);
	tmp2 = qdiv(tmp1, q2);
	qfree(tmp1);
	tmp1 = qexp(tmp2, epsilon);
	qfree(tmp2);
	qfree(epsilon2);
	if (neg) {
		tmp2 = qneg(tmp1);
		qfree(tmp1);
		tmp1 = tmp2;
	}
	return tmp1;
}


/*
 * Calculate the hyperbolic cosine function with a relative accuracy less
 * than epsilon.  This is defined by:
 *	cosh(x) = (exp(x) + exp(-x)) / 2.
 */
NUMBER *
qcosh(q, epsilon)
	NUMBER *q, *epsilon;
{
	long scale;
	FULL n;
	FULL m;
	long bits, bits2;
	NUMBER *sum, *term, *qs, *epsilon2, *tmp;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for exp");
	if (qiszero(q))
		return qlink(&_qone_);
	epsilon = qscale(epsilon, -4L);
	/*
	 * If the argument is larger than one, then divide it by a power of two
	 * so that it is one or less.  This will make the series converge quickly.
	 * We will extrapolate the result for the original argument afterwards.
	 */
	qs = qabs(q);
	scale = zhighbit(q->num) - zhighbit(q->den) + 1;
	if (scale < 0)
		scale = 0;
	if (scale > 0) {
		if (scale > 100000)
			math_error("Very large argument for exp");
		tmp = qscale(qs, -scale);
		qfree(qs);
		qs = tmp;
		tmp = qscale(epsilon, -scale);
		qfree(epsilon);
		epsilon = tmp;
	}
	epsilon2 = qscale(epsilon, -4L);
	bits = qprecision(epsilon) + 1;
	bits2 = bits + 10;
	qfree(epsilon);
	tmp = qsquare(qs);
	qfree(qs);
	qs = tmp;
	/*
	 * Now use the Taylor series expansion to calculate the exponential.
	 * Keep using approximations so that the fractions don't get too large.
	 */
	sum = qlink(&_qone_);
	term = qlink(&_qone_);
	n = 0;
	while (qrel(term, epsilon2) > 0) {
		m = ++n;
		m *= ++n;
		tmp = qmul(term, qs);
		qfree(term);
		term = qdivi(tmp, (long) m);
		qfree(tmp);
		tmp = qbround(term, bits2);
		qfree(term);
		term = tmp;
		tmp = qadd(sum, term);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	qfree(term);
	qfree(qs);
	qfree(epsilon2);
	/*
	 * Now bring the number back up into range to get the final result.
	 * This uses the formula:
	 *	cosh(2 * x) = 2 * cosh(x)^2 - 1.
	 */
	while (--scale >= 0) {
		tmp = qsquare(sum);
		qfree(sum);
		sum = qscale(tmp, 1L);
		qfree(tmp);
		tmp = qdec(sum);
		qfree(sum);
		sum = qbround(tmp, bits2);
		qfree(tmp);
	}
	tmp = qbround(sum, bits);
	qfree(sum);
	return tmp;
}


/*
 * Calculate the hyperbolic sine with an accurary less than epsilon.
 * This is calculated using the formula:
 *	cosh(x)^2 - sinh(x)^2 = 1.
 */
NUMBER *
qsinh(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for sinh");
	if (qiszero(q))
		return qlink(q);
	epsilon = qscale(epsilon, -4L);
	tmp1 = qcosh(q, epsilon);
	tmp2 = qsquare(tmp1);
	qfree(tmp1);
	tmp1 = qdec(tmp2);
	qfree(tmp2);
	tmp2 = qsqrt(tmp1, epsilon);
	qfree(tmp1);
	if (qisneg(q)) {
		tmp1 = qneg(tmp2);
		qfree(tmp2);
		tmp2 = tmp1;
	}
	qfree(epsilon);
	return tmp2;
}


/*
 * Calculate the hyperbolic tangent with an accurary less than epsilon.
 * This is calculated using the formula:
 *	tanh(x) = sinh(x) / cosh(x).
 */
NUMBER *
qtanh(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2, *coshval;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for tanh");
	if (qiszero(q))
		return qlink(q);
	epsilon = qscale(epsilon, -4L);
	coshval = qcosh(q, epsilon);
	tmp2 = qsquare(coshval);
	tmp1 = qdec(tmp2);
	qfree(tmp2);
	tmp2 = qsqrt(tmp1, epsilon);
	qfree(tmp1);
	if (qisneg(q)) {
		tmp1 = qneg(tmp2);
		qfree(tmp2);
		tmp2 = tmp1;
	}
	qfree(epsilon);
	tmp1 = qdiv(tmp2, coshval);
	qfree(tmp2);
	qfree(coshval);
	return tmp1;
}


/*
 * Compute the hyperbolic arccosine within the specified accuracy.
 * This is calculated using the formula:
 *	acosh(x) = ln(x + sqrt(x^2 - 1)).
 */
NUMBER *
qacosh(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2, *epsilon2;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for acosh");
	if (qisone(q))
		return qlink(&_qzero_);
	if (qreli(q, 1L) < 0)
		math_error("Argument less than one for acosh");
	epsilon2 = qscale(epsilon, -8L);
	tmp1 = qsquare(q);
	tmp2 = qdec(tmp1);
	qfree(tmp1);
	tmp1 = qsqrt(tmp2, epsilon2);
	qfree(tmp2);
	tmp2 = qadd(tmp1, q);
	qfree(tmp1);
	tmp1 = qln(tmp2, epsilon);
	qfree(tmp2);
	qfree(epsilon2);
	return tmp1;
}


/*
 * Compute the hyperbolic arcsine within the specified accuracy.
 * This is calculated using the formula:
 *	asinh(x) = ln(x + sqrt(x^2 + 1)).
 */
NUMBER *
qasinh(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2, *epsilon2;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for asinh");
	if (qiszero(q))
		return qlink(&_qzero_);
	epsilon2 = qscale(epsilon, -8L);
	tmp1 = qsquare(q);
	tmp2 = qinc(tmp1);
	qfree(tmp1);
	tmp1 = qsqrt(tmp2, epsilon2);
	qfree(tmp2);
	tmp2 = qadd(tmp1, q);
	qfree(tmp1);
	tmp1 = qln(tmp2, epsilon);
	qfree(tmp2);
	qfree(epsilon2);
	return tmp1;
}


/*
 * Compute the hyperbolic arctangent within the specified accuracy.
 * This is calculated using the formula:
 *	atanh(x) = ln((1 + u) / (1 - u)) / 2.
 */
NUMBER *
qatanh(q, epsilon)
	NUMBER *q, *epsilon;
{
	NUMBER *tmp1, *tmp2, *tmp3;

	if (qisneg(epsilon) || qiszero(epsilon))
		math_error("Illegal epsilon value for atanh");
	if (qiszero(q))
		return qlink(&_qzero_);
	if ((qreli(q, 1L) > 0) || (qreli(q, -1L) < 0))
		math_error("Argument not between -1 and 1 for atanh");
	tmp1 = qinc(q);
	tmp2 = qsub(&_qone_, q);
	tmp3 = qdiv(tmp1, tmp2);
	qfree(tmp1);
	qfree(tmp2);
	tmp1 = qln(tmp3, epsilon);
	qfree(tmp3);
	tmp2 = qscale(tmp1, -1L);
	qfree(tmp1);
	return tmp2;
}

/* END CODE */