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

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