New polynomials manipulation module

[mattpap]

So, after 5 weeks of coding new module is done. Well, almost done because
many tests are missing and some functionalities may still be quite buggy.
However, after minor fixes, SymPy can work with new module very well, so
lets dig into new code and share your opinions.

Sixteen patches are available from: http://github.com/mattpap/sympy-polys/commits/polys Please read commit message of 7f4e5466414302c3106b2ac9e2224eac2f4d75d2 to
get brief idea what was done. Some parts of new code I already showed at
EuroSciPy conference but now all this is really working. Feel free to
experiment.

Is this the end? ... no, sorry, it smells like the beginning of work on
really functional module. Plans for near future? Besides testing and
documentation writing, proper support for algebraic number fields
(isomorphism, inclusion, splitting fields and factorization).

Is this the right approach to implementing polynomials in SymPy? Don't know
but it seems to work, it's fast and does not require too much algebra
knowledge to work with (at least as a user). With some minor fixes on new
polynomials module side (or Cython side -> nested functions) it could be
Cythonized to obtain even more computational power.

Original issue for #4697: http://code.google.com/p/sympy/issues/detail?id=1598 Original author: https://code.google.com/u/101069955704897915480/ Referenced issues: #4746, #3298, #4686, #4568, #3416, #4442, #4448, #3425, #4194, #4196, #4454, #3687, #3882, #4205, #3694, #4531, #4342, #4667, #4221 Original owner: https://code.google.com/u/101069955704897915480/

[asmeurer]

Wow!  This is great.  Simplification is much better with this.  For example, consider the complex expression 
from issue 4661 .  In your branch:
>>> a = wronskian([sin(x), cos(x), sin(x)*x, cos(x)*x, 1], x)
>>> print a
8*cos(x)**2*sin(x)**2 + 4*cos(x)**4 + 4*sin(x)**4
>>> a = factor(a)
>>> print a
4*(cos(x)**2 + sin(x)**2)**2
>>> trigsimp(a)
4

The original wronskian is simpler (I am assuming because of better intermediate simplification), and the 
expression can be factored!  I tried running a more complex expression through simplify, but it hangs.

>>> simplify(-576*cos(x)**26*sin(x)**30/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
6336*cos(x)**24*sin(x)**32/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
31680*cos(x)**22*sin(x)**34/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
95040*cos(x)**20*sin(x)**36/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
190080*cos(x)**18*sin(x)**38/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
266112*cos(x)**16*sin(x)**40/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
266112*cos(x)**14*sin(x)**42/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
190080*cos(x)**12*sin(x)**44/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
95040*cos(x)**10*sin(x)**46/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
31680*cos(x)**8*sin(x)**48/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
6336*cos(x)**6*sin(x)**50/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44) - 
576*cos(x)**4*sin(x)**52/(32*cos(x)**18*sin(x)**30 + 224*cos(x)**16*sin(x)**32 + 
672*cos(x)**14*sin(x)**34 + 1120*cos(x)**12*sin(x)**36 + 1120*cos(x)**10*sin(x)**38 + 
672*cos(x)**8*sin(x)**40 + 224*cos(x)**6*sin(x)**42 + 32*cos(x)**4*sin(x)**44))

It takes 218.32 s to complete.  If I substitute sin(x) and cos(x) for symbols first, that drops to 88.54 s.  On the 
other hand, cancel(together(expr)) takes 0.30 s to produce the same result.

This is from wronskian([sin(x)**2, cos(x)**2, sin(x)*cos(x), sin(x), cos(x)], x) in the current SymPy.  Fortunately, 
if you run the wronskian straight in your branch, it comes out nice:
In [1]: wronskian([sin(x)**2, cos(x)**2, sin(x)*cos(x), sin(x), cos(x)], x)
Out[1]: -72*cos(x)**6*sin(x)**2 - 108*cos(x)**4*sin(x)**4 - 72*cos(x)**2*sin(x)**6 - 18*cos(x)**8 - 
18*sin(x)**8

In [2]: trigsimp(factor(wronskian([sin(x)**2, cos(x)**2, sin(x)*cos(x), sin(x), cos(x)], x)))
Out[2]: -18

I posted some comments on the code in your branch.  I like the philosophy that you outline in your commit 
7f4e5466414302c3106b2ac9e2224eac2f4d75d2.  I don't know much about abstract algebra (yet), so I can't 
say much about the internal part, though the code looks clean.

Referenced issues: #4661 Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c1 Original author: https://code.google.com/u/asmeurer@gmail.com/

[certik]

Wow, Mateusz, you are a genius. 

I tried this:

In [1]: a = (x+y+z)**20

In [2]: b = a.expand()

In [6]: time factor(b)
CPU times: user 3.08 s, sys: 0.00 s, total: 3.08 s
Wall time: 3.11 s
Out[7]: 
           20
(x + y + z)  

And try maxima, it can't do it! :) So we are slowly overtaking maxima, that's exciting.

This fails for me, but it's minor:

$ bin/test sympy/utilities/
============================= test process starts ==============================
executable:   /usr/bin/python  (2.6.2-final-0)

