/* q-Super Catalan number */
new; cls;
n=5; q=2;
print/rd scatalanq(n,q);

/*
**  scatalanq.txt - q-Super Catalan number.
** (C) Copyright 2012 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Computes q-Super Catalan number.
**
**  Format:     s=scatalanq(n,q);
**
**  Input:      n    scalar,    integer 0,1,2, . . .
**
**              q    scalar,    integer 1,2, . . .
**
**  Output:     s    scalar,    positive integer
**
**  Notice:     This program computes coefficients from coefficients of large Schroder 
**              number by factoring out {1,1} or 1+q. When q=1, the result is the same
**              as Super Catalan number. It may be easily exploded for large q. 
**              It may be very useful for applications of prime number and cryptography 
**              when it is a prime.
*/
proc scatalanq(n,q);
   local index,c,j,k,sum,cindex,sum_1,s;
   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;
      sum=factor11(sum);
      j=1;
      do while j<=rows(sum);
         sum[j]=sum[j]*q^(j-1);
         j=j+1;
      endo;
      s=sumc(sum);
      retp(s);
   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;