Computing e

Introduction

One of the things computers are useful for is computing large (and small) numbers. "Euler's number" (the number called e) is the base of natural logarithms, and is one of those computable "transcendental numbers".

This program computes the value of e to 5,000 decimal places, and was inspired by the article "The Impossible Dream: Computing e to 116,000 places with a Personal Computer" by Stephen Wozniak, published in Byte Magazine Vol 6, Issue 6 (June 1981, pp392).

Essentially, the program implements a partial expansion of the Taylor series equation for e¹
    e = 1 + (1/1!) + (1/2!) + (1/3!) + ... + (1/n!)
by inverting the equation, such that it reads...
   e = (((((...((((((1)/n) + 1) / (n-1)) + 1) ... / 3) + 1) / 2) + 1) / 1) + 1

This last equation translates well into computer code. The pseudocode for this equation reads:

n = STARTING_VALUE
result = 1
WHILE n IS GREATER THAN 0
  result = result / n
  result = result + 1
  n = n - 1
END-WHILE

I originally wrote this program in Z80 Assembler code, to run under CP/M (Dec 1986), and subsequently have recoded the program in 370 Assembler to run under MVS and VSE, and C to run under MSDOS and Linux

The natlog program

Program prologue code

First off, I set up some constants and the storage where the program will hold it's computations. I have two constants that define how deep the program will go into the equation that it uses for the computation (P and N), and two constants that define how much space to allocate to the computation workarea (X and I).

Selecting the precision (P)

I selected the value for P as the number of digits I wanted to see to the right of the decimal place (5,000). The program will use this value to control the number of times we loop through the multiplication logic that controls the output.

Selecting the Largest Factorial (N)

I chose the value for N based on the number of digits I wanted to see in N!. The program will use this value to control the number of times it loops through the division logic. In this case, I want 5,000 decimal places, so I found a value of N in the following manner:

Using the above formula, I calculated

and found that the lowest factorial that meets my 5000 digit requirements is 1776!, with 5002 digits.

Selecting the size of the fraction (X)

I set X to the number of 8-bit bytes the program will need to store the binary representation of the 5002 decimal digits in the fraction portion of the computation. I found this value through the following computation:

To obtain the number of whole bytes needed to store our fraction, I rounded this result up to 2078 .

Allocating the space for our computation

The code allocates the computation workarea, e, as an array, large enough to fit the planned computation, and will treat it as a large-precision fixed-point number. I set the size of the integer portion of this fixed-point number to 1 byte (I), as the program only needs to store in it a value of zero or one. I set the size of the fraction portion to the value of X, which I calculated earlier.

Note: The program will use this computational workarea to hold a "fixed point binary" number that will become our result. The program partitions the workarea into two components: an 8 bit integer, and a 16624-bit fraction. You can mentally place the "decimal point" between the first and second byte of the workarea.

The program prelude

/*
** Copyright (c) Lew Pitcher, 1999
*/
#include <stdio.h>

#define P 5000 /* number of multiply iterations     */

#define N 1776 /* number of divide iterations       */

#define X 2078 /* number of bytes in fraction part  */
#define I 1    /* number of bytes in integer part   */

char e[I+X];   /* holds e fixed point binary value  */

main()

The main() function performs the bulk of the work. It will first compute the value of e in fixed-point binary, then it will print the value in characters, converting it on-the-fly from the fixed-point binary number to a decimal display value.

Function initialization

The first part of main() establishes forward references for the fixed-point multiplication and fixed-point division routines, and allocates the variables that the remaining logic will use to control the divisions and multiplications.

Compute phase

The compute phase starts by initializing the e working value to 1.something. Ideally, we should set e to 1.0, but we take a shortcut here, because the uninitialized (and random) fractional part of the initial value will divide down to insignificance (and disappear off the low-order fractional end of the value) as we perform our division loop.

As per our rearranged equation
   e = (((((...(((((1/n) + 1) / (n-1)) + 1) ... / 3) + 1) / 2) + 1) / 1) + 1
we start the computation at n = N, and work our way back to n = 1. We compute in a loop of three steps:

1. divide the current value by n (giving us the ((...)/n) value),
2. add 1 to the value (giving us the ((...)/n)+1) value), and
3. subtract 1 from n (giving us the constantly decreasing divisor

Note:

• On the second last iteration of the loop, we divide by 2, and then add 1. This corresponds to the ( (...)/2 + 1) portion of the equation.
• On the last iteration of the loop, we divide by 1, and then add 1. This corresponds to the ( (...) + 1) portion of the loop, as that can be re-expressed as ( (...)/1 + 1).

Report phase

The report phase works in two parts:

1. convert the integer portion to a printable character and print it, and
2. repeatedly convert the fraction portion to printable characters, and print them

The integer portion is easy: since the integer portion is (from our computations) guaranteed to be a single digit integer (2, actually), we add the character '0' to it, to come up with a printable digit.

The fraction portion is almost as easy: we zero out the integer portion, and then multiply the whole number by 10. This shifts the next (fraction) decimal digit into the integer portion of the result, so that it can be converted to a printable character. To convert that integer (which will be a binary value from 0 to 9), we simply add the character '0' to it to come up with a printable character, just like we did with the integer portion.

We will do this fractional portion in a 5000-iteration loop, zeroing out the integer, multiplying by 10, and converting the resulting integer to a character for each pass. This way, we print all 5000 decimal places of the result.

Of course, the code looks a little more complex than this explanation. The complexity comes from our output formatting requirements. We want each line of the output to hold 50 decimal places, and each line of the result should line up under the previous line. After every 5 decimal places, we want to print a space, to make it easier to read the number.

The main() function

int main(void)
{
    void divide(char *, unsigned, unsigned),
         multiply(char *, unsigned, unsigned);

    unsigned n, /* divisor for calculation loop */
             p; /* precision for print loop     */


    /*********** Compute phase ***********/
    e[0] = 1;  /* start the series at 1.anything */
    for (n = N; n > 0; --n)
    {
        divide(e,n,I+X);
        e[0] += 1;
    }


    /*********** Report phase ***********/
    fputs(" e == ",stdout);
    fputc(e[0]+'0',stdout);
    fputc('.',stdout);

    for (p = 0; p < P; ++p)
    {
        e[0] = 0;
        multiply(e,10,I+X);
        fputc(e[0]+'0',stdout);
        if (p%5 == 4)
            if (p%50 == 49)
                fputs("\n\t",stdout);
            else
                fputc(' ',stdout);
    }
    fputc('\n',stdout);

    return 0;
}

divide()

The function divide() will divide a multibyte fixed-point dividend by an unsigned integer divisor. As division progresses, since we process the dividend from high-order to low-order, the quotient gradually replaces the dividend without interfering with the division. The number of bytes in the multibyte fixed point dividend is specified as a parameter, permitting the function to be used on many different size dividends.

I've written this function as an implementation of "long division", performing the division one "bit" at a time.

Note: There are, of course, much more efficient methods of division that I could have used. I could have, for instance, performed the division one byte at a time, rather than one bit at a time. I could also have scaled the division based on the magnitude of the divisor, knowing that the earlier invocations of this function would compute low-order values that would roll off the end of the result as the later divisions were computed. However, I was not interested in efficiency, but rather in exploring the nature of the algorithm.

Function Initialization

At the beginning of the function, we define three function arguments, and allocate three automatic variables:

• function argument dvd is a pointer to a character, which will point to the first byte of the dividend array,
• function argument dvr is an unsigned integer containing the divisor,
• function argument prc is an unsigned integer indicating the number of bytes in the dividend array, and indicating the "precision" of the division.
• automatic variable Wk_q is a pointer to a character, which the function will use to traverse the dividend array from high-order byte to low-order byte
• automatic variable bit is an integer which the function will use to enumerate the bits within each byte of the dividend, and
• automatic variable rmdr is an integer which the function will use to hold the remainder during the division process

Before any computations, the program initializes the remainder (rmdr) to zero.

Byte-by-byte division loop

Starting at the high-order byte of the dividend array, and working downward toward the low-order byte, the code divides each byte by the divisor. As it goes, it replaces the dividend byte with the computed quotient byte, and keep track of the intra-byte remainders.

Bit-by-bit division loop

Starting at the high-order bit of the selected dividend byte, and working downward toward the low-order bit, the code divides each bit by the divisor. As it proceeds, it replaces the dividend bit with the computed quotient bit, and keep track of the intra-bit remainders.

Single bit division

As each lower bit represents ½ of the previous bit, and the rmd remainder is what was left over after division by that previous bit, the function first doubles the remainder. This gives the remainder the same magnitude as the bit we are now dividing.

The current bit is then added to the remainder, so that the remainder now represents the value (at the current magnitude) to be divided.

The dividend byte is now manipulated so that the next bit to be divided is properly positioned for the next iteration of bit-by-bit division loop. This also initializes the corresponding quotient bit to zero. After 8 successive loops through this logic, the entire 8-bit byte will contain quotient bits, in the same order as the original dividend bits.

If the remainder is now greater than or equal to the divisor

• the quotient bit for this dividend bit is set to 1, and
• the remainder is reduced by the amount of the divisor.

However, if the remainder was less than the divisor, then

• the quotient bit for this dividend bit is left unchanged (at zero), and
• the remainder is not reduced.

The remainder is retained for the next bit in the bit-by-bit division loop, and for the next byte in the byte-by-byte division loop.

The divide() function

void divide(char *dvd, unsigned dvr, unsigned prc)
{
    char *Wk_q;
    unsigned bit, rmdr;

    rmdr = 0;

    for (Wk_q = dvd; Wk_q < dvd+prc; ++Wk_q)
    {
        for(bit = 0; bit < 8; ++bit)
        {
            rmdr *= 2;
            if (*Wk_q & 0200) rmdr += 1;
            *Wk_q <<= 1;
            if (rmdr >= dvr)
            {
                *Wk_q += 1;
                rmdr -= dvr;
            }
        }
    }
}

multiply()

The function multiply() will multiply a multibyte fixed-point multiplicand by an unsigned integer multiplier. As multiplication progresses, since we process the multiplicand from low-order to high-order, the product gradually replaces the multiplicand without interfering with the multiplication. The number of bytes in the multibyte fixed point multiplicand is specified as a parameter, permitting the function to be used on many different size multiplicands.

Note: There are, of course, much more efficient methods of multiplication that I could have used. I could have, for instance, performed the multiplication one byte at a time, rather than one bit at a time. I could also have scaled the multiplication based on the magnitude of the multiplicand, knowing that the later invocations of this function would find low-order zeros that could be ignored as they would not affect the results of the computation. However, I was not interested in efficiency, but rather in exploring the nature of the algorithm.

Function Initialization

At the beginning of the function, we define three function arguments, and allocate three automatic variables:

• function argument mcd is a pointer to a character, which will point to the first byte of the multiplicand array,
• function argument mpy is an unsigned integer containing the multiplier,
• function argument prc is an unsigned integer indicating the number of bytes in the multiplicand array, and indicating the "precision" of the multiplication.
• automatic variable Wk_r is a pointer to a character, which the function will use to traverse the multiplicand array from low-order byte to high-order byte
• automatic variable bit is an integer which the function will use to enumerate the bits within each byte of the multiplicand, and
• automatic variable ovflo is an integer which the function will use to hold the carry during the multiplication process

Before any computations, the program initializes the carry (ovflo) to zero.

Byte-by-byte multiplication loop

Starting at the low-order byte of the multiplicand array, and working upward toward the high-order byte, the code multiplies each byte by the multiplier. As it goes, it replaces the multiplicand byte with the computed product byte, and keep track of the intra-byte carries.

Bit-by-bit multiplication loop

Starting at the low-order bit of the selected multiplicand byte, and working upward toward the high-order bit, the code multiplies each bit by the multiplier. As it proceeds, it replaces the multiplicand bit with the computed product bit, and keep track of the intra-bit carries.

Single bit multiplication loop

If the bit to be multiplied is equal to 1, then the carry is incremented by the amount of the multiplier (effectively multiplying that bit by the multiplier, and adding it to the carry).

The multiplicand byte is now manipulated so that the next bit to be multiplied is properly positioned for the next iteration of the bit-by-bit multiplication loop. This also initializes the corresponding product bit to zero. After 8 successive loops through this logic, the entire 8-bit byte will contain product bits, in the same order as the original multiplicand bits.

The low-order bit of the carry is then placed into the corresponding bit of the product, and the carry is halved. The carry is retained for the next bit in the bit-by-bit multiplication loop, and for the next byte in the byte-by-byte multiplication loop.

The multiply() function

void multiply(char *mcd, unsigned mpy, unsigned prc)
{
    char *Wk_r;
    unsigned bit, ovflo;

    ovflo = 0;
    for (Wk_r = mcd+prc-1; Wk_r >= mcd; --Wk_r)
    {
        for (bit = 8; bit > 0; --bit)
        {
            if (*Wk_r & 1) ovflo += mpy;
            *Wk_r >>= 1; *Wk_r &= 0177;
            if (ovflo & 1) *Wk_r = *Wk_r | 0200;
            ovflo /= 2;
        }
    }
}

The results from natlog

e == 2.71828 18284 59045 23536 02874 71352 66249 77572 47093 69995
	95749 66967 62772 40766 30353 54759 45713 82178 52516 64274
	27466 39193 20030 59921 81741 35966 29043 57290 03342 95260
	59563 07381 32328 62794 34907 63233 82988 07531 95251 01901
	15738 34187 93070 21540 89149 93488 41675 09244 76146 06680
	82264 80016 84774 11853 74234 54424 37107 53907 77449 92069
	55170 27618 38606 26133 13845 83000 75204 49338 26560 29760
	67371 13200 70932 87091 27443 74704 72306 96977 20931 01416
	92836 81902 55151 08657 46377 21112 52389 78442 50569 53696
	77078 54499 69967 94686 44549 05987 93163 68892 30098 79312
	77361 78215 42499 92295 76351 48220 82698 95193 66803 31825
	28869 39849 64651 05820 93923 98294 88793 32036 25094 43117
	30123 81970 68416 14039 70198 37679 32068 32823 76464 80429
	53118 02328 78250 98194 55815 30175 67173 61332 06981 12509
	96181 88159 30416 90351 59888 85193 45807 27386 67385 89422
	87922 84998 92086 80582 57492 79610 48419 84443 63463 24496
	84875 60233 62482 70419 78623 20900 21609 90235 30436 99418
	49146 31409 34317 38143 64054 62531 52096 18369 08887 07016
	76839 64243 78140 59271 45635 49061 30310 72085 10383 75051
	01157 47704 17189 86106 87396 96552 12671 54688 95703 50354
	02123 40784 98193 34321 06817 01210 05627 88023 51930 33224
	74501 58539 04730 41995 77770 93503 66041 69973 29725 08868
	76966 40355 57071 62268 44716 25607 98826 51787 13419 51246
	65201 03059 21236 67719 43252 78675 39855 89448 96970 96409
	75459 18569 56380 23637 01621 12047 74272 28364 89613 42251
	64450 78182 44235 29486 36372 14174 02388 93441 24796 35743
	70263 75529 44483 37998 01612 54922 78509 25778 25620 92622
	64832 62779 33386 56648 16277 25164 01910 59004 91644 99828
	93150 56604 72580 27786 31864 15519 56532 44258 69829 46959
	30801 91529 87211 72556 34754 63964 47910 14590 40905 86298
	49679 12874 06870 50489 58586 71747 98546 67757 57320 56812
	88459 20541 33405 39220 00113 78630 09455 60688 16674 00169
	84205 58040 33637 95376 45203 04024 32256 61352 78369 51177
	88386 38744 39662 53224 98506 54995 88623 42818 99707 73327
	61717 83928 03494 65014 34558 89707 19425 86398 77275 47109
	62953 74152 11151 36835 06275 26023 26484 72870 39207 64310
	05958 41166 12054 52970 30236 47254 92966 69381 15137 32275
	36450 98889 03136 02057 24817 65851 18063 03644 28123 14965
	50704 75102 54465 01172 72115 55194 86685 08003 68532 28183
	15219 60037 35625 27944 95158 28418 82947 87610 85263 98139
	55990 06737 64829 22443 75287 18462 45780 36192 98197 13991
	47564 48826 26039 03381 44182 32625 15097 48279 87779 96437
	30899 70388 86778 22713 83605 77297 88241 25611 90717 66394
	65070 63304 52795 46618 55096 66618 56647 09711 34447 40160
	70462 62156 80717 48187 78443 71436 98821 85596 70959 10259
	68620 02353 71858 87485 69652 20005 03117 34392 07321 13908
	03293 63447 97273 55955 27734 90717 83793 42163 70120 50054
	51326 38354 40001 86323 99149 07054 79778 05669 78533 58048
	96690 62951 19432 47309 95876 55236 81285 90413 83241 16072
	26029 98330 53537 08761 38939 63917 79574 54016 13722 36187
	89365 26053 81558 41587 18692 55386 06164 77983 40254 35128
	43961 29460 35291 33259 42794 90433 72990 85731 58029 09586
	31382 68329 14771 16396 33709 24003 16894 58636 06064 58459
	25126 99465 57248 39186 56420 97526 85082 30754 42545 99376
	91704 19777 80085 36273 09417 10163 43490 76964 23722 29435
	23661 25572 50881 47792 23151 97477 80605 69672 53801 71807
	76360 34624 59278 77846 58506 56050 78084 42115 29697 52189
	08740 19660 90665 18035 16501 79250 46195 01366 58543 66327
	12549 63990 85491 44200 01457 47608 19302 21206 60243 30096
	41270 48943 90397 17719 51806 99086 99860 66365 83232 27870
	93765 02260 14929 10115 17177 63594 46020 23249 30028 04018
	67723 91028 80978 66605 65118 32600 43688 50881 71572 38669
	84224 22010 24950 55188 16948 03221 00251 54264 94639 81287
	36776 58927 68816 35983 12477 88652 01411 74110 91360 11649
	95076 62907 79436 46005 85194 19985 60162 64790 76153 21038
	72755 71269 92518 27568 79893 02761 76114 61625 49356 49590
	37980 45838 18232 33686 12016 24373 65698 46703 78585 33052
	75833 33793 99075 21660 69238 05336 98879 56513 72855 93883
	49989 47074 16181 55012 53970 64648 17194 67083 48197 21448
	88987 90676 50379 59036 69672 49499 25452 79033 72963 61626
	58976 03949 85767 41397 35944 10237 44329 70935 54779 82629
	61459 14429 36451 42861 71585 87339 74679 18975 71211 95618
	73857 83644 75844 84235 55581 05002 56114 92391 51889 30994
	63428 41393 60803 83091 66281 88115 03715 28496 70597 41625
	62823 60921 68075 15017 77253 87402 56425 34708 79089 13729
	17228 28611 51591 56837 25241 63077 22544 06337 87593 10598
	26760 94420 32619 24285 31701 87817 72960 23541 30606 72136
	04600 03896 61093 64709 51414 17185 77701 41806 06443 63681
	54644 40053 31608 77831 43174 44081 19494 22975 59931 40118
	88683 31483 28027 06553 83300 46932 90115 74414 75631 39997
	22170 38046 17092 89457 90962 71662 26074 07187 49975 35921
	27560 84414 73782 33032 70330 16823 71936 48002 17328 57349
	35947 56433 41299 43024 85023 57322 14597 84328 26414 21684
	87872 16733 67010 61509 42434 56984 40187 33128 10107 94512
	72237 37886 12605 81656 68053 71439 61278 88732 52737 38903
	92890 50686 53241 38062 79602 59303 87727 69778 37928 68409
	32536 58807 33988 45721 87460 21005 31148 33513 23850 04782
	71693 76218 00490 47955 97959 29059 16554 70505 77751 43081
	75112 69898 51884 08718 56402 60353 05583 73783 24229 24185
	62564 42550 22672 15598 02740 12617 97192 80471 39600 68916
	38286 65277 00975 27670 69777 03643 92602 24372 84184 08832
	51848 77047 26384 40379 53016 69054 65937 46161 93238 40363
	89313 13643 27137 68884 10268 11219 89127 52230 56256 75625
	47017 25086 34976 53672 88605 96675 27408 68627 40791 28565
	76996 31378 97530 34660 61666 98042 18267 72456 05306 60773
	89962 42183 40859 88207 18646 82623 21508 02882 86359 74683
	96543 58856 68550 37731 31296 58797 58105 01214 91620 76567
	69950 65971 53447 63470 32085 32156 03674 82860 83786 56803
	07306 26576 33469 77429 56346 43716 70939 71930 60876 96349
	53288 46833 61303 88294 31040 80029 68738 69117 06666 61468

References

"The Calculation of e to Many Significant Digits", A. H. J. Sale, 1968

https://academic.oup.com/comjnl/article/11/2/229/378761

"The Impossible Dream: Computing e to 116,000 places with a Personal Computer", Stephen Wozniak, 1981

Byte Magazine Vol 6, Issue 6

Euler's number to 10,000 digits

https://www.math.utah.edu/~pa/math/e.html