sympy/utilities/tests/test_code_quality.py[2] ..                            [OK]
sympy/utilities/tests/test_codegen.py[10] ..........                        [OK]
sympy/utilities/tests/test_decorator.py[1] .                                [OK]
sympy/utilities/tests/test_iterables.py[5] .....                            [OK]
sympy/utilities/tests/test_lambdify.py[25] ...........f.............        [OK]
sympy/utilities/tests/test_pickling.py[28] ..........fffffff...........     [OK]
sympy/utilities/tests/test_pytest.py[1] .                                   [OK]
sympy/utilities/tests/test_source.py[2] ..                                  [OK]
sympy/utilities/tests/test_tests_names.py[1]
['/home/ondrej/repos/sympy/sympy/mpmath/tests/test_rootfinding.py',
'/home/ondrej/repos/sympy/sympy/polys/tests/test_rootfinding.py']
F                            [FAIL]

________________________________________________________________________________
______ sympy/utilities/tests/test_tests_names.py:test_no_duplicate_names _______
  File "/home/ondrej/repos/sympy/sympy/utilities/tests/test_tests_names.py", line 21,
in test_no_duplicate_names
    assert False, message % (fname)
AssertionError: Test 'test_rootfinding.py' has a duplicate name, py.test will not run it!

======= tests finished: 66 passed, 1 failed, 8 xfailed in 20.82 seconds ========
DO *NOT* COMMIT!

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c2 Original author: https://code.google.com/u/104039945248245758823/

[certik]

(the test tests a bug in py.test, that it fails to execute two files in the sympy
tree with the same filenames, e.g. one in mpmath and one in polys)

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c3 Original author: https://code.google.com/u/104039945248245758823/

[certik]

Mateusz, would you manage to use cython to speed up the factorization a bit? Just for
fun, so that I can show it on scipy09, as an example of using cython. Here is the
result in Sage, that uses pari:

sage: var("x y z")
(x, y, z)
sage: a = (x+y+z)**20
sage: b = a.expand()
sage: time factor(b)
CPU times: user 0.19 s, sys: 0.00 s, total: 0.19 s
Wall time: 0.30 s
(x + y + z)^20

e.g. it's only 10x faster. Pari is optimized in C, so --- I don't know, but cython
can speed up things easily by a factor of 10x in many cases. Let's try it?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c4 Original author: https://code.google.com/u/104039945248245758823/

[certik]

And for this:

In [5]: a = (x+y+sin(z))**20

In [6]: b = a.expand()

In [7]: time factor(b)
CPU times: user 3.38 s, sys: 0.00 s, total: 3.38 s
Wall time: 3.46 s
Out[8]: 
                20
(x + y + sin(z))  

Sage is only 4x faster:

sage: a = (x+y+sin(z))**20
sage: b = a.expand()      
sage: time factor(b)      
CPU times: user 0.13 s, sys: 0.01 s, total: 0.14 s
Wall time: 0.80 s
(x + y + sin(z))^20

Could you beat it? That'd be fun. :)

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c5 Original author: https://code.google.com/u/104039945248245758823/

[smichr]

SO NICE to see the changes...this is going to make some of my other work much easier.
Thanks!

I noticed that expressions with NumberSymbols can be made into polys, but the new
module doesn't like reals...can those just be changed internally to constants (like
cse does)?

###
>>> Poly(pi-x)
Poly(pi - x, pi, x, domain='ZZ')
>>> Poly(.1-x)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 43
5, in __new__
    result = _init_poly_from_basic(rep, *gens, **args)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 39
4, in _init_poly_from_basic
    domain, coeffs = _find_out_domain(rep, **args)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 12
5, in _find_out_domain
    raise NotImplementedError('inexact coefficients')
NotImplementedError: inexact coefficients
>>>
###

/c

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c6 Original author: https://code.google.com/u/117933771799683895267/

[mattpap]

Thanks guys for comments. So ...

Ondrej:

What's really slow in your examples is expand(). See the following phenomenon:

In [1]: a = (x+y+sin(z))**20

In [2]: %time b = a.expand()
CPU times: user 5.57 s, sys: 0.00 s, total: 5.57 s
Wall time: 5.68 s

In [4]: %time b = b.expand()
CPU times: user 7.77 s, sys: 0.00 s, total: 7.77 s
Wall time: 7.90 s

Expanding an expanded expression is much slower than the original one. This is slow:

In [6]: %time c = factor(b)
CPU times: user 12.94 s, sys: 0.00 s, total: 12.95 s
Wall time: 13.72 s

However you can turn off expanding in Poly or any function in new module by setting
`expand=False`, provided you know the input expression is already expanded:

In [8]: %time c = factor(b, expand=False)
CPU times: user 5.24 s, sys: 0.00 s, total: 5.24 s
Wall time: 5.37 s

This way you really test efficiency of the factoring algorithm.

btw. Creating a Poly instance without expanding is also fast:

In [12]: %time p = Poly(b, expand=False)
CPU times: user 0.01 s, sys: 0.00 s, total: 0.01 s
Wall time: 0.01 s

