Is there an efficient algorithm to compute the Jacobsthal matrix or quadratic character in GF(q)?

Question

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 ⁿ 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 J^T which should equal q I - U (for p > 2) where U is the "all 1s" matrix.)

Answer 1

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

ex = e0 * x modulo polynomial

                    1                  q = 1 = quotient
          -----------
1 1 1 2 2 | 1 0 0 0 0                  poly | ex
            1 1 1 2 2                  poly * q
            ---------
              2 2 1 1                  remainder

                    2                  q = 2 = quotient
          -----------
1 1 1 2 2 | 2 2 1 1 0                  poly | ex
            2 2 2 1 1                  poly * 2
            ---------
              0 2 0 2                  remainder

Basic math code with initialization:

typedef unsigned char BYTE;

/* GFS(3) */
#define GFS 3
/* GF(3^2) */
#define GF 81
/* alpha = 1x + 0 */
#define ALPHA 3

typedef struct{                         /* element of field */
    int     d;                          /* = dx^3 + cx^2 + bx + a */
    int     c;
    int     b;
    int     a;
}ELEM;

typedef struct{                         /* extended element of field */
    int     e;                          /* = ex^4 + dx^3 + cx^2 +bx + a */
    int     d;
    int     c;
    int     b;
    int     a;
}ELEMX;

/*----------------------------------------------------------------------*/
/*      GFAdd(i0, i1)                                                   */
/*----------------------------------------------------------------------*/
static int GFAdd(int i0, int i1)
{
ELEM e0, e1;
    e0 = aiI2E[i0];
    e1 = aiI2E[i1];
    e0.d = (e0.d + e1.d);
    if(e0.d >= GFS)e0.d -= GFS;
    e0.c = (e0.c + e1.c);
    if(e0.c >= GFS)e0.c -= GFS;
    e0.b = (e0.b + e1.b);
    if(e0.b >= GFS)e0.b -= GFS;
    e0.a = (e0.a + e1.a);
    if(e0.a >= GFS)e0.a -= GFS;
    return (((((e0.d*GFS)+e0.c)*GFS)+e0.b)*GFS)+e0.a;
}

/*----------------------------------------------------------------------*/
/*      GFSub(i0, i1)                                                   */
/*----------------------------------------------------------------------*/
static int GFSub(int i0, int i1)
{
ELEM e0, e1;
    e0 = aiI2E[i0];
    e1 = aiI2E[i1];
    e0.d = (e0.d - e1.d);
    if(e0.d < 0)e0.d += GFS;
    e0.c = (e0.c - e1.c);
    if(e0.c < 0)e0.c += GFS;
    e0.b = (e0.b - e1.b);
    if(e0.b < 0)e0.b += GFS;
    e0.a = (e0.a - e1.a);
    if(e0.a < 0)e0.a += GFS;
    return (((((e0.d*GFS)+e0.c)*GFS)+e0.b)*GFS)+e0.a;
}

/*----------------------------------------------------------------------*/
/*      GFMpy(i0, i1)           i0*i1       using logs                  */
/*----------------------------------------------------------------------*/
static int GFMpy(int i0, int i1)
{
    if(i0 == 0 || i1 == 0)
        return(0);
    return(aiExp[aiLog[i0]+aiLog[i1]]);
}

/*----------------------------------------------------------------------*/
/*      GFDiv(i0, i1)           i0/i1                                   */
/*----------------------------------------------------------------------*/
static int GFDiv(int i0, int i1)
{
    if(i0 == 0)
        return(0);
    return(aiExp[(GF-1)+aiLog[i0]-aiLog[i1]]);
}


/*----------------------------------------------------------------------*/
/*      GFPow(i0, i1)           i0^i1                                   */
/*----------------------------------------------------------------------*/
static int GFPow(int i0, int i1)
{
    if(i1 == 0)
        return (1);
    if(i0 == 0)
        return (0);
    return(aiExp[(aiLog[i0]*i1)%(GF-1)]);
}

