c_cpp/etc/calc/cal/pi.cal

   1 /*
   2  * pi - various routines to calculate pi
   3  *
   4  * Copyright (C) 1999-2004  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  * @(#) $Revision: 30.1 $
  21  * @(#) $Id: pi.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
  22  * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/pi.cal,v $
  23  *
  24  * Under source code control:   1991/05/22 21:56:37
  25  * File existed as early as:    1991
  26  *
  27  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  28  */
  29
  30 /*
  31  * Calculate pi within the specified epsilon using the quartic convergence
  32  * iteration.
  33  */
  34
  35
  36 define qpi(epsilon)
  37 {
  38         local niter, yn, ym, tm, an, am, t, tn, sqrt2, epsilon2, count, digits;
  39         local bits, bits2;
  40
  41         if (isnull(epsilon))
  42                 epsilon = epsilon();
  43         digits = digits(1/epsilon);
  44         if      (digits <=  8) { niter = 1; epsilon =   1e-8; }
  45         else if (digits <= 40) { niter = 2; epsilon =  1e-40; }
  46         else if (digits <= 170) { niter = 3; epsilon = 1e-170; }
  47         else if (digits <= 693) { niter = 4; epsilon = 1e-693; }
  48         else {
  49                 niter = 4;
  50                 t = 693;
  51                 while (t < digits) {
  52                         ++niter;
  53                         t *= 4;
  54                 }
  55         }
  56         epsilon2 = epsilon/(digits/10 + 1);
  57         digits = digits(1/epsilon2);
  58         sqrt2 = sqrt(2, epsilon2);
  59         bits = abs(ilog2(epsilon)) + 1;
  60         bits2 = abs(ilog2(epsilon2)) + 1;
  61         yn = sqrt2 - 1;
  62         an = 6 - 4 * sqrt2;
  63         tn = 2;
  64         for (count = 0; count < niter; ++count) {
  65                 ym = yn;
  66                 am = an;
  67                 tn *= 4;
  68                 t = sqrt(sqrt(1-ym^4, epsilon2), epsilon2);
  69                 yn = (1-t)/(1+t);
  70                 an = (1+yn)^4*am-tn*yn*(1+yn+yn^2);
  71                 yn = bround(yn, bits2);
  72                 an = bround(an, bits2);
  73         }
  74         return (bround(1/an, bits));
  75 }
  76
  77
  78 /*
  79  * Print digits of PI forever, neatly formatted, using calc.
  80  *
  81  * Written by Klaus Alexander Seistrup <klaus@seistrup.dk>
  82  * on a dull Friday evening in November 1999.
  83  *
  84  * Inspired by an algorithm conceived by Lambert Meertens.
  85  *
  86  * See also the ABC Programmer's Handbook, by Geurts, Meertens & Pemberton,
  87  * published by Prentice-Hall (UK) Ltd., 1990.
  88  *
  89  */
  90
  91 define piforever()
  92 {
  93         local k = 2;
  94         local a = 4;
  95         local b = 1;
  96         local a1 = 12;
  97         local b1 = 4;
  98         local a2, b2, p, q, d, d1;
  99         local stdout = files(1);
 100         local first = 1, row = -1, col = 0;
 101
 102         while (1) {
 103                 /*
 104                 * Next approximation
 105                 */
 106                 p = k * k;
 107                 q = k + ++k;
 108
 109                 a2 = a;
 110                 b2 = b;
 111
 112                 a = a1;
 113                 a1 = p * a2 + q * a1;
 114                 b = b1;
 115                 b1 = p * b2 + q * b1;
 116
 117                 /*
 118                 * Print common digits
 119                 */
 120                 d = a // b;
 121                 d1 = a1 // b1;
 122
 123                 while (d == d1) {
 124                         if (first) {
 125                                 printf("%d.", d);
 126                                 first = 0;
 127                         } else {
 128                                 if (!(col % 50)) {
 129                                         printf("\n");
 130                                         col = 0;
 131                                         if (!(++row % 20)) {
 132                                                 printf("\n");
 133                                                 row = 0;
 134                                         }
 135                                 }
 136                                 printf("%d", d);
 137                                 if (!(++col % 10))
 138                                         printf(" ");
 139                         }
 140                         a = 10 * (a % b);
 141                         a1 = 10 * (a1 % b1);
 142                         d = a // b;
 143                         d1 = a1 // b1;
 144                 }
 145                 fflush(stdout);
 146         }
 147 }