Searching factorials to 100! finds no integral solutions, except for 1^20 = 1!., but we can search for the generalized factorial over the reals.

First trying for integral solution, with standard factorials:

for r=20:20:100

  disp([r factorial(r) r^20])

end

narrowed this down to between 20 and 40.

Then

for r=sym(20):40

  f=factorial(r);

  p=r^20;

  disp(r)

  disp(f)

  disp(p)

  disp(' ')

end

showed that the passing point was between 27 and 28.  Neither integral value had an equality.

The generalized factorial is the Gamma function of n+1, so:

low=27; high=28; dx=.01;

found=false; prev=0;

while ~found

for i=low-dx:dx:high

  diff=genFact(i)-i^20;

  if diff==0

    found=true;

    break

  end

  if i>low-dx/2

    if prev*diff<0

      low=i-dx;

      high=i;

      dx=dx/100;

      break

    end

  end

  prev=diff;

end

end

function v=genFact(n)

  v=gamma(n+1);

end

>> rFactorialEqScoreExponent

Operation terminated by user during rFactorialEqScoreExponent 

>> i

i =

          27.5247695102146           this is r

>> gamma(i+1)

ans =

      6.22979589706583e+28          and the equal values

>> i^20

ans =

      6.22979589706581e+28           of the functions

 

 

I haven't explored negative r, but for example, (-2)^20 = 1048576 and (-3)^20 = 3486784401. In between these x values, the Gamma function of x+1 has two arms extending without limit upward from about two and a third. There must be two crossing points in that range. I did just one of these, r=-2.00000095366484 where each function evaluates to about 1048586.