more fun with pi

[Akira Rabelais <akirarabelais@xxxxxxxxxxx>]

> counts off the numbers of pi.


you might like this:

// piextraction.c

#include <math.h>
#include "obviate.h"
#include "Utilitaire.h"
#include "piextraction.h"

#define NHX 16
#define eps 1e-17

/*
This function employs the recently discovered digit piextraction scheme
to produce hex digits of .  This code is valid up to ic = 2^24 on
systems with IEEE arithmetic.	(David H. Bailey     960429)
*/

void piextraction(void)
{
	double pid, s1, s2, s3, s4;
	int ic = 1000000; char chx[NHX];

/*  ic is the hex digit position -- output begins at position ic + 1. */

	s1 = series(1,ic); s2 = series(4,ic);
	s3 = series(5,ic); s4 = series(6,ic);
	pid = 4. * s1 - 2. * s2 - s3 - s4;
	pid = pid - (int) pid + 1.;
	ihex (pid,NHX,chx);

 printf ("Pi hex digit computation\n");
 printf ("position = %i + 1\n %20.15f\n %12.12s\n", ic, pid, chx);

}

/*
This returns, in chx, the first nhx hex
digits of the fraction of x.
*/

void ihex(double x,int nhx,char chx[])
{
	int i; double y;
	char hx[] = "0123456789ABCDEF";

	y = fabs (x);
	for (i = 0;i < nhx;i++){
		y = 16. * (y - floor(y));
		chx[i] = hx[(int)y];
	}
}

/*
This routine evaluates the series  sum_k 16^(ic-k)/(8*k+m)
using the modular exponentiation technique.
*/

double series(int m,int ic)
{
	int k; double ak,p,s,t;
	double expm(double x,double y);

 	s = 0.;

/*  Sum the series up to ic. */
	for(k = 0; k < ic;k++){
		ak = 8 * k + m;
		p = ic - k;
		t = expm (p, ak);
		s = s + t / ak;
		s = s - (int) s;
	}

/*  Compute a few terms where k >= ic. */
	for(k = ic; k <= ic + 100;k++){
		ak = 8 * k + m;
		t = pow (16.,(double)(ic - k))/ak;
		if (t < eps) break;
		s = s + t; s = s - (int) s;
	}
	return s;
}


/*
expm = 16^p mod ak.  This routine uses the left-to-right binary
exponentiation scheme.  It is valid for  ak <= 2^24.
*/

double expm (double p,double ak)
{
	int i, j; double p1, pt, r;
	#define ntp 25
	static double tp[ntp];
	static int tp1 = 0;

/*  If this is the first call to expm,
fill the power of two table tp. */
	if(tp1 == 0) {
		tp1 = 1; tp[0] = 1.;
		for (i = 1;i < ntp;i++)
			tp[i] = 2. * tp[i-1];
	}
	if(ak == 1.) return 0.;

/*  Find the greatest power of two less than or equal to p. */
	for(i = 0;i < ntp;i++) if(tp[i] > p) break;
	pt = tp[i-1]; p1 = p; r = 1.;

/*  Perform binary exponentiation algorithm modulo ak. */
	for (j = 1;j <= i;j++){
		if (p1 >= pt){
			r = 16. * r;
			r = r - (int)(r/ak)*ak;
			p1 = p1 - pt;
		}
		pt = 0.5 * pt;
			if (pt >= 1.){
			r = r * r;
			r = r - (int) (r / ak) * ak;
		}
	}
	return r;
}





// eof



              Akira


_________________________________________________________________
Get your FREE download of MSN Explorer at http://explorer.msn.com