Sympy solve fails to find nullspace of easy equation

I am trying to find the nullspace of various function-sets with pythons module sympy. It managed to find the solution to some sets like
{(x - 1)!, x * (x - 2)!, (x - 2)!}
My code is

from sympy import solve, factorial
from sympy.abc import a, b, c, x
eq = a * factorial(x - 1) + b * x * factorial(x - 2) + c * factorial(x - 2)
print(solve(eq, a, b, c, set=True))  #  output: ([a, b, c], {((-b*x + b - c)/x, b, c)})
eq = -b * x + b - c - a * x
print(solve(eq, a, b, c, set=True)) #  output: ([a, b], {(-c, c)})

However, it struggled with the set
{upper_gamma(x, -1), upper_gamma(x - 1, -1), x * upper_gamma(x - 1, -1), (-1) ^ x}
Which is fair enough, even wolfram-alpha cannot find the solution. However, I was surprised when it also failed on the simplified problem
{(x - 1) * y - z, y, x * y, z}
My code was again the same

eq = a * (x-1) * y - a * z + b * y + c * x * y + d * z
print(solve(eq, a, b, c, d, set=True))  # output: ([a, b, c, d], {((b*y + c*x*y + d*z)/(-x*y + y + z), b, c, d)})

I expected the solution ([a, b, c], {d, d, -d}).
Am I using solve wrongly or is the equation to hard for this solver?

You are depending on an implicit behaviour of solve that given a single equation not in a list that is linear in all unknowns it will assume that you are looking for an undetermined coefficients solution which is semantically something quite different from an ordinary solution.

The difference between these can be seen like:

In [5]: solve(a*x + b*y, a, b, set=True)
Out[5]: ([a, b], {(0, 0)})

In [6]: solve([a*x + b*y], a, b, set=True)
Out[6]: 
⎛        ⎧⎛-b⋅y    ⎞⎫⎞
⎜[a, b], ⎨⎜─────, b⎟⎬⎟
⎝        ⎩⎝  x     ⎠⎭⎠

Here the first case is returning values of a and b that can make the expression identically zero for all possible values of x and y. The second case is treating x and y as being fixed but unspecified parameters and returning the set of possible values for a and b as functions of the symbolic parameters x and y that could make the expression zero.

Conceptually these are two very different things and it would be better if solve did not mix them up in a single end-user function. The difference that you see is because your equation is being treated like the second case rather than the first. That being said I cannot reproduce the output you report with the most recent version of SymPy (1.11.1):

In [1]: a, b, c, d = symbols('a, b, c, d')

In [2]: eq = a * (x-1) * y - a * z + b * y + c * x * y + d * z
   ...: print(solve(eq, a, b, c, d, set=True))  # output: ([a, b, c, d], {((b*y + c*x*y + d*z)/(-x*y + y 
   ...: + z), b, c, d)})
([a, b, c, d], {(d, d, -d, d)})

If I pass eq in a list then I do see it:

In [7]: solve([eq], a, b, c, d, set=True)
Out[7]: 
⎛              ⎧⎛      b⋅y          c⋅x⋅y          d⋅z             ⎞⎫⎞
⎜[a, b, c, d], ⎨⎜- ─────────── - ─────────── - ───────────, b, c, d⎟⎬⎟
⎝              ⎩⎝  x⋅y - y - z   x⋅y - y - z   x⋅y - y - z         ⎠⎭⎠

In any case it is best not to depend on using solve this way. I suggest to use solve only for its primary purpose which is case 2 above and for that it should always be called with a list (even for a single equation).

There is instead a function solve_undetermined_coeffs which is the explicit function for doing what you want to do here:

In [9]: solve_undetermined_coeffs(eq, [a, b, c, d], [x, y, z])
Out[9]: {a: d, b: d, c: -d}

https://docs.sympy.org/latest/modules/solvers/solvers.html#sympy.solvers.solvers.solve_undetermined_coeffs