Looking for a Primitive Roots program

[roadrunner]

Below is a primitive roots program, for prime numbers only, that I was working on a while ago but had to put it aside because real work got in the way. I recall that somewhere above around 1000 it started giving incorrect answers; I'm not sure why.

Code:

#cas
n_to_z_mod_p(n,z,p):=
BEGIN
 return n^z MOD p;
END;
#end

EXPORT smallest_p_root(p)
BEGIN
//debug;
 if p>1000 then return "that hurts my brain";end;
 local a,b,c,j,k;
 a:=p-1;
 b:=ifactors(a); 
 c:=length(b);
 for j from 2 to a do
  for k from 1 to c step 2 do 
   if n_to_z_mod_p(j,a/b(k),p) == 1 then
    break;
   end;
   if k == c-1 then return j;end;
  end;
 end; 
END;

EXPORT p_root(p)
BEGIN
//debug
 if p>1000 then return "that hurts my brain";end;
 local a,b,c,d,j,k;
 a:=p-1;
 b:=ifactors(a); 
 c:=length(b);
 d:={};
 for j from 2 to a do
  print();
  print("checking " + j);
  for k from 1 to c step 2 do 
   if n_to_z_mod_p(j,a/b(k),p) == 1 then
    break;
   end;
   if k == c-1 then
    d(0):=j;
   end;
  end;
 end;
 return d; 
END;