/*----------------------------------------------------------------------*/
/*      GFMpyA(i0)              i0*ALPHA    using low level math        */
/*----------------------------------------------------------------------*/
/* hard coded for elements of size 4 */
static int GFMpyA(int i0)
{
ELEM e0;
ELEMX ex;
int q;                                  /* quotient */
    e0 = aiI2E[i0];                     /* e0 = i0 split up */
    ex.e = e0.d;                        /* ex = e0*x */
    ex.d = e0.c;
    ex.c = e0.b;
    ex.b = e0.a;
    ex.a = 0;
    q = ex.e;
/*  ex.e -= q * pGFPoly.aata[0] % GFS;  ** always == 0 */
/*  if(ex.e < 0)ex.d += GFS;            ** always == 0 */
    ex.d -= q * pGFPoly.data[1] % GFS;
    if(ex.d < 0)ex.d += GFS;
    ex.c -= q * pGFPoly.data[2] % GFS;
    if(ex.c < 0)ex.c += GFS;
    ex.b -= q * pGFPoly.data[3] % GFS;
    if(ex.b < 0)ex.b += GFS;
    ex.a -= q * pGFPoly.data[4] % GFS;
    if(ex.a < 0)ex.a += GFS;
    return (((((ex.d*GFS)+ex.c)*GFS)+ex.b)*GFS)+ex.a;
}

/*----------------------------------------------------------------------*/
/*      InitGF  Initialize Galios Stuff                                 */
/*----------------------------------------------------------------------*/
static void InitGF(void)
{
int i;
int t;

    for(i = 0; i < GF; i++){    /* init index to element table */
        t = i;
        aiI2E[i].a = t%GFS;
        t /= GFS;
        aiI2E[i].b = t%GFS;
        t /= GFS;
        aiI2E[i].c = t%GFS;
        t /= GFS;
        aiI2E[i].d = t;
    }

    pGFPoly.size = 5;           /* init GF() polynomial */
    pGFPoly.data[0] = 1;
    pGFPoly.data[1] = 1;
    pGFPoly.data[2] = 1;
    pGFPoly.data[3] = 2;
    pGFPoly.data[4] = 2;

    t = 1;                      /* init aiExp[] */
    for(i = 0; i < GF*2; i++){
        aiExp[i] = t;
        t = GFMpyA(t);
    }

    aiLog[0] = -1;              /* init aiLog[] */
    for(i = 0; i < GF-1; i++)
        aiLog[aiExp[i]] = i;
}


/*----------------------------------------------------------------------*/
/*      main                                                            */
/*----------------------------------------------------------------------*/
int main()
{
    InitGF();
    return(0);
}

---

Code to display a list of primitive polynomials for GF(3^4)

    pGFPoly.size = 5;           /* display primitive polynomials */
    pGFPoly.data[0] = 1;
    pGFPoly.data[1] = 0;
    pGFPoly.data[2] = 0;
    pGFPoly.data[3] = 0;
    pGFPoly.data[4] = 1;
    while(1){
        i = 0;
        t = 1;
        do{
            i++;
            t = GFMpyA(t);}
        while(t != 1);
        if(i == (GF-1)){
            printf("pGFPoly:    ");
            ShowVector(&pGFPoly);}
        pGFPoly.data[4] += 1;
        if(pGFPoly.data[4] == GFS){
            pGFPoly.data[4]  = 1;
            pGFPoly.data[3] += 1;
            if(pGFPoly.data[3] == GFS){
                pGFPoly.data[3]  = 0;
                pGFPoly.data[2] += 1;
                if(pGFPoly.data[2] == GFS){
                    pGFPoly.data[2]  = 0;
                    pGFPoly.data[1] += 1;
                    if(pGFPoly.data[1] == GFS){
                        break;}}}}}

This produces this list:

1 0 0 1 2   The one normally used x^4 + x + 2
1 0 0 2 2
1 1 0 0 2
1 1 1 2 2   I used this to test all 5 terms
1 1 2 2 2
1 2 0 0 2
1 2 1 1 2
1 2 2 1 2