calculating e

I've written a small code to calculate the constant e with 40 digits precision, but, apparently due to the limitations of long double, it's not accurate. Here's the code:

Code:

---

#include <stdio.h>

long double fact(long double n)
{
        long double sum;
        sum = n;

        if(n == 0)
                return 1;
        else
        {
                while(n != 1)
                {
                        sum = sum * (n - 1);
                        n = n - 1;
                }
       
        return sum;
        }
}

       
main()
{

        int counter;

        long double e, k;

        k = e = 0;

        for(counter = 0; counter < 40; counter++)

            {
                e = e + (1.0/fact(k));
                ++k;
            }

        printf("%.40f\n", e);

}

---

For some reason, the results of 1/21! to 1/40! are calculated but are not added to the rest of the series, resulting in a not-so-precise output. I can't figure out why this happens.


Could anyone suggest anything, perhaps another algorithm, for this calculations which does not rely on any special library?

Floating point is an approximation, and when adding relatively large number with a relatively small number, the precision is not retained. You may only have 15 or 18 digits of precision for your variable, for example.

[edit]So you may need to look into a math library.

Your algorithm would be more accurate if you looped from 39 to 0 instead of from 0 to 39. But you'll still be limited by the long double's precision.

I realise you said you don't want to rely on any special library, but you don't really have much choice, unless you want to reinvent an arbitrary precision arithmetic system.

As said, you are limited by your machine's floating point accuracy, so your options are:

1) Start doing all the maths yourself by storing your values/intermediate results in your own data structure or array.

2) Use a pre-existing arbitrary math library. One example is libgmp.

3) printf("2.7182818284590452353602874713526624977572 47093699959574966967627724076630353547594571382178 525166427427466391932003059921817413596629\n"); :p

Quote:

---

Originally Posted by cwr

1) Start doing all the maths yourself by storing your values/intermediate results in your own data structure or array.

---

Could you show me how? Of course, I'm not allowed to use structures, only arrays.

Your guess is as good as mine. Math was never my strong point. I guess one inefficient way would be to create several arrays of char, having each element represent one digit. You could use a fixed offset for the decimal point. Example for representing pi:

Code:

---

char a[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3, 8, 4, 6 };

---

If you were feeling brave, you could do all the math in binary, then convert it to decimal for the output.

Anyway, once you have such arrays, you just need functions for addition and multiplication of them. One way would be to implement the code just like how people used to add and multiply numbers before they gave all the kids calculators.

I'm sure there are other better ways.

it has taken me a really long time to write this program, to find e (it works can you believe it). give it some time and it will do 10000 decimal places correctly (i just looked at the last few numbers assumming the site i was looking at was correct) with 10010 precision.

and thanks to cwr for providing the value of e to so many decimal places- i used it to check my program.

Code:

---

#include <stdio.h>

#define DPLACES  43                        // number of decimal places
#define PRECISION 50                      // the greater the more accurate (last 3 digits are often incorrect).

#define ctod(x) (x-'0')                      // converts character to digit (thanks to cwr)
#define dtoc(x) (x+'0')                      // converts digit to character (thanks to cwr)

