/* coefficients for q-Super Catalan number */
new; cls;
n=5;
print/rz SCqcoef(n);


/*
**  SCqcoef.txt - coefficients for q-Super Catalan number.
** (C) Copyright 2012 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Computes coefficients for polynomials of q-Super Catalan number.
**
**  Format:     cq=SCqcoef(n);
**
**  Input:      n    scalar,    integer 0,1,2, . . .
**
**  Output:     cq   matrix,    k x 1 vector of coefficients as in
**
**              cq[1]+cq[2]*q+cq[3]*q^2+cq[4]*q^3+. . .
**
**  Notice:     This program computes coefficients from those of large Schroder 
**              number by factoring out {1,1} or 1+q.
**
**  Reference:  M.Ishikawa, H.Tagawa, and J.Zeng, "A q-analogue of Catalan Hankel determinants",
**              RIMS K\^oky\^uroku Bessatsu, B11 (2009), 19--42, arXiv:1009.2004v2 [math.CO]
**              (http://arxiv.org/pdf/1009.2004v2.pdf)
*/
proc SCqcoef(n);
   local index,c,j,k,sum,cindex,sum_1,cq;
   if n==0;
      retp(1);
   else;
      index=zeros(n+1,1);
      index[1]=1; 
      cindex=1;
      c=1; 
      sum_1=c;
      j=0;
      do while j<=n-1;
         sum=zeros(2*(j+1)-1,1)|sum_1;
         k=0;
         do while k<=j+1-1;
            if k==0 and j/=0;
               sum=msum(sum,zeros(2*1*(j+1-1),1)|sum_1);
            elseif k==0 and j==0;
               sum=msum(sum,sum_1);
            elseif k==j;
               sum=msum(sum,c[cindex[k]+1:cindex[k+1]]);
            else;
               sum=msum(sum,zeros(2*(k+1)*(j+1-1-k),1)|mmult(c[cindex[k]+1:cindex[k+1]],c[cindex[j+1-1-k]+1:cindex[j+1-1-k+1]]));
            endif;
            k=k+1;
         endo;
         index[j+2]=rows(sum);
         c=c|sum;
         sum_1=sum;
         cindex=cumsumc(index);
         j=j+1;
      endo;
      cq=factor11(sum);
      retp(cq);
   endif;
endp;

proc mmult(x,y);
   local z,i,j;
   z=zeros(rows(x)+rows(y)-1,1);
   j=1;
   do while j<=rows(y);
      z[j:rows(x)+j-1]=z[j:rows(x)+j-1]+y[j]*x;
      j=j+1;
   endo;
   retp(z);
endp;

proc msum(x,y);
   local z;
   z=zeros(maxc(rows(x)|rows(y)),1);
   z[1:rows(x)]=x;
   z[1:rows(y)]=z[1:rows(y)]+y;
   retp(z);
endp;

proc factor11(x);
   local k,i,c;
   k=rows(x);
   c=zeros(k-1,1);
   c[1]=x[1];
   i=2;
   do while i<=k-2;
      c[i]=x[i]-c[i-1];
      i=i+1; 
   endo;
   c[k-1]=x[k];
   retp(c);
endp;