About Me

Name and Contact Information

Name: Subham Tibra

University: Indian Institute of Technology, Kharagpur

Major: Mathematics and Computing

email: [email protected]

Github Username: shubhamtibra

Time-zone: UTC+5:30

Personal and Programming Background

Hi, I am Subham Tibra a second year undergraduate student at IIT Kharagpur, India. I'm pursuing a degree in Mathematics and Computing.

I use Ubuntu 14.04 LTS with Sublime Text 3 as my primary text editor,I like it mainly because of it's user-friendly GUI and lot's of other cool features.

I'm proficient in C and Python. I started programming in C but then I got introduced to Python and have been using it since then. I like Python because of it's simplicity and readability and the large number of available open source libraries. With Python I can focus more on the actual algorithms and need not to worry about things like pointers and other stuff.

I know how to use Git and Github very well. However, if I am ever stuck at something, I go to Google and always come back with a solution.

Contributions to SymPy

I started contributing to SymPy in late January 2016 by solving Easy-to-Fix bugs and here is a list of all my contributions so far in chronological order:

(Merged) sinc added to the meijerg tables #10480

(Merged) Improper integral does not evaluate #10445 #10515

(Merged) integrate((sign(x - 1) - sign(x - 2))*cos(x), x) raises TypeError #7163 #10530

(Merged) solveset(domain=S.Reals) rewrites the function in the answer. #10425 #10549

(Unmerged) dsolve() converts floats to integers/rationals #10379 #10562

(Merged) dsolve() converts floats to integers/rationals #10379 #10619

(Open) [WIP] A basic representation of Holonomic Functions and Differential Operator #10793.

The last PR is regarding the project, I created classes for holonomic function and Differential Operator and implemented some algorithms.

Theory and Motivation

A function is called holonomic if it satisfies a linear homogeneous ODE with polynomial coefficients. So there must exist polynomials p0, p1, .... pr such that:

p0(x)f(x) + p1(x)f'(x) + p2(x)f''(x) + · · · + pr(x)f^(r)(x) = 0

Holonomic functions are generalization of hypergeometric functions and represents all algebraic, rational, polynomial and hypergeometric functions. (In Discrete cases, a sequence is holonomic if it satisfies a linear homogeneous recurrence relation with polynomial coefficients.)

Examples of Holonomic Functions:

• exp(x): f'(x) − f(x) = 0

• log(1 − x): (x − 1)f''(x) − f'(x) = 0

• Fibonacci numbers, Factorial, Binomial Coefficients (Holonomic sequences)

• Bessel functions, Hankel functions, Struve functions, Airy functions, Polylogarithms, Elliptic integrals, the Error function, Kelvin functions, Mathieu functions etc..

Not all basic functions we know are holonomic e.g. tan(x), sec(x), Riemann Zeta function etc. are not holonomic.

Holonomic functions are closed under addition, multiplication, differentiation, integration, composition with algebraic functions. These closure properties of holonomic functions are algorithmic , so if we are given two holonomic functions then the function constructed by applying these operations can be obtained using an algorithm.

The function is uniquely defined by the differential equation it satisfies and a set of initial conditions. Thus, the function can be easily represented by a matrix for coefficients of the polynomial in the differential equation and a vector of initial conditions.

A very important concept while dealing with holonomic functions is of Ore Algebra. An univariate ore algebra is a set of polynomials with generator Dx and coefficients in a base ring A. The commutation rule for the generator 'Dx' with the elements of base ring is governed by a endomorphism sigma: A->A and a skew derivation delta: A->A.

Dx*a = sigma(a)*Dx + delta(a)

Elements of Ore Algebra called Ore Polynomials works as operators on functions and can also be used to describe a holonomic function. A particular case which we will need is the Differential Operator in which the sigma is identity map and delta is the partial derivative.The generators of this Ore Algebra acts as derivation on the function e.g.

