update:= proc(z,a,b) global p,g,beta; local zz,aa,bb;
     aa:=a: bb:=b:
     if z<p/3 then  zz:=z*beta mod p; bb:=b+1 mod(p-1)
     elif z> 2*p/3 then  zz:=z*g mod p; aa:=a+1 mod(p-1)
     else zz:=z*z mod p; aa:=2*a mod (p-1); bb:=2*b mod (p-1) 
     end if:
     (zz,aa,bb):
end proc:


p:=971; beta:=564 ; g:=43;       # INPUT VALUES; g should be primitive mod p

a0:=0; b0:=0; aX:=a0: bX:=b0:
              aY:=a0: bY:=b0: # INITIALIZATION
x:=g &^aX * beta&^bX mod p:
y:=g &^aY * beta&^bY mod p:
UB:=10^3:                    # UPPER BOUND ON NO. OF ITERATIONS

(x,aX,bX):=update(x,aX,bX):
(y,aY,bY):=update(y,aY,bY):
(y,aY,bY):=update(y,aY,bY):
i:=1:                        # ONE ROUND OUT OF LOOP

while not(x=y) and i<UB do   # ITERATIONS
(x,aX,bX):=update(x,aX,bX):
(y,aY,bY):=update(y,aY,bY):
(y,aY,bY):=update(y,aY,bY):
i:=i+1:
od:

print(`i=`,i);
if i<UB then dB:=bX-bY mod (p-1); dA:=aY-aX mod (p-1) end if;
if igcdex(dB,p-1,'inv','foo')=1 then DL:=(inv*dA) mod (p-1):
                                         print(`DL=`,DL,`beta=`,beta,
                                         `g^DL=`,g &^DL mod p) 
elif igcdex(dB,p-1,'inv','foo')=2 then dA:=dA/2; dB:=dB/2;
                igcdex(dB,(p-1)/2,'inv','foo');
                DL1:=(inv*dA) mod (p-1)/2:
                DL2:= DL1 + (p-1)/2;
                        print(`gcd was 2 and DL1=`,DL1);
                        print(`beta=`,beta, `g^DL1=`,g &^DL1 mod p); 
                        print(`gcd was 2 and DL2=`,DL2);
                        print(`beta=`,beta, `g^DL2=`,g &^DL2 mod p); 
                             
else print(`gcd larger than 2`,gcd(dB,p-1)) end if;


p:=75011; g:=499; beta:= g &^ 45000 mod p; 
beta:= g &^ 45001 mod p; beta:= g &^ 45016 mod p;