> rp:=rand(100000..999999):

> p:=rp():
> u:=[1,[[2,1],[3,1]]];

                      u := [1, [[2, 1], [3, 1]]]

> n:=nops(u[2]);

                                n := 2

> for i from 1 while (isprime(p)=false or u[2][n][1]<9999)do
> p:=rp():
> u:=ifactors(p-1):
> end do:
> 
> 
> p;
> 
> ifactor(p-1);

                                504209


                                4
                             (2)   (31513)

> ifactors(p-1);

                      [1, [[2, 4], [31513, 1]]]

> q:=u[2][n][1];

                              q := 31513

> rg:=rand(p-1):
> g:=3:
> for i from 1 while (modp(g^q,p)<>1) do
> g:=rg():
> end do:
> g;

                                81080

> g^q mod p;

                                  1

> 
>