/* Fractional integration(d) */
new; cls;
y=rndn(100,1);
d=0.4;
print y~fracinteg(y,d);


proc fracinteg(y,d);
   local n,nc,x,j,l,s,k;
   n=rows(y);
   nc=cols(y);
   x=zeros(n,nc);
   l=1;
   do while l<=nc;
      j=n;
      do while j>=1;
         s=c(d,seqa(0,1,j));
         k=1;
         do while k<=j;
            x[j,l]=x[j,l]+s[k]*y[j+1-k,l];
            k=k+1;
         endo;
         j=j-1;
      endo;
      l=l+1;
   endo;
   retp(x);
endp;

proc c(d,j);
   local i,s,k;
   d=-d;
   s=zeros(rows(j),1);
   k=1;
   do while k<=rows(j);
      if j[k]==0;
         s[k]=1;
      elseif j[k]==1;
         s[k]=(-1)^j[k]*(d/1);
      else;
         s[k]=(-1)^j[k]*(d/1);
         i=2;
         do while i<=j[k];
            s[k]=s[k]*((d+1)/i-1);
            i=i+1;
         endo;
      endif;
      k=k+1;
   endo;
   retp(s);
endp;