title: Is there an efficient algorithm to compute the Jacobsthal matrix or quadratic character in GF(q)?
body: <p>Is there an efficient algorithm to compute the <a href="https://en.wikipedia.org/wiki/Jacobsthal_matrix" rel="nofollow noreferrer">Jacobsthal matrix [WP]</a> or equivalently the quadratic character χ in GF(<em>q</em>),</p>
  <p><em><strong><strong>J</strong></strong></em> [ <em>i, j</em> ] = χ ( <em>i - j</em> ) = 0  if  <em>i = j</em>  else  1  if  <em>i - j</em>  is a square in GF(<em>q</em>) else -1,</p>
  <p>where <em>i, j</em> run over the elements of GF(<em>q</em>)?</p>
  <p>The order of the elements &lt;=&gt; rows/columns does not really matter, so it's mainly to know whether an element of GF(<em>q</em>) is a square.
  Unfortunately, when <em>q</em> = <em>p<sup> n</sup></em> with <em>n</em> &gt; 1, one cannot just take <em>i, j</em> ∈ <strong>Z</strong>/<em>q</em><strong>Z</strong> (which works well iff <em>q</em> is a prime &lt;=&gt; <em>n</em> = 1).</p>
  <p>On the other hand, implementing arithmetics in GF(<em>q</em>) appears a nontrivial task to me, at least the naive way (constructing an irreducible polynomial <em>P</em> of degree <em>n</em> over <strong>Z</strong>/<em>p</em><strong>Z</strong> and implementing multiplication through multiplication of polynomials modulo <em>P</em>...).</p>
  <p>The problem is easily solved in Python using the <code>galois</code> package <a href="https://galois.readthedocs.io/" rel="nofollow noreferrer">(see here)</a>, 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 <a href="https://en.wikipedia.org/wiki/Paley_construction" rel="nofollow noreferrer">Paley construction [WP]</a>, 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).</p>
  <p>Since we only need to know which elements are squares, I hoped there might be an efficient way to determine that.</p>
  <p>EDIT: Let me clarify again that the question is whether there exists an efficient way of implementing this function (for arbitrary <em>q = p<sup> k</sup></em>) without implementing general arithmetic in GF(<em>q</em>). It's not difficult to solve the problem using dedicated software: For example, Python's <code>galois</code> package provides the <code>is_quadratic_residue()</code> 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, <em>i.e.</em>, <code>issquare(Mod(i-j, p^k))</code>, will usually yield incorrect results for when <em>k</em> &gt; 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 <em><strong>J J</strong></em><sup>T</sup> which should equal <em>q <strong>I</strong> - U</em> (for <em>p</em> &gt; 2) where <em>U</em> is the &quot;all 1s&quot; matrix.)</p>