Sympythy

Faster arithmetic

2007-08-15T21:10:00.001+01:00

A quick proof-of-concept finite field arithmetic and factorization was done. And the factorization over the finite field was reasonably fast, too, but the big prime reconstruction of the integer factors took some more time. Also, the resulting complete algorithm spent a lot of time with square-free factorization and other preparations. So, it was quite disappointing.

I still started over again, rewriting everything from scratch, now with tests. First the modular integers, with a custom class for each modulus. Then a generic light-weight polynomial class, for one variable, using Python dictionaries. Then some algorithms for polynomials with finite prime field coefficients, with complete factorization, which is also quite fast. Today, I've started to do some more integer polynomial arithmetic, with a gcd using a small primes modular approach. Should be one of the fastest, asymptotically, but probably needs some more tweaking in my implementations. I should now check in which cases these modular methods really are faster than working with rational numbers. All this won't be finished for the Summer of Code deadline, but can be further advanced, while I learn from my new bought book, "Modern Computer Algebra" (Gathen, Gerhard).

Finite field factorization

2007-08-06T08:51:00.000+01:00

I've finished updating the documentation and have started to move to next point, faster univariate factorization. This is done modular, with factorization over a finite field first. Fortunately, this point is a major part of Garthen and Gerhard's "Modern Computer Algebra" and thus well explained. Unfortunately, I've also read there:

Computaionally, (fast) polynomial factorization over a finite field is a much more advanced task than, for example, multiplication or even gcd computation. Before implementing a particular algorithm, one has to implement carefully many other routines for basic polynomial arithmetic, and designing a polynomial factorization package usually requires several woman-years.

I still hope to get something done this week, doesn't have to be optimal, then. Let's see if it is a lot if work to integrate modular (number) arithmetic, or if this is better done "inline".

Refactoring

2007-08-01T15:06:00.000+01:00

Now, that the planned features are roughly present and we have moved to a new core with SymPy, there were some suggestions concerning the structure of the polynomials class. Most of them about the wrapper functions and code duplication, or the internals of the Polynomials class.
The Polynomials class is now immutable, as are all Basic-inherited classes, and tries hard to pretend to be a Basic object itself. We'll see how that will turn out. The wrapper classes essentially only extract the underlying sympy expression of the resulting Polynomial.

The plan for the next days includes:
Documentation (proper docstrings everywhere, algorithms explained)
Testing (more complete tests, with comments)

Then there are some more important features, like fast univariate integer factorization with finite fields. We'll have to look at some number theoretic python package, probably nzmath. Also, the equation system solver could tell us, when we are missing some solutions, that can't be computed but still exist.

Finally, multivariate factorization

2007-07-17T01:04:00.000+01:00

Ok, I gave up with Wang's algorithm, for now. It is not that hard, in general, but the article describing it left quite some questions open for me. So instead of procrastinating the implementation any further, I did another algorithm, just today. It is a lot simpler, also due to Kronecker (20th century!), and not too efficient, of course.

The last feature I want to add now, is solving polynomial equations. It is already possible now to compute a Gröbner base in a convenient order, like lexicographic, to get eliminiation working. Then, for each variable, you have to solve the univariate polynomial in that variable with solve(). You could insert each of these roots in the rest of the equation and proceed with the next variable. The problem with that is, that this approach only works for equations with a finite number of solution, which has to be decided first.

Then there is still some time left for polishing: Documentation and efficiency allow for improvement!

Rational factorization.

2007-07-04T01:04:00.000+01:00

Just committed the new factorization algorithm for polynomials in one variable with rational coefficients. In now factorizes completely, doesn't just look for roots. In theory, factorization of polynomials over the rationals and factorization over the integers are equivalent (by Gauss's lemma). This means, that you can get rational factorization by getting rid of the denominators, then the content (gcd of coefficients), factorize over integers and divide the denominator again. Nothing is lost or gained this way. But in the implementation this situation got me a lot of problems. I always had to assume that the coefficients are integer, checking would occur often and so is expensive. On the other hand, some algorithms, like gcd sometimes introduce non-integer rationals, that need to be corrected afterwards. It seems to be fixed now (all assertions run) but it doesn't look elegant, something simpler should do.

I hoped that the Kronecker algorithm would be fast enough, so that we wouldn't need finite fields for that, but it looks like this is impossible, since the complexity is exponential in (the number of divisors of each of) the coefficients of the polynomials, working fine for simple examples, but probably no real-life problems.

But the next direct step should be multivariate factorization, as a generalization of the univariate case, to complete more of my initial proposal, before worrying too much about efficiency.

Refactoring, Sturm, Sqare-free

2007-06-26T20:09:00.000+01:00

So, most of the refactoring is done, there are now wrapper functions for the common algorithms, that create Polynomial instances and call specialized functions in sub-modules. It's not perfect, the code-reusing and efficiency could be better, but I'm relatively happy with it.

I also implemented Sturm sequences. These can be used for univariate polynomials, to count the real roots in a given interval. This could be extended to isolate the (existing) roots in rational intervals, to give some representation of algebraic numbers, together with an irreducible polynomial.

For that, you need a square-free polynomial (multiple roots are likely to be ignored by the Sturm sequence), so I looked at (and fixed) the square-free decomposition algorithm.

Next things to do include multivariate square-free decomposition, fast factorization with one variable over the integers (Zassenhaus or Berlekamp), some multivariate factorization (found a paper by Wang) and some faster multivariate gcd (using sub-resultants?). All this is needed, before I can attack elimination and finally equation systems.

