I recently ran across an article about calculating the mathematical constant e (the base of natural logarithms) to 116,000 places in the June 1981 issue of Byte magazine, written by some character named Wozniak. (It seems unlikely that he is notable for anything else 6502 related.)

https://archive.org/stream/byte-magazin ... 3/mode/2up

(The Wikipedia page for e also references that article.)

https://en.wikipedia.org/wiki/E_(mathematical_constant)

In the Byte article are program listings for the Apple II; two 6502 assembly programs, and an "Integer" BASIC program. It calculates e by fixed-point arithmetic using the well-known (and more importantly, quickly converging) series:



(Actually, only fractional part is calculated, i.e. excluding the terms to the left of the 1/2! term.) It is an interesting article. A naive implementation would use (roughly) half of available RAM to hold the current (1/n!) term, and the other half of available RAM to hold the current sum, but the article points out that by working from right to left, (roughly) all of available RAM can be used for the sum, and thus twice as many digits can be calculated.

In fact, working from right to left is exactly how the spigot algorithm for e works. (Although they look quite different from one another, the formula in the article that illustrates how the series is calculated from right to left is the same "refactoring" (to borrow a software term) of the above series for e in Rabinowitz and Wagon's classic spigot paper.)

After trying the program in the article, I wrote a short implementation of a spigot algorithm for comparison. Of course, I did not have the patience to wait out a 34500 (decimal) digit calculation (which is the number of digits the programs in article will output). Instead I made NUMPAG 1 (and adjusted elsewhere accordingly) which gives 616 correct digits (after the decimal point). (Note that a 256 byte number is a 616 (decimal) digit number, i.e. 256*log10(256) =~ 616.5)

The spigot program -- which also only calculates (actually it only converts one base to another) the fractional part -- uses a comparable amount of space (254 bytes) for the mixed-radix number, and the largest divisor used is 255 (and divisors are therefore 8 bits wide -- the main reason for cutting off the calculation here). Rabinowitz and Wagon state that n+2 mixed-radix digits suffices for n decimal digits. This is true, but it greatly understates the number of correct digits you will actually get, especially as n gets larger. A divisor of 255 corresponds to the 1/255! term in the series above, so we should expect that there should be about 504 correct digits, since log10(255!) =~ 504.5, and this is exactly how many correct digits there happen to be. It is also double the number of digits (254-2) you would expect according to Rabinowitz and Wagon's criteria.

As you may have noticed, this initial comparison of the two programs is not going to be totally fair (616 digits vs. 504 digits). The idea at this stage was to have them close enough to do an initial rough A-B comparison, rather than attempt to draw any final conclusions.

Woz's program has two passes: the first pass is the fixed point calculation of e, and the second pass converts that result to base 10. In the second pass there is an assembly routine to multiply by 100 (to obtain the next 2 digits); the BASIC program (using division and MOD) splits apart the two digits, and outputs them, formatting them into neat columns.

Because the spigot program is entirely assembly, step one is to convert the BASIC program to assembly, so that we are comparing fixed-point to spigot rather than the BASIC interpreter (and its division, output, etc. routines) to their more streamlined equivalents in assembly.

In addition, there is a routine which will patch the (original) assembly code to use only 1, 2, or 4 (instead of 56 = $38) pages for the fixed-point numbers, and to use only as many terms in the series as needed.

From the Apple II monitor, to patch the program to use 4 pages for the fixed-point number (and 1000 terms) as recommended in the article, and calculate and output the result, type: (the monitor places the (low byte) of hex number before the < in zero page location $42)

8<4403G4400G

1 page (i.e. a page of paper, not a page of memory) of digits is output; and 42 lines of the last (and in this case, only) page are output. There are 60 digits per line (in 12 groups of 5), so that is 42*60 = 2520 digits (after the decimal point).

Likewise, to use 1 (memory) page for the fixed-point number, type:

10<4403G4400G

Note that we want to output all of the correct digits, since we are only outputting entire lines; so we actually output more than the expected number of correct digits, and indeed the last few digits output will be incorrect.

Likewise, to use 2 (memory) pages for the fixed-point number, type:

18<4403G4400G

To patch the assembly back to its original values in the program listings, type:

0<4403G

Here is the assembly listing:



The spigot program is mostly the same as the algorithm described by Rabinowitz and Wagon. It uses base 100 (thus calculating 2 digits as at time, same as Woz's program), i.e. it multiplies by 100 instead of 10. It uses exactly the same divide-by-10 tables as the program above (GENTBLDIV is identical in both programs) to split two base 10 digits. The multiply-by-100 tables are also identical. GENTBLMUL is nearly the same as Woz's INIT (Listing 2), the only difference being that the former uses LDA MUL100L,X (4 cycles, 3 bytes) and the latter uses PHA and PLA (7 cycles, 2 bytes).

The spigot program also formats the output into neat columns of 5 groups of 5 digits per line, with a blank line after every 5 lines of digits. This makes it very easy to see how many digits you have. I chose 5 groups instead of 10, since the former will fit on the Apple II 40 column text screen, and that's much easier to read than the 80 column screen. (Paging BDD to the white courtesy phone for remarks about aging and vision )

Note that the addresses of the tables, and in particular of N (the mixed-radix) number, were chosen to avoid page boundary crossings.

The LDA #DIGITS&1 line looks superfluous, but it means the program will still work if you reduce the value of DIGITS. You could use conditional assembly here instead. DIGITS, by the way, is the total number of digits (before and after the decimal point).

None of these routines (including Woz's assembly routines) are Apple specific, so it's easy to run them on another machine. As you'd expect, OUTPUT ($FDED on the Apple II) outputs the character in the accumulator (on the Apple II (normal video) characters have their high bit set, hence the EOR #$80). OUTCRLF ($FD8E on the Apple II) outputs the appropriate CR,LF sequence for your system (on the Apple II, this is a CR only).

What would we expect in a fair comparison?

We'd expect the advantage of the fixed-point calculation to be that it's looping over a smaller amount of memory (a 504 (decimal) digit can be stored in 210 bytes, whereas the mixed-radix number is stored in 254 bytes).

We'd expect the advantage of the spigot routine to be that it only has one pass, since it doesn't actually calculate e, it just converts to base 10. (In the mixed-radix base, e is exactly 1.111111... (i.e. repeating forever) in the same way that 1/3 is exactly .333333... in base 10). We can also make a spigot calculation speed up as it goes along (Rabinowitz and Wagon's paper doesn't mention this); the spigot program currently doesn't do this.

In the unfair comparison we would expect the spigot routine to be faster for a couple of reasons (thus if it weren't, there would be no need to make the test more fair). First, it's calculating fewer digits (504 vs. 616). Second, it uses 8-bit divisors (the fixed-point calculation uses 16-bit divisors, since it goes beyond the 1/255! term).

Indeed, the spigot is faster in the unfair comparison. I measured 16 seconds for the spigot, and 27 seconds for fixed-point, both at ~1.02 MHz.

Thus the next step (which I have not yet taken) would be to do a fair comparison. First, calculate the same number of digits for both. Second, use 8 bit divisors for both for the terms before 1/255!, and 16 bit divisors for both for the terms beyond 1/255!. Third, add the speed-up-as-you-go-along optimization to the spigot calculation (this is fair, since we want to compare an optimal spigot algorithm to an optimal fixed-point algorithm -- rather than comparing an 8-bit division routine to a 16-bit division routine). (By the way, do you see why the fixed-point calculation can't do the spigot speed-up optimization?)

Here is the spigot program:

Well-timed, as ever!

...and, of course, a very interesting investigation. I had half a mind to construct a URL to run in visual6502 or JSBeeb, but I haven't managed to do it yet.

Right .....


That brings back memories, of doing something similarly silly back in 1984 / secondary school ... Don't have my 6502 assembly program anymore but a yellowed printout of the 20K decimals I calculated with it on my Acorn BBC remains. Those were the days. Exact same thing for pi (first year at university I redid that in some compiled language (Pascal at that time, I think) and ran it on one of the mainframes overnight : way quicker than on the Acorn BBC ).

In 1984 and 85, I had an account on a CDC Cyber 170 at CSU Sacramento, for running the IBM 360 emulator and the Pascal compiler as part of my normal Computer Science curriculum. The Cyber was a cumbersome beast, but a very powerful one for its time. I was able to run some number-crunching experiments on it, but I had to be careful about the length of my line printer output, because the sysop would get annoyed at long printouts that didn't look like they belonged to any recognizable class assignment. I think I still have some 17" wide tractor feed printouts in my attic somewhere, near the Forth issue of Byte that I still need to package up and send to Garth. Summer is over, so I should be able to crawl up there and look around without risking a heat stroke. If the printouts are still readable, I'll try to see if I can do something interesting with some on-line emulators, and post it over at anycpu.org

Mike B.