02-28-2009, 05:13 PM
This question was originally ask by Gene in the following post: http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv014.cgi?read=53988:
Quote:I was not a member of this fine forum at the time of this post. But its never too late (just under 5 years late).
Gene: This is a write-up I just received from Palmer Hanson, who edited the TI PPC Notes for many years in the 1980s - early 1990s. Thought it might be interesting. Looks like the race to 1000 digits on a machine prior to 1990 now belongs to TI? Can we beat the 5.81 hours on the TI95 with the HP71B? :-)
I wrote a spigot program in FORTH and ran it on my 71B this morning. I selected the spigot algorithm for the following reasons:
- Spigot programs are very small. Mine is 2x larger because I unrolled the first outer loop. This avoids a conditional statement in the inner loop. Consider the initial loop a method to initialize the array A.
- This spigot algorithm requires integers only. A benefit for FORTH.
- Spigot programs allow you to watch the action. No need to wait for the final result.
Output:
3141592653589793238462643383279502884197169399375105820974944592307816406286
2089986280348253421170679821480865132823066470938446095505822317253594081284
8111745028410270193852110555964462294895493038196442881097566593344612847564
8233786783165271201909145648566923460348610454326648213393607260249141273724
5870066063155881748815209209628292540917153643678925903600113305305488204665
2138414695194151160943305727036575959195309218611738193261179310511854807446
2379962749567351885752724891227938183011949129833673362440656643086021394946
3952247371907021798609437027705392171762931767523846748184676694051320005681
2714526356082778577134275778960917363717872146844090122495343014654958537105
0792279689258923542019956112129021960864034418159813629774771309960518707211
3499999983729780499510597317328160963185950244594553469083026425223082533446
8503526193118817101000313783875288658753320838142061717766914730359825349042
8755468731159562863882353787593751957781857780532171226806613001927876611195
909216420198
TIME: 11362.03 SEC1000 digits in 3.16 hours. Almost 2x faster than the TI95. Over 2x faster than BASIC versions in the aforementioned thread.
Code (my first attempt, more optimization can be done):
: SPIGOT ;
: 04D. <# # # # # #> TYPE ;
: DUM* DUP ROT ROT UM* SWAP 2SWAP UM* D+ ;
1000 CONSTANT N
N 4 / 14 * CONSTANT M
5 CONSTANT S
CREATE A M S * NALLOT
VARIABLE B
VARIABLE C
VARIABLE E
10000 CONSTANT F
VARIABLE K
76 4 / CONSTANT MAXL
M DUP B ! C ! 1 K !
: SPIGOTM
0,
BEGIN B @ 1 - DUP B ! 0 > WHILE
B @ DUM*
20000000,
D+
B @ 2* 1 -
M/MOD
ROT
A B @ S * + !
REPEAT
F M/MOD 04D.
DUP E ! 0
BEGIN C @ 14 - DUP DUP C ! B ! 0 > WHILE
BEGIN B @ 1 - DUP B ! 0 > WHILE
B @ DUM*
A B @ S * + @ F UM*
D+
B @ 2* 1 -
M/MOD
ROT
A B @ S * + !
REPEAT
F M/MOD E @ 0 D+ 04D.
DUP E ! 0
1 K +! K @ MAXL = IF CR 0 K ! THEN
REPEAT
;
: DOIT
CLOCK
SPIGOTM 2DROP CR CR
." TIME: " CLOCK X<>Y F- STD F. ." SEC" CR
;
Original C program:
#include <stdio.h>
#define NDIGITS 1000
#define LEN (NDIGITS/4+1)*14
int main()
{
int a[LEN];
int b;
int c = LEN;
int d = 0;
int e = 0;
int f = 10000;
int g;
int h = 0;
while ((b = c -= 14) > 0) {
while (--b > 0) {
d *= b;
if (h == 0)
d += 2000 * f;
else
d += a[b] * f;
g = b + b - 1;
a[b] = d % g;
d /= g;
}
printf("%04d", e + d / f);
d = e = d % f;
h = 4;
}
return (0);
}
More spigot info:
http://sense.net/~egan/hpgcc/#Spigot%20Algorithm.
Edited: 28 Feb 2009, 5:35 p.m.