evalPolyMat:=proc(f,x,alpha)
   local n,m,k,BabySteps,GiantSteps,i,j,g,h;

   n:=degree(f,x);
   m:=ceil(sqrt(n));
   k:=ceil(n/m) + 1;

   # divide f into k pieces of degree <m
   for j from 0 to k-1 do
      g[j]:=0;
      for i from 0 to m-1 do
         g[j]:=g[j]+coeff(f,x,j*m+i)*x^i;
      od;
      # print(g[j]);
   od;

   # Compute the baby steps alpha^i, i=0..m-1
   BabySteps[0]:= 1;
   for i from 1 to m do
      BabySteps[i]:=BabySteps[i-1] . alpha;
   od;

   # Compute the giant steps alpha^(m*j), j=0..k-1
   GiantSteps[0]:=1;
   for j from 1 to k-1 do
      GiantSteps[j]:=GiantSteps[j-1].BabySteps[m];
   od;

   # evaluate each of the g[j]'s at alpha, using the baby steps
   for j from 0 to k-1 do
      h[j]:=add(coeff(g[j],x,i)*BabySteps[i],i=0..m-1);
   od;

   # combine the evaluations
   add(h[j].GiantSteps[j],j=0..k-1);

end: