WIP: Incremental progress on sparse polynomials

[mattpap]

The name of this branch (new-polys) may be a bit confusing, because it doesn't go that far, but there is enough work done to merge it soon. In this branch:

• improved support for real and complex domains (previously Improved support for reals in polys #1921)
• switched PolynomialRing and FractionField to use PolyElement and FracElement instead of DMP and DMF
• use code generation to improve speed of monomial ops (instead of looping over items of monomials, unpack monomial tuples and execute every op explicitly, see MonomialOps)
• other optimizations to PolyElement, FracElement and PythonRational (some significant)
• added caching of PolyRing and FracField (faster == or just use is)
• added @public decorator (very rough and limited implementation)
• started cleanup of polys' modules (for now split monomialtools.py into monomials.py and orderings.py)
• create custom PolyElement per PolyRing to simplify isinstance(g, PolyElement) and ring == g.ring to isinstance(g, ring.dtype) (among other benefits) (the same for FracElement/FracField)

What remains to do is to:

• fix two test failures in test_risch.py (test_residue_reduce, test_integrate_nonlinear_no_specials)
• fix test failure in test_solvers.py (test_exclude)
• fix doctest failure in solvers.py (solve(x**2 - y**2/exp(x), y, x))
• fix very weird doctest import failure in Py3 (bin/doctest sympy/polys crashes)
• remove Cython support
• uncomment a test in test_risch.py (test_hermite_reduce)
• investigate Piecewise integral

There is one downside of this PR. I couldn't, yet, remove old PolynomialRing and FractionField (those using DMP and DMF), because I didn't want to dig into agca module. I will rewrite agca at some point, unless someone would like to step up (@ness01?).

[asmeurer]

Regarding caching in general, also see https://code.google.com/p/sympy/issues/detail?id=3222. I think we should use LRU for the cache. It's dead simple to implement, and performs quite well.

[asmeurer]

https://tech.dropbox.com/2012/10/caching-in-theory-and-practice/

[mattpap]

I have a few questions about the new caching

This cache is very small, about 2500 entries per bin/test. I only cache domains (polynomial rings and fraction fields). No polynomials are being cached. Cache control can be added later on. I don't see it very useful at this point (in this context).

[mattpap]

History of support for this integral is interesting. It used to work in 0.7.2. Then it was broken in cd42ae8 (bisected to that commit). Then 8d7d239 fixed the exception, but lead to symbolic, yet correct, result. The following commit allowed to return numerical result, but due to earlier commits the result is 0.

[asmeurer]

So maybe try bisecting while doing a git cherry-pick --no-commit 8d7d239 at each stage to see what broke the numeric behavior.

[certik]

Cache: ok, I would leave it as it is then.

[mattpap]

So here is the story. Last time ratint(1/(0.000871222*t + 0.995)**2, t) used to work was in 0.6.7 (didn't check if this wasn't by coincidence). At that time there was no support for piecewise functions in the integrator. This was added in 0.7.0, but at that time ratint() was already broken. However, the integral in consideration did work. This is because ratint() wasn't used, just simple pattern matching. It worked this way until integration3 was merged. It must have significantly change heuristics in integrate(). Since then ratint() was used instead of pattern matching, but ratint() has been broken for a long time. new-polys partially fixed the problem. If you want a quick fix of this integral, I suggest investigating heuristics in integrate() (those changed by integration3). I can safely say that this integral didn't ever work the way we want it to work now. It worked but in a different setup. I don't expect that the situation with inexact numbers will be fixed in this branch (it's a too deep issue). Sill, I will disable automatic truncation (one more commit --- need to fix two simple test failures), at least for now, until I have time to look more closely on this issue.

[asmeurer]

I don't care if it works or not, just as long as we don't get any wrong results.

[asmeurer]

