proc Sudoku(intmat M)
{
   list E=createE();
   ring R=32003,x(1..81),dp;
   ideal F=createF();
   ideal G=createG(F,E);
   G=addPreass(G,M);
   option(redSB);
   G=std(G);
   intmat N=prepareRes(G);
   return(N);
}

proc createE()
{
   int i,j,k,l,a,b;
   list E;
   for(j=1;j<=9;j++)
   {
      for(k=1;k<=9;k++)
      {
         i=(j-1)*9+k;
         a= j mod 3;if(a==0){a=3;}
         b= k mod 3;if(b==0){b=3;}

         for(l=k+1;l<=9;l++)
         {
             E[size(E)+1]=list(i,i+l-k);
         }
         for(l=j+4-a;l<=9;l++)
         {
             E[size(E)+1]=list(i,(l-1)*9+k);
         }
         if(a!=3)
         {
             E[size(E)+1]=list(i,j*9+k-b+1);
             E[size(E)+1]=list(i,j*9+k-b+2);
             E[size(E)+1]=list(i,j*9+k-b+3);
         }
         if(a==1)
         {
             E[size(E)+1]=list(i,(j+1)*9+k-b+1);
             E[size(E)+1]=list(i,(j+1)*9+k-b+2);
             E[size(E)+1]=list(i,(j+1)*9+k-b+3);
         }

      }
   }
   return(E);
}

proc createF()
{
   int i;
   ideal F;
   poly p=(var(1)-1);
   for(i=2;i<=9;i++)
   {
      p=p*(var(1)-i);
   }
   F[1]=p;
   for(i=2;i<=nvars(basering);i++)
   {
      F[i]=subst(p,var(1),var(i));
   }
   return(F);
}

proc createG(ideal F, list E)
{
   int i;
   ideal G=F;
   for(i=1;i<=size(E);i++)
   {
      G[size(G)+1]=(F[E[i][1]]-F[E[i][2]])/(var(E[i][1])-var(E[i][2]));
   }
   return(G);
}

proc addPreass(ideal G, intmat M)
{
   int i,j,k;
   for(i=1;i<=9;i++)
   {
      for(j=1;j<=9;j++)
      {
         if(M[i,j]!=0)
         {
            k=(i-1)*9+j;
            G[size(G)+1]=var(k)-M[i,j];
         }
      }
   }
   return(G);
}

proc prepareRes(ideal G)
{
   intmat M[9][9];
   int i,j,k;
   for(i=1;i<=size(G);i++)
   {
      j=((i-1) div 9) +1;
      k=i mod 9;if(k==0){k=9;}
      M[j,k]=int(leadcoef(-G[size(G)-i+1][2]));
   }
   return(M);
}

/*
intmat B[9][9] =  0,6,0,1,0,4,0,5,0, 0,0,8,3,0,5,6,0,0, 2,0,0,0,0,0,0,0,1,
                  8,0,0,4,0,7,0,0,6, 0,0,6,0,0,0,3,0,0, 7,0,0,9,0,1,0,0,4,
                  5,0,0,0,0,0,0,0,2, 0,0,7,2,0,6,9,0,0, 0,4,0,5,0,8,0,7,0;
print(Sudoku(B));
*/