Solving Pi
kl at unido.UUCP
kl at unido.UUCP
Tue Jul 30 04:18:00 AEST 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 at unido.UUCP
# Univ. of Dortmund hpfcla!hpbbn!unido!kl
# West Germany mcvax!unido!kl
#
More information about the Comp.sources.unix
mailing list