## C: Sum not converging to e as it should.

A place to discuss the implementation and style of computer programs.

Moderators: phlip, Moderators General, Prelates

Blatm
Posts: 638
Joined: Mon Jun 04, 2007 1:43 am UTC

### C: Sum not converging to e as it should.

I have this formula for e:
$e = \frac{1 + \sqrt{1 + \alpha\beta}}{\alpha}$
where
$\alpha = \sum_{n=0}^{\infty}\sum_{i=n}^{2n}\sum_{k=0}^{2(i-n)} (-1)^{n+i+k}\frac{ (2n)!(2i)!(4n+1)}{16^n (2n-i)!(n!)^2 i! k!}$
and
$\beta = \sum_{n=0}^{\infty}\sum_{i=n}^{2n}\sum_{k=0}^{2(i-n)} (-1)^{n+i}\frac{ (2n)!(2i)!(4n+1)}{16^n (2n-i)!(n!)^2 i! k!}$

I try to implement this using some code. For n = 0..2, I get 2.718281 57... (compared to the correct value of 2.718281 83...), which is about as good as I could have hoped for. However, for n = 0..3, I get 2.835..., and for higher n the number I get quickly converges to 2.83424495... While it's certainly possible that my formula is wrong, or that I implemented it wrong, the fact that I get such a good result for n = 0..2 makes me wonder if this isn't a problem with c handling numbers. Can someone please look at my code and see if there might be such problems somewhere, and if there's anything could do to fix them? Thanks.

Code: Select all

#include <stdio.h>#include <stdlib.h>#include <string.h>#include <time.h>#include <math.h>int fact (int n){   int prod;      if(n!=0)   {      prod = n*fact(n-1);   }      else   {      prod = 1;   }         return prod;}int main (void){   int nmax = 15;   int imin, imax, kmin, kmax;   int n = 0, i = 0, k = 0;   double term, terma, termb;   double alpha = 0, beta = 0;   double e;   double moredigits;      for(n = 0; n < nmax + 1; n++)   {      imin = n;      imax = 2*n;            for(i = imin; i < imax + 1; i++)      {         kmin = 0;         kmax = 2*(i-n);                  for(k = kmin; k < kmax + 1; k++)         {            term = (4*n + 1)*fact(2*n)*fact(2*i) / ( pow(16,n)*pow(fact(n),2)*fact(2*n-i)*fact(i)*fact(k) );            terma = term*pow((-1), (n+i+k)%2);            termb = term*pow((-1), (n+i)%2);                        alpha += terma;            beta += termb;         }      }   }      e = (1 + pow((1+alpha*beta), 0.5))/alpha;   /*printf("%f ", e);*/   moredigits = (e*1000000);   printf("%f", moredigits);}

