Large numbers, floating point number questions

C

Chad

The following program is my attempt to calculate the formula at the
following url...

http://i23.photobucket.com/albums/b363/CompressorX/1.jpg

#include <stdio.h>
#include <stdlib.h>

unsigned long fact(unsigned long);
unsigned long sum_fact(unsigned long);
float harmonic(unsigned long);

unsigned long fact(unsigned long n)
{
unsigned long acc = 1;
unsigned long i;

for (i = 1; i <= n; i++)
acc *= i;

return acc;
}

unsigned long sum_fact(unsigned long n)
{
unsigned long acc = 0;
unsigned long i;

for (i = 1; i <= n; i++) {
acc += fact(i);
}

return acc;
}

float harmonic(unsigned long n)
{
float acc = 0;
float i;

for (i = 1; i <= n; i++) {
acc += ((100) *1/i);
}

return acc;
}

int main(int argc, char **argv)
{
unsigned long s;

if (argc != 2) {
fprintf(stderr, "Invalid number of args\n");
return 1;
}

s = atoi(argv[1]);

printf("The final value is: %f\n", sum_fact(s)/(harmonic(s)/100));

return 0;
}


The problem is that the code seems to break around n=6 or n=7. How
would I go about calculating stuff up to say the number 100.
 
E

Eric Sosman

The following program is my attempt to calculate the formula at the
following url...

http://i23.photobucket.com/albums/b363/CompressorX/1.jpg

#include<stdio.h>
#include<stdlib.h>

unsigned long fact(unsigned long);
unsigned long sum_fact(unsigned long);
float harmonic(unsigned long);

unsigned long fact(unsigned long n)
{
unsigned long acc = 1;
unsigned long i;

for (i = 1; i<= n; i++)
acc *= i;

return acc;
}

You say you're interested in values of n up to "say the
number 100." Are you aware that 100! is approximately 1E158?
You'll need an `unsigned long' of about 525 bits or so to
get that high ...
unsigned long sum_fact(unsigned long n)
{
unsigned long acc = 0;
unsigned long i;

for (i = 1; i<= n; i++) {
acc += fact(i);
}

return acc;
}

... and a few more bits by the time you've added a bunch
of those enormous quantities together ...
[...]
The problem is that the code seems to break around n=6 or n=7. How
would I go about calculating stuff up to say the number 100.

I haven't studied the code to find why it "breaks," because
the fundamental approach simply isn't going to work with the
native types on today's machines. Find a different formula, or
find a multiple-precision arithmetic package.
 
C

Chad

The following program is my attempt to calculate the formula at the
following url...


unsigned long fact(unsigned long);
unsigned long sum_fact(unsigned long);
float harmonic(unsigned long);
unsigned long fact(unsigned long n)
{
   unsigned long acc = 1;
   unsigned long i;
   for (i = 1; i<= n; i++)
     acc *= i;
   return acc;
}

     You say you're interested in values of n up to "say the
number 100."  Are you aware that 100! is approximately 1E158?
You'll need an `unsigned long' of about 525 bits or so to
get that high ...
unsigned long sum_fact(unsigned long n)
{
   unsigned long acc = 0;
   unsigned long i;
   for (i = 1; i<= n; i++) {
     acc += fact(i);
   }
   return acc;
}

     ... and a few more bits by the time you've added a bunch
of those enormous quantities together ...
[...]
The problem is that the code seems to break around n=6 or n=7. How
would I go about calculating stuff up to say the number 100.

     I haven't studied the code to find why it "breaks," because
the fundamental approach simply isn't going to work with the
native types on today's machines.  Find a different formula, or
find a multiple-precision arithmetic package.

What other formula do have in mind?

Chad
 
B

Ben Bacarisse

Chad said:
The following program is my attempt to calculate the formula at the
following url...

http://i23.photobucket.com/albums/b363/CompressorX/1.jpg

#include <stdio.h>
#include <stdlib.h>

unsigned long fact(unsigned long);
unsigned long sum_fact(unsigned long);
float harmonic(unsigned long);

unsigned long fact(unsigned long n)
{
unsigned long acc = 1;
unsigned long i;

for (i = 1; i <= n; i++)
acc *= i;

return acc;
}

unsigned long sum_fact(unsigned long n)
{
unsigned long acc = 0;
unsigned long i;

for (i = 1; i <= n; i++) {
acc += fact(i);
}

return acc;
}

float harmonic(unsigned long n)
{
float acc = 0;
float i;

for (i = 1; i <= n; i++) {
acc += ((100) *1/i);
}

Using a float as a loop index is not a good idea.
return acc;
}

