SUM of two Holonomic Functions using Euclidean Algorithm - sympy/sympy GitHub Wiki

A holonomic function can be represented by the differential equation it satisfies. Ore polynomials can be used to represent such differential equations. The particular case of Ore polynomials which we will deal with are the differential operators. Let's take for example a holonomic function say exp(x): The differential equation it satisfies is f(x) - f'(x) = 0 This can also be written as (Dx - 1).f(x) = 0 Where Dx-1 is the differential operator(also an Ore polynomial). The operator Dx satisfies the multiplication rule Dx*a = a*Dx + a'. For instance Dx*x = x*Dx + 1. An operator L acting on f such that L.f = 0 is called an annihilator of f. Thus the annihilator of exp(x) is Dx-1. Let annihilator of two functions f and g be L and M respectively, then annihilator of f + g is given by LCLM(L, M) where LCLM stands for least common left multiple which can be obtained by using the Extended Euclidean Algorithm.

def lclm(A, B):
    u, v = euclid(A, B)
    return u*A.normalize()

def euclid(A, B):
    a11, a12, a21, a22 = 1, 0, 0, 1
    while not B.is_zero():
        (C, q, alpha, beta) = rem_seq(A, B)
        r = r2; bInv = beta**-1
        a11, a12, a21, a22 = a21, a22, bInv*(alpha*a11 - q*a21), bInv*(alpha*a12 - q*a22)
    c = a21.denominator().lcm(a22.denominator())
    return (c*a21, c*a22)

def rem_seq(A, B):
    orddiff = A.order() - B.order()
    alpha = B.leading_coefficient()
    newRem = (alpha*A).quo_rem(B)
    beta = newRem[1].content()
    r2 = newRem[1].map_coefficients(lambda p: p//beta)
    return (B, r2, newRem[0], alpha, beta)

def quo_rem(A_ORIG, B_ORIG):
    A = A_ORIG
    B = B_ORIG
    if (A.order() < B.order()):
        return (B.zero(), A)
    cfquo = one
    quo = zero
    op = lambda x, y: x/y
    orddiff = A.order() - B.order()
    while(orddiff >= 0):
        currentOrder = p.order()
        cfquo = op(A.leading_coefficient(), B.leading_coefficient()) * Dx**(orddiff)
        quo = quo+cfquo
        A = A - cfquo*B
        if A.order()==currentOrder:
            A = A_ORIG
            B = B_ORIG
            op = lambda x,y:x/y
        orddiff = A.order() - B.order()
    return (quo,A)

where content is greatest common divisor of all the coefficients.