Documentation preview for this PR (built with commit 973dffa; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

@fchapoton @mantepse @tscrim
I will be happy to have your opinion on this PR, in particular concerning the three following points:

• Is it reasonable to set a coercion map as I do? Could it break easily the coercion system?
• I changed the behavior of PowerSeriesRing and LaurentSeriesRing to return lazy series when called with default_prec=Infinity. How does it sound? And did I implement this correctly?
• I changed the behavior of _first_ngens to include generators of the base when the required number of generators exceeds that of the ring.

Of course, any other comment is welcome!

The more I think about it, the more I am convinced that completion deserves their own class. Indeed, we would like some methods available, e.g. polynomial_ring* (to recover the underlying polynomial ring) or place(to recover the ideal at which we completed). Also the methodconstruction` does not reflect the construction.

The reason why I rejected this option at first is that it somehow breaks the unicity assumption for parents, in the sense that I will need to include A[x] and P in the key. In particular, we will have the following:

sage: A.<t> = QQ[]
sage: A1.<u> = A.completion(t - 1)
sage: A2.<u> = A.completion(t - 2)
sage: A1 is A2
False

although A1 and A2 are both power series ring in u over $\mathbb Q$.
Nonetheless, I underline that this does not work currently because it overloads the coercion graph.

So I finally think that having two different copies of PowerSeriesRing(QQ, name='u') is a lesser problem.
Moreover, this also solves the problem with _first_ngens because, here, I can simply redefine _defining_names.

What do you think about this?

okay, makes sense. In that case, I suggest design the API as follows.

internally elements of B are represented by elements of C = PowerSeriesRing(A.quotient(f), 'u').
define a coercion morphism A → B sending x to u+a, where a = A.quotient(f)(x).
then B(f) is an uniformizer of B, although it likely would have infinitely many terms internally.
for c in C, write B._from_internal_repr(x) to be the element of B representing the above.

the above are what you did. The below aren't.

We must not expose the internal implementation detail above to the user.
For all they know, we might have sent $x$ to $-u+a$, and they wouldn't have noticed.
To do so, implement the following behavior:

1. Indexing
   let $S ⊆ A$ be the subset of elements that is the image of the operation a ↦ a % f,
   there's a natural set bijection (but not ring bijection!) between $S$ and A.quotient(f).
   each element $b ∈ B$ can be represented uniquely as $b = c_0 + c_1 B(f) + c_2 B(f)^2 + ⋯$
   then for all integer $i$, b[i] = c_i.

   the % operation actually has a nice math description, namely the image are the polynomials
   of degree less than $f$. For multivariate polynomial ring the choice is somewhat arbitrary
   (based on Gröbner basis).

   On the other hand, this indexing operation cannot be particularly efficiently implemented. We may want to implement .list() in almost-linear time, then tell the user to use that (to avoid quadratic behavior when indexing multiple elements), or for even more performance-critical code, they may access the internal detail in one way or the other.

2. String representation
   Printing out $b$ must print it as c_0 + c_1 * f + c_2 * f^2 + O(f^n), similar to p-adic.

3. More coercion
   there is a coercion A → B as above. There must not be a coercion B → C or vice versa.

4. Generator
   B has a single generator which is the image of x over the coercion map.

Anyway, the above would require making separate parent and element classes.

I am not really in favor to choose $f$ as a uniformizer (i.e., as the variable of the power series ring) because, as I said, doing this transformation requires to solve the equation $f(x) = u$ in $B$ (which is certainly possible by Hensel's lemma but affects the complexity and the ease of implementation).
What's actually the problem with choosing another uniformizer?

I'm not also sure that x is the best choice for a generator. What is your argument for this?
In any case, we will need letters for the uniformizer, and the (Teichmüller representative of the) generator of the residue field, won't we?

I invite @roed314 to the discussion; he might have good suggestions.

I think there was a typo somewhere.

Let $𝔭$ be the maximal ideal, then any element of $𝔭 ∖ 𝔭^2$ is an uniformizer.

You may send x to B._from_internal_repr(u+a). Then B(f) = B._from_internal_repr(f(u+a)) is also an uniformizer, apart from u.

(B(f) would be some polynomial in B._from_internal_repr(u) with zero constant term and finitely many nonzero coefficients then.)

we will need letters for the uniformizer

I guess you can give the user

• .canonical_uniformizer() = B(f), string representation e.g. 1 * (x^2+x+1), somewhat complicated internal representation.
• .efficient_uniformizer() = B._from_internal_repr(u), string representation depends on implementation detail.

No solving system of equation needed either way.

[edit: okay I guess Witt vectors display format may be more convenient in some settings, but mimicking the current behavior of Zp is more elementary, and you don't have to choose a generator name.]

the (Teichmüller representative of the) generator of the residue field

The element B(x) (after mapping to the residue field) is a generator of the residue field. This will just have string representation x.

[sorry, mistake. I guess you can provide a method for this, such thing would have representation like x + something * (x^2+x+1)^2 + ⋯, I guess. Internally it would be efficiently represented as B._from_internal_repr(a).]

I'm not also sure that x is the best choice for a generator. What is your argument for this?

Same thing happens in quotient of polynomial ring, e.g. in A.<x> = QQ[]; S.<xbar> = A.quotient(f), then xbar = S(x).

If the user wants some letter, I guess one can do something similar to print_mode in Zp constructor. Note that by default Zp elements directly use the prime p instead of any new character as well.

Let me summarize because I'm getting lost.

Let $K$ be a field and $A = K[x]$. Let also $\mathfrak p$ be an irreducible polynomial in $A$.
By definition, the completion of $A$ at $\mathfrak p$ is
$A_{\mathfrak p} := \varprojlim_n A/\mathfrak p^n$.
This gives a first way to handle $A_{\mathfrak p}$ and represent elements in there: an element is just a compatible sequence of elements of $A/\mathfrak p^n$ and we can work with them by just retaining the $n$-th element of the sequence (we then work at precision $\mathfrak p^n$).
It's basically what is done for the $p$-adics.

However, it turns out that $A_{\mathfrak p}$ is also isomorphic to a ring of power series with coefficients in the residue field $L := A/\mathfrak p$ (which is, of course, not the case for the $p$-adics). There exist actually many such isomorphisms: basically each choice of a uniformizer of $A_{\mathfrak p}$, i.e. an element of $\mathfrak p A_{\mathfrak p} \backslash \mathfrak p^2 A_{\mathfrak p}$ provides such an isomorphism (actually even several of them).

In what follows, I consider two of them. Denoting by $a$ the image of $x$ in $L$, here are they:
$\alpha : A_{\mathfrak p} \to L[[u]]$, $x \mapsto u + a$,
$\beta : A_{\mathfrak p} \to L[[v]]$, $\mathfrak p \mapsto v$.
It turns out that $\alpha$ is easy to compute since we know the image of $x$. On the other hand, computing $\beta$ requires to solve the equation $\mathfrak p(x) = v$ in $L[[v]]$ and, similarly, computing $\alpha^{-1}$ and $\beta^{-1}$ also requires to solve equations in order to compute the image of $a$ in $A_{\mathfrak p}$.

My proposal is that A.completion(p) just returns $L[[u]]$, so that we just need to compute $\alpha$ which defines the canonical inclusion $A \to L[[u]]$.

I'm not sure what's your proposal. My guess is that you would like to carry computations in $L[[u]]$ and print elements as in $L[[v]]$ (or maybe $A_{\mathfrak p}$, I'm not sure). Is it correct?
In any case, I think you need to compute $\alpha^{-1}$ (and possibly also $\beta$).

Of course, another different option would be to just implement operations in $A_{\mathfrak p}$ but I prefer not doing it because it requires a full new implementation and it will not be possible to take advantage of what is already implemented for power series.

Oh, I just noticed that computing $\alpha \circ \beta^{-1}$ is also easy because it acts trivially on $L$ and maps $v$ to $\mathfrak p(u+a)$. (However, going in the other direction looks more difficult.) Is it your point?

Let me also mention that in the $p$-adic case, although computing the completion of number fields is currently not implemented, it would also probably require at some point to provide names for the "unramified variable" and the "uniformizer" since general $p$-adic fields need them.

I will try to take a look at the code and discussion this evening.

Quick comments:

• I think implementing the completion as a class is good. However, I don't think they need whole new class, but instead could be used on top of the existing (Lazy)PowerSeriesRing by an additional argument of the defining maximal ideal. I think this should be less work to implement.
• Subsequently, having non-equal parents for different (i.e., nonequal) ideals is good. Operations between them don't (immediately) make sense. It also gives you a lot of control over the print representation.
• I don't think _first_ngens should start pulling things from the base ring. It should just return all the generators of the ring at hand. I feel like otherwise it is breaking a (silent) promise.

I just pushed a rough implementation (without doctest so far) using RingExtension.
I think I am rather happy with the result but I'd like to have your feedback.

Here is a short demo:

sage: A.<x> = GF(7)[]
sage: Ap = A.completion(x^2 + 2*x + 2, prec=5)
sage: y = Ap(x)
sage: y
x
sage: y^20
5 + 4*x*(x^2 + 2*x + 2) + (5*x + 4)*(x^2 + 2*x + 2)^2 + (5*x + 5)*(x^2 + 2*x + 2)^3 + (6*x + 5)*(x^2 + 2*x + 2)^4 + ...
sage: y.inverse_of_unit()
(3*x + 6) + (5*x + 3)*(x^2 + 2*x + 2) + (6*x + 5)*(x^2 + 2*x + 2)^2 + (3*x + 6)*(x^2 + 2*x + 2)^3 + (5*x + 3)*(x^2 + 2*x + 2)^4 + O((x^2 + 2*x + 2)^5)
sage: a = y.teichmuller()
sage: a
x + (4*x + 4)*(x^2 + 2*x + 2) + (3*x + 3)*(x^2 + 2*x + 2)^2 + (6*x + 6)*(x^2 + 2*x + 2)^3 + ...
sage: a^2 + 2*a + 2
0

We can also recover the underlying power series ring as follows:

sage: Ap.power_series(names=('u','a'))
Power Series Ring in u over Finite Field in a of size 7^2

However, the construction S.<u,a> = Ap.power_series() does not work because _first_ngens only returns one generator. I continue to think to it would be nice to modify this.

Finally, I used the method quotient instead of extension; the advantage is that it does not require a variable name since it uses xbar as a default.
Hence now the following works:

sage: S.<u> = Ap.power_series()
sage: S
Power Series Ring in u over Univariate Quotient Polynomial Ring in xbar over Finite Field of size 7 with modulus x^2 + 2*x + 2

Looks nice, though I haven't looked too carefully.

I wonder what's a good way out of this. Maybe refactor the code such that

A.<b, c> = f()

desugar to

A, (b, c) = objgens(f(names=('b', 'c')));

then f can return a special object with .objgens()[0] does not return self to "workaround" the issue.

Caveat being currently .<> sugar allows specifying less generator than .gens(), but it's not clear why should anyone want to use this. (e.g. R.<x> = QQ['x', 'y'])

Maybe

A, (b, c) = _objgens_allow_less_generators(f(names=('b', 'c')), 2);

A question: each time, I change one line in a pyx file, sage -b/make build recompiles all the Sage library; do you know why and what should I do to prevent this?

Yes, this is terrible. This happens for everybody. This is supposed to be cured by the switch to meson in #39030

you can already use meson using conda/mamba (I recommend mamba as it's faster to install). See https://doc.sagemath.org/html/en/installation//meson.html

Anyway, I think make build also have incremental build, which is why you (at least, I did) see the incremental compilation time to be 2 minutes instead of 5 hours. But for meson, it could be as little as 5 seconds.

Il faut pas de return dans une methode __init__

The point is that the natural map A → L [ [ u ] ] takes x to a + u and so it is not the composite of the two coercion maps A → L → L [ [ u ] ] which takes x to a . (Here L = A / p , see my previous messages.)

Okay, thank you for the summary. As a reminder $A = K[x]$. While the underlying model for this would be $L[[u]]$, the actual base ring (field) for the parent should be $K$ in this case as we are thinking of this as a $K$-algebra, right? So there would be no automatic coercion from $L$. If instead you wanted to implement it with the base ring being $L$, you can disable the automatic base ring coercion (or set it to a different map) by using _coerce_map_from_base_ring().