% function [degrees,lambda,stats]=CompLeftPoly2c(R,c,n);
%
% Computes degree distributions from the polynomial
% family
%
%  f(x;b,n) =c*x^{2*n}+(1-c)*x^n)
%
% Uses explicit formulas for the Taylor
% coefficients of f and g.
%
% Input: Rate 0<R<1 and parameters 0<c<1,n>0
% 
% Output:
%  degrees=[a_r a_l N M] average and maximal degrees 
%          of right and left degree distributions
%          (N=2*n+1);
%  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,lambda,stats]=CompLeftPoly2c(R,c,n);
%
% Safety bound for M
Mmax=20000;
% Computation of a_r, N
N=2*n+1;Ir=c/(2*n+1)+(1-c)/(n+1);ar=1/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/(2*c);beta=1/c;w01=2*c/(1+c);w02=w01*w01;
A=1-1/n;B=1/2;C=alpha/n*w01;
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 c 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*w02*B.*C 0]+[0 alpha*w01*A.*C];
      A=[A(1) A+1];B=[B(1)+1 B+1/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