> AssertionError: Test 'test_rootfinding.py' has a duplicate name, (...)

OK, so this is py.test bug. I can rename rootfinding to polyroots (or similar).

> would you manage to use cython to speed up the factorization a bit?

I can, I even started to write pxd declarations for dense*.py but didn't yet run it
because of lack of nested defs in Cython. However I can rewrite nested functions, to
make Cython happy and experiment, but I don't guarantee that this can be done in one
day. Another thing is that I tried your Cython usage in sympy example, but I can't
reproduce this great speedup you got. I get only like 30% not an order of magnitude.
I will write more about this on the mailing list in the appropriate thread.

> Could you beat it? That'd be fun. :)

Oh, sure I would be ;) What can be improved? Cython and pure mode is one thing but
also factoring algorithm isn't finished yet. Current implementation is focused on
correctness not speed and one big branch of Wang's algorithm isn't implemented at
all, so improvements are possible. I'm also curious how switching to recursive sparse
polynomial representation would affect performance.

asmeurer:

> Should this rather be "Returns self as if it were an Add instance."

Right ;)

> How is this routine different than radsimp?

This code was a part of Poly.cancel, which was copied from radsimp(). However,
radsimp() tries to do more and is problematic (there are issues opened for this).

> I tried running a more complex expression through simplify, but it hangs.

Yes, but this is once again expand() (or really core problem). Now I see that I
should do some hacking and replace expand() (at least partially) with an algorithm
based only on polynomials, i.e. replace p**n with Poly(p)**n if p is expanded,
replace p_1*p_2*...*p_n with Poly(p_1)*Poly(p_2)*...*Poly(p_n) if p_1,...,p_n are
expanded etc. This should be a lot faster. I will also port fast low-level expanding
algorithm to polynomials module to replace the current powering algorithm. All this
is doable because Basic.as_numer_denom() is fast so I don't have to care about
negative exponents.

smichr:

> but the new module doesn't like reals

I does not and rising this exception is intentional. The problem is that some poly
algorithms don't care about coefficients, mostly arithmetics, some can work with
reals but need special treatment, e.g. div, gcd, and some just won't work, e.g.
factorization. Consider all this a work in progress. Soon I will provide an algebra
for reals and wrap up mpmath to do computations. The idea is to specify tolerance of
computation in this algebra so that zero equivalence could be solved and which would
allow division algorithm to work (this is similar to Mathematica Tolerance argument
to polynomial manipulation functions). This way computations will be fast and safe.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c7 Original author: https://code.google.com/u/101069955704897915480/

[mattpap]

I forgot to mention that you can easily switch between different ground types, just
specify SYMPY_GROUND_TYPES=something before isympy. So lets factor once again, this
time using gmpy's mpz and mpq types for ZZ and QQ algebras respectively:

$ SYMPY_GROUND_TYPES=gmpy isympy 
Python 2.6.2 console for SymPy 0.7.0-git (cache: off)

In [1]: a = (x+y+sin(z))**20

In [2]: b = a.expand()

In [3]: %time c = factor(b, expand=False)
CPU times: user 0.87 s, sys: 0.00 s, total: 0.87 s
Wall time: 0.95 s

and lets compare with standard Python's types:

$ SYMPY_GROUND_TYPES=python isympy 
Python 2.6.2 console for SymPy 0.7.0-git (cache: off)

In [1]: a = (x+y+sin(z))**20

In [2]: b = a.expand()

In [3]: %time c = factor(b, expand=False)
CPU times: user 5.25 s, sys: 0.01 s, total: 5.26 s
Wall time: 5.37 s

If Fraction is not available SymPy can run in a mixed mode, i.e. you can work with
int as ZZ and mpq as QQ. Note, however, this mode is not yet tested well.

btw. I thought, Ondrej, that the conference is this weekend but it starts next
Thursday so I can surely experiment with Cython and try to use it in polynomials
module before then.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c8 Original author: https://code.google.com/u/101069955704897915480/

[certik]

ok, first expand=False only speeds things up a bit for me (from 3.16s to 2.84s),
however, this:

$ SYMPY_GROUND_TYPES=gmpy bin/isympy
Python 2.6.2 console for SymPy 0.7.0-git

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z = symbols('xyz')
>>> k, m, n = symbols('kmn', integer=True)
>>> f, g, h = map(Function, 'fgh')

Documentation can be found at http://sympy.org/ In [1]: a = (x+y+sin(z))**20

In [2]: b = a.expand()

In [3]: %time c = factor(b, expand=False)
CPU times: user 0.44 s, sys: 0.00 s, total: 0.44 s
Wall time: 0.45 s

makes it faster than sage!
sage: var("x y z")        
(x, y, z)
sage: a = (x+y+sin(z))**20
sage: b = a.expand()      
sage: %time c = factor(b) 
CPU times: user 0.15 s, sys: 0.01 s, total: 0.16 s
Wall time: 0.59 s

but that's because sage is using maxima if there is sin(x) in the expression. So
that's really awesome.

