GSoC 2016 Application Srajan Garg: Polynomial manipulation & Multivariate Polys in Cpp - wrat/sympy GitHub Wiki

#Personal Background

###Details

Name : Srajan Garg

University : Indian Institute of Technology, Bombay

Email : [email protected]

GitHub : srajangarg

Blog : srajangarg.github.io

Time-zone : IST (UTC+5:30)

Age : 19

###About Me

I am a second year undergraduate student pursuing Bachelor of Technology in Computer Science and Engineering at Indian Institute of Technology, Bombay. I was introduced to computer programming in grade 10, and I've never stopped since. I have explored a lot of areas related to programming ranging from Django/Flask python web-frameworks to C++ physics libraries like Box2D. I've also actively participated in hacking style CTFs and competitive programming contests with my batch mates, apart from software development.

I was already aware of Sympy, and was introduced to SymEngine by Sumith. I was fascinated by the fact that SymEngine might be the (optional) core of Sympy in the near future. SymEngine was also the first open source project I started contributing to, and it's been a great experience to get to know open source inside out and has been a great learning experience thus far.

###Platform details

OS : Ubuntu 15.10

Hardware Configuration : i7 5600U/ 8GB

Editor : SublimeText 3

###Programming experience I have been programming in C++ for the past two years, and have greatly improved my knowledge from when I began. I am also fairly experienced in Python having used it for about one and a half years. I use C++ in competitive coding exercises. I have also developed an 8 Ball Pool Game using the Box2D physics library. I have also written various python apps/ web-apps which I made mainly for personal use. I am also familiar with git for version control, and also familiar with most of the workflow associated with it. Apart from this I am also experienced with makefiles, PHP, MATLAB/Octave and Latex.

###Contributions to SymEngine/Sympy

  • Implemented as_numerator_denominator function for SymEgnine

  • Implemented a fast parser for text to expression conversions within SymEngine, also made small changes afterwards

  • Added simplifications for trigonometric equations of the form sin(arcsin(x)) = x

  • Implemented the min and max functions for SymEngine

  • Implemented Horner's method for evaluating univariate polynomials

  • Redid the mul_poly function of the Univariate class, and also made smaller changes (pending)

  • Worked on creating creating code quality checks through a python script

  • Fixed a minor bug in the LambertW function

  • Fixed bug related to printing of sets and intervals in Sympy (pending)

#Project Overview

SymEngine currently has two univariate polynomial class, which lack much needed functionality. I intend to amend this by adding the missing functionality by fully furnishing the existing polynomial class. Also, I would like to implement a univariate class for rationals coefficients (discussed below) and also improve the existing general univariate class, while working alongside the UC Davis Team.

Also a good CAS must have fast implementations of multivariate polynomials. I would like to add a efficient multinomial class for SymEngine, following the same pattern as that of the univariate case. I would also like to play with the idea of multivariate polynomials being the super class of univariate polynomials. This would ensure that multivariate polynomials are as fast as our univariate polynomials, and univariate polynomials are just a special case with only one variable.

Implementation of a fast and robust polynomial class forwards the progress in other features like Solvers, Integrals and Series Expansions too.

#Univariate Polynomials

###Overview

Currently in SymEngine, a UnivariateIntPolynomial class exists, which captures integer coefficients. Also, another UnivariatePolynomial has been recently introduced by the UC Davis team. This class deals with polynomials with all other types of coefficients including rationals and expressions. The Int polynomial class will have many efficient algorithms for various functions, which are not possible to implement for a general coefficient class.

###Multiplication

Currently, the UnivariateIntPolynomial class, uses the Kronecker substitution trick to multiply two Univariate Polynomials. To be more precise, the algorithm evaluates the polynomials at a single sufficiently large power of 2. After multiplying the values obtained, it is able to interpolate the coefficients of the new polynomial. These slides call this the KS1 algorithm.

