> restart;
> 
> cheb := proc(n,a,b,fnu)
> 
>   #     Input:    n   -- the number of recursion coefficients desired
>   #               a,b -- arrays of dimension 2*n-1 to be filled with the
>   #                      values of  a(k-1),b(k-1), k=1,2,...,2*n-1
>   #               fnu -- array of dimension  2*n  to be filled with the
>   #                      values of the modified moments  fnu(k-1), k=1,2,
>   #                      ...,2*n
>   #     Output:   alpha,beta-- arrays containing, respectively, the
>   #                      recursion coefficients  alpha(k-1),beta(k-1),
>   #                      k=1,2,...,n, where  beta(0)  is the total mass.
>   #               s  --  array containing the normalization factors
>   #                      s(k)=integral [pi(k)(x)]**2 dlambda(x), k=0,1,
>   #                      2,...,n-1.
> 
>   local alpha, beta, s, s0, s1, s2, nd, l, k, lk;
> 
>   #
>   # The arrays  a,b  are assumed to have dimension  2*n-1, the arrays
>   # fnu,s0,s1,s2  dimension  2*n.
>   #
> 
>   nd := 2*n;
> 
>   #
>   # Initialization
>   #
> 
>   alpha := array(1..n);
>   beta := array(1..n);
>   s := array(1..n);
>   s0 := array(1..nd);
>   s1 := array(1..nd);
>   s2 := array(1..nd);
> 
>   alpha[1] := a[1] + fnu[2]/fnu[1];
>   beta[1]  := fnu[1];
>   if (n=1) then 
>     RETURN([alpha=eval(alpha),beta=eval(beta)]);
>   fi;
> 
>   s[1] := fnu[1];
>   for l from 1 to nd do
>     s0[l] := 0;
>     s1[l] := fnu[l];
>   od;
> 
>   #
>   # Continuation
>   #
> 
>   for k from 2 to n do
>     lk := nd - k + 1;
>     
>     for l from k to lk do
> 
>       s2[l] := s1[l+1] - (alpha[k-1] - a[l])*s1[l] - beta[k-1]*s0[l] + b[l]*s1[l-1];
>       if (l=k) then 
>         s[k] := s2[k]; 
>       fi;
>     od;
> 
>     #
>     # Compute the alpha- and beta-coefficient
>     #
> 
>     alpha[k] := a[k] + s2[k+1]/s2[k] - s1[k]/s1[k-1];
>     beta[k] := s2[k]/s1[k-1];
>     for l from k to lk do
>       s0[l] := s1[l];
>       s1[l] := s2[l];
>     od;
> 
>   od;
> 
>   RETURN(['alpha'=eval(alpha), 'beta'=eval(beta),'s'=eval(s)]);
> end:
> 
> m1 := r -> simplify((1/Pi)*cos(Pi/6)*2^r*GAMMA((r+1+1/3)/2)*GAMMA((r+1-1/3)/2));
> m2 := r -> simplify((2^(2/3)*sqrt(Pi))/(GAMMA(2/3)) * (r!*GAMMA(r+1/3))/((6*r+1)*2^r*GAMMA(r+1/6)));
>
> 
> N := 40:
> a   := [seq(0,l=1..2*N-1)]:
> b   := [seq(0,l=1..2*N-1)]:
> fnu1 := [seq(m1(r),r=0..2*N-1)]:
> fnu2 := [seq(m2(r),r=0..2*N-1)]:
> 
> res:=cheb(N,a,b,fnu1):
> 
> alpha:=op(2,op(1,res)):
> beta:=op(2,op(2,res)):
> if (N = 1) then
>   s:=[]:
> else
>   s:=op(2,op(3,res)):
> fi:
> 
> Digits := 35;
>
> 
> # Alpha values
> print(`Alpha Values`);
> for k from 1 to N do
>   print(evalf(alpha[k]));
> od;
>
>
> # Beta values
> print(`Beta Values`);
> for k from 1 to N do
>   print(evalf(beta[k]));
> od;
>
>
> # s values
> print(`s Values`);
> if (N <> 1) then
>   for k from 1 to N do
>     print(evalf(alpha[k]));
>   od;
> fi;
>