calculating e

This is a discussion on calculating e within the C Programming forums, part of the General Programming Boards category; I've written a small code to calculate the constant e with 40 digits precision, but, apparently due to the limitations ...

  1. #1
    Registered User
    Join Date
    Jan 2006
    Posts
    5

    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?

  2. #2
    Just Lurking Dave_Sinkula's Avatar
    Join Date
    Oct 2002
    Posts
    5,006
    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.
    7. It is easier to write an incorrect program than understand a correct one.
    40. There are two ways to write error-free programs; only the third one works.*

  3. #3
    aoeuhtns
    Join Date
    Jul 2005
    Posts
    581
    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.

  4. #4
    cwr
    cwr is offline
    Registered Luser cwr's Avatar
    Join Date
    Jul 2005
    Location
    Sydney, Australia
    Posts
    869
    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");

  5. #5
    Registered User
    Join Date
    Jan 2006
    Posts
    5
    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.

  6. #6
    cwr
    cwr is offline
    Registered Luser cwr's Avatar
    Join Date
    Jul 2005
    Location
    Sydney, Australia
    Posts
    869
    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.

  7. #7
    Registered User
    Join Date
    Dec 2005
    Location
    Australia - Melbourne
    Posts
    63
    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
    Last edited by peterchen; 01-25-2006 at 10:44 PM. Reason: very good suggestions

  8. #8
    cwr
    cwr is offline
    Registered Luser cwr's Avatar
    Join Date
    Jul 2005
    Location
    Sydney, Australia
    Posts
    869
    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.

  9. #9
    Registered User
    Join Date
    Dec 2005
    Location
    Australia - Melbourne
    Posts
    63
    very good point.

    i guess i was just worried about the other stuff.

    i will edit my code, thanks cwr

  10. #10
    Registered User
    Join Date
    Jan 2006
    Posts
    5
    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?

  11. #11
    cwr
    cwr is offline
    Registered Luser cwr's Avatar
    Join Date
    Jul 2005
    Location
    Sydney, Australia
    Posts
    869
    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

  12. #12
    Registered User
    Join Date
    Jan 2006
    Posts
    5
    I'd already seen that page. Those methods seemed much too complex to implement. Or if not, could you enlighten me?

  13. #13
    Registered User
    Join Date
    Dec 2005
    Location
    Australia - Melbourne
    Posts
    63
    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.
    Last edited by peterchen; 01-27-2006 at 06:33 AM.

  14. #14
    Registered User
    Join Date
    Dec 2005
    Location
    Australia - Melbourne
    Posts
    63
    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 
            }
        }
    }
    Last edited by peterchen; 01-29-2006 at 06:48 AM.

  15. #15
    Registered User
    Join Date
    Jan 2006
    Posts
    5
    Thanks a lot, I will begin deciphering it in the next few days.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Calculating : high numbers
    By MiraX33 in forum C++ Programming
    Replies: 9
    Last Post: 06-09-2006, 12:08 AM
  2. Calculating CPU Usage
    By vitaliy in forum Linux Programming
    Replies: 3
    Last Post: 08-21-2005, 10:38 AM
  3. Recursion
    By Lionmane in forum C Programming
    Replies: 11
    Last Post: 06-04-2005, 01:00 AM
  4. Taking input while calculating
    By Unregistered in forum C Programming
    Replies: 1
    Last Post: 07-12-2002, 05:47 PM
  5. Calculating window sizes
    By Mox in forum Windows Programming
    Replies: 3
    Last Post: 11-08-2001, 09:17 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21