Interpolating a Bivariate Polynomial over a Finite Field

Question

Let $F=GF(p)$ be a finite field, $p$ prime and write $F^\times=\{x_1,\ldots,x_n\}$.

I'm trying to implement an earlier version of Sudan's list-decoding algorithm for Reed Solomon Codes

$$RS(d+1,n+1)=\{(f(x_1),\ldots,f(x_n))|f\in F[x], \deg f \leq d\}$$ I would like to find a way to do the following in Mathematica:

Let $(y_1,\ldots,y_n)\in F^n$.

(1) Is there a way in Mathematica to solve a homogeneous linear system where the number of variables is greater than the number of constraints? Sudan's algorithm involves finding a polynomial $Q(x,y)$ (not identical to 0) with coefficients in $F$ such that $Q(x_i,y_i)=0$ for all $i$. The algorithm specifies that every monomial $cx^jy^k$ in $Q$ satisfy $j+kd\leq m+ld$, where $m=\left\lceil \frac{d}{2} \right\rceil-1$, $l=\left\lceil\sqrt{\frac{2(n+1)}{2}}\right\rceil-1$ and has $(m+1)(l+1)+d_f\binom{l+1}{2}$ unknowns. This entails solving the linear system $Q(x_1,y_1)=0,\ldots,Q(x_n,\ldots,y_n)$. I tried using

Solve[Q(x_1,y_1)==0 && ... && Q(x_n,y_n)==0, <list of the unknowns>, Modulus->n+1]

however, this methods doesn't work if the number of unknowns is greater than the number of constraints ($n$), which is the case here.

(2) Factor $Q(x,y)$ into irreducible factors. Is there a way of doing this given that $Q$ has coefficients over a finite field.

Here's a link to Sudan's paper.

EDITS: As suggested, here is the actual code I used:

The elements of the field are in the vector $x=$ fieldx. The input is the codeword $y=$ rw, and my desired output is a polynomial $Q(x,y)$ such that $Q(x_i,y_i)=0$ for all $i$.

q = 5;
k = 3;
n = q - 1;
d = k - 1;
m = Ceiling[d/2] - 1;
l = Ceiling[Sqrt[2(n + 1)/2]] - 1;
t = d*Ceiling[Sqrt[2(n + 1)/2]] - Floor[d/2];
Var = Array[f, (m + 1)(l + 1) + d*Binomial[l + 1, 2]];
Eqn = Array[b, Num];
Num = Length[Var];
Cons = Array[h, n];
Poly = Array[g, n];
fieldx = {1, 2, 3, 4};
rw = {1, 2, 2, 3};
For[i = 0, i < n, i++; h[i] = {}];
For[i = 0, i < n, i++;
For[j = -1, j < l, j++;
  For[k = -1, k < m + (l - j)*d, k++;
    Cons[[i]] = Join[Cons[[i]], {(fieldx[[i]]^k)*(rw[[i]]^j)}]]
  ]];
Hom = {}; For[i = 0, i < n, i++; Hom = Join[Hom, {{0}}]];
For[i = 0, i < n, i++; Poly[[i]] = 0; Poly[[i]] = Total[Cons[[i]].Eqn]];
Eqn = Array[b, Num];
LinearSolve[Cons, Hom]
Solve[Poly[[1]] == 0 && Poly[[1]] == 2 && Poly[[2]] == 0 && Poly[[
3]] == 0 && Poly[[4]] == 0, Eqn]

This produces the matrix Cons

{{1, 1, 1, 1, 1, 1, 1, 1, 1}, {1, 2, 4, 8, 16, 2, 4, 8, 4}, {1, 3, 9, 27, 81, 
2, 6, 18, 4}, {1, 4, 16, 64, 256, 3, 12, 48, 9}}

MORE EDITS: While the algorithm itself has several other details, I am mainly interested in

FIRST QUESTION: How can I solve a linear system with more unknowns than constraints (which was addressed by Michael E2's answer below). Details on the algorithm are in Sudan's paper. I don't think I should discuss those details here as they would make my post extremely long (for instance, why I took the values of m and l as above).

As to my SECOND QUESTION: How can I factor a bivariate polynomial over F - it turns out that the latest version of Mathematica can do this. Unfortunately, my current copy of Mathematica can only factor polynomials in one variable. When I try to evaluate, say, Factor[6*x*y^2 + 2, Modulus -> q], I get the following error:

Factor::facmm: Factoring multivariate polynomials with respect to a modulus is not yet implemented.

Answer 1

The rest of the solutions are spanned by

NullSpace[Cons]
(* -> {{78, -133, 57, -8, 0, 0, 0, 0, 6},
       {120, -202, 93, -17, 0, 0, 0, 6, 0},
       {48, -82, 33, -5, 0, 0, 6, 0, 0},
       {18, -37, 15, -2, 0, 6, 0, 0, 0},
       {24, -50, 35, -10, 1, 0, 0, 0, 0}} *)

I assumed that was the question. You don't actually state what the question is.

Other remarks

There are several issues in your code, too many for me to address at this moment.

First avoid using capitals to begin names of your own symbols. You may cause conflicts with Mathematica's built-in functions.

Second, Mathematica has efficient ways of doing things that Java or C lack. For instance, your line

 Hom = {}; For[i = 0, i < n, i++; Hom = Join[Hom, {{0}}]];

does the same thing as

 Hom = ConstantArray[{0}, n]

That makes Hom into an n by 1 matrix:

 { {0}, {0}, .... }

Perhaps you meant to have a dimension n vector? If so, that could be done with

 Hom = ConstantArray[0, n]
 (* { 0, 0, ... } *)

Finally, there are other ways to improve your code. There are many wonderful, efficient, built-in functions for you to discover. There's a nice bunch of advice that site users have compiled for people just starting out in Mathematica that point out the more important ones.