There exists another variant of this algorithm, called KS2. The basic idea remains the same, the only difference being that in KS2, polynomial evaluation happens at four carefully chosen points. The coefficients of the resulting polynomial are then interpolated using combinations of products from the previous evaluation. This paper also gives an explanation of multi-point Kronecker substitution.

We can see the relative performance of the KS2 algorithm with respect to KS1 below.

Thus, we see that KS2 starts performing significantly better for polynomials with larger length (> 80). Thus I propose that the mul_poly implementation will switch between KS1 and KS2 depending on the length of the polynomials. This will have to be carefully benchmarked, as to find a suitable switching point value.

The actual function format will probably remain the same. In addition to this the Fast Fourier Algorithm for polynomial multiplication mode can also be added made as an optional core for multiplying UnivariateIntPolynomial. This a very famous algorithm and many implementations and resources can be easily found for it.

###Division

Currently there is no support for univariate polynomial division in SymEngine. I propose to implement the standard remainder division, which outputs q and r when f is divided by g, such that f = q*g + r and deg(r) < deg(f). The prototype function will be called in this fashion.

(void) div_poly(f, g, outArg(q), outArg(r));

eg. div_poly(x^2 + x + 1, x + 1, q, r);
	q = x
	r = 1

As very succinctly explained by this paper, polynomial division takes order polynomial multiplication time. The paper also fully explains the algorithm, and even provides pseudo code for the same. The division is carried out by defining rev_k on every polynomial,

with rev_k(P) = x^k * P(1/x). When k is the degree of the polynomial, rev_k return the polynomial with reversed coefficients. The detailed algorithm is explained in the paper.

This algorithm assumes g to be a monic polynomial. What if an input g isn't monic? This gives rise to UnivariateRationalPolynomial class. I will discuss this later on in the project proposal.

###GCD & LCM

// Need confirmation on this

// ask @bluescarni about the recent implementation of GCD

The GCD of two polynomials is a polynomial, of the highest possible degree, that is a factor of both the two original polynomials. This paper (from page 41 onwards) explains how the simple euclidean algorithm for finding GCD (as learnt in high school) can be extended to polynomials in some field F. Thus the same algorithm also works Univariate integer polynomials. This algorithm simply put is :

GCD(f, g):
	if deg(f) < deg(g) : swap f, g
	if g = 0           : return f
	else               : return GCD(g, f mod g)

Where mod is calculated using our division algorithm. Since this program is tail recursive, we can also define a iterative function which produces the same result! This algorithm is however O(N^2).

The paper extends the algorithm to a 'Fast Extended Euclidean Algorithm' which not only calculates r = GCD(f, g) but also s, t such that r = s*f + t*g. And this is done via a divide and conquer strategy which reduces the time to almost linear. The complete algorithm has be specified on page 60.

The LCM of two polynomials is a polynomial, of the lowest possible degree, of which both the polynomials are a factor. The LCM can be easily calculated once the GCD function is implemented.

LCM(f, g):
	return (f * g)/GCD(f, g)

The function implementation will be of the form :

(UnivariateRatPoly) gcd_poly(f, g);

eg.	gcd_poly(x^2 - 1, x + 1)
	x + 1
	
(UnivariateRatPoly) lcm_poly(f, g);

eg.	lcm_poly(x^2 - 1, x + 2)
	x^3 + 2x^2 - x - 2

###Derivatives and Indefinite Integrals

Helper functions will be made which will return the derivate or indefinite integral of the polynomial with respect to constituent variable. These functions will be relatively easier to implement.

(UnivariateRatPoly) derivate(f)

eg.	derivative(2y^3 + y + 1)
	6y^2 + 1

(UnivariateRatPoly) iintegral(f)

eg.	iintegral(2y^3 + y + 1)
	(1/2)y^4 + (1/2)y^2 + y

###Multipoint evaluation

This method deals with evaluating polynomials at more than one point at a time, and has wide number of applications in maths and physics. Standard evaluation of each of the given points gives rise to a O(N^2) algorithm. An almost linear time algorithm exists, as indicated by the paper here. Another good explanation is provided here. It uses a divide and conquer algorithm. This will speed up the time to O(MlogN) where M is the order of time of multiplication. The same algorithm can be used for both Univariate and Multivariate polynomials.