So let's rather compare this:

$ SYMPY_GROUND_TYPES=gmpy bin/isympy
Python 2.6.2 console for SymPy 0.7.0-git

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z = symbols('xyz')
>>> k, m, n = symbols('kmn', integer=True)
>>> f, g, h = map(Function, 'fgh')

Documentation can be found at http://sympy.org/ In [1]: a = (x+y+z)**20

In [2]: b = a.expand()

In [3]: %time c = factor(b, expand=False)
CPU times: user 0.40 s, sys: 0.00 s, total: 0.40 s
Wall time: 0.40 s

to this:

sage: var("x y z")       
(x, y, z)
sage: a = (x+y+z)**20    
sage: b = a.expand()     
sage: %time c = factor(b)
CPU times: user 0.14 s, sys: 0.00 s, total: 0.15 s
Wall time: 0.15 s

E.g. Sage (pari) is only about 2.6x faster. My lecture is on Thursday at 4pm, so if
you can make it 3x faster with gmpy + Cython, then we can beat Sage.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c9 Original author: https://code.google.com/u/104039945248245758823/

[certik]

Yeah, but once we beat Sage, it's still nothing, since Mathematica can do this instantly:

In[1]:= c = Factor[Expand[(x+y+z)^100]]

                   100
Out[1]= (x + y + z)

so, the journey will be long.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c10 Original author: https://code.google.com/u/104039945248245758823/

[asmeurer]

