pi.cal

Calc Software Package for Number Calc

/*
 * pi - various routines to calculate pi
 *
 * Copyright (C) 1999-2004  David I. Bell
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.5 $
 * @(#) $Id: pi.cal,v 29.5 2004/02/23 14:04:01 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/pi.cal,v $
 *
 * Under source code control:	1991/05/22 21:56:37
 * File existed as early as:	1991
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * Calculate pi within the specified epsilon using the quartic convergence
 * iteration.
 */


define qpi(epsilon)
{
	local niter, yn, ym, tm, an, am, t, tn, sqrt2, epsilon2, count, digits;
	local bits, bits2;

	if (isnull(epsilon))
		epsilon = epsilon();
	digits = digits(1/epsilon);
	if	(digits <=  8) { niter = 1; epsilon =	1e-8; }
	else if (digits <= 40) { niter = 2; epsilon =  1e-40; }
	else if (digits <= 170) { niter = 3; epsilon = 1e-170; }
	else if (digits <= 693) { niter = 4; epsilon = 1e-693; }
	else {
		niter = 4;
		t = 693;
		while (t < digits) {
			++niter;
			t *= 4;
		}
	}
	epsilon2 = epsilon/(digits/10 + 1);
	digits = digits(1/epsilon2);
	sqrt2 = sqrt(2, epsilon2);
	bits = abs(ilog2(epsilon)) + 1;
	bits2 = abs(ilog2(epsilon2)) + 1;
	yn = sqrt2 - 1;
	an = 6 - 4 * sqrt2;
	tn = 2;
	for (count = 0; count < niter; ++count) {
		ym = yn;
		am = an;
		tn *= 4;
		t = sqrt(sqrt(1-ym^4, epsilon2), epsilon2);
		yn = (1-t)/(1+t);
		an = (1+yn)^4*am-tn*yn*(1+yn+yn^2);
		yn = bround(yn, bits2);
		an = bround(an, bits2);
	}
	return (bround(1/an, bits));
}


/*
 * Print digits of PI forever, neatly formatted, using calc.
 *
 * Written by Klaus Alexander Seistrup <[e-mail]>
 * on a dull Friday evening in November 1999.
 *
 * Inspired by an algorithm conceived by Lambert Meertens.
 *
 * See also the ABC Programmer's Handbook, by Geurts, Meertens & Pemberton,
 * published by Prentice-Hall (UK) Ltd., 1990.
 *
 */

define piforever()
{
	local k = 2;
	local a = 4;
	local b = 1;
	local a1 = 12;
	local b1 = 4;
	local a2, b2, p, q, d, d1;
	local stdout = files(1);
	local first = 1, row = -1, col = 0;

	while (1) {
		/*
		* Next approximation
		*/
		p = k * k;
		q = k + ++k;

		a2 = a;
		b2 = b;

		a = a1;
		a1 = p * a2 + q * a1;
		b = b1;
		b1 = p * b2 + q * b1;

		/*
		* Print common digits
		*/
		d = a // b;
		d1 = a1 // b1;

		while (d == d1) {
			if (first) {
				printf("%d.", d);
				first = 0;
			} else {
				if (!(col % 50)) {
					printf("\n");
					col = 0;
					if (!(++row % 20)) {
						printf("\n");
						row = 0;
					}
				}
				printf("%d", d);
				if (!(++col % 10))
					printf(" ");
			}
			a = 10 * (a % b);
			a1 = 10 * (a1 % b1);
			d = a // b;
			d1 = a1 // b1;
		}
		fflush(stdout);
	}
}