MATLAB Answers

Translated by

このページのコンテンツは英語から自動翻訳されています。自動翻訳をオフにする場合は「<a class="turn_off_mt" href="#">ここ</a>」をクリックしてください。

0

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

Jai Mechanical さんによって質問されました 2019 年 2 月 22 日
最新アクティビティ John D'Errico
さんによって コメントされました 2019 年 2 月 22 日
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);

  3 件のコメント

John D'Errico
2019 年 2 月 22 日
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.
Thanks for the reply!
BTW my program worked by replacing
factorial=factorial*i with factorial=sym(factorial*i)
Thanks anyway!
John D'Errico
2019 年 2 月 22 日
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.

サインイン to comment.

3 件の回答

回答者: John D'Errico
2019 年 2 月 22 日
編集済み: John D'Errico
2019 年 2 月 22 日
 採用された回答

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.

  5 件のコメント

Adam Danz
2019 年 2 月 22 日
But the wandering is what's interesting. Good read!
Rik
2019 年 2 月 22 日
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
2019 年 2 月 22 日
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.
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.

サインイン to comment.


回答者: Rik
2019 年 2 月 22 日
編集済み: Rik
2019 年 2 月 22 日

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

  0 件のコメント

サインイン to comment.


回答者: Adam Danz
2019 年 2 月 22 日
編集済み: Adam Danz
2019 年 2 月 22 日

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.

  3 件のコメント

"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
2019 年 2 月 22 日
Or:
f = factorial(n);
if isinf(f)
s=NaN;
else
s = sum(int2str(f)-'0');
end
Adam Danz
2019 年 2 月 22 日
+1, I forgot about that one!

サインイン to comment.



Translated by