Pierre-Étienne Meunier
email: pmeunier@caltech.edu.

The self-assembly of paths and squares at temperature 1

published or not. The documents may come

The self-assembly of paths and squares at temperature 1

Introduction

Self-assembly is the process through which unorganized, simple, components automatically coalesce according to simple local rules to form some kind of target structure. It sounds simple, but the end result can be extraordinary. For example, researchers have been able to self-assemble a wide variety of structures experimentally at the nanoscale, such as regular arrays [START_REF] Winfree | Design and self-assembly of two-dimensional DNA crystals[END_REF], fractal structures [START_REF] Fujibayashi | Toward reliable algorithmic self-assembly of DNA tiles: A fixed-width cellular automaton pattern[END_REF][START_REF] Paul | Algorithmic self-assembly of DNA Sierpinski triangles[END_REF], smiling faces [START_REF] Paul | Folding DNA to create nanoscale shapes and patterns[END_REF][START_REF] Wei | Complex shapes self-assembled from single-stranded DNA tiles[END_REF], DNA tweezers [START_REF] Yurke | A DNA-fuelled molecular machine made of DNA[END_REF], logic circuits [START_REF] Qian | Scaling up digital circuit computation with DNA strand displacement cascades[END_REF][START_REF] Seelig | Enzyme-free nucleic acid logic circuits[END_REF], neural networks [START_REF] Qian | Neural network computation with DNA strand displacement cascades[END_REF], and molecular robots [START_REF] Lund | Molecular robots guided by prescriptive landscapes[END_REF]. These examples are fundamental because they demonstrate that self-assembly can, in principle, be used to manufacture specialized geometrical, mechanical and computational objects at the nanoscale. Potential future applications of nanoscale self-assembly include the production of smaller, more efficient microprocessors and medical technologies that are capable of diagnosing and even treating disease at the cellular level.

Controlling nanoscale self-assembly for the purposes of manufacturing atomically precise components will require a bottom-up, hands-off strategy. In other words, the self-assembling units themselves will have to be "programmed" to direct themselves to assemble efficiently and correctly. Molecular self-assembly is rapidly becoming a ubiquitous engineering paradigm, and we need to develop a theory to inform us of its algorithmic capabilities and ultimate limitations.

In 1998, Erik Winfree [START_REF] Winfree | Algorithmic Self-Assembly of DNA[END_REF] introduced the abstract Tile Assembly Model (aTAM), a simplified discrete mathematical model of algorithmic DNA nanoscale self-assembly pioneered by Seeman [START_REF] Nadrian | Nucleic-acid junctions and lattices[END_REF]. The aTAM is an asynchronous nondeterministic cellular automaton that models crystal growth processes. Put another way, the aTAM essentially augments classical Wang tiling [START_REF] Hao | Proving theorems by pattern recognition -II[END_REF] with a mechanism for sequential growth of a tiling. This contrasts with Wang tiling in which only the existence of a valid mismatch-free tiling is considered, and not the order of tile placement. In the aTAM, the fundamental components are translatable but un-rotatable square or cube tiles whose sides are labeled with colored glues colors, each with an integer strength. Two tiles that are placed next to each other interact if the glue colors on their abutting sides match, and they bind if the strengths on their abutting sides match and sum to at least a certain (integer) temperature. Self-assembly starts from a seed tile type and proceeds nondeterministically and asynchronously as tiles bind to the seed-containing-assembly. Despite its deliberate simplification, the aTAM is a computationally expressive model. For example, by using cooperative binding (that is, by having some of the tiles bind on two or more sides), Winfree [START_REF] Winfree | Algorithmic Self-Assembly of DNA[END_REF] proved that it is Turing universal, which implies that self-assembly can be directed by a computer program. Here, we study noncooperative binding.

Tile self-assembly in which tiles can be placed only in a noncooperative fashion is colloquially referred to as "temperature-1 self-assembly". Despite the esoteric name, this is a fundamental and ubiquitous form of growth: it refers to growth from growing and branching tips where each new tile is added if it can match on at least one side.

It has been known for some time that a more general form of growth where some of the tiles must match on two or more sides, i.e. cooperative growth, leads to highly non-trivial behavior: arbitrary Turing machine simulation [START_REF] Lathrop | Computability and complexity in self-assembly[END_REF][START_REF] Paul | The program-size complexity of self-assembled squares (extended abstract)[END_REF], efficient production of n × n squares and other simple shapes using Θ(log n/ log log n) tile types [START_REF] Adleman | Running time and program size for self-assembled squares[END_REF], efficient production of arbitrary finite connected shapes using a number of tile types that is within a log factor of the Kolmogorov complexity of the shape [START_REF] Soloveichik | Complexity of self-assembled shapes[END_REF], and even intrinsic universality: the existence of a single tile set that simulates arbitrary tile assembly systems [START_REF] Doty | The tile assembly model is intrinsically universal[END_REF].

However, the capabilities of two-dimensional noncooperative self-assembly remain largely unknown: several generalizations and restrictions of this model have been studied, that conjectured in all cases that noncooperative binding could not be as powerful as cooperative binding [START_REF] Chandran | Tile complexity of approximate squares[END_REF][START_REF] Cook | Temperature 1 self-assembly: deterministic assembly in 3D and probabilistic assembly in 2D[END_REF][START_REF] Doty | Limitations of self-assembly at temperature 1[END_REF][START_REF] Maňuch | Two lower bounds for self-assemblies at temperature 1[END_REF][START_REF] Patitz | Exact shapes and Turing universality at temperature 1 with a single negative glue[END_REF][START_REF] Paul | The program-size complexity of self-assembled squares (extended abstract)[END_REF]. The first fully general separation results were only proven recently [START_REF] Meunier | Intrinsic universality in tile self-assembly requires cooperation[END_REF], in the context of intrinsic universality [START_REF] Demaine | The two-handed tile assembly model is not intrinsically universal[END_REF][START_REF] Doty | The tile assembly model is intrinsically universal[END_REF][START_REF] Doty | Intrinsic universality in self-assembly[END_REF]. However, the computational capabilities, in the Turing sense, of this model, remain largely unknown.

The conjecture that we prove in this paper was first stated by Rothemund and Winfree [START_REF] Paul | The program-size complexity of self-assembled squares (extended abstract)[END_REF] 1 : the minimal number of tile types to assemble n × n squares is 2n -1. A restriction of this result, where it is required that all the tiles be assembled with all their neighbors, appeared in the same paper. Moreover, computing the minimal number of tile types required to deterministically assemble a shape from a seed of size one is known to be NP-complete [START_REF] Adleman | Combinatorial optimization problems in selfassembly[END_REF], and Σ p 2 -complete in the non-deterministic case [START_REF] Bryans | The power of nondeterminism in self-assembly[END_REF].

Main results

Although a number of terms have not been formally defined, we give an overview of our two main results now. See the definitions in section 2. Our ultimate goal is to prove Rothemund and Winfree's conjecture [START_REF] Paul | The program-size complexity of self-assembled squares (extended abstract)[END_REF] that a tile assembly system assembling only squares of size n × n, from a seed of size 1, has at least 2n -1 tile types. Our first result disproves a statement stronger than this conjecture; namely, that the tile complexity of a square is the same as the tile complexity of its diagonal. Although widely believed, this statement is false: Theorem 3.1. Let n be an integer. There is a tile assembly system T n = (T n , σ n , 1), and two 1 In Figure 2 of that paper 2 points A, B ∈ Z 2 , such that --→ AB 1 = n, the terminal assemblies of T n are all finite, they all include a path from A to B, and

|T n | = 4n/5 + O(1).
The fact that the constructions of this theorem are possible, even though they are quite elementary, is not obvious at all; indeed, the intuition from words and automata theory is that any attempt to "reuse" tile types will enable us to "pump" the path, as in the pumping lemma of finite automata [START_REF] Sipser | Introduction to the Theory of Computation[END_REF], and thus any tile assembly system that can produce paths repeating tile types will also be able to produce ultimately periodic, infinite paths. This intuition is valid in a restricted setting where two adjacent tiles always agree on their abutting sides [START_REF] Doty | Limitations of self-assembly at temperature 1[END_REF][START_REF] Maňuch | Two lower bounds for self-assemblies at temperature 1[END_REF].

Then, we will prove the following theorem, which gives the optimal lower bound on the tile complexity of squares, when growth starts from a corner: Theorem 4.2. Let T = (T, σ, 1) be a temperature 1 tile assembly system, with σ a single tile at (0, 0), and n an integer. If all terminal assemblies producible by T are of domain {0, . . . , n -1} 2 , then |T | ≥ 2n -1.

Definition and preliminaries

We begin by defining the two-dimensional abstract tile assembly model. A tile type is a unit square with four sides, each consisting of a glue label and a nonnegative integer strength. We call a tile's sides north, east, south, and west, respectively, according to the following picture:

West East

South

North

We assume a finite set T of tile types, but an infinite supply of copies of each type. An assembly is a positioning of the tiles on the discrete plane Z 2 , that is, a partial function α : Z 2 T . We say that two tiles in an assembly interact, or are stably attached, if the glue labels on their abutting side are equal, and have positive strength. An assembly α induces a weighted binding graph G = (V, E), where V = dom(α), and there is an edge (a, b) ∈ E if and only if a and b interact, and this edge is weighted by the glue strength of that interaction. The assembly is said to be τ -stable if any cut of G has weight at least τ .

A tile assembly system is a triple T = (T, σ, τ ), where T is a finite tile set, σ is called the seed, and τ is the temperature. Throughout this paper, we will always have |dom(σ)| = 1 and τ = 1. Therefore, we can make the simplifying assumption that all glues have strength one without changing the behavior of the model.

Given two assemblies α and β, we write α → T 1 β if we can get β from α by the binding of a single tile, dom(α) ⊆ dom(β), and |dom(β) \ dom(α)| = 1. We say that γ is producible from α, and write α → T γ if there is a (possibly empty) sequence α = α 1 , . . . , α n = β such that

α 1 → T 1 . . . → T 1 α n . A sequence of k ∈ Z + ∪ {∞} assemblies α 0 , α 1 , . . . over A T is a T -assembly sequence if, for all 1 ≤ i < k, α i-1 → T 1 α i .
The productions of a tile assembly system T = (T, σ, τ ), written A[T ], is the set of all assemblies producible from σ. An assembly α is called terminal if there is no β such that α → T 1 β. The set of terminal assemblies is written A [T ].

An important fact about temperature 1 tile assembly, that we will use heavily, is that any path of the binding graph can grow immediately from the seed, independently from anything else.

Formally, a path P is a sequence of tile types along with positions, that is, of elements of T × Z 2 , such that no position occurs more than once in P , and for all i, the positions of P i and P i+1 are adjacent in the lattice grid of Z 2 .

For any path P and integer i, we write (x P i , y P i ) the coordinates of P i 's position. Moreover, if i < n, we say that the output side of P i is the side adjacent to P i+1 , and if i > 1, that its input side is the side adjacent to P i-1 . This means that the first tile of a path does not have an input side, and the last one does not have an output side. Remark that this definition of input/output sides is only relative to a path, and not to the tiles themselves; indeed, the tiles, including the first and last ones, may have other glues, not used by the path. Also, for any path P = P 0 , P 1 , . . . , P |P |-1 , and any integer i, we call the right-hand side (respectively left-hand side) of P i the side of that is between its input and output sides, in counterclockwise (respectively clockwise) order. When there is no ambiguity, we will say "the right-hand side (respectively left-hand side)" of P itself to mean "the right-hand side (respectively left-hand side) of some tile of P ".

Finally, for A, B ∈ Z 2 , we use the notation --→ AB to mean "the vector from A to B", and the

Manhattan distance between A = (x A , y A ) and B = (x B , y B ), written --→ AB 1 , is |x A -x B |+|y A -y B |.
We also call O the origin on Z 2 , i.e. the point of coordinates (0, 0).

The known upper bound

The only known way to assemble squares of size n × n at temperature 1 with 2n -1 tile types, is by using the "comb" design of Figure 1, already described in Rothemund and Winfree's paper [START_REF] Paul | The program-size complexity of self-assembled squares (extended abstract)[END_REF].

Building efficient paths

A major obstacle in proving the claimed lower bound for squares, is that building only a square's diagonal can require less tile types than building the whole square. In this section, we show the following result:

Theorem 3.1.
There is a tile assembly system T = (T, σ, 1), such that for all terminal assembly a ∈ A [T ], a is finite and of width 5(|T |+2)

4

-23.

Proof. Since a path of height n that is monotonic in the y-dimension, and has less than n tile types with input side south, can be "pumped", our path will need to have "caves", or non-monotonic subpaths, and reuse them several times. Moreover, in order for all the assemblies of T to be finite, the caves must be exited by a different path every time; however, since we want the caves to "save" tile types, these "exit paths" cannot all be new tile types. Therefore, one possible way to solve these constraints is to grow a regular monotonic path P 0 first, then build a cave C, and reuse a part of P 0 as its exit path. The next time we want to reuse C, we can use another part of P 0 as its exit path.

If the exit paths used in previous instance of a cave are all blocked in new instances, we will get only finite assemblies. Figure 2 is an example of such a path. This figure, however, is not a terminal assembly; since our caves, and their exit paths, are used several times, the same assemblies can grow from all their repetitions. Fortunately, we can arrange the shape of our main path so that no collision ever happens between these repetitions, like on Figure 3. Now, this figure has 38 tile types, and is of width 27; it does not yet save tile types. But by inserting:

• n new tile types in place of glue 6,

• n new tile types in place of glue 14,

• n new tile types in place of glue 24, and

• n new tile types in place of glue 26 Zooming in may be needed to read these numbers on Figures 2 and3. Printable versions are included, in Appendices A and B.

We add only 4n tile types, but the assembly is now 5n wider, and the result follows. An example of path that actually saves tile types is given on Figure 4. In this section, we prove that if a tileset has less than 2n-1 tile types and all its terminal assemblies are of domain {0, 1, . . . , n -1} 2 then one of its productions is a path that does not "save" tile types.

Our technique to prove this will be the following: assuming we are given such a tileset, we will first choose, using Algorithm 4.1, the assembly sequences that "lose" as much information as possible about their past, so as to "confuse" the tileset; then, if these assembly sequences can still build efficient paths, we will prove a necessary condition on these (this is Lemma 4.4), that cannot hold in completely filled squares (Theorem 4.5).

A path-building algorithm

We first define the algorithm we use to find assembly sequences that suit our purposes. There are two possible "priority modes" for this algorithm, namely right-priority and left-priority:

Algorithm 4.1.

Let T = (T, σ, 1) be a tile assembly system, be an ordering on T , and S be a non-empty subset of the set of paths producible by T . Let P 0 be the initial path, with just σ. For any i ≥ 0, if several tiles can bind to P i , let P i+1 be the one such that:

1. If P i is of the same type as a previous tile P h on P , with h < i, whose output side s is distinct from P i 's input side, then grow P h+1 , on side s of P i .

2. Else, if it is possible, let P i+1 be the tile such that:

(a) P i+1 binds to the first side of P i from its input side, in clockwise order if we are building a left-priority path, and in counterclockwise order if we are building a right-priority one. (b) There is at least one path in S, of which P 0 , P 1 , . . . , P i+1 is a prefix.

(c) If there are several such choices, we choose the smallest tile with respect to .

3. Finally, if none of the previous cases is possible to follow (for instance, because following case 1 resulted in a collision with a previous part of the assembly), but a new branch can still grow from the existing assembly, and produce a path of S, then follow it. Else, the algorithm halts.

We call the last tile that was placed before the present case the collision tile. In the special case that the current assembly before applying this case is already in S, we will adopt the convention that the empty branch can start.

Remark that this algorithm does not guarantee that the final path will be in S. If all producible paths end up being paths of S, it will be the case. But else, it means that we can "prevent" paths of S from growing. We will use this property heavily in the rest of our proof: S will be the set of all paths that reach some point in the square, and there is no way, in a filled square, to prevent a point from being covered.

Building filled squares from a seed in a corner

We begin by showing that a "path editing" operation is possible on the paths, built by Algorithm 4.1, that save tile types. Definition 4.2. A left-tentacular (respectively right-tentacular ) path is a path P such that the following conditions all hold:

• at least two tiles of P are of the same type. Let i and j be the indices of two tiles of P of the same type.

• P forked after a collision in case 3 of Algorithm 4.1, at some position k > j, and this branch started on the left-hand side (respectively right-hand side) of P k .

• This same branch can also start and completely grow from tile P k-j+i , after P is itself completely grown.

Moreover, we call --→ P j P i a contraction vector of P , and early branchings of P k,k+1,...|P |-1 are tentacles.

We say that a path is two-way tentacular if it is both left-and right-tentacular.

--→ P i P j Figure 5: A right-tentacular path, with the branch in blue, and the collision tile in red. A contraction vector is drawn on the right-hand drawing.

Definition 4.3. Let T = (T, σ, 1) be a tile assembly system. A path P ∈ A[T ] is said to be fragile if:

• a tile type is repeated on P , at positions i and j (i < j), and P j,j+1,...,|P |-1 , translated by --→ P j P i , can also start growing from P immediately after i.

• but P i,i+1,...,|P |-1 cannot grow completely from the resulting assembly.

In the rest of the proof, we will call breaking P , the choice to start growing P j,j+1,...|P |-1 , translated by --→ P j P i , immediately after i before growing P i,i+1,...,|P |-1 . Moreover, the smallest suffix of P that can branched earlier than its original first point on P is called the breaking branch.

See Figure 6 for an example of a fragile path (at this point, though, the colors of that figure

are not yet defined). Lemma 4.4. Let A = (x A , y A ) be a point of Z 2 , and T = (T, σ, 1) be a tile assembly system such that the following conditions all hold:

• dom(σ) = {O}.

• All assemblies can reach A, that is, ∀α ∈ A [T ], A ∈ dom(α).

• There is a path P ∈ A[T ] from O to A, built using Algorithm 4.1 with parameter S of the algorithm being the set of all producible paths of T from O to A, and P has strictly less than -→ OA 1 -1 tile types.

• P is such that for all i, y P i ∈ {0, . . . , y A }.

Then P is either fragile or two-way tentacular.

Proof. Let P be a path from O to A, built using Algorithm 4.1. Without loss of generality, we assume that y A ≥ 0 and x A ≥ 0 (we get other cases by rotating or flipping the argument).

We first introduce the idea of visible glues: we say that a glue between two tiles of P is visible from the east (respectively from the south) if P i interacts with P i+1 on its north face (respectively on its east face), and the glue between them is the rightmost (respectively lowest) one, on an infinite horizontal (respectively vertical) line between P i and P i+1 . Moreover, we adopt the convention that the rightmost tile on row y = y A has its north glue visible from the east, and its east glue visible from the south (even if this tile has 0-strength glues on these sides).

Figure 6: Tiles with glues visible from the east are in green or yellow; yellow tiles are of the same type, and red tiles are collisions. In this case, P is fragile.

Let us call P E the tiles of P that have their north glue visible from the east, and P S the tiles of P that have their east glue visible from the south. First, there are at least y A tiles in P E , and x A tiles in P S . We prove this now for P S : indeed, for all x ∈ {0, . . . , x A -1}, the visible tile of P on column x has output side east. To see this, suppose, for the sake of contradiction, that some tile P i is visible from the south, and has output side W. Then, draw an infinite horizontal half-line from the east side of P i to the south. Along with P 0,...,i and a horizontal infinite line at y = 0, it partitions the plane into two connected components, by Jordan's curve theorem (see Figure 7).

Therefore, because P does not intersect itself, and for all i, y P i ≥ 0, it must necessarily cross the vertical line at the east of P i again before reaching A. Thus, P i cannot be visible from the south, since P has at least one tile below P i .

A Figure 7: Jordan's curve theorem in action, with the enclosed region in blue

A similar argument, rotated by π/2, shows that P E has at least y A tiles. Furthermore, except possibly for the last tile of P , P E and P S are disjoint, because each tile has at most one output side: P E ∩ P S ⊆ {P |P |-1 }, and therefore |P E ∪ P S | ≥ -→ OA 1 -1. Now, because P was grown using Algorithm 4.1, if a tile type is repeated in P E , say at positions P i and P j , with i < j, then the translation of P i,i+1,...,j by --→ P i P j started growing immediately after P j , by case 1 of Algorithm 4.1, ultimately crashing into something, possibly after several repetitions (because we assumed that all productions of T are finite). We call this collision "collision C 0 ". Immediately after it happened, case 3 of the algorithm necessarily got used, and a new (possibly empty) branch was started from the existing assembly, and reached A (because of our hypothesis that all terminal assemblies have A in their domain).

Therefore, an early restart of this branch can also grow from P i,i+1,...,j . There are two cases:

1. If it can grow completely from P , it means that P is right-tentacular, by definition of righttentacular.

2. Else, this new branch crashes into something. We call this crash "collision C 1 ". There are two cases, according to where it crashes:

(a) Either it crashes into a part of P grown before P i . In this case, first observe that the north glue of P i is still visible from the east after this operation: indeed, the north glue of P j was visible from the east before, so no P k , for k > j, is on the right of P j on a horizontal line at y = y P j . And since we only translated these tiles by --→ P j P i , these translated tiles are not on the right of P i on a horizontal line at y = y P i either (see for instance Figure 6). We now prove that P is fragile. Indeed, grow it until index i. Then, grow a translation by --→ P j P i of the path that grew after P j , as long as it can possibly grow. We claim that the original P cannot grow anymore from this new assembly, that we call Q, and this is the definition of P being fragile. This is because of the definition of P E and P S : we first argue that collision C 1 encloses a portion of the plane: this is indeed an application of Jordan's curve theorem, because a collision between two connected paths closes a curve in the plane. from the south in Q, and the side of P i exposed to the south or the east is its right-hand side (this follows from the fact that P i 's output side is its north side if visible from the east, and its east side if visible from the south), the breaking branch starts from some P k (for k > j), on branches from its right-hand side. Therefore, if we grow this branch early, at position k -j +i, and it crashes into a previous part of the assembly, then by Jordan's curve theorem, we enclose the left-hand sides of its tiles; but since this branch forks from the right-hand side of the path, the enclosed region is thus the region in which P grew, because P is then a "left turn" from the breaking branch B. Therefore, since for all i, y P i ∈ {0, . . . , y A }, B cannot grow strictly higher than y A ; thus, P cannot reach A anymore from this assembly, which means that it cannot grow to its original endpoint, and therefore, it is fragile.

(b) Or it crashes into a part of P grown after P i , in which case we also choose to grow it before that part. Since P is a sequence of points, it will not be able to grow completely after that, since a tile will already be present next to the collision tile of collision C 1 .

Finally, the same argument on the tiles of P W and P N (the tiles with their north glue visible from the left, and their east glue visible from the north, respectively), proves that P is also lefttentacular or fragile. Remark that P W ∪ P N is not necessarily disjoint from P E ∪ P S (in the case where P has at most 2n -1 tile types, these two sets may even be equal).

Remark. Consider the "tentacles" of Figure 3, restarted from the path of Figure 2: The right tentacles are after glue 32, and then glue 36, and the left one after glue 35.

We can now prove the claimed result: Theorem 4.5. Let T = (T, σ, 1) be a tile assembly system whose terminal assemblies' domains are all {0, . . . , n -1} 2 , and such that dom(σ) = {(0, 0)}. Then:

|T | ≥ 2n -1
Proof. Let A, B, C and D be the following points of a square that T can assemble:

A B C D = O
Since T can fill this square, it must contain in particular a path P 0 from O to A, that we can build in the right-priority mode of Algorithm 4.1.

We define a sequence (S i ) i of assemblies, S 0 being the assembly where only P 0 has grown. For any i ≥ 0, let A i be the leftmost point on the right-hand side of S i , and not in S i . If there are several such points, let A i be the highest one. Now, let P i+1 be a right-priority path from O to A i , built using Algorithm 4.1. Moreover, let Q i+1 be the prefix of P i+1 that stops at the last occurrence of P i+1 's highest point.

There are three main cases:

1. If Q i+1 is fragile, we break it. Let R i+1 be the assembly resulting from that operation. In the case where Q i+1 shares parts with some other paths of S i , this operation may also break these paths. Let S i+1 be the assembly containing R i+1 and all the parts of S i that can grow from it.

2. If Q i+1 is right-tentacular, there are two subcases:

