#  Input: u - polynomial over Z/(p)[x]
#         f - the modulos polynomial over Z/(p)[x]
#
# Output: Rem(u^mm,f) mod p  over Z/(p)[x]
#
myPowmod:=proc(u,mm,f,x,p)
   local m,w,v;
   m:=mm;  # a copy of mm which we will pull of the bottom bit
   w:=Rem(u,f,x) mod p;   # w starts at u^1 mod (f,p)
   v:=1;
   while (m <> 0) do
      if irem(m,2) = 1 then
         v:=Rem(v*w,f,x) mod p;
      fi;
      w:=Rem(w*w,f,x) mod p;   # w := w^2 mod (f,p)
      m:=iquo(m,2);
   od;
   v;
end:
