e =2.7182818284590455 at 66th iterations with loss of precision after 16th digit when using unsigned long double.

Actual result should be:

__2.7182818284590452__3536... (Wikipedia)

At 67th iteration, overflow occurs for e.

If you use Windows Calculator:

QUOTE

1/64! = 1/1.2688693218588416410343338933516e+89

1/65! = 1/8.2476505920824706667231703067855e+90

But, when compiled using Visual C++ 2008, this program produces:

QUOTE

+ 1/64! = e[64] = 1/9.2233720368547758000000000000000e+018

+ 1/65! = e[65] = 1/9.2233720368547758000000000000000e+018

VC++ gives an erroneous result starting at 1/21!

QUOTE

1/21! = 1/1.4197454024290337000000000000000e+019

whereas Windows Calculator displays:

QUOTE

1/21! = 1/5.1090942171709440000e+019

The maximum limit in VC++:

QUOTE

unsigned long long (unsigned __int64) -> 18,446,744,073,709,551,615 or 1.8446744073709551615e+019

Interestingly, the binary form of IEEE-754 (1985) for the erroneous 21! = 14197454024290336768 in double-precision 64-bit format for VC++:

CODE

VC++ = 1 10001010000 0111011111010011011010111000110001000000000000000000 (0xC5077D36B8C40000)

The correct 21! = 51090942171709440000 in Windows Calculator:

CODE

WCalc = 1 10001010000 0111011111010011011010111000110001000000000000000000 (0xC5077D36B8C40000)

Limit = 1 11111111111 1111111111111111111111111111111111111111111111111111 (0xFFFFFFFFFFFFFFFF)

~~So, there is a bug in VC++ 2008 when interpreting IEEE-754 for 21! in decimal form, wheres it is correct in binary form!~~Ok, the problem lies when passing variable from unsigned long long to unsigned long double; it causes IEEE-754 truncation!

To correct this, edit the original factorial() function as follows by making p and return call as unsigned long double:

CODE

ud factorial(ui n)

{

ud p = 1;

//note: 0! = 1, 1! = 1

for(ui i = 1; i > 0 && i <= n; i++) { p *= i;}

return (p);

}

With the correction above, it is correct that factorial() overflows at 171st iteration!

Original code:

CODE

#include <iostream>

#include <iomanip>

#define ui unsigned long long

#define ud unsigned long double

using namespace std;

ui factorial(ui n)

{

ui p = 1;

//note: 0! = 1, 1! = 1

for(ui i = 1; i > 0 && i <= n; i++) { p *= i;}

return (p);

}

//e = 1/0! + 1/1! + 1/2! + 1/3! +...+1/k!

int main()

{

ui k = 66;

ud e = 0.00;

cout << setiosflags(ios_base::scientific) << setprecision(31);

for(ui i=0; i < k; i++)

{

e += 1/(static_cast<ud>(factorial(i)));

//debug purpose to verify correctness of factorials

cout << "+ 1/"<< i << "! = e[" << i << "] = 1/";

cout << static_cast<ud>(factorial(i)) << "\n";

}

//e =2.7182818284590455 at 66th iteration

cout << "Final result:" << e <<"\n";

return 0;

}

Original code output:

__» Click to show Spoiler - click again to hide... «__

+ 1/0! = e[0] = 1/1.0000000000000000000000000000000e+000

+ 1/1! = e[1] = 1/1.0000000000000000000000000000000e+000

+ 1/2! = e[2] = 1/2.0000000000000000000000000000000e+000

+ 1/3! = e[3] = 1/6.0000000000000000000000000000000e+000

+ 1/4! = e[4] = 1/2.4000000000000000000000000000000e+001

+ 1/5! = e[5] = 1/1.2000000000000000000000000000000e+002

+ 1/6! = e[6] = 1/7.2000000000000000000000000000000e+002

+ 1/7! = e[7] = 1/5.0400000000000000000000000000000e+003

+ 1/8! = e[8] = 1/4.0320000000000000000000000000000e+004

+ 1/9! = e[9] = 1/3.6288000000000000000000000000000e+005

+ 1/10! = e[10] = 1/3.6288000000000000000000000000000e+006

+ 1/11! = e[11] = 1/3.9916800000000000000000000000000e+007

+ 1/12! = e[12] = 1/4.7900160000000000000000000000000e+008

+ 1/13! = e[13] = 1/6.2270208000000000000000000000000e+009