(a) One of its contraction vectors --→ v i+1 is to the left, i.e. x--→ v i+1 ≤ 0. In this case, we can do the same as in case 1 above. Indeed, since the last point of Q i+1 is next to the leftmost point on the right of S i , contracting it to the left will necessarily end in a collision with S i , before the end of Q i+1 . Indeed, this is clear in the case of Q i+1 = P i+1 , i.e. A i is the highest point of P i+1 . Else, P i+1 is longer than Q i+1 , and the subpath P i+1 \ Q i+1 grows after Q i+1 , by definition of Q i+1 as a prefix of P i+1 . Therefore, P i+1 \ Q i+1 is on the right of Q i+1 : indeed, by Jordan's curve theorem, on the closed curve pictured on Figure 8 (the curve defined by Q i+1 , two horizontal lines, immediately above and below Q i+1 , and a vertical line on the right of the square), P i+1 \ Q i+1 is inside a closed region of the plane, including a point on the right of Q i+1 0 . Therefore, A i is on the right of Q i+1 , because Q i+1 is part of the border of this closed zone.

P i+1 \ Q i+1 is in yellow.
Therefore, Q i is necessarily already in S i , because A i was chosen as the leftmost point on the right of S i . Thus, there is an h < i such that A h is the highest (and last) point of Q i+1 , and we can apply the same argument as in the case where Q i+1 = P i+1 , proving that P i+1 is fragile.

(b) Else, all its contraction vectors are strictly to the right. It is not possible that B ∈ dom(P i+1 ), for else B would also be in dom(Q i+1 ) (indeed, by definition of Q i+1 , it includes the last occurrence of the highest point on P i+1 ). and then one of the tentacles, when grown alone, would grow out of {0, . . . , n -1} 2 , contrarily to our hypothesis that all the assemblies stay within {0, . . . , n -1} 2 . Therefore, in this case, we let the tentacle grow as long as it can, until it reaches its highest point or else crashes into something. We continue this process with S i+1 being the union of S i , P i+1 , and the tentacles. (b) if A i = B, the construction is over. Indeed, in this case P i+1 = Q i+1 , and P i+1 is a path from O to B, with at least --→ OB 1 -1 = 2n -1 distinct tile types.

In order to conclude the proof, we need to argue that this construction halts, even though cases 1 and 2a seem to make the assembly smaller. First, all the paths we grow in this construction are right-priority paths. Therefore, after these two cases happen, the same points A i will appear again as A j , for j > i, but in this case, the path that we can build from O to A j will necessarily be "less right-priority" than the original ones, meaning that they will turn left earlier than P i . Indeed, in both case 1 and 2a, we break P i , and add the resulting assembly to S i .

To justify this, we examine what can possibly happen to these "broken paths" in greater detail now. When a path is broken, and we branch and grow a new one to reach A j , one of two things could happen:

1. Either P j is fragile, in which case nothing can happen to earlier broken parts. The assembly sequence that we use to prove this is the following: First grow all the parts that can grow on the left of P j , then the broken P j . Then, grow the earlier broken parts on its right. Either these parts are on the left of P j 's breaking branch, in which case they are also enclosed, or they are on its right, but in this case, even if a collision happens between the breaking branch of P j , and these parts, this collision is either with the original fragile path -which still keeps it broken -or with an earlier breaking branch -which still keeps the original path enclosed, and therefore broken.

2. Else, if P j is right-tentacular, the construction will continue. If P j 's tentacles crash into something, then this crash does not affect any of the arguments, since we only let them grow until their first collision. Proof. Let Σ be the position of the seed. We can use the technique of Theorem 4.2 on all the rectangles with diagonals (ΣA), (ΣB), (ΣC) and (ΣD). At least one is of width and height at least n/2, and verify the assumptions of Lemma 4.4, and the proof of Theorem 4.5 can be applied, either directly, or by "turning in the other direction", i.e. looking for a non-left-tentacular path.

Future work

The next step, in proving the fully general conjecture, is to extend Lemma 4.4 to the case where the seed can be anywhere in the square. The reason why it does not apply to that case, is that a "lower restart" is not necessarily an "early restart"; indeed, let P be a path from O to D, and Q be a path from O to B. If P i and Q j are of the same type for some i and j, then restarting P i,i+1,...,|P |-1 from Q j does not result in a competition for growth between this "lower branch" and Q, since the tiles of Q, to the north of Q j , were grown before Q j . Therefore, despite Corollary 4.6, the general question remains open:

Open Problem 1. Is there a tile assembly system with less than 2n -1 tile types, that can assemble a filled square at temperature 1, starting from a single tile anywhere in the square?

Even though the recent results about its intrinsic universality [START_REF] Meunier | Intrinsic universality in tile self-assembly requires cooperation[END_REF], and the present paper, have made significant advances in that direction, the exact computational power of temperature 1 systems is still completely unknown. Moreover, the existence of single tileset simulating, at temperature 1, any other temperature 1 tile assembly system, is still open. The following open problem is particularly puzzling, especially in regard of the impressive results of [START_REF] Cook | Temperature 1 self-assembly: deterministic assembly in 3D and probabilistic assembly in 2D[END_REF], showing that threedimensional temperature 1 tile assembly are capable of Turing computation: Open Problem 2. Is is decidable whether two tile assembly systems have the same terminal assemblies? Moreover, Theorem 3.1 is the first two-dimensional construction at temperature 1, with less tile types than its Manhattan diameter, in the general aTAM. On the other hand, our solution to the original conjecture relies heavily on the fact that our squares need to be completely filled. This leaves the following question open:

Open Problem 3. all n, is there a tile assembly system T n = (T n , σ n , 1), with |dom(σ n )| = 1 and T n < 2n -1, whose terminal assemblies all contain at least a square frame, that is, such that for all α ∈ A [T n ], ({0, n -1} × {0, . . . , n -1}) ∪ ({0, . . . , n -1} × {0, n -1}) ⊆ dom(α)?

The initial construction of a square with 2n -1 tile types, by Rothemund and Winfree, used a fairly simple design. Our result shows that this is optimal, but it does not discard the possibility of other designs:

Open Problem 4. Is there a way to self-assemble a square of size n × n at temperature 1 with 2n -1 tile types, that is not a trivial variation of Rothemund and Winfree's "comb" design (Figure 1)? How many are there?

Figure 1 :

 1 Figure 1: The comb design for size 10 × 10

Figure 3 :Figure 4 :

 34 Figure 3: Early branches, not blocking the main assembly 54 53 55 54 56 55 57 56 58 57 59 58 60 59 61 60 62 61 63 62 64 63 65 64 66 65 67 66 68 67 69 68 70 69 71 70 72 71 73 72 74 73 75 74 75

Figure 8 :

 8 Figure 8: Jordan's curve theorem. The enclosed region is colored blue, Q i+1 is in grey, and P i+1 \ Q i+1 is in yellow.

3 .

 3 Else, by Lemma 4.4, Q i+1 has at least OA i 1 -1 tile types. There are two subcases: (a) if A i = B, we let S i+1 = S i ∪ P i+1 , and resume the construction.

Corollary 4 . 6 .

 46 Let T = (T, σ, 1) be a tile assembly system whose terminal assemblies are all of domain {0, . . . , n -1} 2 , and such that |dom(σ)| = 1. Then |T | ≥ n.

Acknowledgments

I thank David Doty and Damien Woods for friendly support and discussions, Erik Winfree and his group for hosting me during this research, and Paul Rothemund for encouragements. I am especially grateful to Damien Woods, for all the comments that helped me unshuffle a previous, particularly tortured, version of this paper, and more generally for taking me with him on the journey [29]. Section 1 of this paper is partially taken from [14], I would like to thank my co-authors on that paper. Also, Tom Hirschowitz and Christophe Raffalli suggested the name "tentacular".

Supported in part by National Science Foundation Grant CCF-1219274.


Andreas Enge
email: andreas.enge@inria.fr

Emmanuel Thomé
email: emmanuel.thome@inria.fr

Computing class polynomials for abelian surfaces

Keywords:

published or not. The documents may come

Introduction

Igusa class polynomials describe the complex multiplication points in the moduli space of principally polarised abelian surfaces, that is, they parameterise abelian varieties of dimension 2 with complex multiplication by a maximal order of a quartic CM field. Such abelian surfaces are Jacobians of hyperelliptic curves of genus 2, so that by computing Igusa class polynomials one may obtain genus 2 curves over finite fields with known Jacobian cardinality.

In the dimension 1 case of elliptic curves, several approaches have been described in the literature. While the output of the algorithms (a large polynomial) is of exponential size in the input (a number field described by a single integer), all of these approaches may lead to an algorithm with a complexity that is quasi-linear (up to logarithmic factors) in its output size: The complex analytic method uses floating point approximations to the roots of the class polynomials [START_REF] Enge | The complexity of class polynomial computation via floating point approximations[END_REF]; the p-adic approach starts from a curve with the given endomorphism ring over a small finite field and lifts its invariants to a p-adic field [START_REF] Couveignes | Action of modular correspondences around CM points[END_REF][START_REF] Bröker | A p-adic algorithm to compute the Hilbert class polynomial[END_REF]; the Chinese remaindering approach combines curves over several small prime fields [START_REF] Belding | Computing hilbert class polynomials[END_REF].

In principle, the same approaches apply to abelian surfaces. A 2-adic algorithm is described in [START_REF] Gaudry | The 2-adic CM method for genus 2 curves with application to cryptography[END_REF], and there are currently attempts at making the Chinese remainder based method more efficient [START_REF] Lauter | Improved CRT algorithm for class polynomials in genus 2[END_REF]. So far, 2-adic lifting appears to have been the most successful approach: The Echidna database maintained by Kohel1 contains Igusa class polynomials, the largest of which is of degree 576 and has been obtained by lifting from a curve over F 2 6 2 .

A detailed description of the complex analytic approach, together with complexity analyses of its different steps, has recently been given in [START_REF] Streng | Computing Igusa class polynomials[END_REF][START_REF] Streng | Complex multiplication of abelian surfaces[END_REF]. Our work pushes the limits for the attainable degrees of Igusa class polynomials: We present an example of degree 20 016. Moreover, relatively small class polynomials (say, below degree 150) can be computed in matters of seconds. The key tool in this approach is the use of a quasi-linear algorithm for the computation of ϑ-constants, initially described in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF].

This article is organised as follows. §2 presents the necessary background material for discussing the complex multiplication theory of abelian surfaces and states the general algorithm. § §3 and 4 show how to explicitly (providing concrete descriptions for the occurring number fields, maps between them and their embeddings) and symbolically compute an appropriate set of reduced period matrices, which form the input of the computationally expensive step of computing ϑ-constants, detailed in §5. The recognition as algebraic numbers of the coefficients of the Igusa class polynomials from their approximations by complex embeddings is described in §6, and experimental results are given in § §7 and 8.

All computations presented in this article have been achieved with the software package Cmh [START_REF] Enge | Cmh -Computation of Igusa Class Polynomials[END_REF], released under the GNU General Public License.

Complex multiplication theory

In this section, we provide a concise introduction to the theory of complex multiplication of principally polarised abelian surfaces or, equivalently, Jacobians of genus 2 hyperelliptic curves over the complex numbers, to the extent needed to describe our algorithms and implementation. The presentation follows [START_REF] Streng | Complex multiplication of abelian surfaces[END_REF], and proofs are given in [START_REF] Shimura | Complex Multiplication of Abelian Varieties and its Applications to Number Theory[END_REF][START_REF] Shimura | Abelian Varieties with Complex Multiplication and Modular Functions[END_REF][START_REF] Streng | Complex multiplication of abelian surfaces[END_REF][START_REF] Streng | Computing Igusa class polynomials[END_REF].

Quartic CM fields and abelian surfaces

A CM field K is an imaginary-quadratic extension of a totally real number field K 0 . We denote by κ indiscriminately the complex conjugation on C and the automorphism generating Gal(K/K 0 ). For any embedding ϕ : K → C, we have κ • ϕ = ϕ • κ, which justifies the notation ϕ = κ • ϕ.

Quartic CM fields K of degree 4 over Q come in three Galois types. Generically, K/Q is not Galois, the Galois closure L/K is of degree 2, and Gal(L/Q) is isomorphic to the dihedral group D 4 . The Galois closure L is itself a CM field, and the complex conjugation of L, which we denote again by κ, restricts to the complex conjugation of K. If K/Q is Galois, it may be either cyclic or biquadratic. We will not consider the biquadratic case in the following, since then the abelian surfaces of which it is the endomorphism algebra are products of elliptic curves; so from now on, all Galois quartic CM fields are tacitly understood to be cyclic.