int main() {
   
    char ans[DPLACES+3]= {'0', '.'};        // a string that stores 'digit'
    char tem[DPLACES+3]= {'1', '.'};        // another one
    int i, x, y=0, z;                        // counters, except for y
   
    // set the elements to '0' of ans and tem (i = position in array) i=0 and i=1 already defined
    for(i=2; i< DPLACES+2; ++i) {
        tem[i] = ans[i] = '0';
    }
   
    /* not simple (it took me very long to debug the code and i can't believe it all worked out).
      basically tem stores the values of 1/i!. and ans adds them all up to get e.
      using the formula e = 1 +1/i! where i starts from one and approaches PRECISION */
    for (i=1; i< PRECISION; ++i) {
       
        // add up tem and ans right to left (just like doing addition on paper)
        for (x=DPLACES+1; x > -1; --x) {                      // x = position of 'digit' in tem and ans starting from right
            if ((ctod(ans[x]) + ctod(tem[x])) > 9) {          // 1. if 'digit's tem add ans at position x is greater than 9       
                if (ctod(ans[x-1]) != 9) ans[x-1] = dtoc((ctod(ans[x-1])+1));    // a. when the 'digit' after x is not 9 increase that number by one e.g 11+9 = 20 (the 1 to 2) 
                else {                                        // b. when the 'digit'/s after x is 9 do appropriate changes e.g 1999+1 = 2000 (the 199 to the 200 in ans)
                    for (z=x; z > -1 ; --z) {               
                        ans[z-1] = '0';
                        if(ans[z-2] != '9') {
                            ans[z-2] = dtoc(ctod(ans[z-2])+1);
                            break;
                        }   
                    }       
                }       
                ans[x] = dtoc(ctod(ans[x]) + ctod(tem[x])-10); // change 'digit' in ans at position x e.g 9+9=18 (the 9 to an 8 in ans)
            }
            else ans[x] = dtoc(ctod(ans[x]) + ctod(tem[x]));  // 2. if tem add ans at position x is less than 9, ans becomes tem add ans e.g 1+2 = 3 (the 1 to the 3 in ans)
            if (x==2) --x;                                    // skip the decimal point
        }   
       
        y=0; // important and stores the 'remainder'
       
        /* divide 'digit' in tem by i left to right and store the answer in tem. method is similar to how division is done
        on paper. */
        for (x=0; x<DPLACES+2; ++x) {        // x = position of 'digit' in tem starting from left         
            y = y*10 + ctod(tem[x]);        // 'remainder' left from before*10 + 'digit' at position x
            if ( y >= i) {                  // 1. 'remainder' is greater than i 
                tem[x] = dtoc(y/i);          // changes the 'digit' at position x to y/i e.g 5/2 =2 (the 5 to 2)
                y = y - ctod(tem[x])*i;
            }
            else {                          // 2. 'remainder' is less than i   
                tem[x] = '0';                // 'digit' at x becomes '0'
            }   
            if (x==0) ++x;                  // skip the decimal point
        }       
    }   
               
    printf("e = %s\n", ans);
    printf("\ncorrect to about %d decimal places, precision %d\n", DPLACES-3, PRECISION);     
   
    while (getchar() != '\n');
    return 0;
}

---

here is the link to check 10000 decimal places of e: http://www-groups.dcs.st-and.ac.uk/~...s/e_10000.html

peterchen,

Cool, but it's interesting you made those ctod() and dtoc() functions. Replacing them with simple macros:

Code:

---

#define ctod(x) (x-'0')
#define dtoc(x) (x+'0')

---

yields identical behaviour (provided the passed values are in range, but makes the program run in excess of thrice the speed on my machine).

The switch/case statements are unnecessary. Standard C guarantees that the characters that represent '0' through to '9' are contiguous, so a simple addition or subtraction of '0' suffices.

very good point.

i guess i was just worried about the other stuff.

i will edit my code, thanks cwr

Thanks everyone.

As a sidenote, could anyone suggest an algorithm to calculate 2^(-1/2) with more than 18 digits precision? In general, is there a formula for calculating the root of 2?

Quote:

---

Originally Posted by warg

In general, is there a formula for calculating the root of 2?

---

http://en.wikipedia.org/wiki/Methods...g_square_roots

I'd already seen that page. Those methods seemed much too complex to implement. Or if not, could you enlighten me?

i will see what i can do.

it might take some time like the other one.

i think i will probably adopt the method described in the website below - it would be could if i could figure one out myself though
http://www.geocities.com/oosterwal/p...uareroots.html

hope somebody else tries to come up one as well.

well it seems no one has posted one before me.

sorry about the long wait. it is because i've been busy with chinese new year celebrations and also because it is so hard to debug when there is something wrong.

i adopted the method of calculating square roots used in the link i provided in my previous post.

here is the link to check 1000 and more decimal places of sqrt two: http://www.rossi.com/sqr2.htm

suggestions are greatly welcomed.

anyway here it is.

Code:

---

#include <stdio.h>

#define LESS 0
#define MORE 1
#define DPLACES 102                              // actual decimal places is DPLACES-2