(My apologies for the lack of documentation. Let me know if there's anything you want clarified.)

Sagekilla
Posts: 382
Joined: Fri Aug 21, 2009 1:02 am UTC
Location: Long Island, NY

### Re: C: Sum not converging to e as it should.

Floating point numbers can't represent real numbers perfectly, so you could be getting round off errors. I have to take a close look at the code though.
It's also very possible you could be getting overflow errors too.

Edit: You realize you're doing integer division right? Do (4 * n + 1.0) not (4 * n + 1). This makes sure your expression is a floating point number.
Last edited by Sagekilla on Wed Nov 17, 2010 3:15 am UTC, edited 1 time in total.
http://en.wikipedia.org/wiki/DSV_Alvin#Sinking wrote:Researchers found a cheese sandwich which exhibited no visible signs of decomposition, and was in fact eaten.

masher
Posts: 821
Joined: Tue Oct 23, 2007 11:07 pm UTC
Location: Melbourne, Australia

### Re: C: Sum not converging to e as it should.

Code: Select all

          terma = term*pow((-1), (n+i+k)%2);            termb = term*pow((-1), (n+i)%2);

Why are the "%2"'s there? They aren't in your original formulae

Sagekilla
Posts: 382
Joined: Fri Aug 21, 2009 1:02 am UTC
Location: Long Island, NY

### Re: C: Sum not converging to e as it should.

masher wrote:

Code: Select all

          terma = term*pow((-1), (n+i+k)%2);            termb = term*pow((-1), (n+i)%2);

Why are the "%2"'s there? They aren't in your original formulae

He's trying to do (-1)^n:

If n is even, then (-1)^n is +1.
If n is odd, then (-1)^n is -1.

If you do (-1)^(n % 2), you're essentially doing the same exact thing. If n is even, you get +1. If n is odd, you get -1.

As a quick check, the smallest factorial you can fit in a signed integer (32-bits) is 12!, with a value of 479,001,600.
You are, without a doubt, overflowing.

Try rewriting it but only using ints for indices, doubles for everything else. Use the following factorial function:

Code: Select all

double fact(int n){   double m = 1;   for (int i = 1; i < n+1; i++)      m *= i;   return m;}

This version is more numerically stable for larger n.
Last edited by Sagekilla on Wed Nov 17, 2010 3:42 pm UTC, edited 1 time in total.
http://en.wikipedia.org/wiki/DSV_Alvin#Sinking wrote:Researchers found a cheese sandwich which exhibited no visible signs of decomposition, and was in fact eaten.

masher
Posts: 821
Joined: Tue Oct 23, 2007 11:07 pm UTC
Location: Melbourne, Australia

### Re: C: Sum not converging to e as it should.

Sagekilla wrote:As a quick check, the smallest factorial you can fit in a signed integer (32-bits) is 12!, with a value of 479,001,600.
You are, without a doubt, overflowing.

Yes.

Just tested with ints and longs (in java) and I get quite a difference, except that with longs, I've started calcing -e....

Sagekilla
Posts: 382
Joined: Fri Aug 21, 2009 1:02 am UTC
Location: Long Island, NY

### Re: C: Sum not converging to e as it should.

masher wrote:
Sagekilla wrote:As a quick check, the smallest factorial you can fit in a signed integer (32-bits) is 12!, with a value of 479,001,600.
You are, without a doubt, overflowing.

Yes.

Just tested with ints and longs (in java) and I get quite a difference, except that with longs, I've started calcing -e....

You're probably overflowing in just the right way that it becomes positive again.
Try using the factorial function I posted above.

Also, try calculating your terms like this:

Code: Select all

double a = fact(2 * n) / fact(2 * n - i);double b = fact(2 * i) / (fact(i) * pow(16, n));double c = (4 * n + 1) / (fact(i) * fact(k));double d = pow(fact(n), 2);double term = (a / d) * (b * c);

Depending on what your compiler does, that may be more stable.
http://en.wikipedia.org/wiki/DSV_Alvin#Sinking wrote:Researchers found a cheese sandwich which exhibited no visible signs of decomposition, and was in fact eaten.

hotaru
Posts: 1045
Joined: Fri Apr 13, 2007 6:54 pm UTC

### Re: C: Sum not converging to e as it should.

Sagekilla wrote:Use the following factorial function:

Code: Select all

double fact(int n){   double m = 1;   for (int i = 0; i < n+1; i++)      m *= i;   return m;}

This version is more numerically stable for larger n.

Code: Select all

double fact(int n){ double a = n + 0.5;  return round(sqrt(2 * M_PI) * exp(a * (log(a) - 1)) *               pow(a /                   (a + 1.0 / 24 /                    (a + 3.0 / 80 /                     (a + 18029.0 / 45360 /                      (a + 6272051.0 / 14869008 / a)))), a)); }

Sagekilla
Posts: 382
Joined: Fri Aug 21, 2009 1:02 am UTC
Location: Long Island, NY

### Re: C: Sum not converging to e as it should.

hotaru wrote:i'd recommend this one instead:

Code: Select all

double fact(int n){ double a = n + 0.5;  return round(sqrt(2 * M_PI) * exp(a * (log(a) - 1)) *               pow(a /                   (a + 1.0 / 24 /                    (a + 3.0 / 80 /                     (a + 18029.0 / 45360 /                      (a + 6272051.0 / 14869008 / a)))), a)); }

I fixed my above snippet (used the wrong lower bound by accident).
I tried out your version and it's actually a little bit less accurate once you get to around n = 17. (read: very very small difference)

If concerned with speed, I would just calculate a table of factorials before hand using the
most accurate version you can use. There's an upper limit to the number of terms that series
would produce any meaningful results in floating point numbers, so it wouldn't be a very large table.
http://en.wikipedia.org/wiki/DSV_Alvin#Sinking wrote:Researchers found a cheese sandwich which exhibited no visible signs of decomposition, and was in fact eaten.

Thesh
Posts: 6598
Joined: Tue Jan 12, 2010 1:55 am UTC

### Re: C: Sum not converging to e as it should.

For the factorial function, I would cache the previous values if you are planning on doing a lot of calculations.

e.g.

Code: Select all

double fact(int n){   static double a[100] = {1,1};   static int i = 2;   if (n<i)      return a[n];   while (i <= n)      a[i++] = a[i-1] * i;   return a[n];}
Summum ius, summa iniuria.

gorcee
Posts: 1501
Joined: Sun Jul 13, 2008 3:14 am UTC

### Re: C: Sum not converging to e as it should.

Thesh wrote:For the factorial function, I would cache the previous values if you are planning on doing a lot of calculations.

To add some detail, this process is called "memoization." There's a old but good page that describes it well: http://perl.plover.com/MiniMemoize/memoize.html

MHD
Posts: 630
Joined: Fri Mar 20, 2009 8:21 pm UTC
Location: Denmark

### Re: C: Sum not converging to e as it should.

One way to get around this is to use some extended precision integer format.
EvanED wrote:be aware that when most people say "regular expression" they really mean "something that is almost, but not quite, entirely unlike a regular expression"

SWGlassPit
Posts: 312
Joined: Mon Feb 18, 2008 9:34 pm UTC
Location: Houston, TX
Contact:

### Re: C: Sum not converging to e as it should.

Check out something like the GNU MultiPrecision library
Up in space is a laboratory the size of a football field zipping along at 7 km/s. It's my job to keep it safe.

Erdös number: 5

Blatm
Posts: 638
Joined: Mon Jun 04, 2007 1:43 am UTC

### Re: C: Sum not converging to e as it should.

Major thanks guys. Using Thesh's factorial function everything is working smoothly. However, I have another problem now: I can't get enough digits to display. With nmax = 4, all 16 digits I can get the program to display as correct, which isn't especially helpful. Is there any way to get more than 16 digits of precision in the output? (The program is the same as in the OP, except using Thesh's factorial function).