A CM type of a quartic CM field K is a set Φ = {ϕ 1 , ϕ 2 } of two embeddings K → C such that ϕ 2 = ϕ 1 ; that is, it contains one out of each pair of complex-conjugate embeddings. Two CM types Φ and Φ ′ are equivalent if there is an automorphism σ of K such that Φ ′ = Φ • σ; in particular, Φ and Φ are equivalent. If K/Q is Galois, there is only one equivalence class of CM types; otherwise, there are two inequivalent classes

Φ = {ϕ 1 , ϕ 2 } and Φ ′ = {ϕ 1 , ϕ 2 }.
For a given CM type Φ = {ϕ 1 , ϕ 2 }, its reflex field is the field K r generated over Q by the type traces, that is, K r = Q {ϕ 1 (x) + ϕ 2 (x) : x ∈ K} ; it is itself a quartic CM field and we denote by K r 0 its real-quadratic subfield. Equivalent CM types yield conjugate reflex fields. In the Galois case, K and K r are isomorphic, while in the dihedral case, they are not isomorphic, but the two reflex fields for the two inequivalent CM types are. In both cases, there is a natural way of defining a dual CM type Φ r = {ϕ r 1 , ϕ r 2 } of K r , and the reflex field of K r is isomorphic to K. Define the (dual) type norm

N Φ r : K r → K by x → ϕ r 1 (x)ϕ r 2 (x), so that N Φ r N Φ r = N; (1)
this map extends to ideals and ideal classes.

In §3, we provide explicit equations for all occurring number fields and consider their embeddings from an effective point of view.

Let a be a fractional ideal of

O K . A CM type Φ = {ϕ 1 , ϕ 2 } induces an embedding K → C 2 , x → (ϕ 1 (x), ϕ 2 (x)
), under which Φ(a) is a lattice of rank 4. Its cokernel C 2 /Φ(a), a complex torus of genus 2, is an abelian surface. Let δ -1 K = {y ∈ K : Tr(xy) ∈ Z ∀x ∈ O K } be the codifferent ideal of K. Assume that (aaδ K ) -1 is principal and generated by some ξ ∈ K such that ϕ 1 (ξ), ϕ 2 (ξ) ∈ iR >0 ; in particular, ξξ ∈ K 0 is totally negative. Then E Φ,ξ : Φ(K) 2 → Q, (Φ(x), Φ(y)) → Tr(ξxy) is a symplectic form over Q which takes integral values on Φ(a) 2 . By tensoring with R, one obtains a symplectic form C 2 → R such that (x, y) → E Φ,ξ (ix, y) is symmetric and positive definite, a principal polarisation on C 2 /Φ(a).

The principally polarised abelian surface A(Φ, a, ξ) = C 2 /Φ(a), E Φ,ξ has complex multiplication by O K ; conversely, any such surface can be obtained up to isomorphism in this way. Two principally polarised abelian surfaces A(Φ, a, ξ) and A(Φ ′ , a ′ , ξ ′ ) are isomorphic if and only if Φ = Φ ′ (up to equivalence) and there is a u ∈ K * such that a ′ = ua and ξ ′ = (uu) -1 ξ. In particular this implies that uu ∈ K 0 is totally positive, and that we may assume a to be an integral ideal of O K .

The Shimura group, its type norm subgroup and cosets

The Igusa invariants to be defined in §2.3 determine the moduli space M of principally polarised complex abelian surfaces, which has a model over Q. Let M K,Φ be the subset of surfaces A(Φ, a, ξ) obtained from an integral ideal of O K and the CM type Φ as

described in §2.1. Then M K,Φ is stable under Gal(Q/K r 0 ). If K is cyclic, then M K,Φ is even stable under Gal(Q/Q).
Otherwise let Φ ′ be inequivalent with Φ. Then M K,Φ and M K,Φ ′ are disjoint and conjugate under Gal(K r 0 /Q) [34, Lemmata 1.1 and 2.1]. Let the Shimura class group C be defined by C = (a, u) : a a fractional ideal of O K , aa = uO K , and u ∈ K 0 totally positive}/ ∼

(2) with component-wise multiplication. The equivalence relation denoted ∼ above is the one induced by principal ideals, more precisely the equivalence modulo the subgroup given by the (vO K , vv) with v ∈ K * and vv ∈ K 0 totally positive.

By the discussion of §2.1, the Shimura class group C acts regularly on

M K,Φ via (b, u) • A(Φ, a, ξ) = A(Φ, b -1 a, uξ). ( 3
)
Consider the dual type norm map N Φ r :

Cl K r → C, b → (N Φ r (b), N(b))
, which is well defined by (1). For any A(Φ, a, ξ), the action induced by N Φ r (Cl K r ) is that of the Galois group of the field of moduli of A(Φ, a, ξ) over K r [34, Theorem 9.1]; otherwise said, the field of moduli is the fixed field of ker(N Φ r ) inside the Hilbert class field of K r . The cokernel of N Φ r is elementary abelian of exponent 1 or 2 [34, Theorem 2.2], so M K,Φ splits into orbits under C of size | im(N Φ r )|, and the number of orbits is a power of 2.

As stated above, these orbits are in fact defined over K r 0 , with the orbits of M K,Φ and M K,Φ ′ being mapped to each other by Gal(K r 0 /Q).

ϑ-functions, Igusa invariants and class polynomials

Given an ideal a and a principal polarisation E Φ,ξ as in §2.1, one may choose a Z-basis

(α 1 , α 2 , α 3 , α 4 ) of a such that v 1 = Φ(α 1 ), v 2 = Φ(α 2 ), w 1 = Φ(α 3 ), w 2 = Φ(α 4 ) form a symplectic basis, for which E Φ,ξ becomes 0 id 2 -id 2 0
. That the change of basis is defined over Z and not only over R follows from the principality of the polarisation; we also call this basis of a symplectic.

Let V = v 1 v 2 , W = w 1 w 2 ∈ C 2×2 .
Rewriting the ambient vector space C 2 and Φ(a) in the basis spanned by w 1 and w 2 , we obtain Φ(a) = Ω Φ,a,ξ id 2 Z 4 with the period matrix

Ω Φ,a,ξ = W -1 V (4) in the Siegel half space H 2 = Ω ∈ C 2×2 : Ω symmetric
(Z). Let a = a 1 a 2 , b = b 1 b 2 ∈ 1 2 Z 2 be two vectors of ϑ-characteristics. Then for Ω ∈ H 2 , ϑ 16a 1 +8a 2 +4b 1 +2b 2 (Ω) = ϑ a,b (Ω) = n∈Z 2 e 2πi( 1 2 (n+a) ⊺ Ω(n+a)+(n+a) ⊺ b) . ( 5
)
Only the even ϑ-constants ϑ i for i ∈ T = {0, 1, 2, 3, 4, 6, 8, 9, 12, 15} are not identically 0.

The following duplication formulae relate the values of the squares of the ten even ϑ-constants in the argument Ω with the values of the four fundamental ϑ-constants ϑ 0 , . . . , ϑ 3 (which have a = 0) in the argument Ω/2 (omitted from the formulae for the sake of conciseness).

4ϑ 2 0 (Ω) = ϑ 2 0 + ϑ 2 1 + ϑ 2 2 + ϑ 2 3 4ϑ 2 1 (Ω) = 2ϑ 0 ϑ 1 + 2ϑ 2 ϑ 3 4ϑ 2 2 (Ω) = 2ϑ 0 ϑ 2 + 2ϑ 1 ϑ 3 4ϑ 2 3 (Ω) = 2ϑ 0 ϑ 3 + 2ϑ 1 ϑ 2 4ϑ 2 4 (Ω) = ϑ 2 0 -ϑ 2 1 + ϑ 2 2 -ϑ 2 3 4ϑ 2 6 (Ω) = 2ϑ 0 ϑ 2 -2ϑ 1 ϑ 3 4ϑ 2 8 (Ω) = ϑ 2 0 + ϑ 2 1 -ϑ 2 2 -ϑ 2 3 4ϑ 2 9 (Ω) = 2ϑ 0 ϑ 1 -2ϑ 2 ϑ 3 4ϑ 2 12 (Ω) = ϑ 2 0 -ϑ 2 1 -ϑ 2 2 + ϑ 2 3 4ϑ 2 15 (Ω) = 2ϑ 0 ϑ 3 -2ϑ 1 ϑ 2 (6)
Denote by h j the following modular forms of weight j:

h 4 = i∈T ϑ 8 i , h 6 = 60 triples (i,j,k)∈T 3 ±(ϑ i ϑ j ϑ k )
j 1 = h 4 h 6 h 10 , j 2 = h 2 4 h 12 h 2 10 , j 3 = h 5 4 h 2 10 . ( 8
)
The principally polarised abelian surfaces A(Φ, a, ξ) are parameterised by the triples of singular values (j 1 (Ω), j 2 (Ω), j 3 (Ω)) in the period matrices Ω = Ω Φ,a,ξ , which may be obtained from the action of the Shimura class group C on a fixed base point β = (Φ, a Φ , ξ Φ ). The singular values lie in the subfield of the Hilbert class field of K r given in §2.2. Following the discussion there, the Igusa class polynomials I i (X) = (Φ,a,ξ) (X -j i (Ω Φ,a,ξ )) are defined over Q. More precisely their irreducible factors, over K r 0 in the dihedral case or Q in the cyclic case, are given by

C∈N Φ r (Cl K r ) X -j i (Ω CC ′ •β ) ,
where Φ is one CM type and

C ′ ∈ C/ N Φ r (Cl K r ).
In the following, we fix a CM type Φ (for its explicit description, see §3) and a base point β = (Φ, a Φ , ξ Φ ) and let

H 1 (X) = C∈N Φ r (Cl K r ) (X -j 1 (Ω C•β )) . ( 9
)
As elements of the same class field, the singular values of j 2 and j 3 are rational expressions in the singular value of j 1 . Computationally, it is preferable to use the Hecke representation in the trace-dual basis to keep denominators small. We thus define polynomials Ĥ2 and Ĥ3 through

j i H ′ 1 (j 1 ) = Ĥi (j 1 ) with Ĥi (X) = C∈N Φ r (Cl K r ) j i (Ω C•β ) D∈N Φ r (Cl K r )\{C} (X -j 1 (Ω D•β )) ( 10
)
for i ∈ {2, 3}, where H 1 , Ĥ2 , Ĥ3 ∈ K r 0 [X] in the dihedral case and ∈ Q[X] in the cyclic case.

Algorithm for Igusa class polynomials

We briefly summarise the algorithm for computing class polynomials.

Algorithm 1

Input: CM field K and CM type

Φ = {ϕ 1 , ϕ 2 } of K Output: Irreducible class polynomials H 1 , Ĥ2 , Ĥ3 ∈ K r 0 [X] in the dihedral case and ∈ Q[X] in the Galois case 1) Compute N Φ r (Cl K r ) = {(b 1 , u 1 ), . . . , (b h , u h )} ⊆ C. 2) Compute a base point β = (Φ, a Φ , ξ Φ ) such that (a Φ a Φ δ K ) -1 = (ξ Φ ), ϕ 1 (ξ Φ ), ϕ 2 (ξ Φ ) ∈ iR >0 . 3) Enumerate {C • β = (Φ, b -1 i a Φ , u i ξ Φ ) : C = (b i , u i ) ∈ N Φ r (Cl K r )} and compute the associated period matrices Ω i = Ω C•β for i = 1, . . . , h. 4) For i = 1, . . . , h, compute the fundamental ϑ-constants ϑ 0 (Ω i /2), . . . , ϑ 3 (Ω i /2);
then deduce the squares of the ten even ϑ-constants ϑ 2 k (Ω i ) by ( 6), the values h k (Ω i ) by [START_REF] Brent | Fast multiple-precision evaluation of elementary functions[END_REF] and finally the triples