int main(int argc, char **argv)
{
unsigned long s;

if (argc != 2) {
fprintf(stderr, "Invalid number of args\n");
return 1;
}

s = atoi(argv[1]);

printf("The final value is: %f\n", sum_fact(s)/(harmonic(s)/100));

return 0;
}


The problem is that the code seems to break around n=6 or n=7. How
would I go about calculating stuff up to say the number 100.

What accuracy do you need? S(100) is exactly

262898806796746281676079930191080203927243930357368966202376910097033877607478900552047729559394146124839310937700464163196753332648214687830620489578849276792933475322073718652663784818718016211136 / 14466636279520351160221518043104131447711

(that was a Haskell program) or about

18172766752207502666566054742592774333334408288875508533667627526462215026676578187355017535365067611072187901276460031983419028727944397925589152640251401379.88848916234112956007944544952276803589241196007806492632486867816494303869269492189114953628059466014191088986813481278364231972207711080807475859000106626980248937302365258665295763512500573893636009

but you can get a result of 1.81727667522075e+157 in plain C by using
double precision and being careful about how you do the sum. If you
need more accurate results you will have to use some sort of extended
precision float or bignum library.
 
C

Chad

Why not use float (or, better, double) for fact() and sum_fact() too? Then
you won't get the overflow problems so early.

Okay, I don't get how using float (or double) for fact() and
sum_fact() would prevent overflow problems early on.
 
C

Chad

Using a float as a loop index is not a good idea.

What's the alternative?
  return acc;
}
int main(int argc, char **argv)
{
  unsigned long s;
  if (argc != 2) {
    fprintf(stderr, "Invalid number of args\n");
    return 1;
  }
  s = atoi(argv[1]);
  printf("The final value is: %f\n", sum_fact(s)/(harmonic(s)/100));
  return 0;
}
The problem is that the code seems to break around n=6 or n=7. How
would I go about calculating stuff up to say the number 100.

What accuracy do you need?  S(100) is exactly

Maybe around S(20).
 
N

Nobody

The following program is my attempt to calculate the formula at the
following url...

http://i23.photobucket.com/albums/b363/CompressorX/1.jpg
The problem is that the code seems to break around n=6 or n=7. How
would I go about calculating stuff up to say the number 100.

If you don't need an exact answer, just use "double" throughout; for
n=100, that gives a result which is accurate to 14 decimal digits. If you
need an exact answer, use a multi-precision arithmetic library which
supports rational arithmatic, e.g. GMP:

http://gmplib.org/

Or use a high-level language with built-in bignum support, e.g. Python,
Haskell.
 
B

bartc

Chad said:
Okay, I don't get how using float (or double) for fact() and
sum_fact() would prevent overflow problems early on.

If your 'unsigned long' is, say, 32-bits, then the largest factorial it can
store is 12!

13! would require about 33 bits I think, and if you try and do it using 32
bits, it will give funny results.

Summing the factorials has similar problems; integers just have a too narrow
range.

Floats have a much wider range: typically a float might go up to around
1e38, and a double up to about 1e300. But a double only gives some 16
digits of precision.

Simplify your program to *only* show factorials from 1 to 100 (forget
command line input, just use a loop), and you will
see where it goes wrong. Then put in doubles, and see the difference.
 
B

Ben Bacarisse

Chad said:
What's the alternative?

An integer type is preferable. I don't think there is anything actually
wrong with the above as far as C is concerned (because of the guarantees
C makes about floating point arithmetic) but you'll make everyone read
the code several times just to be sure.

Maybe around S(20).

That's not an accuracy. Do you need to know what S(20) is exactly
(39750532290178371794884752 / 55835135) or will an approximation do? If
so, what accuracy do you need? Would 7.12e17 be good enough or do you
need 7.119268591394715e17 or better?
 
C

Chad

If your 'unsigned long' is, say, 32-bits, then the largest factorial it can
store is 12!

13! would require about 33 bits I think, and if you try and do it using 32
bits, it will give funny results.

Summing the factorials has similar problems; integers just have a too narrow
range.

Floats have a much wider range: typically a float might go up to around
1e38, and a double up to about 1e300. But a double only gives some 16
digits of precision.

Simplify your program to *only* show factorials from 1 to 100 (forget
command line input, just use a loop), and you will
see where it goes wrong. Then put in doubles, and see the difference.


How many numbers would 1e300 produce? What can I say, my math isn't up
to par.
 
C

Chad

An integer type is preferable.  I don't think there is anything actually
wrong with the above as far as C is concerned (because of the guarantees
C makes about floating point arithmetic) but you'll make everyone read
the code several times just to be sure.




