Calculate the sum of digits in 100 factorial

Question

Edit - Changed title to match the actual problem statement.

I'm programming a function that calculates the sum of digits in 100! but I seem to be having two big issues.

1. The actual result of 100! is only accurate to the first few numbers (actual result is 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000)

2. My method for adding up the digits of the resulting number doesn't output the correct result.

This is my current code:

void factorialSum()
{
    double fact100 = factorial(100);
    double suma = 0;

    printf("100! is equal to: %.0f", fact100);

    while (fact100 > 0)
    {
        double temporal = fmod(fact100, 10);
        suma = suma + temporal;
        fact100 = fact100/10;
    }
    printf("\nThe sum of all digits in 100! is: %.0f", suma);
}

And the function factorial() is defined as:

double factorial (double n)
{
    double mult = 1;
    double i = n;

    while (i>=1)
    {
        mult *= i;
        i = i - 1;
    }
    return mult;
}

The program outputs 93326215443944102188325606108575267240944254854960571509166910400407995064242937148632694030450512898042989296944474898258737204311236641477561877016501813248 as a result for 100! and says the summation of its digits is equal to 666.

Any help is appreciated, thank you.

Answer 1

In C, a double typically has 53 bits of precision, which corresponds to 16 or 17 digits of precision. So once you go beyond 22!, a double can no longer represent the exact result, as demonstrated by the following code. Note that at 23!, the trailing zeros disappear, since the double no longer represents the exact value.

#include <stdio.h>
#include <stdint.h>

int main( void )
{
    double y;

    y = 1;
    for ( int i = 2; i < 30; i++ )
    {
        y *= i;
        printf( "%2d %32.0lf\n", i, y );
    }
}

Here's the output from the program

 2                                2
 3                                6
 4                               24
 5                              120
 6                              720
 7                             5040
 8                            40320
 9                           362880
10                          3628800
11                         39916800
12                        479001600
13                       6227020800
14                      87178291200
15                    1307674368000
16                   20922789888000
17                  355687428096000
18                 6402373705728000
19               121645100408832000
20              2432902008176640000
21             51090942171709440000
22           1124000727777607680000
23          25852016738884978212864
24         620448401733239409999872
25       15511210043330986055303168
26      403291461126605650322784256
27    10888869450418351940239884288
28   304888344611713836734530715648
29  8841761993739700772720181510144

If you want to compute the exact value of 100! you need to use arrays of digits (aka bignums), to do the calculations. You can either find a bignum library to use, or implement bignum multiplication yourself. The wikipedia article on bignums provides the pseudocode for calculating factorials.

Answer 2

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

typedef unsigned char byte;

int main(){
    size_t size =16;//initial size
    size_t temp, len;
    byte *n = malloc(size);
    byte i, carry = 0;
    *n = 1;
    len = 1;

    for(i=2;i<=100; ++i){
        size_t k;
        for(k = 0; k < len; ++k){
            if(k == size){
                n = realloc(n, size*=2);//expand, check omitted
            }
            temp = n[k] * i + carry;//Calculation is performed on the promoted to int
            carry = (temp >= 10) ? temp / 10 : 0;//or simply temp/10
            n[k] = temp % 10;
        }
        while(carry){
            if(k == size){
                n = realloc(n, size*=2);
            }
            temp = carry;
            carry = (temp >= 10) ? temp / 10 : 0;
            n[k++] = temp % 10;
            ++len;
        }
    }
    temp = 0;
    while(len){
        printf("%u", n[--len]);
        temp += n[len];
    }
    puts("");
    printf("sum=%zu\n", temp);
    free(n);

    return 0;
}
#if 0
anticipate in advance the required size
100! = 1*2*...99*100
size > (log(1)+log(2)+...log(99)+log(100))/log(10)
(∫1->100 log(x)dx)/log(10)
f(x)=log(x) => F(x)=x(log(x)-1)+C
(∫1->101 log(x)dx)/log(10) = (101*log(101)-101+1)/log(10) ≒ 159.00..
160 digits are sufficient.
http://en.wikipedia.org/wiki/Stirling%27s_approximation
log10(n!)≒log10(√(2nπ)(n/e)^n)=log10(√(2nπ))+nlog10(n/e)=157.9696..
size=158
#endif

Answer 3

As others have mentioned, you can either write your own bignum library for this (using arrays of bytes), or you can use something like OpenSSL's BIGNUM implementation.

Here is the OpenSSL version of the fact function (compiled with gcc main.c -lcrypto on nearly any distro). You just need to include <openssl/bn.h>.

BIGNUM* fact (unsigned long n, BN_CTX* ctx) {

    BIGNUM *mult, *i, *one;

    mult = BN_new();
    i = BN_new();
    one = BN_new();

    BN_one(one);

    BN_set_word(mult, 1L);
    BN_set_word(i, n);


    while (BN_cmp(i, one) >= 0) {
        BIGNUM *a, *b;

        a = BN_new();
        b = BN_new();

        BN_mul(a, i, mult, ctx);

        BN_sub(b, i, one);

        BN_free(mult);
        BN_free(i);
        mult = a;
        i = b;
    }

    BN_free(one);
    BN_free(i);

    return mult;
}

With that, you can then use BN_bn2bin to generate a string representation of this number, which you can then use to calculate the sum of the digits.

Answer 4

Here is something that might help. Suppose you have a number like 2356. How can you add it's digits. Well, you can extract the least significant digit, add it to result and throw it away (right shift). You extract the least digit by taking number mod 10. You shift by dividing by 10. So 2356 mod 10 = 6 and 2356 / 10 = 235.

Taking the mod is easy, you just go over the numbers in your product, mod them and multiply the mods together. For example 12*6 mod 10 = (12 mod 10) * (6 mod 10) = 2 * 6 = 12 and then take one last mod at the end: 12 mod 10 = 2. If you were to do it by multiplying 12*6 you will get 72 which is equal to 2 mod 10.

The hard part is dividing by 10. Now if the product contains a multiple of 10 (like 100 for example) you divide that number by 10 and you are done. If no such number exists you can still divide by 10 but that will screw up the results i guess. So if you can find how to do this you can solve the problem without implementing a BIGNUM. If not, then it's not that hard either, you only need to implement a function to add two BIGNUMs together (multiplication is just repeated addition).