Can matlab solve this simple factorial equation?

I am interested in knowing whether MATLAB can solve the equation x! == y! to get x == y. I have tried

syms x1 x2
eqns = [factorial(x2) == factorial(x1)];
S = solve(eqns, [x1 x2]);
S.x1
S.x2

But the response is

 “Warning: Unable to solve symbolically.
Returning a numeric approximation instead”.
In solve (line 304) 
ans =
0
ans =
0

I used the gamma function instead of factorial but get the same response. Any suggestions?

[Steven Lord]

Your proposed solution is incomplete. (x = 0, y = 1) and (x = 1, y = 0) are also solutions.

Answer 1

[John D'Errico]

Um, the problem is, this is arguably more complicated than you think.

First of all, lets extend the problem to the real line. We know that

factorial(x) = gamma(x+1)

The gamma function is defined for all values of x, whereas factorial is classically defined for only non-negative integers. You need to do this extension if you will use solve. Solve is not really designed to solve problems with integer variables (thus essentially Diophantine equations.)

If you want to solve that relationship, is it true that x1 == x2? ALWAYS? Yes, it is true that is one solution. But in fact, there are infinitely many solutions to the relationship

gamma(x1+1) == gamma(x2+1)

such that x1 is NOT equal to x2. You can easily convince your self of that fact simply by a plot.

ezplot(@(x) gamma(x+1))
grid on

For example, I'll pick a value of x1.

x1 = 0.2;
gamma(x1+1)
ans =
    0.91817

x2 = fzero(@(x) gamma(x + 1) - gamma(x1 +1),2)
x2 =
    0.74604

gamma(x2+1)
ans =
    0.91817

As you can see, I chose one value for x1, then I found a second value x2, where x1 was not equal to x1, but the two produce the same value.

So there are infinitely many solutions with x1~=x2, if you use the extension of factorial to the real line. vpasolve can give some of them to us, but not solve.

x1 = 0.2;
syms x

x2 = vpasolve(gamma(x1 + 1) == gamma(x+1),x,3)
x2 =
0.74604419294554377739652086343099

x2 = vpasolve(gamma(x1 + 1) == gamma(x+1),x,0)
x2 =
0.19999999999999984718524773398979

x2 = vpasolve(gamma(x1 + 1) == gamma(x+1),x,-2.3)
x2 =
-4.1549616274283604766122794976542

Lets look at it another way. There is no analytical inverse to the gamma function. So a symbolic engine won't understand that it can just take an "inverse" of the gamma function here, getting x1==x2. Symbolic tools don't look at your problem from a high level.

[Rob Burns]

Thanks John. Your comment about solve not being designed to solve integer equations is interesting. Is there another way within matlab to do this? I haven’t seen one but I am new to matlab.

The issue with the gamma function is a bit of a red herring - my fault for not realising the implications that using gamma has on the domain and range of the function. Your point about that is spot on so ignore that I mentioned gamma. I am only interested in integer solutions so let’s stick to the factorial restricted to positive integers. Since factorial is injective on this domain I expected it to be able to simplify x! ==y! but it seems that it can’t unless I am missing something. Just to be certain matlab knows we are dealing with positive integers I tried adding in some assumptions as follows:

#
syms x1 x2
assume(x1,'integer')
assumeAlso(x2,'integer')
assumeAlso(x1>0)
assumeAlso(x2>0)
assumptions
eqns = [factorial(x2) == factorial(x1)];
S = solve(eqns, [x1 x2], 'ReturnConditions', true);
S.x1
S.x2
S.conditions
S.parameters

It didn’t help as the response was

ans = [ in(x1, 'integer'), in(x2, 'integer'), 0 < x1, 0 < x2]
ans = z
ans = z1
ans = 0 < z & 0 < z1 & factorial(z1) - factorial(z) == 0
ans = [ z, z1]

which unfortunately isn’t very helpful. It is interesting that matlab can be helpful in this situation for other functions. E.g. if we look at the square function on the positive integers we have

#
syms x1 x2
assume(x1,'integer')
assumeAlso(x2,'integer')
assumeAlso(x1>0)
assumeAlso(x2>0)
assumptions
eqns = [x2^2 == x1^2];
S = solve(eqns, [x1 x2], 'ReturnConditions', true);
S.x1
S.x2
S.conditions
S.parameters
###
ans = [ in(x1, 'integer'), in(x2, 'integer'), 0 < x1, 0 < x2]
ans = k
ans = k
ans = 1 <= k & in(k, 'integer')
ans = k.

So regarding my original question, I really just want to know one way or another whether there is a way that the equation x! == y! can be solved to give x == y for positive integers. You imply that the answer is no. That’s fine. Thanks.

[Walter Roberson]

S  = solve(eqns,[x1,x2]);

S.x1, S.x2 now both show up as 1 .

This is the not-very-helpful way that MATLAB outputs range() expressions generated by the symbolic engine: it finds a representative value within the range and outputs the representative value. A few months ago I deduced how the representative values were chosen and wrote it up in an Answer. The TL;DR is that 1 is chosen to stand in for (0,inf)

[Rob Burns]

Thanks Walter. Do you have a link to your post that you mentioned. It does seem a rather obscure way of providing a solution. From what you have indicated the solution is in the form of a range rather than a parameterisation. How is this range then related back to the solution set of the variables x1 and x2?

[Walter Roberson]

https://www.mathworks.com/matlabcentral/answers/356969-solving-absolute-value-inequality#answer_281887

[Walter Roberson]

One thing to keep in mind with the version that returns a pair of 1's is that you are asking to solve a single equation for two variables: normally no result would be expected from that.

[Rob Burns]

Ok Walter. I had a look at your analysis. Thanks. I think the correct answer to my original question is that matlab won’t provide a parametrised solution set to the equation x! == y! in the same way it does for the equation x^2 == y^2 for example. That means it isn’t suitable for the purpose I had in mind.

Answer 2

[Torsten]

syms x1 x2
eqn = factorial(x2) == factorial(x1);
S = solve(eqn, x1)

Best wishes

Torsten.

[Rob Burns]

With this coding I get

“Warning: Unable to find explicit solution. For options, see help”.

The help link says

“web(fullfile(docroot, 'symbolic/troubleshoot-equation-solutions-from-solve-function.html”.