(vec_rat) multi_eval_poly(f, {p0, p1, ..})

eg.	f = x^2 - x
	multi_eval_poly(f, {0, 1, 3, 7})
	{0, 0, 6, 42}

###Interpolation

Polynomial interpolation stated simply is : given some points, find a polynomial which goes exactly through these points. The problem of interpolating polynomials is fundamental to many problems for numeric computation and manipulation. Interpolation of a points to produce a polynomial, will most likely have rational coefficients, which pushes the need for a separate class for Rational coefficients even more.

The problem of multipoint evaluation and interpolation are solved very similarly, using divide and conquer. The papers provided in multipoint evaluation also cover polynomial interpolation, and provides succinct explanation.

(UnivariateRatPoly) interpolate_poly({x0, x1, ..}, {y0, y1, ..})

eg.	interpolate_poly({0, 1, -1}, {0, 3, 1})
	2x^2 + x

The function will throw an error, when interpolating such a function isn't possible.

###Factorization

Polynomial factorization is useful in all types of computations especially in a future Solver module. Factorization is done in steps, the first two being Square-Free factorization and Distinct-Degree factorization. This does not completely cover the whole topic of polynomial factorization, as it's not the main aim of my project.

#####Square free factorization

A square-free polynomial is a polynomial that is not a multiple of the square of a non unit factor. In layman's terms a square-free polynomial a polynomial with no repeated roots. A square-free factorization of a polynomial is a factorization into powers of square-free factors.

f = f1 * f2^2 * ... * fn^n

where each fi is square free. Also, each of them are pairwise co-prime. This is a very common polynomial manipulation problem. I will use the fast and popular Yun's algorithm. A small write-up on the algorithm is also given in this paper.

The result is returned as a list in Sympy, where the first element is f1, second element is f2 and so on. But I believe for larger and sparser polynomials, the square free decomposition may also be sparse. So, another alternative to returning a vector of polynomials, is to return a dictionary. The keys in this dictionary are the powers of each polynomial and the values are the polynomials themselves. This will have to be discussed with mentors.

(list/dict) square_free_factorize(f)

eg. square_free_factorize(x^5 - x^4 + 6x^3 - 6x^2 + 9x - 9)
	[1, x - 1, x^2 + 3]
			or
	{1: x - 1, 2: x^2 +3}

Another useful function will be

(bool) is_square_free(f)

which indicates whether a polynomial is square free or not.

#####Distinct degree factorization

This method splits an already square-free polynomial into a product of polynomials whose irreducible factors all have the same degree. This wikipedia page provides some pseudo code, while this lecture discusses the algorithm in detail. The same argument regarding the choice of returned data structure can be made here.

(list/dict) distinct_degree_factorize(f)

#MultiVariate Polynomials

###Overview

SymEngine currently has no class for Multivariate Polynomials. From the discussions on gitter, I found out that the UC Davis team has already started working on the multivariate class. I will work alongside them on this, and contribute to the design of the class. The following are the ideas I had about implementing the class.

I will make this class for either the rational coefficients or integer coefficients depending on how things have shaped up with the univariate counterpart. The class will have the following structure :

  • A dictionary : The keys of the dictionary can be thought of as the tuple formed by the powers of the variables in each term. The values will be the coefficients associated with said term. The representation of tuples will be discussed below.
  • Symbols : The number of variables, and the symbols associated with them. Also, these variables will be indexed, so there is some sort of ordering within them.
  • Constructor methods and other internal methods.

Canonicalization will have a very similar approach as the univariate counterpart. The tuples will be stored in increasing order of sum, ties being decided by which term has a smaller power on the variables in order.

So, 3*x^2*y^3 + 2*x^3*y^2 + x^2*y - x*y^2 + 3 will be stored in the order

3   ->  -x*y^2  ->  x^2*y    ->    3*x^2*y^3   ->    2*x^3*y^2