That's not an accuracy.  Do you need to know what S(20) is exactly
(39750532290178371794884752 / 55835135) or will an approximation do?  If
so, what accuracy do you need?  Would 7.12e17 be good enough or do you
need 7.119268591394715e17 or better?

--


Okay, I misread that. 7.12e17 is good enough.
 
B

Ben Bacarisse

Chad said:
On Apr 20, 4:05 pm, Ben Bacarisse <[email protected]> wrote:

Okay, I misread that. 7.12e17 is good enough.

Then just use double to calculate the factorial and factorial sum and
you should be fine. Even using float for the harmonic sum will give you
more than three significant figures but I'd use double for that as
well. There is no need to be careful about the order in which you sum
the terms.

However, if all you need it S(1) to S(20) it would make sense to
pre-calculate these values and put them in a table.
 
D

Dann Corbit

The following program is my attempt to calculate the formula at the
following url...

http://i23.photobucket.com/albums/b363/CompressorX/1.jpg

#include <stdio.h>
#include <stdlib.h>

unsigned long fact(unsigned long);
unsigned long sum_fact(unsigned long);
float harmonic(unsigned long);
{snip}

Using float gives very little precision and dynamic range.
Probably, you want a toolkit if you are determined to write it yourself,
but a math program would be better. For instance, you could download
pari/gp and calculate things like this with ease.
 
N

Nick Keighley

How many numbers would 1e300 produce? What can I say, my math isn't up
to par.

well I make s(100) to be about 1e157 so it will fit comfortably into
1e300.

S(100) is (assuming I didn't screw up)
26289880679674628167607993019108020392724393035736896620237691009703387760747890
05520477295593941461248393109377004641631967533326482146878306204895788492767929
33475322073718652663784818718016211136/14466636279520351160221518043104131447711

this was written in scheme (Gambit scheme). Note the two posts that
have given exact answers have switched to a different language. Out-of-
the-box C isn't good for this sort of problem. As someone else
mentioned you need a multiple-precision arithmetic package. Google for
it.
 
B

bartc

How many numbers would 1e300 produce? What can I say, my math isn't up
to par.

My maths isn't quite up to it either.

But, modifying your code to use double as necessary, and just making it loop
over a set of numbers, on my machine (double=64-bit) it clearly starts going
after N=170.

#include <stdio.h>
#include <stdlib.h>

double fact(unsigned long);
double sum_fact(unsigned long);
double harmonic(unsigned long);

double fact(unsigned long n)
{
double acc = 1.0;
unsigned long i;

for (i = 1; i <= n; i++)
acc *= i;

return acc;
}

double sum_fact(unsigned long n)
{
double acc = 0;
unsigned long i;

for (i = 1; i <= n; i++) {
acc += fact(i);
}

return acc;
}

double harmonic(unsigned long n)
{
double acc = 0;
int i;

for (i = 1; i <= n; i++) {
acc += 1.0/i;
}

return acc;
}

/* int main(int argc, char **argv) */
int main(void)
{
unsigned long s;

/*
if (argc != 2) {
fprintf(stderr, "Invalid number of args\n");
return 1;
}
s = atoi(argv[1]);
*/

for (s=1; s<=200; ++s)
/* for (s=100; s<=100; ++s) */
printf("S(%d) = %g\n", s, sum_fact(s)/(harmonic(s)));

return 0;
}
 
B

bartc

jacob said:
Chad a écrit :

The lcc-win compiler featrues 105 digits precision floating point. The
dynamic range is +- 10 ^10000

I've tried it, and it sort of works. (Although some of the conversions
between qfloat (this extended type) and ordinary ints for example are
fiddly:
ltoq(&i,&qi);
acc += 1.0q/qi;
)

Using qfloats (and any necessary alterations), the function s() seems to
converge to about 1.5e9863 (luckily just within the limits of qfloat).

However in the neighbourhood of s(5000), the calculation
sum_fact()/harmonic() (see OP), was taking several seconds each, so could
only do a few spot checks.

(Qfloats seem to be about 20 times slower than doubles.)
 
E

Eric Sosman

The following program is my attempt to calculate the formula at the
following url...
[... code that needs >525-bit integers ...]

I haven't studied the code to find why it "breaks," because
the fundamental approach simply isn't going to work with the
native types on today's machines. Find a different formula, or
find a multiple-precision arithmetic package.

What other formula do have in mind?

None. What are you trying to calculate? Can you calculate
its logarithm, say, instead of the quantity itself? (That's how
I got to the ">525-bit" estimate, for example: by calculating
lgamma(101) instead of attempting 100 factorial.)
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

Forum statistics

Threads
473,743
Messages
2,569,477
Members
44,898
Latest member
BlairH7607

Latest Threads

Top