(x*Dx^2 + Dx + 2).f(x) = (x*f''(x) + f'(x) +2*f(x))

For the differential operator multiplication becomes:

Dx*a = a*Dx + a'(w.r.t. to x)

e.g.

Dx*x = x*Dx + 1
(Dx - 1)*(x*Dx + 1) = x*Dx^2 + (-x + 2)*Dx - 1

An annihilator of a holonomic function is a operator L such that L acted on f, L.f = 0 e.g.

annihilator(exp(x)) = Dx - 1
annihilator(sqrt(1-x)) = 2*(1-x)*Dx - 1
annihilator(sin(x)) = Dx^2 + 1

The fact that these functions cover a lot of usual mathematical functions we generally encounter and they have algorithmic closure properties (i.e. given annihilating operators of two holonomic functions f and g, an annihilating operator of a function h obtained defined by different operations can be determined using different algorithms.) and the easiness of representing them on a computer, makes them a powerful tool in computer algebra. As we know that the product of two hypergeometric functions is not generally hypergeometric, converting a SymPy expression into hypergeometric is not possible every time, but to holonomic is much more feasible as they can also be multiplied.

The Project

The Idea

The project is to create a class that can represent a holonomic function and implement different operations:

Given two holonomic functions f(x) and g(x) which satisfies differential equation/annihilators diffeq1/ann1 and diffeq2/ann2 with initial conditions list l1 and l2 respectively. We want to implement methods to find a differential equation/annihilator diffeq3/ann3 and initial conditions l3 for a function h(x) such that:

• Addition h(x) = f(x) + g(x)

• Multiplication h(x) = f(x)*g(x)

• Power h(x) = f(x)**n, for n = 0, 1, 2....

• Application with a algebraic function h(x) = f(a(x)), where a(x) is some algebraic function of x.

• Differentiation h(x) = f'(x)

• Integration h(x) = Integrate(f(x), x)

Other than these basic operation, we want to compute the function's numerical value at any point, convert a SymPy expression tree to holonomic(may be not the complete expression, as much leaves as possible), since all hypergeometric functions are holonomic we want a method to convert hypergeometric function to holonomic and find a formal power series for a given holonomic function and then convert it to hypergeometric if possible or may be linear combination of those.

We also need to create a class to represent Differential Operators as some of the algorithms are heavily dependent on the operators of holonomic function.(ore_algebra package in sage have a good implementation for this.)

Implementation by Other CAS

• Sage: Sages's package ore_algebra supports different operations on Ore polynomials. It supports addition(computing LCLM), product, power and composition with algebraic function. A lot of our algorithms will be inspired from this package.
• Maple: gfun in Maple has support for addition, multiplication, composition,convert to recurrence relation of Taylor coefficients
• Mathematica: HolonomicFunctions package of mathematica deals with multivariate holonomic functions in general and DifferentialRoot has some basic implementation for it.

Proposed Design

Note The presented prototype for the API is tentative and the final API will be decided after discussing it with my mentor(s).

Differential Operator

I already have built a preliminary class to represent a Differential Operator. A lot of internal algorithms will be using the operations on Differential Operators. Addition of these operator will be same as SymPy Expression, so we just need to add the multiplication rule i.e.

Dx*a(x) = a*Dx + Derivative(a(x), x)

I plan to use SymPy expressions as Differential Operators with Dx as a symbol for the generator and override the Expressions's multiplication rule, I have written the method for multilication of two operators. e.g.

    In []: DifferentialOperator(Dx*x + 1)
    Out []: DifferentialOperator(x*Dx + 2)

    In []: DifferentialOperator(Dx*x + 1)*DifferentialOperator(x*Dx -1)
    Out []: DifferentialOperator(x**2*Dx + 2*x*Dx -2)

Recurrence Operator

I will also define a class to represent recurrence operator which will be similar to Differential Operator with respect to internal storage and working.

    >>> RecurrenceOperator(Sn*n)
    RecurrenceOperator((n+1)*Sn)
    >>> RecurrenceOperator(Sn*(n**2)+1)
    RecurrenceOperator((n**2 + 2*n + 1)*Sn + 1)

Holonomic Function

Since a Holonomic Function is defined by the differential equation/annihilator it satisfies and a set of initial conditions, I propose to make a class which takes the differential equation/annihilator as it's first argument and the variable for the function as the second one, optional arguments will be a list representing the initial conditions in the form of [ f(x0), f'(x0), ... f''..(n-1)times(x0) ] (where n is the degree of diff. eq.), and the point for initial condition if the same is also provided(default is assumed x0 = 0).

