new; cls;
k=5;   /* k x k matrix */
NN=3; m=6;            /* 0<=x<=NN^m-1 over GF(NN^m) */
x=floor((NN^m)*rndu(k,k));
irp10q1=findirp10b1(m,NN);
q=dectobary(irp10q1[m],NN);
print/rz "x=" x;
print/rz "x_1'=" sinvGF(x,NN,m,q);
print/rz "x*x_1=" mulmatGF(x,sinvGF(x,NN,m,q),NN,m,q);

/*
**  sinvGF.txt - 
** (C) Copyright 2012 Yosuke Amijima. 
**
**  Purpose:    Calculates inverse matrix on GF(NN^m)
**
**  Format:     y=sinvGF(x,NN,m,q)
**
**
**  Input:      x        matrix,  k x k matrix consisting of integers from 0 to NN^m-1
**
**              NN       scalar,  prime in GF(NN^m)
**             
**              m        scalar,  multiple in GF(NN^m) for m=1,2,3,...
**
**              q        vector,  (m+1) x 1 coefficient vector for irreducible polynomial
**
**
**  Output:     y        matrix,   k x k inverse matrix consisting of integers from 0 to NN^m-1
**
**  Notice:     There is no error routine implemented for singular case. 
**
*/
proc sinvGF(x,NN,m,q);
   local k,i,j,l,y,a;
   k=rows(x);
   y=eye(k);
   j=1;
   do while j<=k;
      i=1;
      do while i<=k;
         if i==j;
            a=Kdiv(1,x[i,j],NN,m,q);
            l=1;
            do while l<=k;
               x[i,l]=Kmul(a,x[i,l],NN,m,q);
               y[i,l]=Kmul(a,y[i,l],NN,m,q);
               l=l+1;
            endo;
         endif;
         i=i+1;
      endo;
      i=1;
      do while i<=k;
         if i/=j;
            a=Ksub(0,Kdiv(x[i,j],x[j,j],NN,m,q),NN,m);
            l=1;
            do while l<=k;
               x[i,l]=Kadd(Kmul(a,x[j,l],NN,m,q),x[i,l],NN,m);
               y[i,l]=Kadd(Kmul(a,y[j,l],NN,m,q),y[i,l],NN,m);  
               l=l+1;
            endo;
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   retp(y);
endp;

proc mulmatGF(x,y,NN,m,q);
   local i,j,k,s,z;
   z=zeros(rows(x),cols(y));
   i=1;
   do while i<=rows(x);
      j=1;
      do while j<=cols(y);
         s=0;
         k=1;
         do while k<=cols(x);
            s=Kadd(s,Kmul(x[i,k],y[k,j],NN,m,q),NN,m);
            k=k+1;
         endo;
         z[i,j]=s;
         j=j+1;
      endo;
      i=i+1;
   endo;
   retp(z);
endp;

proc Kadd(a,b,NN,m);
   local s,add;
   s=a|b|(NN^m-1);
   s=dectobary(s,NN);
   add=barytodec((s[.,1]+s[.,2])%NN,NN);
   retp(add);
endp;

proc Kmul(a,b,NN,m,q);
   local x,y,mul;
   x=dectobary(a,NN);
   y=dectobary(b,NN);
   mul=barytodec(barymod0(polymulb(x,y,NN),q,NN),NN);
   retp(mul);
endp;

proc Ksub(a,b,NN,m);
   local sub,k;
   k=0;
   do while k<=NN^m-1;
      if Kadd(b,k,NN,m)==a;
         sub=k;
         break;
      endif;
      k=k+1;
   endo;
   retp(sub);
endp;

proc Kdiv(a,b,NN,m,q);
   local div,k;
   if b==0;
      div=miss(0,0);
   endif;
   k=0;
   do while k<=NN^m-1;
      if Kmul(b,k,NN,m,q)==a;
         div=k;
         break;
      endif;
      k=k+1;
   endo;
   retp(div);
endp;

proc dectobaryvec(y,NN,nb);
   local x,i;
   x=zeros(nb,1);
   i=1;
   do while i<=nb;
      x[i]=floor(y/(NN^(nb-i)));
      y=y-(NN^(nb-i))*x[i];
      i=i+1;
   endo;
   retp(x);
endp;

proc dectobary(x,NN);           /* dec to b-ary */
   local n,z,y;
   if x==0;
      y=0;
   else;
      n=maxc(log(x)/log(NN))+1;
      z=reshape(NN,n,rows(x));
      y=rev(recserrc(x,z)); 
   endif;
   retp(y);
endp;

proc barytodec(y,NN);           /* b-ary to dec */
   if rows(y)==1;
      retp(y);
   else;
      retp(polyeval(NN,y));
   endif;
endp;

proc findirp10b1(q,NN);
   local irp10,i,irp10q1,np;
   if q<1;
      errorlog "ERROR: The degree must be q=1,2,3,...";
      retp(".");
   endif;
   irp10=seqa(NN,1,NN);
   irp10q1=irp10[1];
   i=2;
   do while i<=q;
      np=rows(irp10);
      irp10=nextirp10b(i,NN,irp10);
      irp10q1=irp10q1|irp10[np+1];
      i=i+1;
   endo;
   retp(irp10q1);
endp;

proc nextirp10b(q,NN,irp10);
   local n,x,i,j,a,b;
   if q<2;
      errorlog "ERROR: The degree is at least 2 to be set.";
      retp(irp10);
   endif;
   n=rows(irp10);
   x=seqa(0,1,NN^q);
   j=1;
   do while j<=NN^q;
      a=1|dectobaryvec(x[j],NN,q);
      i=1;
      do while i<=n;
         b=dectobary(irp10[i],NN);
         if barymod0(a,b,NN)==0;
            break;
         endif;
         i=i+1;
      endo;
      if i==n+1;
         irp10=irp10|barytodec(a,NN);
      endif;
      j=j+1;
   endo;
   retp(irp10);
endp;

proc polymulb(a,b,NN);
   local ra,rb,c,i,j,x;
   ra=rows(a); rb=rows(b);
   c=zeros(ra+rb-1,1);
   i=0;
   do while i<=ra+rb-2;
      x=0;
      j=maxc(0|i-ra+1);
      do while j<=minc(rb-1|i);
         x=(x+a[i-j+1]*b[j+1])%NN;
         j=j+1;
      endo;
      c[i+1]=x;
      i=i+1;
   endo;
   retp(c);
endp;

proc barymod0(a,b,NN);
   local ra,rb,q,r,j;
   ra=rows(a);
   rb=rows(b);
   if ra<rb;
      r=zeros(rb-1,1);
      r[rb-1-(ra-1):rb-1]=a;
      retp(r);
   else;
      q=Kdivp(a[1],b[1],NN);
      r=(NN+(a[2:rb]-q*b[2:rb])%NN)%NN;
      j=2;
      do while j<=ra-rb+1;
         r=r|a[rb+j-1];
         if r[1]>=1;
            q=Kdivp(r[1],b[1],NN);
            r=(NN+(r[2:rb]-q*b[2:rb])%NN)%NN;
         else;
            r=r[2:rb];
         endif;
         j=j+1;
      endo;
   endif;
   retp(r);
endp;

proc Kdivp(a,b,NN);
   local div,k;
   if b==0;
      div=miss(0,0);
   endif;
   k=0;
   do while k<=NN-1;
      if (b*k)%NN==a;
         div=k;
         break;
      endif;
      k=k+1;
   endo;
   retp(div);
endp;