/*---------------------------------------------------------------------- File : euler.c Contents: compute e (Euler's number) up to a given precision using the formulae e = \sum_{n=0}^\infty 1/n! e = \lim_{n\to\infty}(1+1/n)^n e = \lim_{n\to\infty}(1-1/n)^-n Author : Christian Borgelt History : 05.11.1997 file created 07.11.1997 comments on the ++-operator added 12.11.1997 error in comment on e2 removed ----------------------------------------------------------------------*/ #include /*---------------------------------------------------------------------- Preprocessor Definitions ----------------------------------------------------------------------*/ #define PRECISION 0.000001 /* desired precision of e */ /*---------------------------------------------------------------------- The '#define' directive is not evaluated by the compiler, but by another program called the preprocessor. This program is run on the input files before the actual compiler is invoked. It can be used to make certain replacements in the source code. The above directive defines such a replacement. Whenever the preprocessor encounters the expression following the '#define' in the source code (here: if it finds the string 'PRECISION'), it replaces it by the expression on the right (here: by '0.000001'). This feature comes in handy, if a certain constant appears several times in the source code. If you had to write it as a number every time it is needed, and if you then wanted to change that number, you would have to recall carefully all the places in which it is used. The preprocessor frees you from such considerations. Just define a symbolic constant (like 'PRECISION') with the desired value and use this symbolic constant in the source code. If you then want to change the value of constant, you can simply change the definition and leave the rest to the preprocessor. ----------------------------------------------------------------------*/ /*---------------------------------------------------------------------- Auxiliary Functions ----------------------------------------------------------------------*/ double power (double x, int n) { /* --- compute x^n */ double t; /* temporary buffer */ if (n == 1) return x; /* x^1 can be computed directly */ t = power(x, n/2); /* compute x^{n/2} or x^{(n-1)/2} */ if (n % 2 == 0) /* if the exponent is even, */ return t*t; /* then x^n = (x^{n/2})^2 */ else /* if the exponent is odd, */ return t*t*x; /* then x^n = (x^{(n-1)/2})^2 * x */ } /* power() */ /*---------------------------------------------------------------------- Note that it is not necessary to compute power(x, (n-1)/2) in case n is odd, since the quotient n/2 is computed as an integer number. That is, any fraction is discarded from the result. Hence, for odd n, it is n/2 = (n-1)/2. Note also that power(x, n) works properly only for n > 0. ----------------------------------------------------------------------*/ /*---------------------------------------------------------------------- Main Functions ----------------------------------------------------------------------*/ double e1 (void) { /* --- e = \sum_{n=0}^\infty 1/n! */ double e; /* e as computed in current step */ double t; /* next term of the sum */ int n; /* loop variable (summation index) */ e = t = 1.0; /* start summing with the term 1/0! */ n = 1; /* initialize n for the second term */ while (t > PRECISION) { /* while des. precision not reached */ t /= n; /* compute next term of the sum, */ e += t; /* add it to the value of e that */ n++; /* has been computed up to now and */ } /* increment n for the next term */ return e; /* return computed approximation */ } /* e1() */ /*---------------------------------------------------------------------- Note that the computation is made efficient by reusing the value of the term t from the previous execution of the loop. It would be much less efficient to compute n! anew in each step, since computing 1/n! from scratch takes n multiplications and one division (or n divisions, which is usually even worse), whereas computing it from the previous value of the term t = 1/(n-1)! takes only one division. For the experts: It is possible to write the loop body with only two statements, namely 't /= n++;' and 'e += t;'. This is possible since the ++-operator increments the value of the variable it is applied to only after this value has been used in the term it appears in. That is, 't /= n++;' first computes 't /= n;' and then increments the value of n by one. But be careful. Do not write: 't /= ++n;'. Both 'n++' and '++n' increment n, but the first does so after the value of n has been used in the term it appears in, the second before it is used. ----------------------------------------------------------------------*/ double e2 (void) { /* --- e = \lim_{n\to\infty}(1+1/n)^n */ double e_curr; /* e as computed in current step */ double e_prev; /* e as computed in previous step */ int n; /* loop variable */ e_curr = 1.5; /* start with e computed with n = 1 */ n = 2; /* initialize n for second trial */ do { /* repeat computing e for a new n */ e_prev = e_curr; /* note value of previous trial */ e_curr = power(1+1.0/n, n); /* compute new approximation of e */ n++; /* increment n for next trial */ } while (e_curr -e_prev > PRECISION); /* while des. precision not reached */ return e_curr; /* return computed approximation */ } /* e2() */ /*---------------------------------------------------------------------- Note that it is necessary to write 1.0/n instead of 1/n, since (see comment on power(n, x)) 1/n is computed as an integer number and thus any fraction is discarded from the result. Hence 1/n = 0 for all n > 1, and the loop would terminate after only two executions with result 1. The result computed shows, that the stopping criterion is not well suited for this type of formula, since the values of (1+1/n)^n increase very slowly. If two consecutive values differ less than PRECISION, this does not guarantee that the computed value of e differs less than PRECISION from the true value. For the experts: In this case it is dangerous to try to save one statement by writing the loop body as 'e_prev = e_curr' and 'e_curr = power(1+1.0/n, n++);' The reason is that C does not fix the order in which the different arguments of a function are evaluated. If the second argument is evaluated first, then the first argument will be computed with a wrong value of n. ----------------------------------------------------------------------*/ double e3 (void) { /* --- e = \lim_{n\to\infty}(1-1/n)^-n*/ double e_curr; /* e as computed in current step */ double e_prev; /* e as computed in previous step */ int n; /* loop variable */ e_curr = 4; /* start with e computed with n = 2 */ n = 3; /* initialize n for second trial */ do { /* repeat computing e for a new n */ e_prev = e_curr; /* note value of previous trial */ e_curr = 1.0 /power(1-1.0/n, n); /* compute new approximation of e */ n++; /* increment n for next trial */ } while (e_prev -e_curr > PRECISION); /* while des. precision not reached */ return e_curr; /* return computed approximation */ } /* e3() */ /*---------------------------------------------------------------------- The same comments apply as for function e2(). In addition, note that the stopping criterion is now e_prev -e_curr instead of e_curr -e_prev, since the values of (1-1/n)^-n decrease with increasing n, n > 1. This is also the reason, why the loop is started with n = 3 instead of n = 2, since for n = 1, we get e_curr = 0 and the loop would terminate after only one execution. ----------------------------------------------------------------------*/ int main (void) { /* --- compute e (Euler's number) */ printf("e = 2.7182818284590452354\n"); /* correct value */ printf("e1 = %.15f\n", e1()); /* e1 = \sum_{n=0}^\infty 1/n! */ printf("e2 = %.15f\n", e2()); /* e2 = \lim_{n\to\infty}(1+1/n)^n */ printf("e3 = %.15f\n", e3()); /* e3 = \lim_{n\to\infty}(1-1/n)^-n */ return 0; /* return 'ok' */ } /* main() */