n4=500;k=.5;G=50; a=(1+k*(G-1))/G; 
x=0:1:n4-1;z=ones(1,n4);
y=1:1:n4;q=zeros(1,n4);
c=zeros(1,n4);s=1/(1-k);t=a^s;q(1)=t;q(2)=t^2*a*(1-a);q(3)=q(2)*t*(1-a);
for j=4:n4
	q(j)=q(j-1)*t*(1-a);
end
p=zeros(1,n4);c(1)=1;
p(1)=t;c(2)=s;p(2)=s*q(2);
for i=3:n4
	d=y(1,1:i-1);
	e=x(1,1:i-1);
	h=z(1,1:i-1);
	f=(s*i).*h;
	l=(f-e)./d;c(i)=prod(l)/i;
	p(i)=c(i)*q(i);
end
loglog(y,p)