> This code was a part of Poly.cancel, which was copied from radsimp(). However,
> radsimp() tries to do more and is problematic (there are issues opened for this
OK.  I didn't know that was from Poly.cancel().  Looking at the radsimp source, it looks like the only difference 
between your code and radsimp is that radsimp tries to handle the case a + b*sqrt(c) where b, c == 0 and it tries 
to collect terms.  Maybe you could add a keyword argument to radsimp that turns this stuff off (like radsimp(expr, 
simple=False)), and then use that argument in the call in simplify().  That way, you are not duplicating code.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c11 Original author: https://code.google.com/u/asmeurer@gmail.com/

[asmeurer]

It looks like factor also chokes on reals, and you have a bug:
>>> factor(x*(y + z)**(1/2)*(1 + y))
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "./sympy/polys/polytools.py", line 1839, in factor
    F = Poly(f, *gens, **args)
  File "./sympy/polys/polytools.py", line 435, in __new__
    result = _init_poly_from_basic(rep, *gens, **args)
  File "./sympy/polys/polytools.py", line 367, in _init_poly_from_basic
    rep, gens = dict_from_basic(ex, **args)
  File "./sympy/polys/polyutils.py", line 168, in dict_from_basic
    return _dict_from_basic_no_gens(ex, **args)
  File "./sympy/polys/polyutils.py", line 123, in _dict_from_basic_no_gens
    base, exp = _analyze_power(*factor.as_Pow())
  File "./sympy/polys/polyutils.py", line 65, in _analyze_power
    raise PolynomialError("%s is not a valid polynomial term" % term)
NameError: global name 'term' is not defined

It would be nice if factor just returned the expression if it can't do it, instead of raising PolynomialError.  That 
way, you don't have to have a try block every time you want to factor something.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c12 Original author: https://code.google.com/u/asmeurer@gmail.com/

[vks]

Chris, could you please try to rebase this branch on master? There are some
non-trivial conflicts with your improvements, I think it's easier for you to resolve
them than for me.

**Cc:** smichr  

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c13 Original author: https://code.google.com/u/Vinzent.Steinberg@gmail.com/

[smichr]

Do you mean the polys branch?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c14 Original author: https://code.google.com/u/117933771799683895267/

[vks]

Yeah, please try to rebase polys on master, I think you know the code better.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c15 Original author: https://code.google.com/u/Vinzent.Steinberg@gmail.com/

[smichr]

OK, give me a day or two...

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c16 Original author: https://code.google.com/u/117933771799683895267/

[vks]

Great, looking forward!

Maybe we should release 0.6.6 before without new polynomials?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c17 Original author: https://code.google.com/u/Vinzent.Steinberg@gmail.com/

[smichr]

ok, there are some issues to attend to as the test will show and there are a couple
of things discussed in the commit. I think I've got it rebased properly. See the
polys2 branch at smichr's github.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c18 Original author: https://code.google.com/u/117933771799683895267/

[smichr]

Was I right in deleting trim and cancel from simplify, Matteusz?

Also, factor still should be doing something different. It's failing to factor this
function:

>>> from  sympy import *
>>> var('a x')
(a, x)
>>> f21 = (1 - x)**2*(4 + x)/((2 + x)*(a + x))
>>> f=diff(f21,x,3)
>>> type(factor(f))
<class 'sympy.core.add.Add'>  <=== we should have gotten a Mul

>>> n,d=f.as_numer_denom()
>>> type(factor(n))
<class 'sympy.core.mul.Mul'>  <=== it *is* factorable

>>> factor(f)==f.expand()
True  <=== but factor(f) is not factoring it

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c19 Original author: https://code.google.com/u/117933771799683895267/

[smichr]

I posted a rebased-on-master version of polys named polys2 to my smichr github
account. How shall we keep track of work that is to be done on this? Will each person
that has worked on the upgrade eventually push their versions to github and then each
version will be merged? How do we keep track of what each person is doing and the
current state of this module?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c20 Original author: https://code.google.com/u/117933771799683895267/

[vks]

I think everyone can post his progress here. git is decentralized, so everyone can
work on his own branch and pull changes from others. So if Mateusz wanted to get your
changes to rebase on master, he could do something like

$ git pull chris polys2

Could you please rebase your branch on master? I'm getting conflicts when trying to
do this. Maybe your forgot to push your recent changes to github?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c21 Original author: https://code.google.com/u/Vinzent.Steinberg@gmail.com/

[smichr]

But who merges all the changes. I can keep pushing my changes here, but unless (like
for master) someone will oversee the updating as changes are made, you can end up
with a merge headache later.

My most recent re-rebase is uploaded.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c22 Original author: https://code.google.com/u/117933771799683895267/

[asmeurer]

I think you need to try again.  You left a merge line in test_basic.py.  Also, a ton of tests are failing now, so it is 
definitely a good idea I think to postpone the merge for this until 0.7.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c23 Original author: https://code.google.com/u/asmeurer@gmail.com/

[mattpap]

> But who merges all the changes. I can keep pushing my changes here, but unless
> (like for master) someone will oversee the updating as changes are made, you can
> end up with a merge headache later

That's right and we would like to avoid unnecessary headaches. So, the official polys
branch is polysn (where n is an integer) in my github repo. The latest is: http://github.com/mattpap/sympy-polys/commits/polys3 (each rebase will increase n).

Recent changes:

1. All tests, doctests, examples are OK.
2. Simplified polytools.py and improved coverage (75%).
3. Lots of bugfixes here and there (in polys of course).

Plan for this week (actually weekend) is to test 100% of polytools.py, finish
implementing algebraic factorization routines and write finally some documentation.

> Also, factor still should be doing something different. It's failing to factor
> this function:

Just use factor(f, frac=True). This will use cancel() to simplify the fraction and
factor both numerator and denominator, and return a factored fraction. The same was
implemented in sqf().

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c24 Original author: https://code.google.com/u/101069955704897915480/

[mattpap]

**Labels:** Polynomial  

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c25 Original author: https://code.google.com/u/101069955704897915480/

[smichr]

in light of Mateusz's comment above I will stop trying to rebase this on master.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c26 Original author: https://code.google.com/u/117933771799683895267/

[certik]

Let's improve this as Mateusz has said, then put it in (probably after the release).

**Labels:** -NeedsReview NeedsBetterPatch  

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c27 Original author: https://code.google.com/u/104039945248245758823/

[smichr]

I have a patch for consideration...factor() in polytools should not be returning two
arguments. 

The commit with memo to that effect is pushed to the polys2 branch. 

Also, I made a gcdfactor function for simplify and showed it being used in
Add.as_numer_denom(). That is a separate commit in the same branch. What do you think?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c28 Original author: https://code.google.com/u/117933771799683895267/

[mattpap]

> I have a patch for consideration...factor() in polytools should not be returning
> two arguments

I already fixed that in polys3 and written tests.

> Also, I made a gcdfactor function for simplify and showed it being used in
> Add.as_numer_denom(). That is a separate commit in the same branch. What do
> you think?

terms_gcd() function (equivalent of gcdfactor()) is available in polys3, so you don't
need to reimplement its functionality. The rest of the patch seems OK. The only
question is if using non-core functionality in Add.as_numer_denom() is the right
thing to do. I expect as_numer_denom() to be as fast as possible and not necessarily
returning optimal results. We have fraction() for "magic". But if terms_gcd() (or
gcdfactor()) is not slowing down as_numer_denom() too much then why not, at least we
cut down on cancel() in this case.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c29 Original author: https://code.google.com/u/101069955704897915480/

[smichr]

The functionality for gcdfactor is there, but a simple wrapper for someone who
doesn't want to get into Polys is not there. Compare

    expr = gcdfactor(expr) 

with

    if (expr.is_Add or expr.is_Mul):
        expo, poly = Poly(expr, expand=False).terms_gcd()
        if not all(i is S.Zero for i in expo):
            terms = [x**expo[i] for i,x in enumerate(poly.gens)]
            terms.append(poly.as_basic())
            expr = Mul(*terms)

That's a lot of overhead for someone that just wants to factor out the gcd from an
expression. gcdfactor isn't duplicating terms_gcd, it is just putting a non-Poly
wrapper on it. I put it in the simplify module since separatevars should be using
gcdfactor (and not pure factor).

As for fraction and magic...it's a little underwhelming. You have to do
as_numer_denom() to get an expression to have a numer and denom or else fraction
doesn't do anything for you...and even then, a simple removal of a gcd is not done. 

>>> eq = 1/x + 1/x**2
>>> fraction(eq)
(1/x + x**(-2), 1)
>>> n,d=(1/x**2+1/x).as_numer_denom()
>>> fraction(n/d)
(x + x**2, x**3)

I (and root finders) would really like to see

(1 + x, x**2) 

come back from eq.as_numer_denom(). Whether this slows things down too much, that's
an issue to consider. If it is, then having gcdfactor in simplify is a nice tool
since before doing as_numer_denom() one can pull out a gcd and run as_numer_denom on
the residual and build the result on their own:

### simulated output
>> eq = 1/x+1/x**2
>> gcdfactor(eq)
1/x*(1+1/x)
>> Mul(*[m.as_numer_denom() for m in make_list(_, Mul)]) #existing as_numer_denom()
(1+x)/x**2
###

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c30 Original author: https://code.google.com/u/117933771799683895267/

[mattpap]

1. I'm talking about this commit: http://github.com/mattpap/sympy-polys/commit/8e1c2fd03f29128eefb584bce28c3ddde42dcf2c Example:

In [3]: terms_gcd(x**3*y-x*y**3)
Out[3]: 
    ⎛ 2    2⎞
x⋅y⋅⎝x  - y ⎠

In [4]: factor(_)
Out[4]: x⋅y⋅(x - y)⋅(x + y)

The only difference is that you still have to add expand=False.

2. OK. I thought fraction() is something better than as_numer_denom(). So, is there
any reason to have both?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c31 Original author: https://code.google.com/u/101069955704897915480/

[vks]

I think we can drop fraction(), if we don't improve it. Currently as_numer_denom()
seems to be more powerful.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c32 Original author: https://code.google.com/u/Vinzent.Steinberg@gmail.com/

[smichr]

see issue 4830 .

Also, Poly.terms() is giving out monom, coeff whereas the old Poly.iter_terms() gave
out coeff, monom...is this intentional?

Referenced issues: #4830 Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c33 Original author: https://code.google.com/u/117933771799683895267/

[smichr]

I'm running on polys3 and observe the following:

>>> n
6*(-a + x**2)**8*(-a - 3*x**2) - 48*x**3*(-a + x**2)**6*(a*x + x**3) + 24*x*(-a +
x**2)**7*(a*x + x**3) + 24*x**2*(-a + x**2)**7*(a + 3*x**2) - 36*x**2*(-a + x**2)**8
+ 6*(-a + x**2)**9
>>> factor(n)
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 1854,
in factor
    coeff, factors = _factor(f)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 1846,
in _factor
    (coeff, factors), result = F.factor_list(**args), S.One
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 1099,
in factor_list
    result = f.rep.factor_list(**args)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyclasses.py", line 1419,
in factor_list
    coeff, factors = dmp_factor_list(f.rep, f.lev, f.dom, **args)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\factortools.py", line 1031,
in dmp_factor_list
    coeff, factors = _dmp_inner_factor(f, u, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\factortools.py", line 999,
in _dmp_inner_factor
    coeff, factors = dmp_zz_factor(f, v, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\factortools.py", line 923,
in dmp_zz_factor
    sign, H = dmp_zz_wang(g, u, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\factortools.py", line 848,
in dmp_zz_wang
    raise NotImplementedError("if this happened we need an extra loop here")
NotImplementedError: if this happened we need an extra loop here
>>>

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c34 Original author: https://code.google.com/u/117933771799683895267/

[asmeurer]

You are right.  fraction() is like as_numer_denom() except for a little "magic", but I think we should move all of the code there into 
as_numer_denom() as it is duplicated (actually the only "magic" it performs that as_numer_denom doesn't as far as I can tell is that 
it handles assumptions on variable exponents, which as_numer_denom should maybe be doing the same with a similar keyword 
argument.  See issue 4830 ).  

To Chris, separatevars() does want to use normal factor.  That way, it can separate 1 + x + y + x*y.  If you are doing it right, you 
should see a nice xpass test for it with the new polys module (actually it needs to be fixed to not bother with PolynomialError 
anymore).

Referenced issues: #4830 Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c35 Original author: https://code.google.com/u/asmeurer@gmail.com/

[smichr]

Also, is there a way to make monic work even when coefficients don't reduce to
multiples of the leading coefficient?

###
>>> Poly(3*x**3+5*x**2+7*x+11,x).monic()
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 1014,
in monic
    result = f.rep.monic()
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyclasses.py", line 1370,
in monic
    return f.per(dmp_ground_monic(f.rep, f.lev, f.dom))
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densetools.py", line 1391,
in dmp_ground_monic
    return dup_monic(f, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densetools.py", line 1385,
in dup_monic
    return dup_quo_ground(f, lc, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densearith.py", line 143,
in dup_quo_ground
    return [ K.quo(cf, c) for cf in f ]
  File "C:\Documents and Settings\chris\sympy\sympy\polys\algebratools.py", line 387,
in quo
    raise ExactQuotientFailed('%s does not divide %s in %s' % (b, a, self))
ExactQuotientFailed: 3 does not divide 5 in ZZ
###

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c36 Original author: https://code.google.com/u/117933771799683895267/

[smichr]

re comment 35:

Yes, factor should be used whenever there is a term in the resultant Mul that has
more than one symbol in it, otherwise it doesn't need to overwork itself. If one of
the terms is 1-x**2, there is no need to factor it into its two factors is there?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c37 Original author: https://code.google.com/u/117933771799683895267/

[mattpap]

> Also, Poly.terms() is giving out monom, coeff whereas the old Poly.iter_terms()
> gave out coeff, monom...is this intentional?

Yes, because this is coherent with dictionary representation, e.g. when iterating:

for monom, coeff in poly.terms():
    pass

for monom, coeff in poly.as_dict().iteritems():
    pass

Also now you can easily do dict(poly.terms()). 

> NotImplementedError: if this happened we need an extra loop here

This means that parameters for Wang's algorithm were set too low and the algorithm
wasn't able to find a suitable substitution integers for all but the main variable.
If you run factor() several times for this example, you'll see that it fails only
from time to time. In my local repo I fixed this already. I just run the algorithm
recursively, improving parameters if needed, until it returns without error.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c38 Original author: https://code.google.com/u/101069955704897915480/

[mattpap]

Forgot to say, that there always exists a suitable evaluation set. Our task is to
find the one which results in a univariate polynomial which has exactly the same
number of factors as the initial multivariate polynomial (correctness) and those
univariate factors are of minimal max-norm (speed).

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c39 Original author: https://code.google.com/u/101069955704897915480/

[smichr]

One other observation on polys3...when sending roots a Poly, even if the heuristics
fail, the whole process shouldn't fail if the Poly is 4th order or less--a call to
the roots_foo (e.g. roots_cubic) should be called:

in the following, poly was the 4th order Poly((-42 + 6*a**2)*x**4 + (96 +
48*a**2)*x**3 + (288 + 648*a + 144*a**2)*x**2 + (384 + 864*a + 624*a**2)*x + 192 +
432*a + 312*a**2 + 108*a**3, x, domain='ZZ[a]')

Traceback (most recent call last):
  File "C:\Documents and Settings\chris\sympy\my\___my_util\symp4.py", line 1036, in
solv1
    sols=roots(poly).keys()
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyroots.py", line 427, in
roots
    result = _try_decompose(f)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyroots.py", line 343, in
_try_decompose
    factors = f.decompose()
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 1052,
in decompose
    result = f.rep.decompose()
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyclasses.py", line 1393,
in decompose
    return map(f.per, dup_decompose(f.rep, f.dom))
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densetools.py", line 1849,
in dup_decompose
    result = _dup_decompose(f, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densetools.py", line 1809,
in _dup_decompose
    h = _dup_right_decompose(f, s, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densetools.py", line 1779,
in _dup_right_decompose
    g[s-i] = K.quo(coeff, i*r*lc)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\algebratools.py", line 387,
in quo
    raise ExactQuotientFailed('%s does not divide %s in %s' % (b, a, self))
ExactQuotientFailed: DMP([12, 0, -84], ZZ) does not divide DMP([48, 0, 96], ZZ) in ZZ[a]

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c40 Original author: https://code.google.com/u/117933771799683895267/

[mattpap]

> Also, is there a way to make monic work even when coefficients don't reduce to
> multiples of the leading coefficient?

Yes, there is only one issue:

In [1]: f = 3*x**3+5*x**2+7*x+11

In [2]: g = Poly(f)

In [3]: g.get_domain()
Out[3]: ZZ

By default Poly tries to find a minimal domain for coefficients. In this case it's ZZ
which supports only exact division. There are several solutions to fix this,
depending what you are actually doing.

You can use monic() function:

In [4]: monic(f)
Out[4]: 
                2     
       7⋅x   5⋅x     3
11/3 + ─── + ──── + x 
        3     3       

If it is given basic expression it will default to a field (QQ in this case). Or you
can force field usage in Poly directly:

In [5]: g = Poly(f, field=True)

In [6]: g.get_domain()
Out[6]: QQ

In [7]: g.monic()
Out[7]: Poly(x**3 + 5/3*x**2 + 7/3*x + 11/3, x, domain='QQ')

You can always set an explicit domain:

In [8]: g = Poly(f, domain='QQ')

In [9]: g.monic()
Out[9]: Poly(x**3 + 5/3*x**2 + 7/3*x + 11/3, x, domain='QQ')

Most functions in polytools.py are strict and they won't change initial ground
domain. There are however some which do this intentionally, e.g. Poly.integrate(),
for users' convenience.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c41 Original author: https://code.google.com/u/101069955704897915480/

[mattpap]

> ExactQuotientFailed: DMP([12, 0, -84], ZZ) does not divide DMP([48, 0, 96], ZZ)
> in ZZ[a]

That's a bug. There should be K.exquo instead of K.quo in:

    g[s-i] = K.quo(coeff, i*r*lc)

In new polys exquo() does what quo() did previously. Now quo() computes also the
remainder and checks if it's zero. If not it raises ExactQuotientFailed exception.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c42 Original author: https://code.google.com/u/101069955704897915480/

[smichr]

PolynomialError should be imported in polyclasses.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c44 Original author: https://code.google.com/u/117933771799683895267/

[mattpap]

I went ahead and finished refactoring polytools.py and writing tests for it. Now
coverage for this file is 100%. Also added tests for monomialtools.py and fixed all
bugs mentioned in this issue so far (e.g. ExtraneousFactors exception). All this is
in polys3 branch. One more day like today and polys will have near 100% coverage.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c45 Original author: https://code.google.com/u/101069955704897915480/

[smichr]

Also, perhaps roots_cubic and roots_quartic should use the monic function rather than
the monic method for reasons discussed in comment 41:

###
>>> var('x');p=Poly(3*x**3+5*x**2+7*x+11)
x
>>> roots(p)
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyroots.py", line 417, in
roots
    result = _try_decompose(f)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyroots.py", line 337, in
_try_decompose
    for r in _try_heuristics(h):
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyroots.py", line 380, in
_try_heuristics
    result += roots_cubic(f)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyroots.py", line 36, in
roots_cubic
    one, a, b, c = f.monic().all_coeffs()
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polytools.py", line 1053,
in monic
    result = f.rep.monic()
  File "C:\Documents and Settings\chris\sympy\sympy\polys\polyclasses.py", line 1400,
in monic
    return f.per(dmp_ground_monic(f.rep, f.lev, f.dom))
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densetools.py", line 1391,
in dmp_ground_monic
    return dup_monic(f, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densetools.py", line 1385,
in dup_monic
    return dup_quo_ground(f, lc, K)
  File "C:\Documents and Settings\chris\sympy\sympy\polys\densearith.py", line 143,
in dup_quo_ground
    return [ K.quo(cf, c) for cf in f ]
  File "C:\Documents and Settings\chris\sympy\sympy\polys\algebratools.py", line 400,
in quo
    raise ExactQuotientFailed('%s does not divide %s in %s' % (b, a, self))
ExactQuotientFailed: 3 does not divide 5 in ZZ
>>>

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c46 Original author: https://code.google.com/u/117933771799683895267/

[smichr]

or should that be exquo in line 143 of densearith.py? 

return [ K.exquo(cf, c) for cf in f ]

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c47 Original author: https://code.google.com/u/117933771799683895267/

[smichr]

Yes, making the change in densearith.py solves the problem:

>>> var('x');p=Poly(3*x**3+5*x**2+7*x+11)
>>> roots(p)
{-1/3 + (65/54 + 5*21**(1/2)/18)**(1/3)*(1/2 + I*3**(1/2)/2) - 5/(9*(1/2 +
I*3**(1/2)/2)*(65/54 + 5*21**(1/2)/18)**(1/3)): 1, -1/3 + 5/(9*(65/54 +
5*21**(1/2)/18)**(1/3)) - (65/54 + 5*21**(1/2)/18)**(1/3): 1, -1/3 + (65/54 +
5*21**(1/2)/18)**(1/3)*(1/2 - I*3**(1/2)/2) - 5/(9*(1/2 - I*3**(1/2)/2)*(65/54 +
5*21**(1/2)/18)**(1/3)): 1}
>>>

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c48 Original author: https://code.google.com/u/117933771799683895267/

[mattpap]

> or should that be exquo in line 143 of densearith.py? 

No, because you are changing dup_quo_ground(). See a few lines below that, there is a
function dup_exquo_ground() which uses K.exquo (actually its replacement). In a field
there is no difference between quo and exquo, and in a ring there is an assumption
that K.exquo is equivalent to // operator and that operator really does its job (this
is not true for / operator for certain types). So you can use both // and K.exquo.

As to c#46 solution is to avoid using Poly explicitly if not sure. If you pass Basic
expression instead, roots() will force field=True and there will be no error. It
would be however useful it didn't fail in this simple case.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c49 Original author: https://code.google.com/u/101069955704897915480/

[smichr]

I wonder if roots(basic) should accept a basic expression or if it should perhaps
give a tuple as a result. With the expanded capabilities of polys, you can now get a
solution but may not know what generator it is given in terms of:

>>> roots(exp(x)+exp(2*x)-2)
{1: 1, -2: 1}

You will have to make a Poly of the expression to see what the generator was:

>>> Poly(exp(x)+exp(2*x)-2)
Poly(exp(x)**2 + exp(x) - 2, exp(x), domain='ZZ')

Do you think roots should return (gen, dict) as a result?

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c50 Original author: https://code.google.com/u/117933771799683895267/

[mattpap]

> Do you think roots should return (gen, dict) as a result?

It might be useful. Actually I didn't care much about root-finding code yet, when
working on new polys. For example there should be a method roots() in Poly and the
function roots() should be implemented in polytools.py, whereas root-finding logic
should be placed in polyroots.py (e.g. _try_heuristics()). I'm going to work on this
in near future.

btw. Today I implemented minpoly() function (based on Groebner bases). Available in
polys3 branch.

Original comment: http://code.google.com/p/sympy/issues/detail?id=1598#c51 Original author: https://code.google.com/u/101069955704897915480/