Also, after taking derivative, evaluating and addition/subtraction of multivariate polynomials, one or more variables may disappear from it, and this too needs to be updated. A helper method can also be made for this purpose, and can be used elsewhere.

###Goal

My main objective is to implement a fast sparse multivariate polynomial class. Additionally, I believe that the only a single multivariate class should exist. All univariate polynomials are in fact multivariate polynomials with only one variable. This is also an ideal design structure for the class.

Thus, when only one variable will be declared, it itself uses the whole integer to store the powers (acting like the dictionary in the univariate case). Similarly, all functions in this class are defined for n number of variables, and should work when n=1. It may be difficult to achieve the speeds of just the univariate class in the multivariate parent class, but the results must be benchmarked, and a trade-off between design and speeds.

If speeds are severly affected by using the multivariate class as a univariate one, I will have to implement separate functions for both the classes, and some functions which combine the two classes together too (discussed later).

###Implementation of monomial tuples

// Could use more inputs on this

// as mentioned by the ideas page, we can substitute the tuple by packing it into a single integer. But that has it's limitations (we have to use the machine ints, for a fixed size)

The trivial way to implement the dictionary is to construct a map from a tuple of integers (corresponding to the respective powers of each variable in the monomial) to another integer, which corresponds to the coefficient of this monomial. This tuple can be implemented as a vector in C++.

For example, 3*x^2*y^4 - x^3 will be represented as {(2, 4): 3, (3, 0): -1}, where the tuples are being stored in a std::vector.

###Addition & Subtraction Because of this sparse approach, addition and subtraction of multivariate polynomials are trivial. We just need to take care of disappearing variables, as mentioned above.

(MultiRatPoly) add_poly(f, g)
(MultiRatPoly) sub_poly(f, g)

eg.	add_poly(x^2*y - 1, x^2*y + x)
	2*x^2*y + x - 1

eg.	sub_poly(x^2 - y, x^3 - y)	// decrease in symbols
	x^2 - x^3

The above methods may also be overloaded for differentiating between univariate and multivariate addition/subtraction.

###Partial derivatives and integrals As in the univariate case, helper functions will be made which will return the derivate or indefinite integral of the polynomial with respect to a given variable. These functions can then be called by the diff/integrate method externally.

(MultivariateRatPoly) derivate(f, wrt_symbol)

eg.	derivative(x^2*y - x*y, 'x')
	2*x*y - y
	
(MultivariateRatPoly) iintegral(f, wrt_symbol)

eg.	iintegral(x^2*y - x*y, 'y')
	(1/2)*x^2*y^2 - (1/2)*x*y^2

###Multiplication

The most common algorithm for multivariate polynomial multiplication is by using Kronecker's trick. On a high level, what it basically does is convert the two multivariate polynomials to their univariate counterpart (by using a suitable power of one of the variables) and then use a fast univariate multiplication technique, to get the resulting univariate polynomial. Now the Kronecker inverse of this new polynomial, gives us back the required polynomials. This paper and these slides cover up on this technique and provide algorithms for the same.

This will require a new method which converts multivariate polynomials to univariate polynomials efficiently (using the given encoding to the variables). Also, I would need evaluation methods ready before multiplication can be done.

(MultivariateRat) mul_poly(f, g)

eg.	mul_poly(x*y, x^2 - y + z)
	x^3*y - x*y^2 + x*y*z

Care must be taken while implementing these methods, to check all the different cases. A univariate polynomial with different variables x and y, should give rise to a multivariate polynomial. Also different multivariate polynomials (with mismatch of variables) getting multiplied should also be handled correctly.

Multivariate division requires Gröbner basis, which I'm not proposing to implement.

###Evaluation

Multivariate polynomial evaluation can be done using a generalization of Horner's method proposed by Pena & Sauer here (page 3). These keeps in memory more than one registers to compute the terms of the polynomial iteratively. The algorithm works like this :

Chooses a variable and apply the Horner's evaluation scheme in that variable. Thereby we will be treating the other variables as constants. Afterwards another variable is chosen and the same process is applied again. This is repeated until all variables are processed. The function returns an integer (or a pair for the rational case).

