2940.sagews

# celine's method as per "A=B" Ch.4
# input: F(n,k) is doubly hypergeometric
# That is, F(n+1,k)/F(n,k) and F(n,k+1)/F(n,k) are rational functions
# output: recurrence formula for F(n,k) in the form
# \sum_{i=0}^I \sum_{j=0}^J a_i_j(n)*F(n-j,k-i)=0

var('n k')

F=binomial(n,k)*binomial(2*k,k)*(-2)^(n-k)

IMAX=5
JMAX=5

celine(F)

I=0,J=1 I=1,J=0 I=0,J=2 I=1,J=1 I=2,J=0 I=0,J=3 I=1,J=2 I=2,J=1 Recurrence found: F(n,k) + 0 *F(n-0,k-1) + 2*(2*n - 1)/n *F(n-1,k-0) + -2*(2*n - 1)/n *F(n-1,k-1) + 4*(n - 1)/n *F(n-2,k-0) + -8*(n - 1)/n *F(n-2,k-1) =0. (1, 0, 2*(2*n - 1)/n, -2*(2*n - 1)/n, 4*(n - 1)/n, -8*(n - 1)/n)

# celine's method

def celine(F):
    
    for (I,J) in modified_cartesian_product(IMAX,JMAX):
        
        if (I,J)==(0,0): continue
        
        print "I=%d,J=%d"%(I,J)
    
        
        # finding rational functions F(n-i,k-j)/F(n,k)
        
        Rational=[[n for j in xrange(J+1)] for i in xrange(I+1)]

        for (i,j) in cartesian_product_iterator([xrange(I+1),xrange(J+1)]):
            
            Rational[i][j]=(F.subs(n=(n-i),k=(k-j))/F).factorial_simplify()
        
        #print Rational
        
        
        # determining the sum_{i=0}^I sum_{j=0}^J a_i_j F(n-i,k-j)
        
        a = [[var('a_%d_%d'%(i,j)) for j in xrange(J+1)] for i in xrange(I+1)]
        
        SUM=0

        for i in xrange(I+1):
            for j in xrange(J+1):
                SUM+=a[i][j]*Rational[i][j]
    
        # finding the numerator and trying to solve it to be 0
        
        numerator=(SUM.full_simplify()).numerator()
        coeff_k=[numerator.coefficients(k)[i][0] for i in xrange(numerator.degree(k)+1)]
    
        
        M=matrix(SR, len(coeff_k), (I+1)*(J+1), 0)
        

        for i in xrange(len(coeff_k)):
            
            for j in xrange((I+1)*(J+1)):
        
                j1=int(j / (J+1))
                j2=j % (J+1)
        
                M[i,j]=coeff_k[i].coefficient(a[j1][j2])
        
        
        Ker=M.right_kernel()
        
        if Ker.dimension()==0: continue
        
        t=Ker.random_element()

        while t==0:
            t=Ker.random_element()
        
        t0=t[0]
        
        for i in xrange(len(t)):
            t[i]=(t[i]/t0).full_simplify()
        
        #print t
        
        print
        print "Recurrence found:"
        
        for i in xrange(len(t)):
            i1=int(i/(J+1))
            i2=i% (J+1)
            if i==0:
                print "F(n,k)"
            else:
                print "+",t[i],"*F(n-%d,k-%d)"%(i1,i2)
        
        print "=0."
        print
        
        return t

# modified order for the elements in the cartesian product of xrange(IMAX),xrange(JMAX)

def modified_cartesian_product(A,B):
    
    sumMAX = A+B
    
    list=[]
    
    for sum in xrange(sumMAX+1):
        
        for i in xrange(sum+1):
            
            if i>A: continue
            
            j=sum-i
            
            if j>B: continue
            
            list.append((i,j))
    
    return list