Facing a problem in calculating the sum of digits of 100!.

This program is working fine ( calculate the sum of digits of a factorial of a number.) for smaller numbers but is giving wrong answer for 100!.

The answer should be 648 but mine is coming out to be 659.

clc;
clear;
n=input('enter any number');
i=1;
factorial=1;
sum=0;
while i<=n
    factorial=factorial*i;
   i=i+1;
end
a=factorial;
 while a~=0  
       remainder=rem(a,10);
        sum = sum+remainder;
       a=fix(a/10);
 end
   fprintf('%d',sum);

[John D'Errico]

You did make a credible attempt at this. Compare what you did to the code in my answer. You will see that mine is not really that far off. So you did have the right general approach. Spend some time digging into what I wrote to see why it works. Put the two codes next to each other.

[Jai Mechanical]

Thanks for the reply!

BTW my program worked by replacing

factorial=factorial*i with factorial=sym(factorial*i)

Thanks anyway!

[John D'Errico]

Well, yes. But if you are going to do that, using syms, then the one line solution I showed seems as good:

sum(char(factorial(sym(100))) - '0')

and is far more efficient. But as long as you have something that works, it is best if you wrote the code.

Answer 1

[John D'Errico]

It looks like this question has gotten enough interest that it deserves a complete answer on how to do it using MATLAB. And since you have done something reasonable...

You simply cannot compute factorial of a number as a double, if n is greater than 18.

factorial(1:20) < flintmax
ans =
  1×20 logical array
   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   0   0

Actually, it looks like you might get to 20, if you use uint64. But uint64 overflows beyond that point. So using a double will just fail for the sum of digits of factorial(100).

factorial(100)
ans =
  9.3326e+157

That is WAY past flintmax. Personally, I'd do it elegantly as:

sum(char(factorial(sym(100))) - '0')
ans =
   648

However nothing stops you from doing the multiplication the hard way, using carries.

factdigs = zeros(1,200);
factdigs(end) = 1;
N = 100;
for n = 2:N
    % multiply. note that you multiply ALL digits of
    % the number by n.
    factdigs = factdigs * n;
    % carry op
    carry = 1;
    % a while loop, since there may be multiple carries.
    % for example, a multiple carry would occur when you
    % multiply the number 6 times 100.
    while any(carry)
        % strip off the units digit
        R = rem(factdigs,10);
        % this divide will always be safe, resulting
        % in an exact integer.
        carry = (factdigs - R)/10;
        % The carry is just an add, with a shift.
        % if I knew the highest order digit will
        % always be zero, I could have used circshift.
        factdigs = R + [carry(2:end),0];
    end
end
sum(factdigs)
ans =
   648

This will work as long as I make the vector of digits long enough that factorial(100) does not overflow a 200 digit number. And since we know that factorial(100) is on the order of 9e157, 200 digits is fine.

reshape(factdigs,10,20)'
ans =
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     9     3     3     2     6     2     1     5
     4     4     3     9     4     4     1     5     2     6
     8     1     6     9     9     2     3     8     8     5
     6     2     6     6     7     0     0     4     9     0
     7     1     5     9     6     8     2     6     4     3
     8     1     6     2     1     4     6     8     5     9
     2     9     6     3     8     9     5     2     1     7
     5     9     9     9     9     3     2     2     9     9
     1     5     6     0     8     9     4     1     4     6
     3     9     7     6     1     5     6     5     1     8
     2     8     6     2     5     3     6     9     7     9
     2     0     8     2     7     2     2     3     7     5
     8     2     5     1     1     8     5     2     1     0
     9     1     6     8     6     4     0     0     0     0
     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0

Of course, the code as I wote it would fail around n = 122.

factorial(122)
ans =
   9.875e+202

That would be fixable merely by preallocating factdigs to be a longer vector. Or, with a little more care in the code, I could have written it to dynamically expand the vector of digits, so factorials of any size would then theoretically work.

[Rik]

Wouldn't this work? At some point Matlab will respond inf of course, which would be quite a while before this stop working, but that would be a separate fix.

factdigs = zeros(1,1+floor(log10(factorial(N))));
%2+ might be needed instead?

[John D'Errico]

Good try, but no. You are lost as soon as you compute factorial(n). n is a double. so factorial(n) is still a double. It does not carry 157 significant digits. So only the first 16 significant digits of that have any meaning at all.

Computing the log of that number does not help, because the log was computed too late. But you cannot even get tricky and throw logs at it sooner. (I know someone will want to ask this.)

So next, I might think to try computing log(1:n). Would that help? (Asked rhetorically, becaue I know the answer. No it does not help.) That is, in theory, one could compute factorial(n) as

exp(sum(log(1:n)))
ans =
  9.3326e+157

It even looks like it is correct. The problem is still that each of those numbers in the log is a double. And I need EVERY digit of 157 of them to be exact, because I want to compute the sum of the decimal digits.

They would need to carry each roughly 157 significant digits for the result to work. So, for example, let me use my own HPF. First, in 200 decimal digits of precision.

exp(sum(log(hpf(1:n,200))))
ans =
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

But, if I used less than roughly 158 digits of precision (realistically, I would need a few more digits to be safe) then HPF would also fail here.

sum(mantissa(round(exp(sum(log(hpf(1:n,160)))))))
ans =
   648

What you need to recognize is that doubles go only so far when you want to know the digits of an integer. And they fall apart completely as soon as you go past 2^53 (flintmax). A double only carries the information content to give you roughly 16 decimal digits, no matter how you try to twist that information around.

So unless you use high precision in SOME incarnation, you cannot recover all 157 digits in that sum. At best, you might as well use an approximate result, thus, since the average of the decimal digits in a huge number will approach 4.5 as a limit, I might hope to see this be close:

157*4.5
ans =
        706.5

And realistically, 706.5 is not that far off, though a little high. One reason why factorial(100) has a lower sum of digits than my expected value there is because factorial(100) has a lot of zeros at the end. You will have ONE trailing zero in factorial(n) for every power of 5 in that product. There are 24 trailing zeros, (don't forget 25, 50, 75, 100), which each yield TWO factors of 5 in the product.

Finally, is 4.5 a good approximation? It should be. For example, are the decimal digits of pi uniformly distributed?

sum(mantissa(hpf('pi',500000)))/500000
ans =
       4.5002

While I will not argue the digits of pi are truly random, 4.5 for a mean there suggests they are pretty well uniformly distributed. Oh well, I have obviously wandered off way too deeply into the woods here. ;-)

[Adam Danz]

But the wandering is what's interesting. Good read!

[Rik]

I just meant it as a step to get a rough approximate for the number of digits (RE "Or, with a little more care in the code, I could have written it to dynamically expand the vector of digits, so factorials of any size would then theoretically work."). Since in that step we only care about the order of magnitude, it doesn't matter at all if the answer is numerically correct, as long as the power is a good approximate.

But thank you for thoughts, an interesting read.

[John D'Errico]

Ok, yes. This tells me how many decimal digits you have in a big factorial.

sum(log10(1:100000))
ans =
          456573.450899971

So there are 456573 decimal digits in factorial(100000).

Another nice solution would use Stirlings approximation for the factorial.

https://en.wikipedia.org/wiki/Stirling%27s_approximation

Which tells us that

factorial(n) ~ sqrt(2*pi*n)*n^n/exp(n)

This is quite accurate for even small n.

n = 5;
sqrt(2*pi*n)*n^n/exp(n)
ans =
           118.01916795759

As I recall, there are additional terms in the series.

sqrt(2*pi*n)*n^n/exp(n)*(1+ 1/(12*n) + 1/(288*n^2))
ans =
          120.002545641322

And it looks like

sum(floor(100000* 5 .^(-1:-1:-7)))
ans =
       24998

trailing zeros in factorial(100000). So I might expect to see a digit sum for factorial(100000) of

456573*4.5 - 24998*4.5
ans =
                 1942087.5

Which is reasonably close to truth, as measured by my vpij.

sum(digits(factorial(vpij(100000))))
ans =
     1938780

That computation for n=100000 using the loop approach I gave in my answer is a fairly long one, because the sum of decimal digits in factorial(n) should be an O(n^2) operation to compute.

It is always fun playing with numbers. Time to start crunching again though. My search for the second prime in an interesting modified Woodal prime class is now crunching through numbers with over 50000 decimal digits. Lets see if I can melt down my CPU.

Answer 2

[Rik]

This is due to how computers store digits. You will need to rethink your method for big numbers:

clc
%don't use clear, use functions to keep your workspaces clean
n=input('enter any number\n');

if flintmax<factorial(n)
    err_txt={'This value cannot be uniquely represented as a floating point integer.',...
        'This means that the last few digits contain rounding errors.'};
    error('%s\n',err_txt{:})
else
    fact=1;%don't shadow the builtin factorial function
    s=0;%don't shadow the builtin sum funtion
    for ii=1:n%don't use i (or j) as a variable name
        fact=fact*ii;
    end
    a=fact;
    while a~=0
        remainder=rem(a,10);
        s = s+remainder;
        a=fix(a/10);
    end
    fprintf('The sum of the digits of %d! is %d\n',n,s);
end

Answer 3

[Adam Danz]

This solution is fast and simple but only works for numbers <=21 (see notes below) so it's not the best approach.

It computes the factorial, converts it into a character vector and then into a numerical vector so that each number in the integer is separate. Then it just adds up the elements stored in the vector.

% Ask user for number
n=input('enter any number: ');

% Calculate factorial
f = factorial(n); 
% Convert to string
fs = int2str(f); 
% separate string into single numbers (still strings)
fsc = cellstr(fs')'; 
% Convert to numbers, vector
fn = str2double(fsc); 
% Add the numbers in the vector
s = sum(fn); 

% return sum
fprintf('Sum is %d\n', s)

For 100! the sum of each individual number comes to 683 which doesn't match your expectation due to loss of precision.

Note that if matlab returns 'Inf' for the factorial (eg, 1000!), this code will return "NaN'.

Also note that according to Matlab's documentation, for double precision inputs the factorial result is exact when n <= 21. Larger values are accurate up to the first 15 digits.

[Stephen23]

"For 100! the sum of each individual number comes to 683 which doesn't match your expectation "

100 factorial is larger than flintmax, so it cannot be stored as a double without loss of precision. Whatever you do with it after that is irrelevant.

[Rik]

Or:

f = factorial(n); 
if isinf(f)
    s=NaN;
else
    s = sum(int2str(f)-'0');
end

[Adam Danz]

+1, I forgot about that one!