An Exercicse from A=B

edit

Function $r(n,k)$ and thus $g(n,k)$ can be computed from $f(n,k)$ as follows:

var('k n')
f(n,k) = factorial(n)^4 / factorial(k)^2 / factorial(n-k)^2 / factorial(2*n)
print( f.WZ_certificate(n, k) )

which prints 1/2*(2*k - 3*n - 3)*k^2/((k - n - 1)^2*(2*n + 1)) matching given $r(n,k)$.

A tad less esoteric and much more lazier than @Max Alekseyev's excellent answer ; on Sage 10.4.rc3 :

sage: var("n, k")
(n, k)
sage: f(n, k)=factorial(n)^4/(factorial(k)^2*factorial(n-k)^2*factorial(2*n))
sage: r(n, k)=-k^2*(3*n+3-2*k)/(2*(n+1-k)^2*(2*n+1))
sage: g(n,k)=r(n,k)*f(n,k)
sage: %time (f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k)).simplify_factorial().factor()
CPU times: user 5.87 s, sys: 12.5 ms, total: 5.88 s
Wall time: 4.76 s
0

And, BTW, one can be even lazier :

sage: %time (f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k)).is_zero()
CPU times: user 8.09 s, sys: 56.1 ms, total: 8.15 s
Wall time: 6.75 s
True

To understand why, try :

sage: view((f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k)).numerator())
sage: view((f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k)).numerator().simplify_factorial())
sage: view((f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k)).numerator().simplify_factorial().expand())

EDIT : Of course, this method is valid if and only if the denominator(s) aren't zero, which is possible :

sage: (f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k)).denominator().simplify_factorial().factor().solve(k)
[k == n + 1, k == n, k == -1, factorial(-k + n - 1) == 0, factorial(k) == 0]

EDIT 2 : Other ways :

sage: (f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k)).simplify_full() # simplify() isn't enough...
0
sage: (f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k))._sympy_().simplify()
0
sage: (f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k))._giac_().simplify()
0
sage: # Uses Wolfram Engine
sage: (f(n+1,k)-f(n,k)-g(n,k+1)+g(n,k))._mathematica_().FullSimplify() # Simplify() isn't enough
0

HTH,