0

I'm trying to represent the following mathematical expression in C:

P(n) = (n!)(6^n)

The program should compute the answer to expression when n = 156. I have attempted to create the program in C and it fails to produce an answer. The answer is approximately 10^397. The program utilises 2 logarithmic identities. It also utilises Stirling's approximation to calculate the large factorial.

How can I make it produce the correct answer and do you have any suggestions as to how I could improve the code? (I'm fairly new to programming):

#include <math.h>
typedef unsigned int uint;

int main()
{
uint n=156;                               // Declare variables
double F,pi=3.14159265359,L,e=exp(1),P;
F = sqrt(2*pi*n) * pow((n/e),n);          // Stirling's Approximation Formula
L = log(F) + n*log(6);                    // Transform P(n) using logarithms - log(xy) = log(x) + log(y) and log(y^n) = n*log(y)
P = pow(e,L);                             // Transform the resultant logarithm back to a normal numbers
}

Thank you! :)

Sneftel
  • 40,271
  • 12
  • 71
  • 104
thelastpanda
  • 1
  • 2
  • Hey I saw this question on http://codereview.stackexchange.com/ (in Python instead of C, which is why it's not likely to work out for you as easily as it worked there). – barak manos Dec 02 '14 at 18:50
  • 10^397 is not representable as a double. See DBL_MAX. – Sneftel Dec 02 '14 at 18:51
  • Yep, that was posted by me. I'm now trying to represent the same program but in C. I think there's an issue with the precision I am using though.I posted it on Overflow instead of Code review as the program needs fixing. – thelastpanda Dec 02 '14 at 18:53
  • possible duplicate of [I would like to add 2 arbitrarily sized integers in C++. How can I go about doing this?](http://stackoverflow.com/questions/2926219/i-would-like-to-add-2-arbitrarily-sized-integers-in-c-how-can-i-go-about-doin) – Cloud Dec 02 '14 at 19:12
  • Changing your code to use `long double` I get 1.83969e+397 with gcc... – Dmitri Dec 02 '14 at 19:32
  • `long double` *may or may not*` be able to represent the result. On some systems, `long double` is substantially larger in range and precision than `double`. On others, it's the same (though it's still a distinct type). – Keith Thompson Dec 02 '14 at 20:37

4 Answers4

4

Neither integer nor floating point variables in most C implementations can support numbers of that magnitude. Typical 64-bit doubles go up to something like 10308 with substantial loss of precision at that magnitude.

You'll need what's called a 'bignum library' to compute this, which is not part of standard C.

Russell Borogove
  • 18,516
  • 4
  • 43
  • 50
1

One idea is to use the long double type. Its precision isn't guaranteed, so it may or may not be big enough for your needs, depending on what compiler you're using.

Replace double with long double. Add an 'l' (lower case L) suffix to all math functions (expl, logl, powl, sqrtl). Compile with C99 enabled, since the long double math functions are provided in C99. It worked for me using GCC 4.8.1.

#include <math.h>
#include <stdio.h>
typedef unsigned int uint;

int main()
{
    uint n=156;                               // Declare variables
    long double F,pi=3.14159265359,L,e=expl(1),P;
    F = sqrtl(2*pi*n) * powl((n/e),n);          // Stirling's Approximation Formula
    L = logl(F) + n*logl(6);                    // Transform P(n) using logarithms - log(xy) = log(x) + log(y) and log(y^n) = n*log(y)
    P = powl(e,L);                             // Transform the resultant logarithm back to a normal numbers
    printf("%Lg\n", P);
}

I get 1.83969e+397.

Fred Larson
  • 60,987
  • 18
  • 112
  • 174
  • I get approximately -4.68e+66. :( – thelastpanda Dec 02 '14 at 19:40
  • Ouch. What compiler are you using? You may have to tweak your compile options. – Fred Larson Dec 02 '14 at 19:41
  • Try adding an `#include ` and `printf("%d\n", LDBL_MAX_10_EXP);`. That will show you the maximum base-10 exponent `long double` can represent. You might be out of luck with that compiler. – Fred Larson Dec 02 '14 at 19:48
  • Interesting. That should be big enough. Did you miss any of the math functions? – Fred Larson Dec 02 '14 at 19:53
  • Posted what I did below: – thelastpanda Dec 02 '14 at 19:58
  • Please don't post that as an answer. Looks like you got them all. What specifier did you use for `printf`? – Fred Larson Dec 02 '14 at 20:00
  • I meant to print the `long double` result. You need an 'L' prefix, like `"%Lf"` or `"%Lg"`. – Fred Larson Dec 02 '14 at 20:03
  • I tried %Ld, %Lf and %Lg. None of them printed the correct answer. – thelastpanda Dec 02 '14 at 20:12
  • There seems to be a common misconception that `long double` was introduced by C99 ((I've even seen it in published books). In fact, the `long double` type was introduced by the ANSI C standard in 1989. The `` functions with names ending in `l` that operate on `long double` *were* introduced by C99. (I know you said "provided by", not "introduced by".) – Keith Thompson Dec 02 '14 at 20:20
  • @KeithThompson: You're correct. C99 apparently increased support for `long double`, but it existed in C89. http://en.wikipedia.org/wiki/Long_double#History Do you recommend any rewording of my answer to avoid perpetuating the myth? – Fred Larson Dec 02 '14 at 20:25
  • @FredLarson: You could just drop the phrase "provided in C99", and perhaps add something to say that the `*l()` library functions were added by C99. – Keith Thompson Dec 02 '14 at 20:36
0

Loosely speaking, in C a double is represented as a base number raised to a power. As already mentioned, the maximum is roughly 1E308, but as you get to larger and larger numbers (or smaller and smaller), you lose precision because the base number has a finite number of digits and cannot always be accurately represented in this way.

See http://en.wikipedia.org/wiki/Double-precision_floating-point_format for more information

SpikeMF
  • 40
  • 5
  • I don't think this answers the question. The problem is that the mathematical result is outside the range of type `double`, not that there's insufficient precision. – Keith Thompson Dec 02 '14 at 20:29
  • Well, yes. That answer was already given, but I think that if someone is going to be working with large numbers like this, it's important to understand, at least on some level, how it works. – SpikeMF Dec 02 '14 at 20:32
  • Then it would be more appropriate as a comment. – Keith Thompson Dec 02 '14 at 20:36
  • I would have, but it requires 50 reputation to comment on other people's answers. – SpikeMF Dec 03 '14 at 15:53
0
#include <math.h>
#include <float.h>
typedef unsigned int uint;

int main()
{
uint n=156;                               // Declare variables
long double F,pi=3.14159265359,L,e=expl(1),P;
F = sqrtl(2*pi*n) * powl((n/e),n);          // Stirling's Approximation Formula
L = logl(F) + n*logl(6);                    // Transform P(n) using logarithms - log(xy) = log(x) + log(y) and log(y^n) = n*log(y)
P = powl(e,L);                             // Transform the resultant logarithm back to a normal numbers
printf("%d\n", LDBL_MAX_10_EXP);
}
thelastpanda
  • 1
  • 2