[with patch, not ready for review] bug with power series and/or p-adics

This is strange:

sage: (derivative(x)).parent()
Power Series Ring in s over Laurent Series Ring in t over Eisenstein Extension of 11-adic Field with capped relative precision 10 in a defined by (1 + O(11^10))*x^2 + (10*11 + 10*11^2 + 10*11^3 + 10*11^4 + 10*11^5 + 10*11^6 + 10*11^7 + 10*11^8 + 10*11^9 + 10*11^10 + O(11^11))
sage: y.parent()
Power Series Ring in s over Laurent Series Ring in t over Eisenstein Extension of 11-adic Field with capped relative precision 10 in a defined by (1 + O(11^10))*x^2 + (10*11 + 10*11^2 + 10*11^3 + 10*11^4 + 10*11^5 + 10*11^6 + 10*11^7 + 10*11^8 + 10*11^9 + 10*11^10 + O(11^11))
sage: (derivative(x)).parent == y.parent()
False

Here is a much smaller example that triggers the bug:

sage: R.<u> = QQ[]
sage: p = 11
sage: K = Qp(p,10)
sage: J.<a> = K.extension(u^2-p)
sage: C.<t> = LaurentSeriesRing(J)
sage: D.<s> = PowerSeriesRing(C)

sage: y = t + (t^2 + O(t^3))*s
sage: h = y / y
[boom]

Interestingly, if we use the vanilla p-adics instead of a ramified extension, we get a different traceback, which seems to indicate a different (possibly related?) bug:

sage: R.<u> = QQ[]
sage: p = 11
sage: K = Qp(p,10)
sage: C.<t> = LaurentSeriesRing(K)
sage: D.<s> = PowerSeriesRing(C)

sage: y = t + (t^2 + O(t^3))*s

sage: h = y / y
Traceback (most recent call last):
...
File "/Users/david/sage-2.11/local/lib/python2.5/site-packages/sage/rings/polynomial/padics/polynomial_padic_capped_relative_dense.py", line 379, in __getslice__
    raise IndexError, "list index out of range"
IndexError: list index out of range

Even series inversion gives the list index breakage:

sage: p = 11
sage: K = Qp(p,10)
sage: C.<t> = LaurentSeriesRing(K)
sage: D.<s> = PowerSeriesRing(C)
sage: z = 1 + (t^3 + O(t^13))*s
sage: 1/z

Actually even just squaring something really simple can break:

sage: K = Qp(p,10)
sage: C.<t> = LaurentSeriesRing(K)
sage: D.<s> = PowerSeriesRing(C)
sage: z = (1 + O(t)) + t*s^2
sage: z * z
[boom]

ok, the bug is in the polynomial mult instead of the power series mult:

sage: K = Qp(p,10)
sage: C.<t> = LaurentSeriesRing(K)
sage: D.<s> = PolynomialRing(C)
sage: z = (1 + O(t)) + t*s^2
sage: z * z
[boom]

At least part of this bug is caused by some broken code in p-adic polynomial getslice. Here is the example from the docstring:

sage: K = Qp(13,7)
sage: R.<t> = K[]       
sage: a = 13^7*t^3 + K(169,4)*t - 13^4
sage: a[1:2]
(13^2 + O(13^4))*t

but this one dies:

sage: t[0:1]
[boom]

A smaller example of the original bug (this example involves only multiplication, no division):

sage: R.<u> = QQ[]
sage: p = 11
sage: K = Qp(p,10)
sage: J.<a> = K.extension(u^2-p)
sage: C.<t> = LaurentSeriesRing(J)
sage: D.<s> = PowerSeriesRing(C)
sage: x = t + (t^2 + O(t^3))*s
sage: y = t^(-1) + (-1 + O(t))*s + O(s^2)
sage: x * y
[boom]