Another evaluation method should be partial evaluation. Which evaluates the multivariate polynomial at some given values of variables, i.e. not all symbols may have a value. This will return another multivariate/univariate polynomial.

(MultivariateRatPoly) eval(f, dict_values)

eg.	eval(x*y + 1, {x:1, y:2})
    3

	eval(x*y + 1, {x:1})
    y + 1

	eval(x*y + 1, {x:1, z:2})
    ERROR: Symbol not present

#Additional Goals

###Class for Rational coefficients (as a replacement)

According to the current polynomial class structure in SymEngine, polynomials with rational coefficients will be treated under the general class of coefficients UnivariatePolynomial. I propose, that a new class should be implemented for Rational coefficient polynomials. Rational coefficient polynomials are very common in all types of computations, and hence their function implementations must be as fast as integer coefficients polynomials. If they fall under the general class of polynomials (as they do now), they will use the general methods of the UnivariatePolynomial class and will not be able to benefit from the fast methods available to them.

The only difference between the UnivariateInt and the UnivariateRat class will be that the rational class, will store a map from degree to a pair of numbers.

typedef map<unsigned, pair<intclass, intclass> > map_uint_pair_mpz

These numbers represent the numerator and denominator respectively. Maybe the UnivariateInt class, can altogether be removed, (as integers are just rationals with denominator 1). This will have to be decided by benchmarks of how the rational class performs as an exclusive integer class only.

All the above mentioned algorithms, and the already existing ones will work even for rational coefficients (except for the Kronecker substitution trick for multiplication). The Fast Fourier Transform algorithm can be used as a substitute (which works perfectly for rational coefficients). This can also be used as an optional substitute for multiplication in the integer class as well. (if the integer class stays). Also, slight modifications would have to be done to evaluation functions, which will return a pair of integers as a result.

The same idea can be used for the Multivariate class as well. Ideally, we should look forward to implementing only a single class which captures both Univariate and Multivariate polynomials, and both Integer and Rational coefficient. This superclass should be made as fast as the individual classes.

###Seamless operation with the rest of SymEngine

As in Sympy, the polynomial class and other symbolics must seamless operate with one another. Whenever basic operations occur between a Basic and a Polynomial, it must be checked wether the resulting expression can be represented using Polynomial. For eg. add(x, Poly(x**2 + x)) should give Poly(x**2 +2x) instead of creating an instance of the Add class. For this to happen, we need a helper function can_be_poly which will tell us if the other operand can be converted to a polynomial or not. If it can, the result should always be a polynomial for basic operations.

The can_be_poly function can be even given a form of convert_to_poly for constructing polynomial instances out of Basic instances consisting of Mul,Pow and Add. Another functionality that can be implemented is a poly_to_basic which does the reverse of the above mentioned function. It converts polynomial instances to Basic.

This will mainly induce changes in add, sub, mul, div, pow. These are basically all the operations that can be performed on polynomials, which may still result in a polynomial. Only after these seamless operations are implemented, the SymEngine polynomial module can be said to be complete. I will need guidance from mentors on how this integration has to be done, and also have to read up the Sympy code base to see how it's implemented there.

###Implement an alternate dense core for the class Another feature is to implement a dense core for all the classes implemented above. This simply involves an array for Univariate polynomials. Sympy also has both the internal implementations for it's polynomial class. It also chooses which implementation to use, depending on the type of polynomial it is dealing with it.

I will need feedback and ideas from mentors about how this core can be switched from within the class. Also, new functions will have to be defined for dense polynomial representation too. Apart from this, I'll require a function to inter convert the two forms, and also a classifier function which decides which core is better to use for which polynomial.

###SymEngine's own Integer Class

