new; cls; /* Calculates polynomial a % polynomial b in base NN */
a={1,2,3};
b={1,1,2};
NN=3;
print/rz barymod0(a,b,NN)';


proc barymod0(a,b,NN);   /* base NN is prime and b[1]/=1 is allowed */
   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;