When looking for algorithms today, I've noticed that most books only treat the univariate cases, but the articles are not accessible for free/everybody. Where do you usually get the ideas/help?

Back to work - New polynomial representation.

2007-06-13T21:53:00.000+01:00

Sorry for the delay, folks, I am finally back in Heidelberg, Germany and have also got over last week's procrastination problems.

I was about to refactor the function code, to write wrapper functions, looking for variables and coefficient type and then choosing the best specialized algorithm to do the job. But writing the functions, e.g. the division algorithm, I remarked the cumbersome way of using (and keeping in sync) two different representations, the sympy object and the coefficient list. Since the latter is necessary for the leading term computation, I thought it was best to stick with it for all arithmetic like addition and multiplication. This way, it should also be way faster, but we'll see about that.

The way it is implemented does not satisfy me at all, the code that handles variables and order before adding two such objects is ugly and (perhaps) contra-intuitive. Also, much more code-reusing could be made. But I hope, that the use of this class will look cleaner in the basic algorithms.

class `Ideal', next steps

2007-05-30T22:51:00.000+01:00

So, the Summer of Code officially started this week. Unfortunately, I'm spending my last days in Paris right now and will have to spend some time on moving back to Germany this/next week.

On the other hand, some functionality (especially the groebner() function) are already available, even if in some slow, unpolished way. I've also added a new class `Ideal' that represents polynomial ideals and supports simple compositions, like I+J, I*J, I**3 and I.intersect(J), can test for membership, as in f in I --> True, or equality, I == J (equality in terms of the set, not the attributes like variables or order). I hoped to add the functionality of computing the radical of an ideal (the multivariate, multipolynomial equivalent of square-free part) but this turns out to be rather difficult, and not of any practical use outside of (theoretic) algebraic geometry.

The groebner() function also has a nice application in computing the lcm and gcd of two multivariate polynomials with real (and even symbolic) coefficients, which could be used for rational simplifications, if not for speed issues. I know that there are faster ways to compute the gcd and also the whole factorization of polynomials, using modular algorithms, that is, factorizing mod some prime number p, that trying to `lift' the result. Unfortunately these algorithms can only be applied to polynomials in integer (and thus rational) coefficients, but could never handle some factor sin(z) for an unrelated z, as the Gröbner bases can. Nevertheless, these modular algorithms are on the to-do list.

Also on this list is a function generating the Sturm sequence (I'll write about that later) that can count real roots of a polynomial and than perhaps approximate them with rational intervals. Then there is room for those standard formulas to get the real and complex roots of polynomials of degree < style="font-weight: bold;">not be such a formula for higher degree)

These will be the requirements for the whole elimination process, which could be somewhat limited with approximated non-rational roots, but we'll see.

2007-05-24T12:19:00.000+01:00

Nice literature

I forgot to mention some interesting books I use.

The most important one is Cox, Little, O'Shea: Ideals, Varieties, and Algorithms, which is an introduction to algebraic geometry (and commutative algebra) using a relatively constructive approach. That is, for most of the objects, there are not only proofs of existence (and uniqueness) but also algorithms (or at least references). I like it very much and have learned a lot so far; it was already a plan for the summer to revise some algebraic geometry, this fits very well.

Yesterday, I found two other books in the university, more directly concerned with (general) computer algebra:

Davenport, Siret, Tournier: Computer Algebra is a gentle introduction to different features of CAS, problems of data representation and some advanced algorithms. It could be a bit outdated on current implementation details and running times, but gives a quick overview and offers confrontation with problems otherwise ignored. (Already half-way through, not looking at all details)

Gathen, Gerhard: Modern Computer Algebra. Have not looked at it in detail, but looks really nice. One of the authors is an employee of MapleSoft and all methods are presented with a huge list of applications and examples. The content is presented with some information on historic development and stuffed with nice quotations (also in the original language, way to go!)

2007-05-23T21:29:00.000+01:00

Gröbner basis

Hey, my first blog entry so far.

My name is Robert Schwarz, I am a student of mathematics and computer science at the university of Heidelberg, Germany. For the next two weeks, I'll still be staying in Paris, where I passed a year with the Erasmus student exchange program.

The GSoC project is mostly about Gröbner bases and application. The euclidean algorithm for finding the gcd of 2 integers is one of the most popular examples in programming. The same thing also works for polynomials, but unfortunately only in one variable. The Gröbner basis is a set of polynomials that means to replace the role of the gcd.

When looking for the solutions of some system of polynomial equations f(x,y)=0, g(x,y)=0 it is clear that we could add more equations like f+g or f*h for some polynomial h without affecting the set of solutions. The set of such equations/polynomials is called the ideal I = , generated by f and g. If I want to know if some polynomial h(x,y) lies in I, I have to check if it can be written as h = a*f + b*g, that is, I try to divide h by f and g and look at the remainder, similar to the division with integers. Unfortunately, the choice of f,g and even their order can influence the remainder here, and it is the special characterization of the Gröbner bases, that it works with them. Hope, it is a bit clearer for those who have never heard of them.

The natural application lies of course in the solving of polynomial equations, with simplification of the equations and systematic elimination of variables. I think I've read somewhere that Gröbner bases are also needed for some other feature in sympy, I'm very interested in that, please tell me, so that I can start to help.

Yesterday, I finished the implementation of Buchberger's algorithm that computes Gröbner bases, so in theory, it is already available, even if it's still a bit slow.

Items on my to-do list: Speed improvements in groebner(), real-life examples for polynomial equations, elimination, factorization of polynomials, solving for (real) roots.