This is how the API should look like:

    In []: f = Function('f')(x) ; x, Dx = symbols('x, Dx')
    In []: g = x**2*diff(f) - f + diff(diff(f))
    In []: holonomic(g, x) ### or holonomic(x**2*Dx - 1 + Dx**2, x)
    Out []: holonomic(Dx**2 + x**2*Dx -1, x)

    In []: holonomic(g, x, [1, 0], 1) #### Providing f(1) = 1 and f'(1) = 0
    Out []: holonomic(Dx**2 + x**2*Dx -1, x, f(1) = 1, f'(1) = 0)

Addition, Multiplication and Power

These will be implemented by operator overloading i.e. defining __add__, __mul__ and __pow__. Once the addition and multiplication are done, power is just the multiplication of same function n times. Subtraction will be easy as f(x) - g(x) = f(x) + -1*g(x) and multiplication of holonomic with -1 is quite trivial, one just reverses the sign of initial conditions and the differential equation remains the same.

    In []: p = holonomic(diff(f) - f, x, [1])
    In []: q = holonomic(diff(diff(f)) + f, x, [0, 1])
    In []: p+q
    Out []: holonomic(Dx**3 - Dx*2 + Dx - 1, x, f(0) = 1, f'(0) = 2, f''(0) = 1)

Differentiation and Integration

Integration is pretty straightforward, the annihilator is multiplied by 'Dx' from the right. I have planned to define instance methods diff and integrate in the holonomic class that takes an instance of holonomic function class as it's first parameter and an optional parameter is the variable for which we differentiating/integrating.

    In []: p = holonomic(diff(f) - f, x, [1])
    In []: p.diff(x)
    Out []: holonomic(Dx - 1, x, [1])

    In []: p.integrate(x)
    Out []: holonomic(Dx - 1, x, [1])

Application with an algebraic function

Sage have a nice implementation of this annihilator_of_composition. I plan to define a function applyalg which returns the holonomic function after application with the algebraic function. The algorithm will be quite similar to sage's annihilator_of_composition.

    In []: p = holonomic(Dx - 1, x, [1])
    In []: applyalg(p, -x, x)
    Out []: holonomic(Dx + 1, x, [1])

Finding the recurrence relation

This case will be easy and a

Numerical Computation at a point

Conversion from hypergeometric to holonomic function

This will be quite easy after implementation of applyalg. I will define hypertoholo that takes a hypergeometric function and finds the corresponding holonomic representation

    >>> p = hyper([1,0],[0,1],x)
    >>> hypertoholo(p,x)
    holonomic(Dx - 1, x)

Find formal power series and convert to hypergeometric or a linear combination of them

One function anntorec will be defined to find the recurrence relation of taylor coefficients corresponding to the given holonomic function. Then this recurrence relation can be used to find the coefficents and thus the ratio to eventually convert to hypergeometric functions.

The Implementation

Differential Operator

I already have written a basic class that stores differential operator. I have planned to use Dx as a symbol to represent the generator of polynomial. The operators behaves as normal expressions except for the case of multiplication of 'Dx' with an element of the base ring. Internally a list of coefficent polynomials will be stored along with the expression. Multiplication rule can be easily implemented, see this PR regarding this for the implementation. I haven't yet implemented a planned operation shifttoright to convert all Dx*x to x*Dx + 1, but it's pretty trivial so I'll do it anytime soon.

Recurrence Operator

This will be very similar to that of DifferentialOperator just the multiplication rule needs to be changed. Internally it will also work by storing the polynomial coefficients and the expression.

Multiplication rule:

Sn*a(n) = a(n+1)*Sn

The shifttoright operation will be changed to substitute all Sn*n to (n+1)*Sn.

Holonomic Function

The internal data structure for the annihilator will be the expression itself and a list of polynomial coefficients for each power of 'Dx'. I planned to store both and I think it will make the algorithms more efficient as one can use whichever one is better for a particular algorithm. For initial conditions, a list will be stored inside e.g. [1, 0] x0 = 0 means f(0) = 1 and f'(0) = 0

During initialisation __init__ will set the parameters as the attributes of the holonomic function, In case when diff. eq. is provided, A function that converts differential equations to annihilators will be called and the annihilator will be stored inside. A list of polynomials will also be stored as the coefficients of Dx, and a function to convert the list of polynomials to annihilator listtoann will also be defined due to it's frequent use.

    def  __init__(self, diffeq, x, *args):
        if len(args) is 2:
            self.cond = args[0]
            self.cond_point = args[1]
        elif len(args) is 1:
            self.cond = args[0]
            self.cond_point = 0
        if diffeq isinstance(DifferentialOperator):
            ann = diffeq
        else:
            ann = difftoann(diffeq, x)
        self.annihilator = ann
        listofpoly = []
        for i in range(poly(ann.eq, Dx).degree() + 1):
           listofpoly.append(poly(ann.eq, Dx).nth(i))
        self.listofpoly = listofpoly

Addition

I plan to implement the addition similar to what is implemented in sage's ore_algebra lclm_linalg. Given two holonomic functions f, g and their annihilators F, G, with orders d1, d2 respectively. The LCLM(least common left multiple) of F and G is same as the annihilator of f+g. So one just needs to find the LCLM of F and G. To compute the LCLM, we need a0, a1, ... a(t-d1) and b0, b1, ... b(t-d2) such that:

(a0 + a1*Dx + ... a(t-d1)*Dx^(t-d1))*F = (b0 + b1*Dx + ... b(t-d2)*Dx^(t-d2))*G   -(1)
(...) + (...)*Dx + ..... + (...)Dx^t = (...) + (...)*Dx + ..... + (...)*Dx^t

The coefficients represented by (...) are linear combination of a0,a1.. and b0,b1.... We start by taking t= max(d1, d2) and keep increasing it by 1 in each iteration until a nontrivial solution is found. The LCLM then would be either the L.H.S or R.H.S in equation(1). Now, one just multiplies (a0 + a1*Dx + ... a(t-d1)*Dx^(t-d1)) with F.

F = f.annihilator; G = g.annihilator
frows = [F, Dx*F, ... Dx**(t-d1)*F]
grows = [G, Dx*G, ... Dx**(t-d2)*G]
r = [poly(expr).nth(i) for i in range(t+1)] for expr in (frow+grow)] + [[0 for q in range(t+1)](/sympy/sympy/wiki/poly(expr).nth(i)-for-i-in-range(t+1)]-for-expr-in-(frow+grow)]-+-[[0-for-q-in-range(t+1))
solve_linear_system(r)

First t + 1 - d1 solutions will be a0, a1, ..a(t-d1).

lclm = (a0 + a1*Dx + ... a(t-d1)*Dx^(t-d1))*F

I have already create a basic instance method of addition and is similar to what I have explained here. I'll be working on it to make it more accurate and robust as explained the timeline section.

Multiplication

This algorithm is also inspired from ore_algebra's symmetric_product. To find the annihilator of fg, we need a linear dependency in fg, Dx.(fg), Dx^2.(fg), ... Dx^(k).(fg)(where k is the smallest value for which we get a linear dependent system), the product rule is as follows:

Dx.(fg) = (Dx.f)g + f(Dx.g)
Dx^2.(fg) = (Dx^2.f)g + 2*(Dx.f)(Dx.g) + (Dx^2.g)f

and so on. So very term Dx^k.(fg) can written as a linear combination of (Dx^i.f)(Dx^j.g). The coefficient of (Dx^i.f)*(Dx^j.g) in Dx^k.(fg) can be found by this rule and the terms having Dx^d1.f and Dx^d2.g can be replaced by the terms of lower degree from the annihilator of f and g. Initial value of k is taken to be zero and then it will be increased by 1 in each iteration up to a point where the system has a nonzero solution.

    # coeff[i][j] is the coefficient of (Dx^i.f)*(Dx^j.g) in Dx^k.(fg)
    # initial case for `fg`, only coeff[0][0] = 1

    coeff = [[0 for i in range(d2 + 1)] for j in range(d1 + 1)]; coeff[0][0] = 1
    k = 0
    # finding (Dx^d1).f and (Dx^d2).g in terms of lower powers
    list1 = f.listofpoly; list2 = g.listofpoly

    coefflist1 = [-list1[i]/list1[d1] for i in range(d1)]
    coefflist2 = [-list2[j]/list2[d2] for j in range(d2)]

    linsys = [coeff[i][j] for i in range(d1) for j in range(d2)](/sympy/sympy/wiki/coeff[i][j]-for-i-in-range(d1)-for-j-in-range(d2))
    sol = (Matrix(mat).transpose()).gauss_jordan_solve(b) # b is the column matrix filled with zeroes

    while sol.is_zero:
        k += 1
        coeff = [[0 for i in range(d2 + 1)] for j in range(d1 + 1)]
        for i in range(d1+1):
            for j in range(d2+1):
                if i+j == k:
                    coeff[i][j] = binomial(1, i) # coefficent of (Dx^i.f)(Dx^j.g)
        # replace the terms Dx^d1 and Dx^d2 with lower powers
        for i in range(d1 + 1):
                if not c[i][d2] == zero:
                    for j in range(d2):
                        c[i][j] += coefflist2[j]*c[i][d2]
                    c[i][d2] = 0

            for j in range(d2+1): # not b + 1
                if not c[d1][j] == 0:
                    for i in range(d1):
                        c[i][j] += coefflist1[i]*c[d1][j]
                    c[d1][j] = 0

        mat.append([c[i][j] for i in xrange(d1) for j in range(d2)])
        sol = (Matrix(mat).transpose()).gauss_jordan_solve(b) # b is a column matrix filled with zeroes

Integration

Annihilator will be multiplied by 'Dx' from the right to get the integral.

    In []: integrate(holonomic(Dx**2 + x*Dx -2, x), x)
    Out []: holonomic( Dx**3 + x*Dx**2 - 2*Dx, x)

Application with a Algebraic Function

This is similar to that of multiplication and addition, Let the algebraic function be a(x), then we need a linear dependency between f(a), Dx.(f(a)) ....

Dx^k.(f(a)) can be written as a linear combination of (Dx^i.f)(a) with coefficients as derivatives of a, e.g. Dx.(f(a)) = a'*(Dx.f)(a) (because Dx.(f(a)) = a'*f'(a) = a'*(Dx.f)(a) using the chain rule). Let an annihilator F of function f is given then replacing x by a in the annihilator gives us a linear relation in (Dx^i.f)(a) for i = 0,1...d1 and thus the linear system in (Dx^i.f)(a) can be solved.

    Da = Derivative(a, x)
    subs = [-f.listofpoly(i).subs({x:a})/f.listofpoly(d1).subs({x:a}) for i in range(d1)]
    coeffincomp = [0 for i in range(d1)] # coeffincomp[i] is coeff. of (D^i f)(a) in D^k (f(a))
    coeffincomp[0] = 1 # initially f(a)

The matrix for solving the system will be formed by adding coeffincomp as rows in each iteration. Now, in the next part we increase the length of the linear system by 1. Compute the coefficients of (Dx^i.f)(a) in Dx^k.(f(a))( sage does it) and reduce the power greater than d1 using the list subs and solve to find a linear relation. I have looked at the sage's implementation for this and I'm very sure that I'll be able to implement this.

Initial Conditions

To find out the initial conditions of the resultant holonomic functions when doing these operations. I planned if the conditions are already given at the same point then we can directly find them for the resulting function. If both the conditions are not at the same point, then use sympy's dsolve to see if the given holonomic function can be converted to a elementary function, if yes then compute the initial conditions of both the functions at a same nonsingular point symbolically and add/multiply them. If it isn't possible to get the initial conditions symbolically we'll use our function to calculate numerical value of the holonomic function at a same point and use them.

Find the corresponding recurrence relation in the taylor coefficents

This is just a pretty trivial rewriting, convert each x^k.Dx^j in the annihilator to u(n-k+j)*factorial(n - k +j)/factorial(n - k), where u(n) is the formal power series coefficient about origin. The recurrence relation thus obtained can also be converted to the annihilator in shift operator u(n -k + j) = Sn^(j-k).

The algorithm will iterate throughout the annihilator for each x^k.Dx^j converting it into a list of polynomial coefficents in n of the recurrence relation. The resulting list of polynomials in n thus formed will be converted to the annihilator using listtoann by specifying the generator to be Sn.

Differentiation

When the coefficient of f in the diff. eq. or that of 1 in the annihilator is zero the derivative can be found out by just decreasing each power of Dx by 1. When coefficient of f is not zero, take f to L.H.S with R.H.S as a linear cobination of Dx.f, Dx^2.f .... Now differentiate the equation and repeat the first case.

    if f.listofpoly[0] == 0:
        annofderivative = listtoann(f.listofpoly[1:])
    else:
        list = [f.listofpoly[i])/f.listofpoly[0] for i in range(1, d1)] # starting from 1 not zero
        #sol is the annihilator of the derivative
        sol = (DifferentialOperator(Dx)*DifferentialOperator(listtoann(list))).eq # multiply the annihilator by
                                                                                     #Dx from the LEFT
        sol = sol + 1 #add the term from L.H.S in the annihilator 
    return holonomic(sol, x)

