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

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

C: Sum not converging to e as it should.

Postby Blatm » Wed Nov 17, 2010 3:00 am UTC

I have this formula for e:
[math]e = \frac{1 + \sqrt{1 + \alpha\beta}}{\alpha}[/math]
where
[math]\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!}[/math]
and
[math]\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!}[/math]

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.

Postby Sagekilla » Wed Nov 17, 2010 3:13 am UTC

Have you thought perhaps about your limited precision?

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.

Postby masher » Wed Nov 17, 2010 3:15 am UTC

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.

Postby Sagekilla » Wed Nov 17, 2010 3:17 am UTC

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.

Postby masher » Wed Nov 17, 2010 3:40 am UTC

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.

Postby Sagekilla » Wed Nov 17, 2010 3:44 am UTC

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.

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

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

Postby hotaru » Wed Nov 17, 2010 6:51 am UTC

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.

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)); }

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.

Postby Sagekilla » Wed Nov 17, 2010 3:50 pm UTC

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.

User avatar
Thesh
Made to Fuck Dinosaurs
Posts: 6598
Joined: Tue Jan 12, 2010 1:55 am UTC
Location: Colorado

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

Postby Thesh » Wed Nov 17, 2010 4:36 pm UTC

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.

Postby gorcee » Wed Nov 17, 2010 4:56 pm UTC

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

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

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

Postby MHD » Wed Nov 17, 2010 7:40 pm UTC

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"

User avatar
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.

Postby SWGlassPit » Thu Nov 18, 2010 12:21 am UTC

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.
Image
Erdös number: 5

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

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

Postby Blatm » Thu Nov 18, 2010 2:00 am UTC

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.

Postby Sagekilla » Thu Nov 18, 2010 2:03 am UTC

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.

User avatar
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.

Postby phlip » Thu Nov 18, 2010 4:12 am UTC

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]

User avatar
Thesh
Made to Fuck Dinosaurs
Posts: 6598
Joined: Tue Jan 12, 2010 1:55 am UTC
Location: Colorado

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

Postby Thesh » Thu Nov 18, 2010 6:27 am UTC

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.

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

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

Postby hotaru » Thu Nov 18, 2010 6:53 am UTC

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 1
isPrime n 
factorial (1) `mod== 1

User avatar
Thesh
Made to Fuck Dinosaurs
Posts: 6598
Joined: Tue Jan 12, 2010 1:55 am UTC
Location: Colorado

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

Postby Thesh » Thu Nov 18, 2010 6:57 am UTC

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.

Postby Goplat » Thu Nov 18, 2010 6:28 pm UTC

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 whatever
void 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.

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

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

Postby hotaru » Fri Nov 19, 2010 8:51 am UTC

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 1
isPrime n 
factorial (1) `mod== 1

User avatar
Thesh
Made to Fuck Dinosaurs
Posts: 6598
Joined: Tue Jan 12, 2010 1:55 am UTC
Location: Colorado

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

Postby Thesh » Fri Nov 19, 2010 9:25 am UTC

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.

User avatar
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.

Postby thoughtfully » Fri Nov 19, 2010 9:30 pm UTC

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.
Image
Perfection is achieved, not when there is nothing more to add, but when there is nothing left to take away.
-- Antoine de Saint-Exupery

User avatar
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.

Postby Berengal » Sun Nov 21, 2010 2:12 am UTC

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.

User avatar
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.

Postby thoughtfully » Thu Nov 25, 2010 3:22 am UTC

Never misses a thing, Berengal.

Code: Select all

import operator

def 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 25
15511210043330985984000000
15511210043330986055303168.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.
Image
Perfection is achieved, not when there is nothing more to add, but when there is nothing left to take away.
-- Antoine de Saint-Exupery


Return to “Coding”

Who is online

Users browsing this forum: No registered users and 7 guests