Path: newsspool2.news.pas.earthlink.net!stamper.news.pas.earthlink.net!elnk-nf2-pas!newsfeed.earthlink.net!newshub.sdsu.edu!logbridge.uoregon.edu!xmission!news.cc.utah.edu!newsfeed.cs.utexas.edu!news.cs.utexas.edu!not-for-mail
Message-ID: <bjat8b$bca$1@zorkmid.cs.utexas.edu>
From: djimenez@cs.utexas.edu (Daniel A. Jimenez)
Newsgroups: sci.math
Subject: Re: An Informal Call For Program Submissions To Compute pi(n)
Date: 5 Sep 2003 15:55:39 -0500
References: <v7ml21-ujr.ln1@lexi2.athghost7038suus.net>
Lines: 90
Organization: CS Dept, University of Texas at Austin
Sender: djimenez@zorkmid.cs.utexas.edu
NNTP-Posting-Host: zorkmid.cs.utexas.edu
X-Trace: news.cs.utexas.edu 1062795341 7723 128.83.120.248 (5 Sep 2003 20:55:41 GMT)
X-Complaints-To: usenet@cs.utexas.edu
NNTP-Posting-Date: Fri, 5 Sep 2003 20:55:41 +0000 (UTC)
X-Received-Date: Fri, 05 Sep 2003 13:55:46 PDT (newsspool2.news.pas.earthlink.net)
Xref: lexi2.athghost7038suus.net sci.math:36015

In article <v7ml21-ujr.ln1@lexi2.athghost7038suus.net>,
The Ghost In The Machine  <ewill@sirius.athghost7038suus.net> wrote:
>I seem to be collecting these things anyway, so I might as well
>throw out a call for slightly more formal submissions.
>
>[snippage]
>
>[7] Submissions will be picked up from Usenet posts in sci.math.
>    Ideally they'd be in this thread. :-)  (For those who have
>    posted code snippets that I've replied to in other threads,
>    thank you. :-) )

Here you are.  Enjoy.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

unsigned find (unsigned b, unsigned i, unsigned a, unsigned c) {
	unsigned j, k, n;
	n = 6 * i + c;
	b *= 6;
	b += (n-(b%n));
	for (j=b; ; j+=n) {
		k = (j-a)/6;
		if (k * 6 + a == j) if (k != i) return k;
	}
}

unsigned long long pi (unsigned long long v) {
	char A0[65536], A1[65536];
	unsigned b, s, r, i, j, k, l, w, f1, f2, n;
	char	*p0, *p1;

	n = 0;
	w = 1 + v/6;
	if (v >= 2) n++;
	if (v >= 3) n++;
	if (v >= 5) n++;
	if (v >= 7) n++;
	if (v >= 11) n++;
	for (b=0; b<=w; b+=65536) {
		r = b + 65536;
		if (r > w) r = w;
		l = 1 + (unsigned) sqrt (r);
		memset (A0, 1, r-b);
		memset (A1, 1, r-b);
		for (i=0; i<l; i++) {
			f1 = i * 6 + 1;
			f2 = i * 6 + 5;
			if (f1 != 1) {
			for (j=find(b,i,1,1); j<r; j+=f1) A0[j-b] = false;
			for (j=find(b,i,5,1); j<r; j+=f1) A1[j-b] = false;
			}
			for (j=find(b,i,5,5); j<r; j+=f2) A1[j-b] = false;
			for (j=find(b,i,1,5); j<r; j+=f2) A0[j-b] = false;
		}
		s = b;
		p0 = &A0[0];
		p1 = &A1[0];
		if (b == 0) {
			s = 2;
			p0 += 2;
			p1 += 2;
		}
		for (i=s; i<r; i++) {
			k = i * 6 + 1;
			if (k > v) break;
			n += *p0++;
			k = i * 6 + 5;
			if (k > v) break;
			n += *p1++;
		}
	}
	return n;
}

int main (int argc, char *argv[]) {
	unsigned n;
	if (argc != 2) exit (1);
	sscanf (argv[1], "%u", &n);
	printf ("%llu\n", pi (n));
	exit (0);
}
-- 
Daniel A. Jiménez                 djimenez@cs.utexas.edu
Assistant Professor               djimenez@cs.rutgers.edu
Department of Computer Science
Rutgers University