% script to generate statistics on gcd's and number of steps required to
% find them by the Euclidean algorithm, using M randomly generated pairs
% (a,b) between 1 and N
gcdcount=zeros(1,1);
stepcount=zeros(1,1);
for j=1:M
    a=floor(N*rand)+1;
    b=floor(N*rand)+1;
    
    if(a<b)  % switch a and b, if necessary, so that b<a
        c=a; % hang onto the value of a
        a=b; % even while replacing a with b
        b=c; % now replace b with a
    end
    count=0; % reset counter
    while(b>0)
        q=floor(a/b);
        r=a-q*b;
        a=b;
        b=r;
        count=count+1;
    end
    if(length(gcdcount) >=a)
        gcdcount(a)=gcdcount(a)+1; %increment appropriate counter
    else
        gcdcount(a)=1;
    end
    if(length(stepcount) >=count)
        stepcount(count)=stepcount(count)+1; %increment appropriate counter
    else
        stepcount(count)=1;
    end
end
subplot(2,1,1)
plot(gcdcount/M)
title('Distribution of gcds')
subplot(2,1,2)
plot(stepcount/M)
title('Distribution of algorithm steps')
subplot(111) % means the figure window returns to normal single-graph behavior