c_cpp/etc/calc/cal/pix.cal

   1 /*
   2  * pix - iterative method of finding the number of primes less than x
   3  *
   4  * Copyright (C) 1999  Landon Curt Noll
   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: pix.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
  22  * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/pix.cal,v $
  23  *
  24  * Under source code control:   1996/07/09 03:14:14
  25  * File existed as early as:    1996
  26  *
  27  * chongo <was here> /\oo/\     http://www.isthe.com/chongo/
  28  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  29  */
  30
  31 /*
  32  * Here is an iterative method of finding the number of primes less than
  33  * or equal to a given number.  This method is from "Computer Recreations"
  34  * June 1996 issue of Scientific American.
  35  *
  36  * NOTE: For reasonable values of x, the builtin function pix(x) is
  37  *       much faster.  This code is provided because the method
  38  *       is interesting.
  39  */
  40
  41
  42 define pi_of_x(x)
  43 {
  44         local An;       /* A(n) */
  45         local An1;      /* A(n-1) */
  46         local An2;      /* A(n-2) */
  47         local An3;      /* A(n-3) */
  48         local primes;   /* number of primes found */
  49         local n;        /* loop counter */
  50
  51         /*
  52          * setup
  53          */
  54         An1 = 2;
  55         An2 = 0;
  56         An3 = 3;
  57         primes = 1;
  58
  59         /*
  60          * main A(n+1)=A(n-1)+A(n-2) sequence loop
  61          */
  62         for (n = 3; n < x; ++n) {
  63                 An = An2 + An3;
  64                 An3 = An2;
  65                 An2 = An1;
  66                 An1 = An;
  67                 if (An % n == 0)
  68                         ++primes;
  69         }
  70         return primes;
  71 }