Faster arithmetic

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

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".



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

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.

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

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.

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.