Even simpler (doesn't need the outer PowerSeries? ring):

R.<u> = QQ[]
p = 11
K = Qp(p,10)
J.<a> = K.extension(u^2-p)
C.<t> = LaurentSeriesRing(J)
D.<s> = PolynomialRing(C)
x = t + (t^2 + O(t^3))*s
y = t^(-1) + (-1 + O(t))*s
h = x * y
[boom]

I just did some more debugging with David Roe and I'm just going to write down our discoveries since I'm exhausted now.

The problem is, at some point a Laurent series object is coming into existence whose unit part is a power series whose coefficients are all indistinguishable from zero (but not actually zero). Then when calling valuation on that laurent series, something blows up. It only happens when the base ring is unramified extensions of Qp, since something in the underlying polynomial class hasn't been implemented properly. (Maybe.)

The actual exception generated is a bit misleading. It comes up in the last line of this code in polynomial_element.pyx:

        if p is None:
            for k from 0 <= k <= self.degree():
                if self[k]:
                    return ZZ(k)
                
        if not isinstance(p, Polynomial) or not p.parent() is self.parent():
            raise TypeError, "The polynomial, p, must have the same parent as self."

It's entering the loop with p == None, but it should never drop out of the bottom of the loop (it expects to find a nonzero-coefficient).

Getting closer now.

Set up some polynomial rings over Qp and an unramified extension of Qp. Notice that the polynomial ring over Qp itself uses a specialised type, whereas over the ramified extension uses the generic polynomial ring:

sage: R.<u> = QQ[]
sage: p = 7
sage: K = Qp(p, 10)
sage: J.<a> = K.extension(u^2 - p)
sage: 
sage: S.<x> = PolynomialRing(K)
sage: T.<y> = PolynomialRing(J)
sage: 
sage: type(S)
<class 'sage.rings.polynomial.polynomial_ring.PolynomialRing_dense_padic_field_capped_relative'>
sage: type(T)
<class 'sage.rings.polynomial.polynomial_ring.PolynomialRing_field'>

The specialised type lets you have polynomials whose coefficients are indistinguishable from zero, but it still knows that the whole polynomial is "zero":

sage: f = S([O(p^5), O(p^5)])
sage: f
(O(7^5))*x + (O(7^5))
sage: f.__nonzero__()
False

But the generic polynomial type just strips them off, since it thinks they are zero:

sage: g = T([O(a^5), O(a^5)])
sage: g
0
sage: g.coeffs()
[]
sage: g.__nonzero__()
False

Now we can try forcing the generic code to let us have coefficients indistinguishable from zero:

sage: g = T([O(a^5), O(a^5)], check=False)
sage: g
0
sage: g.coeffs()
[O(a^5), O(a^5)]
sage: g.__nonzero__()
True

And that's basically the problem. The generic polynomial code thinks that the polynomial is nonzero because its coefficient list is non-empty.

After discussion with broune and cwitty on IRC, the proposed solution is as follows:

• in the long run, implement the specialised class for polynomials over ramified extension rings; but we don't have time to do this right now.
• change the generic polynomial's __nonzero__ method to first check the base ring's is_exact method. If the base ring is not exact, we need to check whether each coefficient of the polynomial is zero.
• for performance reasons, change is_exact to a cpdef method (also, probably move it up the hierarchy from rings to parents).

Attached patch should fix this bug.

(Note: this does not fix the "list index" exception issue, which also came up in the above discussion. That is now a separate ticket: #3184)

After this patch, there is still a bug (thanks Kiran Kedlaya for pointing this out). Here is an example:

sage: R.<u> = QQ[]
sage: p = 11
sage: K = Qp(p,10)
sage: J.<a> = K.extension(u^2-p)
sage: C.<t> = LaurentSeriesRing(J)
sage: D.<s> = PowerSeriesRing(C)
sage: y = 1 - (t + O(t^2))*s + O(s^5)
sage: 1/y
(1 + O(a^20)) + ((1 + O(a^20))*t + O(t^2))*s + O(s^5)

The output is wrong; there should be s², s³, s⁴ terms.

Kiran's example was:

sage: R.<u> = QQ[]
sage: p = 11
sage: K = Qp(p,10)
sage: J.<a> = K.extension(u^2-p)
sage: C.<t> = LaurentSeriesRing(J)
sage: D.<s> = PowerSeriesRing(C)

sage: y = t + (t^2 + O(t^3))*s
sage: h = y / y
sage: h
(1 + O(a^20)) + O(t^2)*s + O(t^2)*s^2 + O(t^3)*s^3 + O(t^4)*s^4 +
O(t^5)*s^5 + ((7 + 10*a^2 + 10*a^4 + 10*a^6 + 10*a^8 + 10*a^10 + 10*a^12
+ 10*a^14 + 10*a^16 + 10*a^18 + O(a^20))*t^6 + O(t^7))*s^6 + O(s^20)

which is distinguishable from 1.

Someone else seems to have hit a related bug, see the following thread:

http://groups.google.com/group/sage-devel/browse_thread/thread/56a73f331e6b2977

Any movement here? This ticket seems to have gone stale.

Cheers,

Michael