For[i = RandomPrime[5000], PrimeQ[2*i + 1] == False,
  i = RandomPrime[300]];
q = i
p = 2*i + 1
For[g = RandomInteger[{1, p}];, (PowerMod[g, q, p] != 1 &&
     PowerMod[g, 2, p] != 1) == False,
  g = RandomInteger[{1, p}];];(*teorema dave paulius*)
g