J i = j 1 (Ω i ), j 2 (Ω i ), j 3 (Ω i ) by (8). 5) Let H 1 = h i=1 (X -J i,1 ), Ĥk = h i=1 J i,k l =i (X -J l,1 ) ∈ C[X] for k ∈ {2, 3}. 6) Recognise the coefficients of H 1 , Ĥ2 , Ĥ3 as elements of K r 0 or Q, respectively.
The different steps of the algorithm and our implementation are detailed in the following chapters. The symbolic computations related to number fields in Steps 1) and 2) and to the period matrices Ω i in Step 3) are described in §3. Step 1) is treated in §4.2, Step 4) in §5 and Step 6) in §6.

Explicit equations and symbolic period matrices

While Algorithm 1 in fine works with complex approximations obtained via CM types, it starts from an algebraic setting. In this section, we examine how to carry out the computations as far as possible symbolically with algebraic numbers, which relieves us from the need to decide on the necessary precision early on. In particular, in §3.1 we replace the complex embeddings forming a CM type by algebraic embeddings into the compositum L of all involved fields, followed by a "universal" embedding ψ of L into C. Taking preimages under ψ, the entries of the period matrices Ω ∈ C 2×2 may then be interpreted as elements of the reflex field and may be handled symbolically. We then fix a model for the CM field K in §3.2 and derive explicit equations for all considered fields and embeddings.

Recall the notation of §2: K is a quartic CM field, K 0 its real quadratic subfield and L its Galois closure with Galois group G. We consider only the dihedral case [L : K] = 2 and G = D 4 and the cyclic case L = K and G = C 4 . Let Φ = (ϕ 1 , ϕ 2 ) be a CM type, where ϕ 1 , ϕ 2 : K → C are two complex embeddings of K with ϕ 2 = ϕ 1 , and let K r be the reflex field of K with respect to Φ.

Galois theory, embeddings and period matrices

The dihedral case

Galois theory. Let K = Q(y) be a non-Galois quartic CM field. The following statements are easily seen to be true when choosing a generator y such that z = y 2 belongs to the real subfield K 0 , so that

K = Q[Y ]/(Y 4 + AY 2 + B) for some A, B ∈ Q. The Galois closure of K is then L = K(y ′ ) = Q(y, y ′ )
, where the roots of the minimal polynomial of y in L are ±y and ±y ′ (the former could be identified with ϕ 1 (y) in ( 12), the latter with ϕ 2 (y)). The automorphisms in G = Gal(L/Q) are uniquely determined by their images on y and y ′ , and we obtain the following diagram of fields and Galois groups:

Q K 0 K K r 0 K r * L = KK r κ| K ρ κ| K r σ κ
Here the automorphisms are given by ρ : y → y, y ′ → -y ′ of order 2 σ : y → y ′ , y ′ → y of order 2, which fixes the generator y + y ′ of K r τ : y → y ′ , y ′ → -y of order 4 κ = τ 2 : y → -y, y ′ → -y is the complex conjugation.

So G is the dihedral group D 4 with generators τ of order 4 and ρ (or σ) of order 2 and additional relation ρτ ρ = τ 3 , and with κ as its centre.

Embeddings and CM types.

There is a unique embedding ψ : L → C such that ϕ 1 = ψ| K and ϕ 2 = (ψσ)| K (where multiplication denotes composition), which can be seen as follows. First of all, there are two embeddings which, restricted to K, yield ϕ 1 ; we denote them by ψ 1 and ψ ′ 1 = ψ 1 ρ. Now there is s ∈ G, uniquely defined up to multiplication by ρ from the right, such that

ϕ 2 = (ψ 1 s)| K . Since ϕ 2 = ϕ 1 and ϕ 2 = ϕ 1 , the automorphism s is neither 1, ρ, κ = τ 2 nor κρ = τ 2 ρ. This leaves s as one of τ = ρσ, τ ρ = ρσρ, τ 3 = σρ or τ 3 ρ = σ. If s| K = σ| K = (σρ)| K , we may choose ψ = ψ 1 . Otherwise, s| K = ρσ, and (ψ ′ 1 σ)| K = (ψ 1 ρσ)| K = (ψ 1 s)| K = ϕ 2 , so we choose ψ = ψ ′ 1 .
Period matrices. Let (α 1 , . . . , α 4 ) be a symplectic basis for the ideal a of K with respect to E Φ,ξ as defined in §2.3. Then

V = ϕ 1 (α 1 ) ϕ 1 (α 2 ) ϕ 2 (α 1 ) ϕ 2 (α 2 ) = ψ α 1 α 2 α σ 1 α σ 2 , W = ϕ 1 (α 3 ) ϕ 1 (α 4 ) ϕ 2 (α 3 ) ϕ 2 (α 4 ) = ψ α 3 α 4 α σ 3 α σ 4 and Ω Φ,a,ξ = W -1 V = ψ(M ) with M = 1 α 3 α σ 4 -α 4 α σ 3 α 4 α σ 1 -α 1 α σ 4 α 4 α σ 2 -α 2 α σ 4 α 3 α σ 1 -α 1 α σ 3 α 3 α σ 2 -α 3 α σ 2 (11
) by ( 4). The entries of M are invariant under σ and thus elements of K r .

Remark. It is crucial to choose out of the two embeddings ψ : L → C that extend ϕ 1 the one compatible with ϕ 2 . The other one corresponds to the second CM type Φ ′ = (ϕ 1 , ϕ 2 ) with reflex field (K r ) ′ and Gal (L/(K r ) ′ ) = κσ = ρσρ .

The cyclic case

Here we have the much simpler situation

Q K 0 K κ = σ 2
We may choose ψ = ϕ 1 . Then there is a uniquely determined σ ∈ Gal(K/Q) such that ϕ 2 = ϕ 1 σ, and trivially M of (11) has entries in K r . In general, they will not lie in a subfield: Since σ is neither the identity nor complex conjugation, it is of order 4.

Number field computations

In this section we show how to express the elements of the reflex field K r and the normal closure L in consistent ways, so as to be able to compute type norms and entries of period matrices as given by ( 4). We use the same notation for elements of the Galois group G of L/Q as in §3.1.

The dihedral case

Field equations. By choosing generating elements as in §3.1.1 we may assume that

K 0 = Q(z) = Q[Z]/ Z 2 + AZ + B with A, B ∈ Z >0 , A 2 -4B > 0; K = Q(y) = Q[Y ]/ Y 4 + AY 2 + B .
We then select the CM type Φ = (ϕ 1 , ϕ 2 ) with

ϕ 1 (y) = i A + √ A 2 -4B 2 , ϕ 2 (y) = i A - √ A 2 -4B 2 , ( 12
)
where all the real roots are taken to be positive; the other CM type is Φ ′ = (ϕ 1 , ϕ 2 ) with ϕ 2 (y) = -ϕ 2 (y). Recall from §3.1.1 the notations Gal(L/K) = ρ , Gal(L/K r ) = σ , and let ψ : L → C be such that ϕ 1 = ψ| K and ϕ 2 = (ψσ)| K . The reflex field K r is generated by the type traces of K; letting y r = y + y σ , the equality

ψ(y r ) = ψ(y) + (ψσ)(y) = ϕ 1 (y) + ϕ 2 (y) (13)
shows that we may consider y r as a generator of K r . This gives the equations

K r 0 = Q(z r ) = Q[Z r ]/ (Z r ) 2 + A r Z r + B r with A r = 2A, B r = A 2 -4B; K r = Q(y r ) = Q[Y r ]/ (Y r ) 4 + A r (Y r ) 2 + B r .
The minimal polynomials of y r over K and y over K r follow:

(y r ) 2 -2yy r + (2y 2 + A), (y) 2 -y r y + ((y r ) 2 + A)/2.
We write the Galois closure L = KK r as the compositum generated over K or K r by t = y + y r . The minimal polynomial of t is the resultant

h(T ) = Res Y Y 4 + AY 2 + B, (T -Y ) 2 -2Y (T -Y ) + (2Y 2 + A) = Res Y r (Y r ) 4 + A r (Y r ) 2 + B r , (T -Y r ) 2 -Y r (T -Y r ) + ((y r ) 2 + A)/2 = T 8 + 10AT 6 + (33A 2 -14B)T 4 + (40A 3 -70AB)T 2 + 16A 4 -200A 2 B + 625B 2 .
Conversions and Galois actions. We are interested in the action of ρ, the generator of Gal(L/K), on K r , and in the action of σ, the generator of Gal(L/K r ), on K. The defining equations give:

y r + (y r ) ρ = 2y, y r (y r ) ρ = y r (y r ) ρ = 2y 2 + A, y ρ = y, y + y σ = y r , yy σ = (y r ) 2 + A /2, (y r ) σ = y r .
An element of K is converted to an element of L, as a relative extension of K r , using the identity y = t -y r ; in the opposite direction we use y r = t -y. The entries of the matrix M of ( 11) are obtained from elements of K and their images under σ, and need to be expressed as elements of K r . For this we use the identity y σ = y r -y. This allows us to work in the relative extension L/K r and to easily identify elements of K r .

Dual type norms.

For an ideal b of K r , we have

N Φ r (b) = N L/K (bO L ), see [9, §3.1].
Computing dual type norms thus reduces to conversions in relative extensions as described above.

The cyclic case

We may use the same type of equations for K and K 0 as in the dihedral case, and may fix ψ = ϕ 1 as in [START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF]. Fixing an arbitrary element σ ∈ Gal(K/Q) of order 4, we obtain ϕ 2 = ϕ 1 σ. Then the dual type norm for an ideal b of K is computed as

N Φ r (b) = bb σ , see [9, §3.1].

Symbolic reduction of period matrices

Gottschling in [START_REF] Gottschling | Explizite Bestimmung der Randflächen des Fundamentalbereiches der Modulgruppe zweiten Grades[END_REF] for which the condition is violated). The process needs to be iterated and terminates eventually.

In the light of [START_REF] Cohen | A course in algorithmic algebraic number theory[END_REF], Ω = ψ(M ) with M ∈ (K r ) 2×2 and an explicitly given ψ : K r → C, see ( 13) and [START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF]. Letting K r = K r 0 + y r K r 0 as before, we have ψ| K r 0 : Hence the period matrices may be transformed symbolically into the fundamental domain F 2 without computing complex approximations of their entries, which precludes rounding errors: The test whether the matrix is reduced and, if not, the decision which transformation to apply depend on the sign of ψ(α) for some α ∈ K r 0 , that is, on the sign of some explicitly known a The first step of Algorithm 1 requires to enumerate the Shimura group C of (2), or more precisely, its type norm subgroup N Φ r (Cl K r ). We need the following exact sequence, a proof of which can be found in [START_REF] Bröker | Explicit CM theory for level 2-structures on abelian surfaces[END_REF]:

K r 0 = Q( √ D r ) → R and ψ| y r K r 0 : y r K r 0 → iR. So ℜ(M )
+ b √ D r ∈ R,
1 -→ O + K 0 /N K/K 0 (O * K ) u →(O K ,u) -------→ C (a,α) →a -----→ Cl K N K/K 0 ----→ Cl + K 0 -→ 1, (14)
where O + K 0 is the subgroup of totally positive units in O K 0 and Cl + K 0 is the narrow class group of K 0 .

We have algorithms at hand for the basic arithmetic of C. For a finite abelian group, decomposed as a direct product of cyclic groups G i of order d i with d i | d i+1 , we call the d i the elementary divisors and a system of generators of the G i a (cyclic) basis of the group. Such a basis can be computed for the class group Cl K (quickly under GRH) using the function bnfinit in Pari/Gp. Equality testing of (a, α) and (b, β) amounts to testing whether ab -1 is principal (either using bnfisprincipal in Pari/Gp, or by a direct comparison if each ideal is stored together with its generalised discrete logarithm, its coefficient vector with respect to the basis of the class group), and whether

α/β = 1 in O + K 0 / N K/K 0 (O * K ).
Let ε 0 and ε be the fundamental units of K 0 and K, respectively.

If N(ε 0 ) = -1, then O + K 0 = ε 2 0 = N K/K 0 ( ε 0 ) ⊆ N K/K 0 (O * K )
, and the quotient group

is trivial. If N(ε 0 ) = +1, then O + K 0 = ε 0 , and since ε 2 0 = N K/K 0 (ε 0 ) ∈ N K/K 0 (O * K )
, the quotient group is either trivial or ε 0 / ε 2 0 , in which case bnfisunit of Pari/Gp can be used to compute the exponent of the unit.