EDIT: Derp, I completely missed the post right above me. I'll check that out first.

Sagekilla
Posts: 382
Joined: Fri Aug 21, 2009 1:02 am UTC
Location: Long Island, NY

### Re: C: Sum not converging to e as it should.

Do you know the details of floating point numbers? You should read up a bit on them, it's pretty enlightening stuff.

The reason why you can only get 16 decimal places is because that's the maximum a double precision number can represent.

Any more than that and you'll need to use higher precision numbers (through GMP, for example).
http://en.wikipedia.org/wiki/DSV_Alvin#Sinking wrote:Researchers found a cheese sandwich which exhibited no visible signs of decomposition, and was in fact eaten.

phlip
Restorer of Worlds
Posts: 7573
Joined: Sat Sep 23, 2006 3:56 am UTC
Location: Australia
Contact:

### Re: C: Sum not converging to e as it should.

Thesh wrote:

Code: Select all

a[i++] = a[i-1] * i;

This may not do what you want... the standard doesn't specify what order the parts of this are evaluated in, so it could evaluate the i++ first, in which case it would all fall over. Run it on a DeathStation 9000 and you'll probably kill somebody.

Code: Select all

enum ಠ_ಠ {°□°╰=1, °Д°╰, ಠ益ಠ╰};void ┻━┻︵​╰(ಠ_ಠ ⚠) {exit((int)⚠);}
[he/him/his]

Thesh
Posts: 6598
Joined: Tue Jan 12, 2010 1:55 am UTC

### Re: C: Sum not converging to e as it should.

Hmm... It seems to work in the compilers I have used, but looking at the standard it appears you are right. Oh well, guess I shouldn't try so hard to save keystrokes.

Code: Select all

double fact(int n){   static double a[100] = {1,1};   static int i = 2;   if (n<i)      return a[n];   while (i <= n) {      a[i] = a[i-1] * i;      i++;   }   return a[n];}
Summum ius, summa iniuria.

hotaru
Posts: 1045
Joined: Fri Apr 13, 2007 6:54 pm UTC

### Re: C: Sum not converging to e as it should.

Thesh wrote:

Code: Select all

double fact(int n){   static double a[100] = {1,1};   static int i = 2;   if (n<i)      return a[n];   while (i <= n) {      a[i] = a[i-1] * i;      i++;   }   return a[n];}

just make sure no one ever does fact(100)...

Code: Select all

factorial = product . enumFromTo 1isPrime n = factorial (n - 1) mod n == n - 1

Thesh
Posts: 6598
Joined: Tue Jan 12, 2010 1:55 am UTC

### Re: C: Sum not converging to e as it should.

hotaru wrote:just make sure no one ever does fact(100)...

Yes, it isn't a generic factorial function and has a bounds limit, but you can do this if it makes you happy:

Code: Select all

double fact(int n){   static double a[171] = {1,1};   static int i = 2;   if (n<i)      return a[n];   while (i <= n) {      a[i] = a[i-1] * i;      i++;   }   return a[n];}

You exceed the range of the double beyond that.
Summum ius, summa iniuria.