void temtsub(int tem[], int sub[]);              // temtsub = tem take sub
int cifless(int sub[], int tem[]);              // cifless = see if less
void mbyx(int sub[], int x, int pos);            // mbyx = multiply by x
void mbyt(int sub[],int ans[], int pos) ;        // mbyt = multiply by twenty

int main() {
    int ans[DPLACES], sub[DPLACES*2], tem[DPLACES*2];
    int i, x, pos; // counter
   
    //preset values
    for (i=0; i< DPLACES; ++i) ans[i] = 0;
    for (i=0; i< DPLACES*2; ++i) sub[i] = tem[i] = 0;
   
    tem[0] = 0;              // finding square root of tem[0]tem[1]
    tem[1] = 2;

    for (pos=0; pos< DPLACES; ++pos) {
        for (x=9; x>-1; --x) {
            mbyt(sub, ans, pos);                // multiply by twenty
            sub[pos*2+1] = x;                    // add x
            mbyx(sub, x, pos);                  // multiply by x
            if (cifless(sub, tem) == MORE) {
                temtsub(tem, sub);
                ans[pos+1] = x;                  // for some reason ans[DPLACES] = x does not produce error
                break;
            }
        }
    }
   
    for (i=1; i < DPLACES; ++i) {
        if(i != 2) printf("%d", ans[i]);
        else {
            putchar('.');
            printf("%d", ans[i]);
        }   
    }   
    printf("\n");
   
    while (getchar() != '\n');
    return 0;
}   

void mbyt(int sub[], int ans[], int pos) {
    int i; // counter
   
    // clear sub
    for (i=0; i < DPLACES*2; ++i) sub[i] = 0;
   
    // multiply ans by 20 and store in sub
    for (i=pos; i > -1; --i) { 
        if (ans[i] * 2 + sub[i+pos] > 9) {
            sub[i+pos] = ans[i] *2 -10 + sub[i+pos];
            sub[i+pos-1] = 1; 
        }
        else {
            sub[i+pos] = ans[i] * 2 + sub[i+pos];
        }       
    }   
} 

void mbyx(int sub[], int x, int pos) {
    int i, carry=0, carry2=0;
   
    for (i=pos*2+1; i>-1; --i) {
        if (sub[i] * x + carry > 9) {
            carry = x*sub[i]/10;            // get number that needs to be carried e.g 9*9 = 81 the 8
            if (sub[i]*x-carry*10+ carry2 < 10) // carried number is less than 10
                sub[i] = sub[i]*x-carry*10 + carry2;
            else {        // carried number is more than ten set sub[i] to appropriate value and next carry goes up
                sub[i] = sub[i]*x-carry*10+ carry2-10;
                ++carry;
            }
            carry2 = carry;
        }
        else {
            sub[i] = sub[i] * x + carry;
            carry = carry2 = 0;
        }   
    }
}       

int cifless(int sub[], int tem[]) {
    int i;
    for (i=0; i< DPLACES*2; ++i) {        // returns 1 when tem is greater then sub
        if ((tem[i] > sub[i])) return MORE;
        else if (tem[i] < sub[i]) return LESS;
        else if (i==DPLACES*2-1) return MORE;
    }     
}
 
void temtsub(int tem[], int sub[]) {
    int i, borrow;

    for (i= DPLACES*2-1; i>-1; --i) {
        // tem - sub equals > = 0
        if (tem[i]- sub[i] >= 0) tem[i] = tem[i] - sub[i];      // do not need to borrow
        else {                                                  // need to 'borrow' from number after(before) it
            if (tem[i-1] != 0) tem[i-1] = tem[i-1] - 1;        // 'digit' borrowed from is not zer0
            else {
                for (borrow = i-1; borrow > -1; --borrow) {    // when it is borrow from next one and so on
                    tem[borrow] = 9;
                    if (tem[borrow-1] != 0) {
                        tem[borrow-1] = tem[borrow-1] -1;
                        break;
                    }
                }
            }
            tem[i] = 10 + tem[i] - sub[i];          // set digit
        }
    }
}

---

Thanks a lot, I will begin deciphering it in the next few days.