Numerical computation at a point

After getting the recurrence formula, the taylor coefficients can be calculated and the numerical value can be found out at any point using the taylor polynomial. I will implement a method to recursively find the coefficient of the series using matrix multiplication in SymPy and the coefficients thus obtained will be used to create a taylor polynomial up to required order for computing the numerical value.

Power

Once the multiplication is implemented, powers can be computed using recursion in O(log(n)).

    def __pow__(self, power):
        if power == 1:
            return self
        elif power == 0:
            return 1:
        elif power%2 == 1:
            powreduce = self**(power-1)
            return powerreduce*self
        elif power%2 == 0:
            powreduce = self**(power/2)
            return powreduce*powreduce

Hypergeometric to Holonomic

This will use the applyalg to do the job, first the differential operator is formulated from the hypergeometric function and then applyalg finds the annihilator after composition with the z. This is pretty well described by Ondrej here and I have planned to do the same.

Holonomic to Hypergeometric

First the algorithm will try to find the ratio of two consecutive terms from the recurrence relation and get the hypergeometric function. If this isn't possible then try to solve the recurrence using SymPy's rsolve and then try to get the ratio to convert back to hypergeometric and if the above both doesn't work then I will implement a guessing technique, i.e. first guess the result and then test it.

Execution and Timeline