Goplat
Posts: 490
Joined: Sun Mar 04, 2007 11:41 pm UTC

### Re: C: Sum not converging to e as it should.

Memoization is overkill when you can compute every single value in under a microsecond anyway. If you're actually doing a lot of factorials, I'd say just precompute the whole table.

Code: Select all

static double a[100]; // or 171 or whatevervoid fact_init() {   a[0] = 1;   for (int i = 1; i < sizeof(a) / sizeof(*a); i++)      a[i] = a[i-1] * i;}double fact(int n) {   return a[n];}

fact() doesn't get any faster than that, plus you don't have to worry about race conditions if the program is multithreaded.

hotaru
Posts: 1045
Joined: Fri Apr 13, 2007 6:54 pm UTC

### Re: C: Sum not converging to e as it should.

Thesh wrote:

Code: Select all

double fact(int n){   static double a[171] = {1,1};   static int i = 2;   if (n<i)      return a[n];   while (i <= n) {      a[i] = a[i-1] * i;      i++;   }   return a[n];}

ISO/IEC 9899:1999 wrote:1 The behavior is undefined in the following circumstances:
...
— An array subscript is out of range, even if an object is apparently accessible with the
given subscript (as in the lvalue expression a[1][7] given the declaration int
a[4][5]) (6.5.6).

that means that if someone calls your fact with any value less than 0 or greater than 170, anything can happen.

Thesh wrote:You exceed the range of the double beyond that.

that's not guaranteed. the c standard just says that a double has to be able to hold numbers up to at least 1E+37.

Goplat wrote:fact() doesn't get any faster than that,

only if you're computing a lot of factorials... and especially if you're going to want all the factorials up to a certain number. but if you only want a few large ones, it can be made a lot faster than that.

Code: Select all

factorial = product . enumFromTo 1isPrime n = factorial (n - 1) mod n == n - 1

Thesh
Posts: 6598
Joined: Tue Jan 12, 2010 1:55 am UTC

### Re: C: Sum not converging to e as it should.

If you were making an actual application, I would agree with you. If you are writing a program just to do a calculation, it's a waste of time to really bother with the bounds checking on something like that.
Summum ius, summa iniuria.

thoughtfully
Posts: 2253
Joined: Thu Nov 01, 2007 12:25 am UTC
Location: Minneapolis, MN
Contact:

### Re: C: Sum not converging to e as it should.

hotaru wrote:
Thesh wrote:You exceed the range of the double beyond that.

that's not guaranteed. the c standard just says that a double has to be able to hold numbers up to at least 1E+37.

The actual limit of a IEEE 754 double precision value is quite a lot higher, about 1.7976931348623157e308, which will fit 171! but not 172!, so he wasn't using the value from the C standard.

Perfection is achieved, not when there is nothing more to add, but when there is nothing left to take away.
-- Antoine de Saint-Exupery

Berengal
Superabacus Mystic of the First Rank
Posts: 2707
Joined: Thu May 24, 2007 5:51 am UTC
Location: Bergen, Norway
Contact:

### Re: C: Sum not converging to e as it should.

You'll have issues with precision long before that. You can't represent 25! using a double. If you try to calculate it, the result will depend on how you do that calculation.
It is practically impossible to teach good programming to students who are motivated by money: As potential programmers they are mentally mutilated beyond hope of regeneration.

thoughtfully
Posts: 2253
Joined: Thu Nov 01, 2007 12:25 am UTC
Location: Minneapolis, MN
Contact:

### Re: C: Sum not converging to e as it should.

Never misses a thing, Berengal.

Code: Select all

import operatordef factorial(n):   return reduce(operator.mul, range(1,n+1))if __name__ == '__main__':   import sys   print factorial(int(sys.argv[1]))

Code: Select all

#include <stdio.h>#include <stdlib.h>int main(int argc, char **argv){   double answer=1.0;   int i;   for (i=1; i<=strtol(argv[1], NULL, 10); ++i) answer *= i;   printf("%lf\n", answer);   return 0;}

Code: Select all

cfuller@tinyfugue:~/tmp\$ python factorial.py 25;./factorial 251551121004333098598400000015511210043330986055303168.000000

I get matches up to 22! using these codes.

I feel a little silly for not realizing this, having read the specification I linked to. There's about 16 decimal digits of precision in a double, whereas 25! uses 19. I was interested in this fact for another application too, so I paid attention to it

Another thing to consider is how well the printf output routines reflect the internal representation. Doesn't seem to be any problem in this case, though.

Perfection is achieved, not when there is nothing more to add, but when there is nothing left to take away.
-- Antoine de Saint-Exupery