Multiplication is straightforward and can be made more efficient by a reduction step that outputs a smaller (not necessarily unique) representative. To reduce (a, α), one computes an LLL-reduced ideal a ′ = µa (using idealred in Pari/Gp) and lets α ′ = µµα One then tries to reduce the unit contribution in the size of the algebraic number α ′ by multiplying it with an appropriate power of N K/K 0 (ε).

The Shimura group C and its subgroup N Φ r (Cl Kr ) can be enumerated directly; but the map N Φ r : Cl Kr → C being in general non-injective, this can require a large number of expensive principality tests in C to avoid duplicates. More elegantly, we may consider the groups in [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF] as given by cyclic bases or, more generally, generators and relations, and complete the sequence from known data using tools of linear algebra for Z-modules, in particular the Hermite (HNF) and Smith normal forms (SNF), see [11, §2.4].

Algorithm 2

Input: Cyclic bases for Cl K and Cl

+ K 0 Output: Cyclic basis for C 1) Compute a matrix M for N K/K 0 : Cl K → Cl + K 0 . 2) Compute generators a 1 , . . . , a r of the kernel of M .
3) Lift a 1 , . . . , a r to C: Pick arbitrary totally positive

α i ∈ K 0 such that a i a i = α i O K 0 .
4) Compute a basis for the lattice L 0 of relations such that the subgroup of Cl K generated by a 1 , . . . , a r is isomorphic to Z r /L 0 .

5) If O + K 0 / N K/K 0 (O * ) = 1, let r ′ = r; otherwise let r ′ = r+1 and (a r ′ , α r ′ ) = (O K , ε 0 ). 6)
Expand the basis of 4) into a basis for the lattice L of relations between the generators (a 1 , α 1 ), . . . , (a r ′ , α r ′ ) such that C ≃ Z r ′ /L.

7) Determine a cyclic basis of C.

Step 1) requires to apply the generalised discrete logarithm map in Cl + K 0 to the small number of relative norms of the basis elements of Cl K . Step 3) is possible since the a i a i = N K/K 0 (a i ) are trivial in Cl + K 0 . Steps 5) and 6) rely on the exactness of the sequence [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]. If r ′ = r, there is nothing to do. Otherwise, we first add the relation (O K , ε 0 ) 2 = 1. Lifts of relations from L 0 are then in the image of ε 0 / ε 2 0 , and if the unit exponent is odd in the lift, we need to add (O K , ε 0 ) into the relation. Steps 2) and 4) require an HNF, Step 7) an SNF.

The type norm subgroup

Algorithm 2 also provides an algorithm for generalised discrete logarithms in C, which can be used to determine the subgroup N Φ r (Cl K r ) in a similar way: For each generator of Cl K r , we compute the generalised discrete logarithm of its image in C, then the relations between the images using an HNF and a cyclic basis using an SNF. The enumeration of the subgroup is then trivial. In the same vein, it is possible to compute all the cosets C/N Φ r (Cl K r ) if the complete Igusa class polynomial is desired and not only its irreducible factor H 1 , see §2.3.

Computing ϑ-constants

As explained in Step 4) of Algorithm 1, it suffices to compute the fundamental ϑconstants ϑ 0 , . . . , ϑ 3 in the argument Ω/2 to obtain the class invariants for the period

matrix Ω = ω 0 ω 1 ω 2 ω 0 ∈ F 2 .
In §5.1 we describe an algorithm to compute the ϑ-constants directly from their q-expansions, using a lower number of multiplications than approaches described previously in the literature.

As the coefficients of the Igusa class polynomials grow rather quickly, a high floating point precision is needed for evaluating the ϑ-constants. In § §5.2-5.4 we describe an algorithm with a quasi-linear (up to logarithmic factors) complexity in the desired precision, using Newton iterations on a function involving the Borchardt mean. The algorithm is described essentially in Dupont's PhD thesis [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]; for the corresponding algorithm in dimension 1, using the arithmetic-geometric mean instead of the Borchardt mean, see [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]. We provide a streamlined presentation in dimension 2, together with improved algorithms and justifications.

Naive approach

For the fundamental ϑ-constants, (5) specialises as

ϑ 4b 1 +2b 2 (Ω/2) = m,n∈Z (-1) 2(mb 1 +nb 2 ) q m 2 0 q 2mn 1 q n 2 2 ( 15
)
with q k = exp(iπω k /2).

Positive definiteness and reducedness of the binary quadratic form attached to ℑ(Ω) show that the sum converges when taken over, for instance, a square [-R, R] 2 with R → ∞; [15, p. 210 following the proof of Lemma 10.1, with typos] establishes that for R √ 1.02N + 5.43, the truncated sum is accurate to N bits. Better bounds may be reached using summation areas related to the eigenvalues of ℑ(Ω), but using a square allows to organise and reuse computations so as to reduce the number of multiplications of complex numbers. where Õ(N ) = O N (log N ) O (1) , and M(N ) ∈ Õ(N ) is the time complexity of multiplying two numbers of N bits.

Proof. Using symmetries with respect to the signs of m and n, we may write

-R m,n R (-1) 2(mb 1 +nb 2 ) q m 2 0 q 2mn 1 q n 2 2 = 1 + 2 R m=1 (-1) 2mb 1 q m 2 0 + 2 R n=1 (-1) 2nb 2 q n 2 2 + 2 R m=1 (-1) 2mb 1 q m 2 0 R n=1 (-1) 2nb 2 q n 2 2 q 2mn 1 + q -2mn 1 .
We first compute and store the q n 2 2 with 2R + O(1) multiplications via q 2n-1 2 = q 2(n-1)-1 2

•q 2 2 and q n 2 2 = q (n-1) 2 2

•q 2n-1 2

. After computing the inverse q -1 1 , a similar scheme yields the q m 2 0 and q 2m 1 +q -2m 1 without storing them. At the same time, we may compute for any given m the sum over n inside the double sum: The term q 2mn 1 +q -2mn

1 is the n-th element v n of the Lucas sequence v 0 = 2, v 1 = q 2m 1 + q -2m 1 , v n = v 1 • v n-1 -v n-2
, each element of which is computed with one multiplication. Together with the multiplication by q n 2 2 , each term of the innermost sum is thus obtained with two multiplications. For the time complexity, recall that complex inversions can be computed in time O(M(N )), and exponentials in time O(M(N ) log N ), see [START_REF] Brent | Fast multiple-precision evaluation of elementary functions[END_REF]. This algorithm gains an asymptotic factor of 2/3 over [15, Algorithme 15].

Borchardt mean of complex numbers

The key tool in the asymptotically fast evaluation of ϑ-constants is the Borchardt mean, which generalises Lagrange's and Gauß's arithmetic-geometric mean of two numbers to four. The Borchardt mean of four positive real numbers has been introduced in [START_REF] Borchardt | Das arithmetisch-geometrische Mittel aus vier Elementen[END_REF][START_REF] Borchardt | Theorie des arithmetisch-geometrischen Mittels aus vier Elementen[END_REF]. The complex case is treated in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF], where proofs of most (but not all) propositions below may be found. It is made complicated by the presence of several square roots in the formulae, each of which is defined only up to sign.

Definition 4 Let

H = z ∈ C : arg(z) ∈ - π 2 , π 2 ∪ {0} = {z ∈ C : ℜ(z) > 0, or ℜ(z) = 0 and ℑ(z) 0}
be the complex half-plane defining the standard branch of the complex square root function. For a number in H, its square root in H lies in fact in the complex quarter-plane

Q = z ∈ C : arg(z) ∈ - π 4 , π 4 ∪ {0}.
′ 0 = 1 4 (b 0 + b 1 + b 2 + b 3 ) b ′ 1 = 1 2 ( b 0 b 1 + b 2 b 3 ) b ′ 2 = 1 2 ( b 0 b 2 + b 1 b 3 ) b ′ 3 = 1 2 ( b 0 b 3 + b 1 b 2 )
There are up to eight different Borchardt The following result is proved in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]Chapter 7].

Proposition 6 Any Borchardt sequence converges to a limit (z, z, z, z).

When the elements of b are contained in the same half-plane, the Borchardt mean B 2 (b) of b is the limit of the standard Borchardt sequence starting with b (0) = b. The function B 2 is obviously homogeneous.

A standard Borchardt sequence converges quadratically:

b (n) -B 2 (b) = 2 -O(2 n ) .

This implies that the Borchardt mean is computed to a precision of N bits with O(log N ) multiplications in time O(M(N ) log N ).

Comparison of the formulae in Definition 5 and [START_REF] Brent | Modern Computer Arithmetic[END_REF] shows that for any period matrix Ω ∈ H 2 , the sequence (ϑ 2 j (2 n Ω)) j=0,...,3 n 0 is a Borchardt sequence. This fact alone does not solve the sign issue, however. One would hope for the ϑ-sequence to be the standard Borchardt sequence, which would allow it to be computed with the standard choice of complex square roots. This assumption does not hold in general; however, it is true for the fundamental ϑ-constants and Ω ∈ F 2 .

Proposition 7

For Ω ∈ F 2 , n 0 and j = 0, . . . , 3 we have

ϑ j (2 n Ω) ∈ Q. Hence (ϑ 2 j (2 n Ω)) j=0,...,3 n 0 is the standard Borchardt sequence associated to (ϑ 2 j (Ω)) j=0,...,3 . It converges to 1.
The result follows from [15, Propositions 6.1 and 9.1].

Period matrix coefficients from ϑ-constants

For the time being, we consider the inverse of the function we are interested in and describe an algorithm that upon input of the values of the four fundamental ϑ-quotients in a period matrix returns the coefficients of the period matrix. Newton iterations can then be used to invert this function.

By the modularity of the squares of the ϑ-constants, applying a matrix γ ∈ Sp 4 (Z) to their argument Ω permutes the functions and multiplies them by a common projective factor, which depends on γ and Ω. In this way, information on Ω can be gathered; informally, three matrices suffice to obtain the three different coefficients of Ω. We consider three particular matrices, as suggested in [15, §9.2.3], which lead to well-behaved Borchardt means, see Conjecture 9.

Proposition 8 Let J = 0 -id 2 id 2 0 and M j = id 2 m j 0 id 2 with m 0 = 1 0 0 0 , m 1 = 0 1 1 0 , m 2 = 0 0 0 1 . Let Ω ∈ H 2 . Then ϑ 2 j ((JM 0 ) 2 Ω) j=0,1,2,3 = -iω 0 ϑ 2 j (Ω) j=4,0,6,2 , ϑ 2 j ((JM 1 ) 2 Ω) j=0,1,2,3 = (ω 2 1 -ω 0 ω 2 ) ϑ 2 j (Ω) j=0,8,4,12 , ϑ 2 j ((JM 2 ) 2 Ω) j=0,1,2,3 = -iω 2 ϑ 2 j (Ω) j=8,9,0,1 .
A more general statement with the action on the ϑ-constants (not squared) is given in [START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF]Propriété 3.1.24], following [START_REF] Igusa | Die Grundlehren der mathematischen Wissenschaften[END_REF]Chapter 5,Theorem 2]. The explicit form restricted to squares of ϑ-constants, as given here, is found in [15, §6.3.1].

The idea of the algorithm is now to apply the Borchardt mean function B 2 to both sides of the above equations. Conjecturally, the left hand side becomes 1, so that each Borchardt mean of a right hand side yields a coefficient of Ω. So we rely on the following conjecture, for which we have overwhelming numerical evidence, but no complete proof. Notice that it is a priori not even clear if the Borchardt means are well-defined, that is, if the squares of the various four ϑ-values always lie in the same half-plane. [START_REF] Borchardt | Das arithmetisch-geometrische Mittel aus vier Elementen[END_REF] is defined and equal to 1 .

Conjecture 9 Let

U = Ω ∈ H 2 : B 2 (ϑ 2 j (Ω)) j=0,...,
For k ∈ {0, 1, 2} we have (JM k ) 2 F 2 ⊆ U.
Under Conjecture 9, we can now formulate an algorithm to obtain Ω from four values of ϑ-constants. To make the following Newton iterations more efficient, we dehomogenise all modular functions by dividing by appropriate powers of ϑ 0 , which allows to work with only three inputs.

Algorithm 11

Input: Floating point approximations y (n) of Ω ∈ F 2 , to precision 2N bits, and x (n) of f (Ω), to precision N bits. Output: Floating point approximation x (n+1) of f (Ω), to precision 2N bits (see Theorem 12).

