well, you could try a brute force approach since you have φ(p^k - 1) generators.

so you write a subroutine that calculates the order of a (mod p^k - 1), by just checking a, a^2, etc., and test if that equals p^k - 1.

(something like:

a = 2

(#)n = 2

(*)while n < p^k -1, calculate a^n (mod p^k - 1),

(**)....while a^n (mod p^k - 1) > 1,

...............n-->n+1, goto (*)

..........else order(a) = n, goto (***)

....else, output (a)

(***)a-->a+1, goto (#)

or some such equivalent)

it's not pretty, but it should work.