(14:20) gp > p=nextprime(random(2^100))
%1 = 897329727452362169370115077469
(14:20) gp > q=nextprime(random(2^100))
%2 = 538129143948171430946490118399
(14:20) gp > N=p*q
%3 = 482879278073185639238091091438664461654691221595888867252131
(14:20) gp > Phi=(p-1)*(q-1)
%4 = 482879278073185639238091091437229002783290687995572262056264
(14:20) gp > u=3
%5 = 3
(14:21) gp > gcd(u,Phi)
%6 = 3
(14:21) gp > u=5
%7 = 5
(14:21) gp > gcd(u,Phi)
%8 = 1
(14:21) gp > e=1/u%Phi
%9 = 96575855614637127847618218287445800556658137599114452411253
(14:21) gp > contfrac(e/N+0.)
%10 = [0, 5]
(14:21) gp > \p100
   realprecision = 105 significant digits (100 digits displayed)
(14:21) gp > contfrac(e/N+0.)
%12 = [0, 5, 67278734026291537184226560459, 1, 1, 3, 4, 4, 1, 5, 1, 2, 1, 2, 2, 4, 2, 1, 24, 3, 21, 10, 1, 31, 6, 1, 1, 1, 5, 3, 1, 3, 4, 1, 16, 1, 8, 1, 17, 2, 7, 2, 2, 1, 1, 12, 2]
(14:21) gp > contfracpnqn([0,5])
%13 = 
[1 0]

[5 1]

(14:22) gp > u=random(2^10)
%14 = 86
(14:22) gp > u=random(2^30)
%15 = 877903165
(14:22) gp > gcd(u,Phi)
%16 = 1
(14:23) gp > e=1/u%Phi
%17 = 351744164240070134141071444890729698619859099845349328600941
(14:23) gp > e/N+0.
%18 = 0.7284308526214277858310261360068100199288943277889566448659264375495246814535825681669090802246095740
(14:23) gp > contfrac(%)
%19 = [0, 1, 2, 1, 2, 6, 1, 3, 2, 1, 1, 6, 2, 3, 14, 3, 1, 1, 7, 4, 2, 1, 599192341237, 1, 2, 1, 10, 2, 1, 5, 5, 6, 2, 13, 1, 4, 1, 1, 1, 3, 1, 29, 2, 3, 1, 1, 8, 3, 1, 1, 1, 1, 3, 3, 1, 2, 1, 2, 29, 2, 8, 4, 1, 3, 1, 1, 112, 5, 1, 34, 4, 1, 2, 1, 8, 2, 6, 6, 1, 1, 3, 3, 7, 3, 6, 2]
(14:23) gp > contfracpnqn([0,1,2,1,2,6,1,3,2,1,1,6,2,3,14,3,1,1,7,4,2,1])
%20 = 
[639491751 442239983]

[877903165 607113196]

