1/*
2 * pi - various routines to calculate pi
3 *
4 * Copyright (C) 1999-2004,2021  David I. Bell
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	1991/05/22 21:56:37
21 * File existed as early as:	1991
22 *
23 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
24 */
25
26/*
27 * Calculate pi within the specified epsilon using the quartic convergence
28 * iteration.
29 */
30
31
32define qpi(epsilon)
33{
34	local niter, yn, ym, tm, an, am, t, tn, sqrt2, epsilon2, count, digits;
35	local bits, bits2;
36
37	if (isnull(epsilon))
38		epsilon = epsilon();
39	digits = digits(1/epsilon);
40	if	(digits <=  8) { niter = 1; epsilon =	1e-8; }
41	else if (digits <= 40) { niter = 2; epsilon =  1e-40; }
42	else if (digits <= 170) { niter = 3; epsilon = 1e-170; }
43	else if (digits <= 693) { niter = 4; epsilon = 1e-693; }
44	else {
45		niter = 4;
46		t = 693;
47		while (t < digits) {
48			++niter;
49			t *= 4;
50		}
51	}
52	epsilon2 = epsilon/(digits/10 + 1);
53	digits = digits(1/epsilon2);
54	sqrt2 = sqrt(2, epsilon2);
55	bits = abs(ilog2(epsilon)) + 1;
56	bits2 = abs(ilog2(epsilon2)) + 1;
57	yn = sqrt2 - 1;
58	an = 6 - 4 * sqrt2;
59	tn = 2;
60	for (count = 0; count < niter; ++count) {
61		ym = yn;
62		am = an;
63		tn *= 4;
64		t = sqrt(sqrt(1-ym^4, epsilon2), epsilon2);
65		yn = (1-t)/(1+t);
66		an = (1+yn)^4*am-tn*yn*(1+yn+yn^2);
67		yn = bround(yn, bits2);
68		an = bround(an, bits2);
69	}
70	return (bround(1/an, bits));
71}
72
73
74/*
75 * Print digits of PI forever, neatly formatted, using calc.
76 *
77 * Written by Klaus Alexander Seistrup <klaus at seistrup dot dk>
78 * on a dull Friday evening in November 1999.
79 *
80 * Inspired by an algorithm conceived by Lambert Meertens.
81 *
82 * See also the ABC Programmer's Handbook, by Geurts, Meertens & Pemberton,
83 * published by Prentice-Hall (UK) Ltd., 1990.
84 *
85 */
86
87define piforever()
88{
89	local k = 2;
90	local a = 4;
91	local b = 1;
92	local a1 = 12;
93	local b1 = 4;
94	local a2, b2, p, q, d, d1;
95	local stdout = files(1);
96	local first = 1, row = -1, col = 0;
97
98	while (1) {
99		/*
100		* Next approximation
101		*/
102		p = k * k;
103		q = k + ++k;
104
105		a2 = a;
106		b2 = b;
107
108		a = a1;
109		a1 = p * a2 + q * a1;
110		b = b1;
111		b1 = p * b2 + q * b1;
112
113		/*
114		* Print common digits
115		*/
116		d = a // b;
117		d1 = a1 // b1;
118
119		while (d == d1) {
120			if (first) {
121				printf("%d.", d);
122				first = 0;
123			} else {
124				if (!(col % 50)) {
125					printf("\n");
126					col = 0;
127					if (!(++row % 20)) {
128						printf("\n");
129						row = 0;
130					}
131				}
132				printf("%d", d);
133				if (!(++col % 10))
134					printf(" ");
135			}
136			a = 10 * (a % b);
137			a1 = 10 * (a1 % b1);
138			d = a // b;
139			d1 = a1 // b1;
140		}
141		fflush(stdout);
142	}
143}
144