SUM of two Holonomic Functions using Euclidean Algorithm

Shubham Tibra edited this page Mar 1, 2016 · 10 revisions

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.

Clone this wiki locally
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.
Press h to open a hovercard with more details.