#
#    Shank's baby step giant step in GAP
#wdj, 11-2001
##########################################

DL_shanks:=function(a,b,n)
#solves for log_a(b) mod m:
#x where a^x=b mod n 
local m,g,i,y,x,j,S2;
 m:=RootInt(n,2);
 g:=[];
 g[1]:=b;
 for i in [1..m] do
  g[i+1]:=g[i]*a^-m mod n;
 od;
 S2:=List([0..m-1],i->a^i mod n);
 for i in [1..m] do
  if g[i] in S2 then x:=Position(S2,g[i]); return m*(i-1)+x-1; fi;
 od;
end;

###Example
#n:=NextPrimeInt(10000);;
#a:=PrimitiveRootMod(n);;
#b:=Random([2..n-1]);;
#x:=DL_shanks(a,b,n);;
#b=a^x mod n;
### should return true 


dlp_solver:=function(a,b,n)
local m,g,i,y,x,j,S2;
 m:=RootInt(n,2);
 g:=[];
 g[1]:=b;
 for i in [1..m] do
  g[i+1]:=g[i]*a^-m mod n;
 od;
 S2:=List([0..m-1],i->a^i mod n);
 for i in [1..m] do
  if g[i] in S2 then x:=Position(S2,g[i]); return m*(i-1)+x-1; fi;
 od;
end;


#n:=NextPrimeInt(10000);
#10007
#a:=PrimitiveRootMod(n);
#5
#b:=Random([2..n-1]);
#3572
#y:=dlp_solver(a,b,n);
#5688
#a^y mod n;
#3572


dlp_stat:=function(a,b,n)
#enter same a,b,n as in DL_shanks
#returns number of steps the algorithm
#performs
local m,g,i,y,x,j,S2;
 m:=RootInt(n,2);
 g:=[];
 g[1]:=b;
 for i in [1..m] do
  g[i+1]:=g[i]*a^-m mod n;
 od;
 S2:=List([0..m-1],i->a^i mod n);
 for i in [1..m] do
  if g[i] in S2 then x:=Position(S2,g[i]); return i; fi;
 od;
end;


dlp_statistics:=function(a,n,N)
#enter same a,n. returns the average number of
#steps algorithm performs in the
#range 1 < b < N 
local i,b,sum;
 sum:=0;
 for i in [1..N] do
  b:=Random([2..n-1]);
  sum:=sum+(dlp_stat(a,b,n)/N);
 od;
 return sum;
end;


#dlp_statistics(3,113,10);
#63/10