Is there an efficient algorithm to compute the Jacobsthal matrix [WP] or equivalently the quadratic character χ in GF(q),
J [ i, j ] = χ ( i - j ) = 0 if i = j else 1 if i - j is a square in GF(q) else -1,
where i, j run over the elements of GF(q)?
The order of the elements <=> rows/columns does not really matter, so it's mainly to know whether an element of GF(q) is a square. Unfortunately, when q = p n with n > 1, one cannot just take i, j ∈ Z/qZ (which works well iff q is a prime <=> n = 1).
On the other hand, implementing arithmetics in GF(q) appears a nontrivial task to me, at least the naive way (constructing an irreducible polynomial P of degree n over Z/pZ and implementing multiplication through multiplication of polynomials modulo P...).
The problem is easily solved in Python using the galois package (see here), but this is quite heavy artillery which I'd like to avoid to deploy.
Of course dedicated number theory software may also have GF arithmetics implemented. But I needed this just to produce Hadamard matrices through the Paley construction [WP], so I'd like to be able to compute this without using sophisticated software (and anyway I think it would be interesting to know whether there's a simple algorithm to do this).
Since we only need to know which elements are squares, I hoped there might be an efficient way to determine that.
EDIT: Let me clarify again that the question is whether there exists an efficient way of implementing this function (for arbitrary q = p k) without implementing general arithmetic in GF(q). It's not difficult to solve the problem using dedicated software: For example, Python's galois package provides the is_quadratic_residue() function which immediately gives the matrix elements - in spite of its name, since quadratic residues (mod p^k) aren't the same as squares in GF(p^k): Indeed, default modular arithmetic, i.e., issquare(Mod(i-j, p^k)), will usually yield incorrect results for when k > 1. For example, in G(2^k) every element is a square, but 2 and 3 aren't squares mod 2^2). A crude check is to compute J JT which should equal q I - U (for p > 2) where U is the "all 1s" matrix.)
Here is a basic math table setup for GF(3^4) based on 1 x^4 + 1 x^3 + 1 x^2 + 2 x + 2. At the end of this answer is a brute force search for any primitive polynomial (where powers of 3 map all non-zero elements). Numbers are stored as integer equivalents, for example, x^3 + 2 x + 1 = 1*(3^3) + 2*(3) + 1 = 16, so I store this as 16. Add and subtract just map from integer to vector and back. Multiply and divide use exp and log tables. Exp table is generated by taking powers of 3 (multiplying by x). Log table is the reverse mapped exp table. InitGF initializes the exp table using GFMpyA (multiply by alpha == multiply by x). Showing the math starting at integer 27 = 1 x^3 * x, showing the long hand division of
Basic math code with initialization:
Code to display a list of primitive polynomials for GF(3^4)
This produces this list: