\\% Gaussian normal bases
\\ Author: Joerg Arndt
\\ online at http://www.jjj.de/pari/
\\ version: 2009-September-11 (02:36)

read("vecmanip.gpi"); \\ vecrotl()

gauss_test(n, t)=
{ /* test whether a type-t Gaussian normal basis exists for GF(2^n) */
    local( p, r2, g, d );
    p = t*n + 1;
    if ( !isprime(p), return( 0 ) );
    if ( p<=2, return( 0 ) );
    r2 = znorder( Mod(2, p) );
    d = (p-1)/r2;
\\    print1(" d=", d, "  ");
    g = gcd(d, n);
\\    print1( "  t=",t, "  r2=", r2, "  d=", d, "  g=", g);
    return ( if ( 1==g, 1, 0) );
} /* ------------ */

gauss_nb(n, t)=
{ /* return multiplier matrix for type-t Gaussian normal basis */
  /* returned matrix is over Z and has to be multiplied by Mod(1,2) */
    local(p, r, F, w, x, nh, m, ir, ic);
\\    if ( 0==gauss_test(n, t), print("No type-",t, " GNB for n=",n); return(0));
    p = t*n + 1;
    r = znprimroot(p);  r = r^(n);  /* r has order t */
    F = vector(p-1);
    w = Mod(1, p);
    for (k=0, t-1,
        j = lift(w);
        for (i=0, n-1,
            F[j] = i;
            j+=j;  if (j>=p, j-=p);  /* 2*j mod p */
        );
        w *= r;
    );

    m = matrix(n, n);
    for (i=1, p-2,
        ir = F[p-i];  ic = F[i+1];
        m[ ir+1, ic+1 ] += 1;
    );

    if ( 1==(t%2),
        nh = n/2; /* odd t ==> even n */
        for (i=0, nh-1,
            ir = i;  ic = nh + i;
            ir += 1;  ic += 1;
            m[ir, ic] += 1;
            m[ic, ir] += 1;
        );
    );

    return ( m );
} /* ------------ */


vexp(p, t)=
{
    local( ve, ph, c, s, al, be,  cp, sp, tt );
    tt = 2.0*Pi/p;  \\ angle increment
    c = 1.0;  s = 0.0;  ga = ph;  al = 2.0*(sin(0.5*tt))^2; be = sin(tt);
    ve = vector(p);  ve[1] = 1.0;
    if ( t&1,  /* odd t, need complex values */
        for (j=1, (p-1)>>1,
            tt = c;
            c -= (al*tt+be*s);
            s -= (al*s -be*tt);
            ve[j+1]   = c + I*s;
            ve[p-j+1] = c - I*s;
        );
    , /* even t: can use real values */
        for (j=1, (p-1)>>1,
            tt = c;
            c -= (al*tt+be*s);
            s -= (al*s -be*tt);
            ve[j+1]   = c + s;
            ve[p-j+1] = c - s;
        );
    );
    return( ve );
} /* ----- */


gauss_vpoly(n, t)=
{ /* return field polynomial for type-t Gaussian normal basis
     as polynomial over the complex numbers */

    local(p, r, wk, tk1, tk, a, zp, ve);
    p = n*t + 1;
    r = znprimroot(p)^n;  \\ r has order t (mod p)

    ve = vexp(p, t);  \\ precompute trigonometric values

    zp = 1;
    tk1 = Mod(2,p);  tk = Mod(1,p);
    if ( 0==bitand(t,1),  /* even type */
        for (k=1, n,
            tk *= tk1;  \\ == Mod(2,p)^k;
            wk = 0;
            a = tk;
            for (j=0, t-1,
\\                print("a=", a , " incr=", ve[lift(a)+1]);
                wk += ve[lift(a)+1];
                a *= r;
            );
\\            printp("  a(",k,")=", lift(a), "  w(",k,")=", wk );
            zp *= (x-wk);
        );
    , /* else (odd type) */
        for (k=1, n\2,  \\ note: n/2 times
            tk *= tk1;  \\ == Mod(2,p)^k;
            wk = 0;
            a = tk;
            for (j=0, t-1,
\\                print("a=", a , " incr=", ve[lift(a)+1]);
                wk += ve[lift(a)+1];
                a *= r;
            );
\\            printp("  a(",k,")=", lift(a), "  w(",k,")=", wk );
            \\ use (x-(a+I*b))*(x-(a-I*b)) == x^2 - 2*a*x + (a^2 + b^2):
            zp *= (x^2-2*real(wk)*x+norm(wk));
        );
    );

    return ( zp );
} /* ------------ */



gauss_zpoly(n, t)=
{ /* return field polynomial for type-t Gaussian normal basis
     as polynomial over the complex numbers */
    local(p, r, wk, tk1, tk, a, zp);
    p = n*t + 1;
    r = znprimroot(p)^n;  \\ r has order t (mod p)

    zp = 1;
    tk1 = Mod(2,p);  tk = Mod(1,p);
    for (k=1, n,
        tk *= tk1;  \\ == Mod(2,p)^k;
\\        print("tk=", tk );
        wk = 0;
        a = tk;
        for (j=0, t-1,
\\            print("a=", a );
            wk += exp(2.0*I*Pi*lift(a)/p);
            a *= r;
        );
\\        printp("  a(",k,")=", a, "  w(",k,")=", wk );
        zp *= (x-wk);
    );

    return ( zp );
} /* ------------ */

gauss_poly(n, t)=
{ /* return field polynomial for type-t Gaussian normal basis */
    local(pp, zp);
    local(re, im);
    zp = gauss_zpoly(n, t);

    for (k=0, poldegree(zp),
        re = real(polcoeff(zp,k));
        if ( abs(re-round(re)) > 1.0e-6, print("zpol-real-OUCH:", k, ":  ", re); );
        im = imag(polcoeff(zp,k));
        if ( abs(im) > 1.0e-6, print("zpol-imag-OUCH:", k, ":  ", im); );
    );

    pp = round(real(zp));  /* rounds all coefficients */

\\    print("z(x)=", pp );
\\    printp("roots=", polroots(pp) );
    pp *= Mod(1,2);  /* coefficients modulo 2 */

\\    print("pp=", pp );
\\    forstep (k=poldegree(pp), 1, -1, if ( polcoeff(pp,k),  print1("x^",k" + ")); );
\\    if ( polcoeff(pp,0),  print(1) );

    return( pp );
} /* ------------ */




gauss_poly2(n, t)=
{  /* return field polynomial for type-t Gaussian normal basis */
  /* All computations over GF(2) */
  /* NOTE: slower than the algorithm using complex numbers */
    local(p, M, r, F, t21, t2, Z);
\\    if ( 0==gauss_test(n,t),  return(0) );
    p = t*n + 1;
\\    print("  n=", n, "  t=", t, ":  p=", p);
    r = znprimroot(p)^n;  \\ element of order t mod p
\\    print("  :: r=", r, "  ord(r)=", znorder(r), " ==  t=", t);

\\    M = sum(k=0, p-1, 'x^k);  \\ The polynomial modulus
    M = 'x^p - 1;  \\ Use redundant modulus
    M *= Mod(1,2);  \\ ... over GF(2)


    if ( 1==t,  return( sum(k=0, p-1, 'x^k) ) );  \\ for type 1

\\    print("  :: M =", lift(M));
    F = Mod(1, M);
    t21 = Mod(2,p);  t2 = Mod(1,p);
    for (k=1, n,
\\        print("  :: ------- k=", k);

\\        for(j=0, t-1, print( lift(Mod('x^lift(t2*r^j), M) ) ) );

\\        Z = sum(j=0, t-1, 'x^lift(t2*r^j));  \\ faster (but unclean?)
        Z = Mod( sum(j=0, t-1, 'x^lift(t2*r^j)), M);  \\ faster
\\        Z = sum(j=0, t-1, Mod('x^lift(t2*r^j), M) );  \\ fast
\\        Z = sum(j=0, t-1, Mod('x, M)^lift(t2*r^j) );  \\ slower
\\        print("  :: Z =", lift(lift(Z)));

        F = ('x+Z)*F;

\\        print(" :: w=", sum(j=0, t-1, 'x^lift(t2*r^j) ) );
\\        print(" :: w%=", lift(Mod(sum(j=0, t-1, 'x^lift(t2*r^j) ), sum(k=0, p-1, 'x^k) ) ) );
\\        print("  :: F =", lift(lift(F)) );

        t2 *= t21;
    );
    \\ final reduction for redundant modulus:
\\    M = sum(k=0, p-1, 'x^k);  \\ The polynomial modulus
\\    F = lift( Mod( lift(F), M) );

    \\ final reduction for redundant modulus (simplified):
    F = lift(F);
    if ( 0==polcoeff(F,0), F=sum(k=0, n, (1-polcoeff(F,k))*'x^k) );

    return ( F );
} /* ----- */




normal_mult(a, m, b)=
{ /* multiply a,b in GF(2^n) using multiplier matrix m */
    local( n, t, p );
    n = length(a);
    p = vector(n);
    for (i=1, n,
        p[i] = a*(m*b~);
        a = vecrotl( a );
        b = vecrotl( b );
    );
    return( p );
} /* ------------ */

\\ ==== end of file ====