1) Let ε = 2 -N max j x (n) j .
2) Let (e j ) j=1,2,3 be the standard basis of C 3 . By Algorithm 10, compute F (x (n) ) and

∆F ∆x j = 1 ε F (x (n) + εe j ) -F (x (n) ) . 3) Let J = (J ij ) i,j=1,2,3 with J ij = ∆F i ∆x j . 4) Let x (n+1) = x (n) -F (x (n) ) -y (n) J -1
(where all vectors are seen as row vectors).

All computations in the algorithm are carried out at a precision of 2N bits. But even without taking rounding errors into account, the approximation of the Jacobian matrix by finite differences as well as the Newton method itself introduce some inaccuracy in the result, so that the accuracy improves to less than 2N bits. The following proposition addresses this issue.

Theorem 12 Assume the validity of Conjecture 9. Let Ω ∈ F 2 be such that ϑ 0 (Ω/2) = 0, x = f (Ω), and let x (0) be an initial floating point approximation to x. Not taking rounding errors into account, there exist two real numbers ε 0 > 0 and δ > 0, depending on x, such that for x (0) -x < ε 0 , the sequence x (n) defined by successive applications of Algorithm 11 converges to x, with accuracy increasing in each step from N to 2N -δ.

To reach a given accuracy N , the total complexity is dominated by the complexity of the last lifting step, that is:

O(M(N ) log N ).
Proof. By assumption, F is defined and analytic in a neighbourhood of x. In particular, its second partial derivatives are bounded in a close enough neighbourhood of x, so that the Jacobian matrix of F in x (n) is approximated to accuracy 2N -δ 0 bits by the matrix J computed in Steps 2) and 3), where δ 0 depends on x and on the bound on the second partial derivatives. The remaining assertion, with some δ δ 0 , is the standard result for Newton's method (see [START_REF] Von | Modern Computer Algebra[END_REF]Chapter 9 and §15.4] and [6, §4.2]).

The complexity is derived from the superlinearity of multiplication, which makes the last of the O(log N ) Newton steps dominate the whole computation; the logarithmic factor stems from the complexity of computing the Borchardt mean given in Proposition 6.

Notice that for our application of computing class polynomials for primitive quartic CM fields, the assumption of Theorem 12 is satisfied: As Ω id 2 Z 4 is of rank 4, we have ω 1 = 0, and therefore none of the ϑ j (Ω/2) vanish (see [27, Chapter 9, Proposition 2]).

In practice, we use a fixed initial precision for x (0) , computed according to Proposition 3, which determines ε and implicitly δ. The lack of information about δ can be worked around as follows: If x (n-2) and x (n-1) agree to k bits, and x (n-1) and x (n) agree to k ′ bits, we set δ = 2k -k ′ . This value of δ accounts at the same time for bits lost due to rounding errors induced by the floating point computations.

Remark. It is possible to modify Algorithm 10 and consequently Algorithm 11 so as not to rely on Conjecture 9. The conjecture states that the choices of square roots inside the Borchardt mean computations correspond to doubling the argument of the ϑconstants. So by computing very low precision approximations of the ϑ-constants in 2 n Ω as described in §5.1, one can make sure to choose the correct sign. These computations do not deteriorate the asymptotic complexity; moreover, as Algorithm 11 requires the Borchardt means of the same arguments over and over again (albeit with increasing precision), the sign choices may be fixed once and for all in a precomputation step.

In practice, however, we did not come upon any counterexample to Conjecture 9 with tens of thousands of arguments.

Reconstruction of class polynomial coefficients and reduction modulo prime ideals 6.1 The dihedral case

The class polynomials H 1 , Ĥ2 , Ĥ3 of ( 9) and ( 10) for a fixed CM type Φ are defined over K r 0 , but Steps 1) to 5) of Algorithm 1 compute floating point approximations, precisely of the images of the polynomials under an embedding ψ : K r 0 → C. To realise Step 6) of Algorithm 1, we need to invert ψ: Given ψ(α) to sufficient precision, we wish to reconstruct α symbolically as an element of K r 0 = Q(z r ) = Q[Z r ]/ (Z r ) 2 + A r Z r + B r , cf. §3.2. We may limit the discussion to the CM type Φ and the embedding ψ given by ( 13) and [START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF]; the second CM type Φ ′ leads to class polynomials that are conjugate under Gal(K r 0 /Q). Let D r be the discriminant of K r 0 ; as the discriminant of the minimal polynomial of z r is 16B, we have

K r 0 = Q( √ B), and D r B is a rational square. Let w ∈ K r 0 with w 2 = D r satisfy ψ(w) = √ D r > 0. Write α = a+bw c
with coprime a, b, c ∈ Z. Knowing an approximation β to ψ(α) ∈ R at our working precision of n bits, we wish to recover a, b, c, for which there is hope if 2 n > |abc|.

Let e be the exponent of β in the sense that 2 e-1 |β| < 2 e , and let e + = max(e, 0) and e -= max(-e, 0), so that e = e + -e -, and at most one of e + , e -is non-zero. We expect |a| ≈ √ D r |b| (whereas c is usually smaller), so that 2 e -|a| ≈ 2 e -√ D r |b| ≈ 2 e -|cβ| ≈ 2 e + |c|. On the other hand, the floating point approximation β satisfies

β -a+b √ D r c
≈ 2 e-n (up to a small factor accounting for digits lost to rounding errors), whence 2 n+e -cβ -(a + b √ D r ) ≈ 2 e + c is comparative in size to the previous quantities. Consider the integral matrix

    0 0 2 n+e + 2 e -√ D r 0 2 n+e + √ D r 0 2 e + β2 n+e +     .
Using LLL, we find a short vector -b 2 e -√ D r , c2 e + , r in the lattice spanned by the rows of the matrix; the scaling of the last column was chosen, following the arguments above, such that all entries in the vector have comparable sizes. This determines b and c, and we let a =

c β2 n+e + -b 2 n+e + √ D r 2 n+e + ∈ Z.
To get back to our standard representation of K r 0 , we need to relate w and z r . By ( 13) and [START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF],

ψ(z r ) = ψ(y r ) 2 = -A -2 √ B < 0, so that w = D r B • -z r -A 2 . ( 16
)
To obtain abelian varieties over finite fields, we need to reduce the class polynomials modulo certain prime ideals p 1 of K r 0 . Let p be a rational prime that splits as p 1 p 2 in K r 0 . Assume that p 1 splits in K r , so that p 1 = q 1 q 1 , and that the type norm of q 1 is a principal ideal of K. Then the class polynomial splits totally modulo p 1 , and its reduction may be computed as follows: If p 1 = pO K r 0 + (a + bw)O K r 0 with a, b ∈ Z, replace each occurrence of w bya b and reduce modulo p.

The cyclic case

Here the class polynomials are defined over Q, its coefficients may be obtained by a 2-dimensional lattice reduction, and reduction modulo primes is trivial.

Implementation and parallelisation

Our implementation of the algorithms defined here is available in the software package Cmh [START_REF] Enge | Cmh -Computation of Igusa Class Polynomials[END_REF], which can be downloaded from http://cmh.gforge.inria.fr/.

The current version of the Cmh software package is still in development, and will be named Cmh-1.0 once some packaging improvements, alongside with minor bug corrections, are checked in.

The software implements the different steps of Algorithm 1 as follows:

• Steps 1) to 3) of Algorithm 1 are performed by a script in Pari/Gp [1], which does all computations symbolically, and the running time of which is essentially negligible.

• The computation of ϑ-constants in Step 4) of Algorithm 1 is done by a C program, based on the library Gnu Mpc [START_REF] Enge | Gnu Mpc -A library for multiprecision complex arithmetic with exact rounding[END_REF], itself using the Gnu Mpfr [START_REF] Hanrot | Gnu Mpfr -A library for multiple-precision floating-point computations with exact rounding[END_REF] and Gnu Mp [START_REF] Granlund | Gmp -The Gnu Multiple Precision Arithmetic Library[END_REF] libraries. Newton lifting is used for this step from a base precision of 2 000 bits, and it is parallelised through MPI.

• Reconstruction of the class polynomials from the numerical values of the Igusa invariants is done inside the same C program, relying on the library Mpfrcx [START_REF] Enge | Mpfrcx -A library for univariate polynomials over arbitrary precision real or complex numbers[END_REF] for basic operations on polynomials using the FFT and asymptotically fast algorithms on trees of polynomials. In a preparatory step, the leaves of the tree for H 1 are filled with the linear factors of the class polynomial, those for Ĥk , k = 2, 3, are filled with the values of j k . Let the subscripts l and r denote the left and the right descendant, respectively, of a given node. Then an inner node n (1) in the tree for H 1 is computed as n (1) = n

(1)

l • n (1)
r , while an inner node n (k) in the tree for Ĥk , k = 2, 3, is obtained as n

(k) l • n (1) r + n (1) l • n (k)
r , where n (1) denotes the node at the same position in the tree for H 1 ; for details, see [START_REF] Von | Modern Computer Algebra[END_REF]Algorithms 10.3 and 10.9]. By first combining pairs of complex-conjugate leaves in a preprocessing step, all computations are in fact carried out with real floating point polynomials, see [START_REF] Enge | Fast decomposition of polynomials with known Galois group[END_REF]. So if at a given level the tree for H 1 contains m nodes, all nodes at this level of the three trees can be obtained with 5m independent multiplications, which are parallelised using MPI.

• Recognition of the polynomial coefficients as elements in K r 0 is also embedded in the C program, using Fplll [START_REF] Cadé | Fplll[END_REF] for the LLL step.

• Assuming that K = Q(ζ 5 ), for which the result is known, validation of the obtained class polynomials is performed by computing a Weil number π above a prime p ≈ 2 128 , constructing a curve over F p having as endomorphism ring the ring of integers of K using Mestre's algorithm [START_REF] Mestre | Construction de courbes de genre 2 à partir de leurs modules[END_REF], and verifying that the cardinality of the Jacobian matches N K/Q (1 ± π). This step is done in Pari/Gp and also has a negligible cost.

In the following we report on the performance of these different steps, illustrated by both small and large examples.

Computation of ϑ-constants

We report timing results for the computation of fundamental ϑ-constants for two arbitrary period matrices. Table 1 shows that already our implementation of the relatively simple naive algorithm presented in §5.1 may be several orders of magnitude faster 3 than Magma-2.19.4, the performance improvement ratio depending on the period matrix. Newton lifting is preferable above some cut-off value for the precision, here 16 000 and 4 000 bits, respectively. The naive algorithm is rather sensitive to the period matrix; generally speaking, it converges the faster the larger the imaginary parts in Ω are, which correspond to smaller q 0 , q 1 , q 2 . A noticeable difference between our naive algorithm from §5.1 and the implementation in Magma is that the favorable cases are not the same. This is most likely due do different choices of summation regions, as briefly discussed in §5.1. We note that the timings of Newton lifting depend much less on the period matrix entries than those for the naive method. Notice that the running times for Newton lifts are consistent with the theoretical complexity of O(M(N ) log N ). The code in Cmh implements the approach using finite differences for estimating the Jacobian matrix as described in §5.4, as well as an algorithm which computes the exact Jacobian matrix along with the Borchardt mean as given in [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]Algorithme 16]. Both converge equally well, but the latter approach is computationally more expensive by roughly 45 %, accounted for by a larger number of multiplications. Our code distinguishes orbits of the roots of the Igusa class polynomials under complex conjugation. For instance, there are four real roots and 58 pairs of complexconjugate roots in the second example, so that altogether we need to carry out 62 lifts of ϑ-constants. Instead of targeting a given precision based on arguments as developed The timings indicated as "preparation" and "check" in Table 2 correspond to the number theoretic calculations performed in Pari/Gp. The preparation time covers the enumeration of N Φ r (Cl K r ) ⊆ C, and the creation of the relevant set of reduced period matrices. Checking means finding a Weil number over a 128-bit prime and generating a genus 2 curve whose Jacobian has complex multiplication by the maximal order of K.

Ω = -1+5i 2 i 6 i 6 -1+7i 2 Ω = 2+10i

Breakdown of timings for small class polynomial examples

K = Q[X]/(X 4 + 144X 2 + 3500) C = N Φ r (Cl K r ) = Z/2Z × Z/
H 1 , Ĥ2 , Ĥ3 ∈ C[X] 0.1 H 1 , Ĥ2 , Ĥ3 ∈ K r 0 [X] 3×0.3 check 0.8 Total (incl. I/O) 12.4 K = Q[X]/(X 4 + 134X 2 + 712) C = N Φ r (Cl K r ) = Z/2Z × Z

