Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Posting-Version: notesnews 0.1 (unido 12/05/84); site unido.UUCP
Path: utzoo!linus!philabs!cmcl2!seismo!mcvax!unido!kl
From: kl@unido.UUCP
Newsgroups: net.sources
Subject: Re: Solving Pi
Message-ID: <11800005@unido.UUCP>
Date: Mon, 29-Jul-85 14:18:00 EDT
Article-I.D.: unido.11800005
Posted: Mon Jul 29 14:18:00 1985
Date-Received: Thu, 1-Aug-85 05:05:00 EDT
References: <187@ski.UUCP>
Sender: notes@unido.UUCP
Lines: 327
Nf-ID: #R:ski:-18700:unido:11800005:000:6551
Nf-From: unido!kl    Jul 29 16:18:00 1985


# This is another program to calculate a few decimals of pi, using the
# "standard" series for computing. There are two versions given:
# pi, which uses a normal C-program (div.c) to do multi-precision division,
# and pix, which uses an assembler routine (exdiv.s) instead. Exdiv.s contains
# the special "extended division" opcode for a VAX, thereby allowing twice
# the number of digits per integer (see def.h).
# To use:
#   1. Unpack this file
#   2. Edit def.h
#           It should only be necessary to define the number of DIGITS
#           to compute, or the default LOG file to use.
#   3. If you don't have rusage-stuff (4.2BSD), edit status() in pi.c
#   4. do make
#   5. run pi or pix, an argument is the logfile to use.
#
# Note: a status() is given at the beginning and end of execution.
#       Since pi (pix) will probably be a long running program :-), a
#       HUP (#1) signal sent to pi (pix) will cause a status entry in
#       the LOG file. For really long running progs, see ./clock.
#
# Some sample (cpu) times (on a VAX/750):
#                 pix                     pi
# DIGITS
# 1000            about 12 sec.           about 25 sec.
# 10000           about 20 min.           about 44 min.
# 100050          about 35 hrs.           about 78 hrs.
#
# This is a shell archive.  Unpack with "sh <file".
echo x - Makefile
cat > "Makefile" << '!Funky!Stuff!'
CFLAGS = -O

all:    pi pix print

pi:     def.h pi.o div.o print
	cc $(CFLAGS) pi.o div.o -o pi -lm

pix:    def.h pix.o exdiv.o print
	cc $(CFLAGS) pix.o exdiv.o -o pix -lm

pix.o:  pi.c def.h
	cc -DEXDIV -S pi.c
	as -o pix.o pi.s
	rm pi.s

exdiv.o: exdiv.s def.h
	/lib/cpp exdiv.s tmp.s -DEXDIV
	as -o exdiv.o tmp.s
	rm tmp.s

print:  print.o
	cc $(CFLAGS) -o print print.o

clean:  
	-rm *.o core pi pix print
!Funky!Stuff!
echo x - clock
cat > "clock" << '!Funky!Stuff!'
if [ x$1 = x ]
then
	echo "Usage: $0 <pid of pi-program>"
	exit
fi

while true
do
	kill -1 $1
	sleep 3600
done
!Funky!Stuff!
echo x - def.h
cat > "def.h" << '!Funky!Stuff!'
#define DIGITS 10000        /* no. of digits of pi we wish to compute */
#define LOG     "log"       /* the logfile */

#ifdef EXDIV
#  define SIZE 100000000    /* max. power of ten <= maxint */
#  define NSIZE 8           /* no. of digits in SIZE, ie 10^NSIZE = SIZE */
#else
#  define SIZE 10000        /* max. power of ten <= sqrt(maxint) */
#  define NSIZE 4           /* no. of digits in SIZE, ie 10^NSIZE = SIZE */
#endif EXDIV

#define N (DIGITS/NSIZE+2)  /* +2 for rounding */
!Funky!Stuff!
echo x - div.c
cat > "div.c" << '!Funky!Stuff!'
#include "def.h"

div (a,n)
int *a, n;
{
	register int *ra, *ul, q;
	ul = &a[N];
	q = 0;
	for (ra=a; ra<ul;) {
		q = q*SIZE + *ra; 
		*ra++ = q/n;
		q = q % n;
	}
}
!Funky!Stuff!
echo x - exdiv.s
cat > "exdiv.s" << '!Funky!Stuff!'
#include "def.h"
LL0:
	.data
	.comm   _quad,8
	.text
	.align  1
	.globl  _div
_div:
	.word   L12
	jbr     L14
L15:
	clrq    _quad
	clrl    r11
L18:
	cmpl    r11,$N
	jgeq    L17
	emul    $SIZE,_quad,*4(ap)[r11],_quad
	ediv    8(ap),_quad,*4(ap)[r11],_quad
	clrl    _quad+4
L16:
	incl    r11
	jbr     L18
L17:
	ret
	.set    L12,0xc00
L14:
	jbr     L15
	.data
!Funky!Stuff!
echo x - pi.c
cat > "pi.c" << '!Funky!Stuff!'
#include "def.h"