+ 1/14! = e[14] = 1/8.7178291200000000000000000000000e+010

+ 1/15! = e[15] = 1/1.3076743680000000000000000000000e+012

+ 1/16! = e[16] = 1/2.0922789888000000000000000000000e+013

+ 1/17! = e[17] = 1/3.5568742809600000000000000000000e+014

+ 1/18! = e[18] = 1/6.4023737057280000000000000000000e+015

+ 1/19! = e[19] = 1/1.2164510040883200000000000000000e+017

+ 1/20! = e[20] = 1/2.4329020081766400000000000000000e+018

+ 1/21! = e[21] = 1/1.4197454024290337000000000000000e+019

+ 1/22! = e[22] = 1/1.7196083355034583000000000000000e+019

+ 1/23! = e[23] = 1/8.1282916178948260000000000000000e+018

+ 1/24! = e[24] = 1/1.0611558092380307000000000000000e+019

+ 1/25! = e[25] = 1/7.0345352775739638000000000000000e+018

+ 1/26! = e[26] = 1/1.6877220553537094000000000000000e+019

+ 1/27! = e[27] = 1/1.2963097176472289000000000000000e+019

+ 1/28! = e[28] = 1/1.2478583540742619000000000000000e+019

+ 1/29! = e[29] = 1/1.1390785281054474000000000000000e+019

+ 1/30! = e[30] = 1/9.6821651048622981000000000000000e+018

+ 1/31! = e[31] = 1/4.9992130713784156000000000000000e+018

+ 1/32! = e[32] = 1/1.2400865694432887000000000000000e+019

+ 1/33! = e[33] = 1/3.4001982946751283000000000000000e+018

+ 1/34! = e[34] = 1/4.9262775766970532000000000000000e+018

+ 1/35! = e[35] = 1/6.3990185210108969000000000000000e+018

+ 1/36! = e[36] = 1/9.0037378718776689000000000000000e+018

+ 1/37! = e[37] = 1/1.0969079327018189000000000000000e+018

+ 1/38! = e[38] = 1/4.7890132952500142000000000000000e+018

+ 1/39! = e[39] = 1/2.3040777776550380000000000000000e+018

+ 1/40! = e[40] = 1/1.8376134811363312000000000000000e+019

+ 1/41! = e[41] = 1/1.5551764317513712000000000000000e+019

+ 1/42! = e[42] = 1/7.5380587557415813000000000000000e+018

+ 1/43! = e[43] = 1/1.0541877243825619000000000000000e+019

+ 1/44! = e[44] = 1/2.6739968855884431000000000000000e+018

+ 1/45! = e[45] = 1/9.6493954092226314000000000000000e+018

+ 1/46! = e[46] = 1/1.1503310552118067000000000000000e+018

+ 1/47! = e[47] = 1/1.7172071447535813000000000000000e+019

+ 1/48! = e[48] = 1/1.2602690238498734000000000000000e+019

+ 1/49! = e[49] = 1/8.7892672540227666000000000000000e+018

+ 1/50! = e[50] = 1/1.5188249005818642000000000000000e+019

+ 1/51! = e[51] = 1/1.8284192274659148000000000000000e+019

+ 1/52! = e[52] = 1/9.9940505230885519000000000000000e+018

+ 1/53! = e[53] = 1/1.3175843659825807000000000000000e+019

+ 1/54! = e[54] = 1/1.0519282829630636000000000000000e+019

+ 1/55! = e[55] = 1/6.7114893446888817000000000000000e+018

+ 1/56! = e[56] = 1/6.9085218283863409000000000000000e+018

+ 1/57! = e[57] = 1/6.4041186701208453000000000000000e+018

+ 1/58! = e[58] = 1/2.5040013928179958000000000000000e+018

+ 1/59! = e[59] = 1/1.6212958658533786000000000000000e+017

+ 1/60! = e[60] = 1/9.7277751951202714000000000000000e+018

+ 1/61! = e[61] = 1/3.0984765436309012000000000000000e+018

+ 1/62! = e[62] = 1/7.6381049680203612000000000000000e+018

+ 1/63! = e[63] = 1/1.5852670688344146000000000000000e+018

+ 1/64! = e[64] = 1/9.2233720368547758000000000000000e+018

+ 1/65! = e[65] = 1/9.2233720368547758000000000000000e+018

Final result:2.7182818284590455000000000000000e+000

This post has been edited by **eclectice**: Aug 31 2009, 05:21 PM