Note I have on other plans for this summer and I am quite sure that I'll be able to give 40-50 hours a week for the project. I am familiar with SymPy's codebase and development workflow. Note2 My classes starts in mid July but that won't be an issue as it will be beginning of the semester and academics load will be less and I'll be working with same pace. Note3 I'll be writing tests and documentation along with implementation of the algorithms.

Community Bonding Period and Week 1

Goal: Discuss more about the project with mentors, know the community. Implement Recurrence Operator and Holonomic Function with addition and integration.

As I already have created a PR for holonomic functions I will work on that to make it more functional and robust and will also create integration along with addition. Meanwhile I'll also the recurrence operator.

I will work on the DifferentialOperator class to solve the issues mentioned in the PR, I will work on normalization of the final answer in the addition algorithm, and will also add tests and full documentation. After that I'll work on the class for recurrence operators, add it's multiplication rules and all. Integration is easy thus I have also planned to implement it during this period. I will also define some small functions which will be frequently used such as difftoann( convert a differential equation to annihilator, will be needed in preprocessing) and listtoann (convert a list of polynomials to annihilators).

Week 2 and 3

Goal: Implement multiplication, power and derivative

I will implement multiplication, power and differentiation algorithms during this period. The major part will be implementation of multiplication, once that is done power is trivial. Derivatives are also quite easy, as I have explained the algorithm in implementation.(support for initial condition will be done after implementation of numerical computation) I will submit a PR for this at the end of this period