add (a,b)
int *a, *b;
{ 
	register int *ap, *bp;

	for (ap = &a[N], bp = &b[N]; ap>a; ) {
		*--bp += *--ap;
		if (*bp > SIZE-1) { 
			*bp -= SIZE;
			bp[-1]++;
		}
	}
}

sub (a,b)
int *a, *b;
{ 
	register int *ap, *bp;

	for (ap = &a[N], bp = &b[N]; ap>a; ) {
		*--bp -= *--ap;
		if (*bp < 0) {
			*bp += SIZE;
			bp[-1]--;
		}
	}
}

copy(x,y)
int *x, *y;
{ 
	register int *rx, *ry;
	for (rx = &x[N], ry = &y[N]; rx>x; )
		*--ry = *--rx;
}

#include <sys/time.h>
#include <sys/resource.h>
#include <stdio.h>
#define PRINT(x)        fprintf(Log,"x: %d\t", ru.x)

char    *ctime();
struct rusage   ru;
FILE    *Log;
int     cnt5, cnt239;
int     lim5, lim239;

status()
{
	long now;

	time(&now);
	fprintf (Log,"at %s", ctime(&now));
	fprintf(Log,"cnt5= %d(%d)\tcnt239= %d(%d)\n",4*cnt5,lim5,
						    4*cnt239,lim239);

	getrusage(RUSAGE_SELF,&ru);
	fprintf(Log,"usr time: %d.%06d [sec]\n", ru.ru_utime.tv_sec,
						 ru.ru_utime.tv_usec);
	fprintf(Log,"sys time: %d.%06d [sec]\n", ru.ru_stime.tv_sec,
						 ru.ru_stime.tv_usec);

	/* edit these PRINT's according to taste */

	PRINT(ru_maxrss);       PRINT(ru_minflt);       PRINT(ru_majflt);
	fprintf(Log,"\n"); 
	PRINT(ru_nswap);        PRINT(ru_inblock);      PRINT(ru_oublock);     
	fprintf(Log,"\n");
	PRINT(ru_nsignals);     PRINT(ru_nvcsw);        PRINT(ru_nivcsw);
	fprintf(Log,"\n\n");
	fflush(Log);
};

#include <math.h>
#include <signal.h>

FILE    *popen();

main(argc, argv)
int argc;
char *argv[];
{
	int     term[N], pterm[N], sum[N];
	register int i, j;
	char    aux[20];
	FILE    *Print;

	cnt5 = cnt239 = 0;
	signal(SIGHUP,status);

	lim5   = (int)((double)DIGITS/log10(5.0)) + 10; /* +10 for safety */
	lim239 = (int)((double)DIGITS/log10(239.0)) + 10;

	if (argv[1] != NULL)
		Log = fopen(argv[1],"w");
	else
		Log = fopen(LOG,"w");
	if (Log == NULL) {
		perror("fopen");
		exit(1);
	}

	fprintf (Log,"DIGITS = %d\nlim5   = %d\nlim239 = %d\n",
				DIGITS,lim5,lim239);
	fprintf (Log,"start of exec. ");
	status();

	term[0] = 32*(SIZE/100);
	for (j=1;j<N;j++)
		term[j] = 0;
	copy (term,sum);

	i = 1;
	while (i < lim5) {
		i += 2;
		div(term,25);
		copy (term,pterm); 
		div(pterm,i);
		sub(pterm,sum);
		i += 2;
		div(term,25); 
		copy(term,pterm); 
		div(pterm,i); 
		add(pterm,sum);
		cnt5++;
	}

	term[0] = 4*(SIZE/10);
	for (j=1;j<N;j++) term[j] = 0;
	div(term,239);
	sub(term,sum);

	i = 1;
	while (i < lim239) {
		i += 2;
		div(term,57121);
		copy (term,pterm); 
		div(pterm,i);
		add(pterm,sum);
		i += 2;
		div(term,57121); 
		copy(term,pterm); 
		div(pterm,i); 
		sub(pterm,sum);
		cnt239++;
	}

	Print = popen("print","w");
	for (j=0;j<N-2;j++) {
		sprintf(aux,"%0*d",NSIZE,sum[j]);
		fprintf(Print,"%s",aux);
	}
	putchar('\n');
	fclose(Print);

	fprintf (Log,"end of exec. ");
	status();
	fclose(Log);
}
!Funky!Stuff!
echo x - print.c
cat > "print.c" << '!Funky!Stuff!'
#include <stdio.h>
main() {
	char c;
	int i;

	i = 0;
	putchar('\n');
	c = getchar();
	putchar (c);
	putchar(',');
	putchar('\n');
	while ( (c=getchar()) != EOF) {
			putchar(c); 
			i++; 
			if (i%5 == 0 ) putchar(' ');
			if (i%50 == 0) putchar('\n');
	}
	putchar('\n');
}
!Funky!Stuff!
#
#
# Karlheinz Lehner                kl@unido.UUCP
# Univ. of Dortmund               hpfcla!hpbbn!unido!kl
# West Germany                    mcvax!unido!kl
#