The difference from that risch commit is this: The integrand is of the type handleable by risch_integrate, so it passes through that algorithm. But the way the recursive algorithm works, when it is given a rational function, it just passes it off to ratint. Now, actually, it doesn't call ratint directly, because it doesn't work on polynomials. It calls integrate. But the issue is that the pattern matching looks for (a + b*x)**c, but risch canonicalizes the integrand (having first converted it to a Poly), so it is expanded. Hence the match doesn't work any more.

Unfortunately, calling factor wouldn't be an easy fix, because that doesn't seem to work

In [2]: print factor(7.59027773284e-7*t**2 + 0.00173373178*t + 0.990025)
1.0*(7.59027773284e-7*t**2 + 0.00173373178*t + 0.990025)

But long story short, take a look at the history of integrate(expand(1/(0.000871222*t + 0.995)**2), t). That has given CoersionFailed in master and 0.7.2, but 0 in this branch. That started giving 0 in df33365, but, quite interestingly, in 8d7d239, it gave

In [1]: integrate(expand(1/(0.000871222*t + 0.995)**2), t)
Out[1]:
                                                                                     ⎛      __________________________________________________________
                       __________________________________________________________    ⎜    ╲╱ 45773156866350624606286237784751595254442599505912800645    801
8.62950215186415e-18⋅╲╱ 45773156866350624606286237784751595254442599505912800645 ⋅log⎜t - ──────────────────────────────────────────────────────────── + ───
                                                                                     ⎝                 599631266385525559905071765212450                   7

                    ⎞                                                                                        ⎛      ________________________________________
26175720734332347747⎟                          __________________________________________________________    ⎜    ╲╱ 457731568663506246062862377847515952544
────────────────────⎟ - 8.62950215186415e-18⋅╲╱ 45773156866350624606286237784751595254442599505912800645 ⋅log⎜t + ──────────────────────────────────────────
0158479461074971170 ⎠                                                                                        ⎝                 59963126638552555990507176521

__________________                          ⎞
42599505912800645    80126175720734332347747⎟
────────────────── + ───────────────────────⎟
2450                   70158479461074971170 ⎠

I would rather have CoersionFailed than 0. Can we revert whatever part of df33365 is incorrect (I guess the tolerance stuff, based on your previous answer)?

I guess much of this is just reproducing what you already said.

[mattpap]

I made division algorithms not fail silently or hang when they can't reduce to zero. Now PolynomialDivisionError is raised with a detailed explanation. Disabled automatic reduction to zero (at least until I investigate this issue more). Added an XFAIL-ed test for the integral. ratint() now raises AssertionError when resultant is zero. During another round of rewrites of sympy.polys I will try to figure out how to handle inexact numbers properly.

[mattpap]

Work on this failing integral let me to a conclusion that we should have a better XFAIL, such that I could say how exactly it fails, so when behavior changes (still fails) we can catch that. Also (based on earlier work) I think there should be a way to measure (and report) how correct are test results. It often happens that a function returns a correct result but not in the form we want it. There is (usually) only one, currently passing, version. If things start to work the way they should, instead of some sort of XPASS, we get a test failure.

[asmeurer]

https://github.com/asmeurer/sympy/compare/xfail

[asmeurer]

Ok, here goes...

[asmeurer]

This branch removed the minpoly shortcut to minimial_polynomial.

[certik]

Awesome, thanks Mateusz!

[asmeurer]

SymPy Bot Summary: 🔴 Failed after merging mattpap/new-polys (74ecb37) into master (0ef7411).
@mattpap: Please fix the test failures.
🔴 Python 2.5.6-final-0: fail
✅ Python 2.6.8-final-0: pass
✅ Python 2.7.3-final-0: pass
🔴 PyPy 1.8.0-final-0; 2.7.2-final-42: fail
✅ Python 3.2.3-final-0: pass
✅ Sphinx 1.1.3: pass
Doc Coverage: unchanged
36.065% of functions have doctests (compared to 36.065% in master)
41.570% of functions are imported into Sphinx (compared to 41.570% in master)