% function [degrees,rho,lambda,stats]=CompLeftSelf(R,b,n);
%
% Computes degree distributions from the self-similar
% family
%
% f(x;b,n) =(1-b)*x^n/(1-b*x^n)=(1-b)*sum_{k>=1} b^{k-1}x^{k*n}
%
% Uses explicit formulas for the Taylor
% coefficients of f and g.
%
% Input: Rate 0<R<1 and parameters 0<b<1,n>0
%
% Output:
% degrees=[a_r a_l N M] average and maximal degrees
% of right and left degree distributions
% (N=n*N0+1);
% rho 1xN0 (compacted) array of right degrees,
% rho(k) corresponds to x^{k*n}, coefficients
% in between are set to 0
% lambda 1xM array of left degrees
% stats [delta,epsilon,mu,Delta,Lambda2]
% distances to capacity 1-R, comparison to
% lower bounds, and fraction of degree 2 variable
% nodes Lambda2=lambda2*a_l/2
%
% Author: P. Oswald
% Last changed: 10/4/2002
function [degrees,rho,lambda,stats]=CompLeftSelf(R,b,n);
%
% Safety bounds for degrees
Nmax=20000;Mmax=20000;
% Computation of a_r, N, and rho
% by taking a long enough piece of the Taylor series of f
N0=max([1 floor(log(n*0.00000001)/log(b))]);N=N0*n+1;
if N0>Nmax
disp('N exceeds Nmax, set N=Nmax (to avoid, decrease b and/or n)');
N0=n*floor(Nmax/n)+1;
end
rho=[];c=(1-b);Ir=0;
for k=1:N0
rho=[rho c];Ir=Ir+c/(n*k+1);c=c*b;
end
c1=sum(rho);rho=rho/c1;ar=c1/Ir;
%
% Computation of a_l, M, lambda, and delta
% by implementing Step (2) and (3) of Algorithm 1.
% This is based on hand-derived recursions for the Taylor coefficients
% g_k of g, which lead to recursions for s_k and sigma_k and finally
% to finding t such that hat(I)(t)=1/a_l where a_l=a_r*(1-R)
% is the average degree of the left (variable) nodes.
al=ar*(1-R);alpha=(1-b);beta=b;w01=1;w02=w01*w01;
A=1-1/n;B=2;C=alpha/n*w02;
c=C;sk=c;C=C/2;c1=C;sigmak=c1;lambda=[];
if al<2
disp('Average degree of left nodes too small (decrease R or increase b or n)');
M=0;lambda=[];stats=[];
else
if al==2
lambda=1;M=2;delta=c;
else
m=2;q=sigmak/sk;
while q*al>1&m<Mmax
lambda=[lambda c];
C=[beta*w01*B.*C 0]+[0 alpha*w02*A.*C];
% [m;C']
A=[A(1) A+1];B=[B(1)+1 B+2];
c=sum(C);sk=sk+c;
m=m+1;C=C/m;c1=c/m;sigmak=sigmak+c1;q=sigmak/sk;
end
if q*al>1&m==Mmax
disp('M too big (reduce a or increase R or change Mmax)');
M=m;lambda=[];stats=[];
else
M=m;t1=(al*sigmak-sk)/(al*c1-c);
delta=(sk-t1*c);lambda=[lambda (1-t1)*c]/delta;
end
end
%
% Computation of output
epsilon=1-delta/(1-R);
a=ar*log(R);mu=a/log(epsilon);Delta=epsilon/exp(a);
Lambda2=lambda(1)*al/2;
stats=[delta epsilon mu Delta Lambda2];
end
degrees=[ar al N M];
% End CompLeftSelf.m