Couple of discussions over the polynomial module, and also the GSoC Idea page indicate the need for SymEngine's own Integer class. The need arises because, machine integers though fast, are limited to a max of 128 bits. While integer libraries like GMP don't have a restriction on size but are not as fast as machine ints. SymEngine's Integer class should ideally use machine ints till the numbers being dealt with are within limit, and switch to other libraries (like GMP) when the numbers exceed the limit, automatically.

Though I do not know how this class will be implemented, I plan to help in it's development and work alongside my mentors to get it finished, and ready to use.

#Timeline

I also plan to do an internship at a local startup here in Mumbai, along with the Summer of Code. I will still be able to devote 38 hours a week to this project. I plan to spend 4 hours per weekday and 18 hours on the weekend on getting the project complete. My summer break starts at the end of April. While I believe my internship starts in the first week of May. I can use this time to get a head start into the tasks at hand too. Also, I can efficiently use the time before the summers to gather all the necessary information and details on my project, so I can begin work without any delays. My internship will end by July second week, so I dedicatedly work from there onwards too.

I will also maintain a blog so that mentors can regulate my work, and also monitor my progress.

####Pre-GSoC

  • Examine the Sympy source code to look at how the above mentioned structures and algorithms are implemented
  • Discuss the choice and implementation of the structures and algorithms with mentors
  • Chalk out all design details of implementation and decide on undecided areas of implementation
  • Get all the pseudo code for the complex algorithms ready before summers starts, so there is no delay in implementing them

####Week 1

  • Implement Fast Fourier Transform for multiplication and the above mentioned KS2 algorithm
  • Introduce the UnivariateRatPolynomial class and model the functions and methods

####Week 2

  • Benchmark the UnivariateRatPolynomial class against UnivariateIntPolynomial on processing only integer coefficient polynomials. This will determine how the class shapes up in the future. From this point on, the functions will be implemented for both the classes (i.e. if Int class still exists)
  • Implement the fast division algorithm for the class

####Week 3

  • Discuss with @certik and @bluescarni about the GCD algorithm, and implement it for the class. The GCD algorithm was recently developed for Pirahna
  • Implement the derivative and integral functions

####Week 4

  • Implement interpolation and multi-point evaluation for the class, also write tests and increase overall coverage of the whole module

####Week 5

  • Start working on the factorization of univariate polynomials. Discuss the final implementation scheme for it
  • Parallely, start working on the optional core for the class, i.e. the dense polynomial representation of the class

####Week 6

  • Wind up all the leftover work for the the whole of the univariate class
  • Finish up basic algorithms for the optional core (dense implementation)

####Week 7

  • Work on integrating the univariate polynomial class with the rest of the SymEngine symbolics
  • Start working on the multivariate polynomial class (or refactor it if already merged in by the UC Davis team)

####Week 8

  • Benchmark MultivariateUniPoly against UnivariateIntPoly extensively, to see if the former can be used as a superclass or not
  • Work on a MultivaraiteRatPoly. Again, the progress will depend on how the rational class fared in the univariate counterpart

####Week 9

  • Work on optimizing the internal structure of the multivariate polynomial
  • Add basic functionality like addition and subtraction
  • Improve test coverage

####Week 10

  • Write function for integrals and derivatives
  • Finish up all required helper functions, so the class is complete as is

####Week 11

  • Write up the remaining functions for the multivariate class including Evaluation, Multiplication and Division

####Week 12

  • Finish off the module, by adding functions to interface multivariate polynomials with the rest of SymEngine
  • If time permits, start writing the optional core (dense) for the multivariate counterpart and also corresponding functions

####Week 13

  • Wrap up the whole module, and finish off any leftover work
  • Benchmark the performance of the module, and compare statistics.
  • Noting down where the shortcomings and scopes of improvement for the module

####Post GSoC

I will always stay a part of SymEngine, and keep contributing for years to come. It still excites me to know SymEngine is possibly the fastest CAS right now, and I would always want it to have as many functionalities as possible. Google Summer of Code has given me a great opportunity to get started with open source contributions, and be a part of the open source community.

I hope to stay part of SymEngine when (hopefully) it moves ahead of Sympy in the CAS competition and continues to be the fastest CAS ever!