e-Day

[Valentin Albillo]

.
Hi, Patrick:

First of all, welcome to the forum. I fully expect you'll enjoy the place and enrich us all with your own contributions. That said:

Quote:Today I'm celebrating e-Day. At least in the rest of the date format world a look at the calendar shows 2.718!

Indeed, it shows in that order where I live. To contribute to your celebration, here's a very short program I wrote many many years ago to compute e to in excess of 200 digits in an HP-15C. Here's an extract taken from my article "Long Live the HP-15C !":

This 64-step program will compute from 8 to 208 decimal places of Euler’s constant, the well-known transcendental number e = 2.71828+ . It is by no means optimized for performance but tries instead to be as short and straightforward as possible. Although you can compute more decimal places in an HP-15C, for the purposes of this article this simpler program will do nicely.

Program listing:

Very Important: steps 30 and 43 must be entered in USER mode, while all the rest must be entered out of USER mode. An u-like character must appear next to the step number only for steps 30 and 43, and no others.

01  LBL A         23  RCL A         45  FRAC
02  MATRIX 0      24  X=0?          46  CHS
03  MATRIX 1      25  ISG 2         47  RTN
04    1           26  LBL 0         48  LBL 5
05  DIM A         27  RCL A         49  FRAC
06  DIM B         28  RCL÷ I        50  RCL B
07  RESULT B      29  INT           51  INT
08  STO 2         30u STO A         52  RCL RAN#
09  STO I         31  GTO 2         53    *
10   EEX          32  RCL I         54    +
11    8           33  PSE           55   R/S
12  STO A         34  RCL MATRIX A  56  GTO 4
13   1/X          35  MATRIX 7      57  LBL 2
14  STO RAN#      36  TEST 0        58  RCL* I
15  LBL 1         37  GTO 1         59    -
16  RCL MATRIX A  38  RCL MATRIX B  60  RCL RAN#
17  RCL MATRIX B  39  RCL RAN#      61    ÷
18    +           40    *           62  STO+ A
19  RCL 2         41  FIX 8         63  RCL A
20  STO O         42  LBL 4         64  GTO 0
21  ISG I         43u RCL B        
22  LBL 3         44  GTO 5    

The constant e is computed using the well-known formula:

      e = 1 + 1/1! +1/2! + 1/3! + 1/4! + ...

using enough terms to achieve the desired precision. We use the powerful matrix capabilities of the HP-15C, holding the current term in a vector (matrix A), and the running sum in another vector (matrix B). Steps 01-14 dimension and initialize both matrices, as well as some indexes and ancillary constants.

Starting with 1, each term is computed by dividing the previous one by the corresponding divisor up to a maximum of 208 decimal digits, until we arrive at a term that is zero to the specified precision. This multiprecision arithmetic is done by considering each term as composed of N blocks (1-26), each holding 8 digits, the extra digits in each register being carried out to the next block. Since we may use divisors up to 125, each block is limited to 8 digits, lest the carry would make the next block larger than the 10 digits an HP-15C storage register can hold.

The addition of each multi-block term to the running total is done using matrix arithmetic in steps 16-18. Steps 19-25 increment the divisor and keep track of the first non-zero block of each term to optimize speed by avoiding unnecesary operations. Steps 26-31 and 57-64 perform the multiprecision division.

The process ends when the current term has all its blocks equal to zero. As all blocks always hold positive values, we can use an advanced matrix operation, the Row Norm, to test for finalization, because in this case the Row Norm will equal zero if and only if all block values are zero. This is tested in steps 34-37.

Once this condition is met, the result matrix which holds the running sum is scaled down for display using matrix-scalar arithmetic (steps 38-40). The computed answer is then displayed by recalling each block in turn, adding the carry from the previous block, and marking the last one negative (steps 41-56).

As you can see, though simple, this program does use some of the HP-15C’s advanced programming capabilities (as well as a trick or two), including basic matrix operations (MATRIX 1), advanced matrix functions (MATRIX 7), recall arithmetic (RCL÷ I), using registers as indexes including increment-and-test operations (ISG 2) , auto-increment-test-and-loop matrix element access (uSTO A, uRCL B), matrix arithmetic (steps 18 and 40), etc.

Usage instructions:

1) After keying in the program, make sure you’re not in complex mode (press CF 8 if in doubt), then you must commit enough storage registers to the common pool for the matrix operations. To that effect, press:

      2, f DIM (i)

2) Now, enter the number of 8-digit blocks you want to use, from 1 (8 digits) to 26 (208 digits). For example, if you want to compute 200 decimal digits of e, you must specify 200/8 = 25 blocks. Start the program by pressing either

      f PRGM, R/S or GSB A or (in User mode) A

While running, it will briefly show each succesive divisor used (2, 3, ...), then once the computation is over, it will display each block of 8 decimal digits (with an initial “0.”, in order to preserve leading zeros at the left end of the block), starting with decimals 1st-8th. The very last block will be marked negative, to signal the end of the output.

Let’s see an example: to compute the first 24 decimal digits of e, we specify 24/8 = 3 blocks, and proceed like this (in USER mode):

     3, A -> (2.00000000) [divisor = 2]
          ... [after 2’26”]
          -> (25.00000000) [divisor = 25]
          -> 0.71828182 [decimals 1st- 8th]
      R/S -> 0.84590452 [decimals 9th-16th]
      R/S -> -0.35360274 [decimals 17th-24th]

so, after adding a “2.” at the front and writing down all 8-digit blocks (that is, minus the initial “0.” or “-0.”) in their proper order, we finally get:

      e = 2.71828182 84590452 35360274

where, due to the accumulation of rounding errors during the process, the last block of the computed answer comes out as “3536 0274” while correct is “3536 0287”, so we have an error of 13 ulps (units in the last place).

3) Should you need to display the output again, press: GSB 4

For 208 decimal digits, the final result is:

      e = 2.71828182 84590452 35360287 47135266 24977572
            47093699 95957496 69676277 24076630 35354759
            45713821 78525166 42742746 63919320 03059921
            81741359 66290435 72900334 29526059 56307381
            32328627 94349076 32338298 80753195 25101901
            15738281

The last block of the computed answer comes out as “1573 8281” while correct is “1573 8341”, so the error is 60 ulps, and thus after rounding the last two places we’ve got 206 decimals fully correct.

Regards.
V.
.