(14:23) gp > u
%21 = 877903165
(14:23) gp > N
%22 = 482879278073185639238091091438664461654691221595888867252131
(14:25) gp > Phi
%23 = 482879278073185639238091091437229002783290687995572262056264
(14:25) gp > u=random(2^150)
%24 = 557542599667929055063851694278237497597796257
(14:28) gp > gcd(u,Phi)
%25 = 1
(14:28) gp > e=1/u%Phi
%26 = 115334167036184581887113937769632238105008341067821954552513
(14:28) gp > e/N+0.
%27 = 0.2388467931289120811102737029773090547224695730317507823185833315658962795237733169821183941133347824
(14:28) gp > contfrac(%)
%28 = [0, 4, 5, 2, 1, 4, 1, 3, 2, 1, 4, 1, 2, 3, 1, 1, 5, 1, 1, 4, 1, 1, 19, 6, 2, 1, 3, 3, 23, 2, 1, 1, 1, 4, 15, 1, 2, 4, 1, 1, 2, 18, 1, 1, 3, 1, 1, 1, 1, 3, 3, 1, 1, 1, 1, 1, 2, 3, 3, 2, 1, 1, 1, 1, 6, 4, 1, 24, 14, 2, 1, 1, 1, 2, 2, 1, 1, 28, 1, 1, 37, 1, 2, 34, 1, 1, 1, 1, 19, 2, 10, 1, 4, 1, 29, 1, 2, 3, 1, 64, 1, 3, 3, 2, 1, 2, 1, 9, 2, 41, 1, 9, 3]
(14:28) gp > e=3
%29 = 3
(14:31) gp > gcd(e,Phi)
%30 = 3
(14:31) gp > e=5
%31 = 5
(14:31) gp > gcd(e,Phi)
%32 = 1
(14:31) gp > u=1/e%Phi
%33 = 96575855614637127847618218287445800556658137599114452411253
(14:31) gp > up=ceil((N+1)/e)
%34 = 96575855614637127847618218287732892330938244319177773450427
findprime(size)=p=nextprime(size):while(1,if(isprime(2*p+1),break,p=nextprime(p+1))):2*p+1:
(15:13) gp > findprime(100)
%35 = 227
(15:13) gp > N=227
%36 = 227
(15:13) gp > znstar(227)
%37 = [226, [226], [Mod(2, 227)]]
(15:14) gp > g=Mod(2,227)
%38 = Mod(2, 227)
rhoCard(N,g)=w=g:z=g:a=1:b=1:while(1,if(lift(w)<N/2,w=w^2:a=2*a):if(lift(w)>=N/2,w=w*g:a=a+1):if(lift(z)<N/2,z=z^2:b=2*b):if(lift(z)>=N/2,z=z*g:b=b+1):if(lift(z)<N/2,z=z^2:b=2*b):if(lift(z)>=N/2,z=z*g:b=b+1):if(w==z,break)):return(a-b):
(15:14) gp > rhoCard(N,g)
%39 = -144987136
(15:14) gp > %%(N-1)
%40 = 0
(15:14) gp > rhoCard(N,g)
%41 = -144987136
(15:14) gp > factor(%)
%42 = 
[-1 1]

[2 10]

[7 1]

[113 1]

[179 1]

rhoDL(N,g,h)=w=g:ai=1:bi=0:z=g:aj=1:bj=0:nb=0:while(1,nb++:if(lift(w)<N/3,w=w*h:bi=bi+1):if(lift(w)>=N/3&&lift(w)<2*N/3,w=w^2:ai=2*ai:bi=2*bi):if(lift(w)>=2*N/3,w=w*g:ai=ai+1):if(lift(z)<N/3,z=z*h:bj=bj+1):if(lift(z)>=N/3&&lift(z)<2*N/3,z=z^2:aj=2*aj:bj=2*bj):if(lift(z)>=2*N/3,z=z*g:aj=aj+1):if(lift(z)<N/3,z=z*h:bj=bj+1):if(lift(z)>=N/3&&lift(z)<2*N/3,z=z^2:aj=2*aj:bj=2*bj):if(lift(z)>=2*N/3,z=z*g:aj=aj+1):if(w==z,break)):return((ai-aj)/(bj-bi)):
(15:24) gp > N
%43 = 227
(15:24) gp > g
%44 = Mod(2, 227)
(15:24) gp > r=random(N):
(15:25) gp > h=g^r
%46 = Mod(96, 227)
(15:25) gp > rhoDL(N,g,h)
%47 = -3463/1020
(15:25) gp > %%(N-1)
  ***   impossible inverse modulo: Mod(2, 226).
(15:26) gp > %%((N-1)/2)
%48 = 51
(15:26) gp > g^51%N
%49 = Mod(96, 227)
(15:26) gp > p=findprime(1000)
%50 = 2027
(15:27) gp > g=Mod(2,N)
%54 = Mod(2, 2027)
(15:27) gp > r=random(N):
(15:27) gp > h=g^r
%56 = Mod(330, 2027)
(15:27) gp > rhoDL(N,g,h)
%57 = -24080057/45797757
(15:27) gp > %%(N-1)
%58 = 1537
(15:28) gp > g^1537
%59 = Mod(330, 2027)
(15:28) gp > h
%60 = Mod(330, 2027)
(15:28) gp > nb
%61 = 42
(15:28) gp > sqrt(N)
%62 = 45.02221673796171489440354679688639044581140079623630379011764209389486990614939451472468535474054532