The problem statement can be rephrased as follows: Find all N*M matrices with entries 
from {0,...,C-1} such that no two rows are the same, and no two columns are the same.

First of all, let's forget about rows. How many matrices are there such that all columns
are distinct? There are C^N possibilities for the first column: it can look any way we like.
Then there are (C^N)-1 possibilities for the second column: it can look any way we like,
except it cannot match the first column. And so on.

Now let's count all matrices that have both distinct rows and distinct columns. We will
use dynamic programming on the number of rows. Let tables[n] be the number of n*M tables 
with distinct rows and distinct columns. We will compute the values tables[1] through
tables[N], which is the answer we seek.

How can we compute tables[n]? We will take all n*M tables with distinct columns, and 
we will subtract those that have some identical rows.

How many of those are there? Each such table has between 1 and n-1 distinct rows.
Let's call the number of distinct rows k. How many n*M tables with distinct columns
and exactly k distinct rows are there? The answer has a simple combinatorial form:
stirling2[n][k] * tables[k]. This is because if you choose any of the tables[k] k*M
tables with distinct rows and columns, you know how the k row types look like.
Now you can produce multiple n*M tables: one for each partition of its n rows into
k equivalence classes, which is precisely what Stirling numbers of the second kind
count.

Pseudocode:

for n = 1 to N:
    // count all tables with distinct columns: (C^n) * (C^n -1) * ... * (C^n -M+1)
    tables[n] = 1
    for m = 1 to M:
        tables[n] *= (C^n + 1 - m)

    // for each k<n subtract all tables with exactly k different row types
    for k = 1 to n-1:
        tables[n] -= tables[k] * stirling2[n][k]

return tables[N]

Other ways of counting are possible as well. For instance the following one:
Let f(N,M,C) be the answer. You can convince yourself that the answer has to
be independently polynomial in N and M (because of counting via inclusion
and exclusion). In other words, f(N,M,C) must be equal to the following sum:
sum_{x,y} C^{xy} p_x p_y. How can we determine the unknown coefficients p_i?
By taking M=1 and noticing that f(N,1,C) = N*(N-1)*...*(N-C+1). It turns out 
that the coefficients are precisely the (signed) Stirling numbers of the first kind.