Week 4 and 5

Goal: Composition with Algebraic function and Recurrence relation.

In this period I will start by implementing the composition with algebraic functions and then the algorithm to find the recurrence relation of taylor coefficients and then submit a PR for this.

Week 6 and 7

Goal: Coefficients of expansion(formal power series), numerical Computation and support initial conditions

During this period, First I'll implement the method for finding the taylor coefficients by the recursive matrix multiplication method using the recurrence relation. After that is done I'll formulate the function to compute the numerical values using the taylor polynomial. As determining the coefficients would already have been implemented, it will be an easy task. Now I'll add the support for initial conditions to the operations implemented so far.

Week 8

Goal: Fix bugs and complete any left out task

I have planned to keep this week for solving any bugs that might have been generated so far and complete any tasks if left out due to may be any reason.

Week 9 and 10

Goal: Convert Hypergeometric to Holonomic and Holonomic to Hypergeometric

The conversion of hypergeometric to holonomic is pretty trivial and must completed in a short time. I then will focus on the implementation of converting Holonomic function to a hypergeometric one. This part may take some time as converting holonomic to a linear combination of hypergeometric will require some serious thinking and work.

Week 11 and 12

Goal: Converting an expression tree to Holonomic

This part may take some time so I have allotted two weeks to it. I will be working on this method TODO