A large example

Our currently largest example is K = Q[X]/(X 4 + 1357X 2 + 3299), containing K 0 = Q( √ 1828253) of class number 2. Its Shimura class group is C = N Φ r (Cl K r ) ≃ Z/2Z × Z/2Z × Z/5004Z of size 20016. On one core of an Intel Core i5-4570 clocked at 3.2 GHz, the structure of the class group is obtained with our Pari/Gp script in roughly one second, while the computation of the period matrices and their symbolic reduction into the fundamental domain F 2 takes 388 s. The associated ϑ-constants consist of 10008 pairs of complex-conjugate values. For the first ten Newton iterations up to a precision of about 2 000 000 bits, we used 640 cores Intel Xeon X5675 at 3.07 GHz; for the last two iterations, we switched to a machine with only 160 cores Intel Xeon E7-8837 at 2.67 GHz, but with 640 GB of main memory. Table 3 gives the timings (in seconds) for the Newton lifts of one particular period matrix. The small value of δ, estimated as explained at the end of §5.4, shows that the effective precision indeed almost doubles in each step as predicted by Theorem The lifting step accounts for a total of about 420 CPU days, but thanks to its easy parallelisation on 160 to 640 cores, it was finished in less than 4 days wall-clock time (including additional overhead for writing intermediate results to disk).

The computation of the floating point polynomials H 1 , Ĥ2 and Ĥ3 was carried out at a precision of 7 850 071 bits (the lowest lifting precision reached for one of the period matrices). After regrouping complex conjugates, the first step consists of 5 • 10008/2 = 25020 multiplications of monic polynomials of degree 2 with real coefficients, which can be arbitrarily parallelised; we used a machine with 24 Intel Xeon E7540 cores at 2.0 GHz and 512 GB of memory. At degrees 2 048 to 8192, the FFT multiplications required too much memory to be executed in parallel on all cores, so we reduced the number of simultaneous multiplications to the maximum possible, as indicated in Table 4. In the last step, we needed to multiply a degree 16 384 polynomial with a degree 3 632 polynomial. The wall-clock time of this polynomial reconstruction step, counting inputoutput time and the cost of resuming computations, was almost exactly 3 days.

Recognising one coefficient of the floating point polynomials as an element of K r 0 took on average 980 s per coefficient on one Intel Xeon E5-2650 core at 2 GHz. The total CPU time for the 60 045 coefficients was thus about 680 CPU days; with up to 480 cores working in parallel, this took less than 2 wall-clock days.

The uncompressed storage size of the three resulting polynomials in base 10 is about 90 GB. The common denominator of the coefficients of H 1 has 8 884 distinct prime factors, the largest one being 1 506 803 839. It occurs to powers 2 in H 1 and 4 in Ĥ2 and Ĥ3 , consistent with the fact that the power of h 10 in the denominator of j 2 and j 3 is 2 instead of 1 for j

  and ℑ(M ) are the images under ψ of matrices with entries in K r 0 . The condition | det(CΩ + D)| 1 can be rewritten as det(Cψ(M ) + D) det(Cψ(M ) + D) 1 and thus also depends only on the images under ψ of elements of K r 0 .

4 . 1

 41 where √ D r is the positive root of D r and a, b ∈ Q. This sign can be determined from the signs of a and b and the relative magnitudes of a 2 and b 2 D r . 4 Computing the Shimura group and its type norm subgroup Structure of the Shimura group C

Proposition 3

 3 The truncated sum over (m, n) ∈ [-R, R] 2 for the fundamental ϑconstants (15) may be computed with 2R 2 + O(R) multiplications and one inversion using storage for R + O(1) elements.Letting R = ⌈ √ 1.02N + 5.43⌉ and using complex numbers of precision O(N ), we obtain a time complexity of O(N M(N )) or Õ(N 2 ),

Definition and Properties 5

 5 Given a complex quadruple b = (b 0 , . . . , b 3 ) ∈ C 4 , a Borchardt iterate is a quadruple b ′ = (b ′ 0 , . . . , b ′ 4 ) such that there are four choices of square roots ( b j ) j=0,...,3 yielding b

  has determined a finite set of inequalities describing a fundamental domain F 2 for Sp 4 (Z) \H2, which directly translate into an algorithm for reducing an element of H 2 into F 2 . As the Igusa functions introduced in §2.3 are modular for Sp 4 (Z), we may transform all period matrices occurring in Algorithm 1 into F 2 . A period matrix Ω is reduced if ℜ(Ω) has coefficients between -1 2 and 1 2 (which may be obtained by reducing modulo Z), if the binary quadratic form defined by ℑ(Ω) is reduced (which may be obtained using Gauß's algorithm) and if | det(CΩ + D)| 1 for each of

	19 matrices	A B C D	∈ Sp 4 (Z) (which may be obtained by applying to Ω a matrix

  iterates of a given quadruple. If b ∈ H 4 , the standard Borchardt iterate is obtained by choosing square roots in Q, so that b ′ ∈ H 4 again. More generally, if all entries of b lie in the same half-plane, that is, b ∈ (zH)4 for some z ∈ C, choosing all square roots in the same quarter-plane If all entries of b (0) lie in the same half-plane, its standard Borchardt sequence is defined by taking only standard Borchardt iterates.

	choice of sign for	√ zQ (with either √ z) yields the standard Borchardt iterate in the same half-plane.
	A Borchardt sequence is a sequence b (n)	n 0	such that b (n+1) is a Borchardt iterate
	of b (n) for all n	0.

Table 1 :

 1 Computation of ϑ

		1+2i
	7 1+2i 6	6 10 + 8i 4

0 (τ ) (Intel i5-2500, 3.3GHz; Magma-2.19.4; Cmh-1.0)

Table 2

 2 illustrates our class polynomial computations on relatively small examples.

Table 2 :

 2 Timings in seconds for two examples (on one Intel i5-2500, 3.3GHz) in[START_REF] Streng | Computing Igusa class polynomials[END_REF], we simply carry out successive lifting steps until the polynomial reconstruction succeeds. This explains the time needed for failed reconstruction attempts in Table2, which could be avoided if we had a sharper bound on the required precisions. It regularly occurs, even though this is not illustrated by the examples here, that the reconstruction of the class polynomialH 1 ∈ K r 0 [X] succeeds one lifting step before that of Ĥ2 , Ĥ3 ∈ K r 0 [X]. This can be explained by the relative size of the invariants considered by Streng, see[START_REF] Streng | Complex multiplication of abelian surfaces[END_REF] Appendix 3].

		/60Z
	preparation	0.3
	base, 2 000 bits	1.1
	lift, 3 988 bits	1.6
	lift, 7 958 bits	4.4
	reconstruction attempt	0.1
	lift, 15 886 bits	13.1
	reconstruction attempt	0.2
	lift, 31 744 bits H 1 , Ĥ2 , Ĥ3 ∈ C[X] H 1 , Ĥ2 , Ĥ3 ∈ K r 0 [X] check	38.7 0.6 1.8 + 2×1.4 0.7
	Total (incl. I/O)	69.2

Table 3 :

 3 [START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF]. Time for lifting steps for example with #C = 20 016.

	precision	δ	time
	2 000 -	0.03
	3 986 14	0.01
	7 970	2	0.1
	15 932	8	0.3
	31 862	2	1.0
	63 718	6	2.9
	127 434	2	8.2
	254 858 10	24
	509 714	2	61
	1 019 416 12	150
	2 038 832	0	360
	4 077 652 12	940
	8 155 302	2 2 100

Table 4 :

 4 1 . Polynomial reconstruction timings for example with #C = 20 016.

	input degree # multiplications wall-clock time (s) #threads
	2	25 020	380	24
	4	12 510	560	24
	8	6 225	780	24
	16	3 125	1 300	24
	32	1 565	1 800	24
	64	780	2 700	24
	128	390	3 700	24
	256	195	6 300	24
	512	100	11 000	24
	1 024	50	13 000	24
	2 048	25	25 000	10
	4 096	10	20 000	10
	8 192	5	71 000	3
	16 384	5	67 000	5

http://echidna.maths.usyd.edu.au/echidna/dbs/complex_multiplication2.html

Personal communication

Such a quadratic, yet efficient implementation was used by T. Houtmann to compute class polynomials of degree up to 500 (personal communication, no reference exists).

Acknowledgements

The authors would like to acknowledge the work of Régis Dupont, whose thesis has been extensively used as a basis from the outset of this work. We are also grateful to Damien Robert and Marco Streng for fruitful discussions and to David Gruenewald for his comments on an earlier version.

Computer experiments have used a variety of computing resources funded from different projects. We thus acknowledge the support of the Région Lorraine CPER MISN TALC project; the PlaFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS (see https://plafrim.bordeaux.inria.fr/); the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour.

This research was partially funded by ERC Starting Grant ANTICS 278537 and by Agence Nationale de la Recherche grants ANR-09-BLAN-0020-01 and ANR-12-BS01-0010-01.

Algorithm 10

Input: Floating point approximations of (ϑ j (Ω/2)/ϑ 0 (Ω/2)) j=1,2,3 for some Ω ∈ F 2 , and as auxiliary data the sign of ω 1 . Output: Floating point approximations of the coefficients ω 0 , ω 1 , ω 2 of Ω ∈ F 2 1) Use the duplication formulae (6) to compute (ϑ 2 j (Ω)/ϑ 2 0 (Ω/2)) j=0, 1,[START_REF] Belding | Computing hilbert class polynomials[END_REF][START_REF] Borchardt | Das arithmetisch-geometrische Mittel aus vier Elementen[END_REF][START_REF] Borchardt | Theorie des arithmetisch-geometrischen Mittels aus vier Elementen[END_REF][START_REF] Brent | Modern Computer Arithmetic[END_REF][START_REF] Bröker | A p-adic algorithm to compute the Hilbert class polynomial[END_REF][START_REF] Bröker | Explicit CM theory for level 2-structures on abelian surfaces[END_REF][START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF][START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF] .

2 1,[START_REF] Belding | Computing hilbert class polynomials[END_REF][START_REF] Borchardt | Das arithmetisch-geometrische Mittel aus vier Elementen[END_REF][START_REF] Borchardt | Theorie des arithmetisch-geometrischen Mittels aus vier Elementen[END_REF][START_REF] Brent | Modern Computer Arithmetic[END_REF][START_REF] Bröker | A p-adic algorithm to compute the Hilbert class polynomial[END_REF][START_REF] Bröker | Explicit CM theory for level 2-structures on abelian surfaces[END_REF][START_REF] Cosset | Applications des fonctions thêta à la cryptographie sur courbes hyperelliptiques[END_REF][START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF] .

4) Compute

The correctness of Algorithm 10 under Conjecture 9 follows from the discussions above.

Step 1) uses the homogeneity of (6), Step 2) the homogeneity of the Borchardt mean and Proposition 7. The validity of Step 5) follows from Proposition 8 under Conjecture 9, using again the homogeneity of the Borchardt mean.

Notice that Ω is only well-defined up to the subgroup of Sp 4 (Z) for which the ϑconstants are modular. Assuming Ω ∈ F 2 , the fundamental domain for all of Sp 4 (Z), it is necessarily unique; Conjecture 9 implies that this particular representative for Ω is indeed returned by the algorithm.

Newton lift for fundamental ϑ-constants

Denote by

the function computed by Algorithm 10, and by

its inverse on F 2 (where Ω is interpreted as the three-element vector (ω 0 , ω 1 , ω 2 ) and not as a four-element matrix). We use Newton iterations on F to compute f . The standard Newton approach requires to compute the Jacobian matrix J F of F , that is, its partial derivatives with respect to its different coordinates. Heuristically, Algorithm 10 may be modified accordingly to also output J F , see [START_REF] Dupont | Moyenne arithmético-géométrique, suites de Borchardt et applications[END_REF]Algorithme 16], generalising the dimension 1 approach of [5, §2.4] and [START_REF] Dupont | Fast evaluation of modular functions using Newton iterations and the AGM[END_REF]. The description and justification of this algorithm are rather technical. Instead, we opt for using finite differences, which moreover turn out to yield a more efficient algorithm (see §7.1).