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