# Computing e

## Introduction

One of the things computers are useful for is computing large (and small) numbers. The number called *e* is the base of natural logarithms, and is one of those "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^{1}

* 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

`), and two constants that define how much space to allocate to the computation workarea (`

**N**`and`

**X**`).`

**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:

The number of digits in a number X = | CEILING[log(X)] |

and N! = |
(N * N-1 * N-2 ... * 1) |

Thus, the number of digits in N! = | CEILING[log(N!)] |

= | CEILING[log(N * (N-1) * ... * 1)] |

= | CEILING[log(N) +log (N-1) ... +log(1)] |

Using the above formula, I calculated

N | # digits in N! |
---|---|

1775 |
4999 |

1776 |
5002 |

1777 |
5005 |

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:

# values that an 8-bit byte can hold | = 256 |

# digits we want to compute | = 5002 |

Then | 10^{5002} = 256^{X} |

then | 5002 log(10) = X log(256) |

and | X = 5002 (log(10)/log(256)) |

This last can be computed as | X = 5002 * (1/2.408) |

or | X = 2077.04 |

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 (

`), 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`

**I**`, which I calculated earlier.`

**X****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:

- divide the current value by
`n`(giving us the ((...)/n) value), - add 1 to the value (giving us the ((...)/n)+1) value), and
- 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:

- convert the integer portion to a printable character and print it, and
- 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
is a pointer to a character, which will point to the first byte of the`dvd`*dividend*array, - function argument
is an unsigned integer containing the`dvr`*divisor*, - function argument
is an unsigned integer indicating the number of bytes in the dividend array, and indicating the "precision" of the division.`prc` - automatic variable
is a pointer to a character, which the function will use to traverse the dividend array from high-order byte to low-order byte`Wk_q` - automatic variable
is an integer which the function will use to enumerate the bits within each byte of the dividend, and`bit` - automatic variable
is an integer which the function will use to hold the`rmdr`*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
is a pointer to a character, which will point to the first byte of the`mcd`*multiplicand*array, - function argument
is an unsigned integer containing the`mpy`*multiplier*, - function argument
is an unsigned integer indicating the number of bytes in the multiplicand array, and indicating the "precision" of the multiplication.`prc` - automatic variable
is a pointer to a character, which the function will use to traverse the multiplicand array from low-order byte to high-order byte`Wk_r` - automatic variable
is an integer which the function will use to enumerate the bits within each byte of the multiplicand, and`bit` - automatic variable
is an integer which the function will use to hold the`ovflo`*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

Attachment | Size |
---|---|

natlog.c | 2.2 KB |

natlog.results.txt | 6.65 KB |