# Code involving factorials and a given value

3 ビュー (過去 30 日間)
Ryan Johnson 2022 年 5 月 1 日

Hi, I am trying to come up with a code that is capable of displaying the smallest factorial above an input value. An example would be, x=5000 and the output would be 5040. Can anyone help with this?
Thanks
##### 5 件のコメント表示非表示 4 件の古いコメント
John D'Errico 2022 年 5 月 1 日
THINK ABOUT THE HINTS I GAVE YOU. For example, is there a good approximation for the factorial function? It would probably be a sterling way to solve this problem, IF you think. Maybe I spelled sterling wrong. (HINT.) Could you use fzero on the log of that approximation?

サインインしてコメントする。

### 採用された回答

Voss 2022 年 5 月 1 日
@Ryan Johnson Using a while loop to solve this problem is a perfectly good approach. The main thing you need to change is the while condition. You have while i<=x, but the variable i is the "loop iteration counter", if you will (i.e., it's incremented by one on each iteration), so that condition is saying to iterate a certain fixed number of times, but if you knew how many iterations it would take beforehand, you'd use a for loop rather than a while loop. Instead, the condition should be based on the running product of the integers, which is n, since the objective is to iterate until that product reaches a certain fixed value. So just change that condition to while n<=x and you're good to go:
% prompt= 'Enter value of x: ';
x = 5000;
n=1; i=1;
% while i<=x
while n<=x
n=n*i;
i=i+1;
fprintf('n=%d\n',n)
end
n=1 n=2 n=6 n=24 n=120 n=720 n=5040
(I also fixed up the fprintf statement. And you could use the input function to take the value of x from the user.)
##### 2 件のコメント表示非表示 1 件の古いコメント
Voss 2022 年 5 月 1 日
You're welcome!

サインインしてコメントする。

### その他の回答 (2 件)

John D'Errico 2022 年 5 月 1 日
Now that you have a looped solution, let me show how to solve the problem more directly. There are actually several solutions I can think of.
The Stirling approximation for a factorial is a simple one.
In there, we would find a good approximation for the factorial function.
factorial(n) ~~ sqrt(2*pi*n)*n^n/exp(n)
This is a pretty good approximation. Thus
n = 0:10;
stir = @(n) sqrt(2*pi*n).*n.^n.*exp(-n);
stir(n)./factorial(n)
ans = 1×11
0 0.9221 0.9595 0.9727 0.9794 0.9835 0.9862 0.9882 0.9896 0.9908 0.9917
In fact, you should see the Stirling approxiamtion will consistently under-predict the factorial function. You can probably learn that from a second order approximation to the factorial, by seeing that the next term in a series would be always positive. But stir is itself pretty good as an approximation.
So take the log of that. That is pretty simple, as long as we are careful.
log(factorial(n)) ~~ log(pi)/2 + log(2)/2 + log(n)/2 + n*log(n) - n
How can we use this?
logfact = @(n) log(pi)/2 + log(2)/2 + log(n)/2 + n.*log(n) - n;
Now use fzero. Suppose you want to compute the smallest n, such that factorial(n) is larger than 5000? We need only a call to fzero now.
facttarget = 5000;
ntarget = floor(fzero(@(n) logfact(n) - log(facttarget),[eps,100]))
ntarget = 7
factorial(ntarget)
ans = 5040
Sigh. Is this overkill? Perhaps. A simpler solution can be found from a table lookup. If we recognize that any factorial above factorial(170) will overflow double precision arithmetic, then we could just compute the logs of all the factorials, as:
N = (1:170)';
logfactlist = cumsum(log(N));
Now, suppose we wanted to find the smallest value of n such that factorial(n) is at least some target? This should work:
nsolve = @(target) ceil(interp1(logfactlist,N,log(target),'linear'));
nsolve([5000 5040 5041])
ans = 1×3
7 7 8
In fact, this function will work as high up as factorial(170). It would fail beyond that point, since I only generated a lookup table that went that high.
Can you solve the problem more simply, using a while loop? Well, you already got that answer.
##### 0 件のコメント表示非表示 -1 件の古いコメント

サインインしてコメントする。

Walter Roberson 2022 年 5 月 1 日
Just generate factorial(1:18) and discretize() against it. Beyond 18, you would be dealing with integers larger than double precision can typically represent exactly.

サインインしてコメントする。

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by