Adam Shimi
email: adam.shimi@irit.fr

Armando Castañeda
email: armando.castaneda@im.unam.m

K-set agreement bounds in round-based models through combinatorial topology

Keywords: Distributed Computability, Set-agreeement, Round-based models, Combinatorial Topology, Lower bounds, Upper Bounds

Round-based models are very common message-passing models; combinatorial topology applied to distributed computing provides sweeping results like general lower bounds. We combine both to study the computability of -set agreement.

Among all the possible round-based models, we consider oblivious ones, where the constraints are given only round per round by a set of allowed graphs. And among oblivious models, we focus on closed-above ones, that is models where the set of possible graphs contains all graphs with more edges than some starting graphs. These capture intuitively the underlying structure required by some communication model, like containing a ring.

We then derive lower bounds and upper bounds in one round for -set agreement, such that these bounds are proved using combinatorial topology but stated only in terms of graph properties. These bounds extend to multiple rounds when limiting our algorithms to be oblivious -recalling only pairs of processes and initial value, not who send what and when.

INTRODUCTION 1.Motivation

Rounds structure many models of distributed computing: they simplify algorithms, capture the distributed equivalent of time complexity [START_REF] Fraigniaud | Towards a complexity theory for local distributed computing[END_REF], and underly many fault-tolerant algorithms, like Paxos [START_REF] Lamport | The part-time parliament[END_REF].

https://doi.org/10.1145/3382734.3405752 A recent trend, with parallel results by Charron-Bost and Schiper [START_REF] Charron-Bost | The heard-of model: computing in distributed systems with benign faults[END_REF] on one hand, and Afek and Gafni [START_REF] Afek | Asynchrony from synchrony[END_REF] and Raynal and Stainer [START_REF] Raynal | Synchrony weakened by message adversaries vs asynchrony restricted by failure detectors[END_REF] on the other hand, is using this concept of round for formalizing many different models within a common framework. But the techniques used for proving results in these models tend to be ad-hoc, very specific to some model or setting. What is required going forward is a general approach to proving impossibility results and bounds on round-based models.

Actually, there is at least one example of a general mathematical technique used in this context: the characterization of consensus solvability through point-set topology by Nowak et al. [START_REF] Nowak | Topological characterization of consensus under general message adversaries[END_REF]. We propose what might be seen as an extension to higher dimension of this intuition, by applying combinatorial topology (instead of point-set topology) to bear on -set agreement (instead of just consensus).

Combinatorial topology abstracts the reasoning around knowledge and indistinguishability behind many impossibility results in distributed computing. It thus provide generic mathematical tools and methods for deriving such results [START_REF] Herlihy | Distributed Computing Through Combinatorial Topology[END_REF]. Moreover, this approach is the only one that managed to prove impossibility results and characterization of solvability of the -set agreement [START_REF] Chaudhuri | More choices allow more faults: Set consensus problems in totally asynchronous systems[END_REF], our focus problem.

Concretely, we look at closed-above round-models, that is models where constraints happens round per round, and the set of communication graphs allowed is the closure-above of a set of graphs. These models capture some safety properties, where we require some underlying structure in communication, like having an underlying star, ring or tree. This is a strict generalization of the models with a fixed communication graph considered by Castañeda et al. [START_REF] Castañeda | A topological perspective on distributed network algorithms[END_REF].

For our models, we derive upper bound and lower bounds on the for which -set agreement is solvable. And although the proofs of the bounds use combinatorial topology, they are stated in terms of variants of the domination number, a well-known and used combinatorial number on graphs.

Overview

• We start by defining closed-above models in Section 2.

• Then we give various upper bounds for -set agreement in one round on those models in Section 3. These have the advantage of not requiring any combinatorial topology. • Next, we introduce in Section 4 the combinatorial topology necessary for our lower bounds, both the basic definitions and our main technical lemma. • We then go to lower bounds on round-based models forset agreement in one round in Section 5. Recall that these bounds use combinatorial topology, but are stated in terms of graph properties. • Finally, Section 6 generalize both upper and lower bounds to the case of multiple rounds.

Due to size constraints, most of the proof can be found in the full version [START_REF] Shimi | K set-agreement bounds in round-based models through combinatorial topology[END_REF]

Related Works

Round-based models. The idea of using rounds for abstracting many different models is classical in message-passing. This includes the synchronous adversary models of Afek and Gafni [START_REF] Afek | Asynchrony from synchrony[END_REF] and Raynal and Stainer [START_REF] Raynal | Synchrony weakened by message adversaries vs asynchrony restricted by failure detectors[END_REF]; the Heard-Of model of Charron-Bost and Schiper [START_REF] Charron-Bost | The heard-of model: computing in distributed systems with benign faults[END_REF]; and the dynamic networks of Kuhn et al [START_REF] Kuhn | Distributed computation in dynamic networks[END_REF].

Rounds are also used for building a distributed theory of time complexity [START_REF] Fraigniaud | Towards a complexity theory for local distributed computing[END_REF] and for structuring fault-tolerant algorithms like Paxos [START_REF] Lamport | The part-time parliament[END_REF].

Previous work on the solvability of consensus and -set agreement include the characterization of consensus solvability for oblivious round-based models of Coulouma et al. [START_REF] Coulouma | A characterization of oblivious message adversaries for which consensus is solvable[END_REF], the failure-detectorbased approach of Jeanneau et al. [START_REF] Jeanneau | Solving k-set agreement using failure detectors in unknown dynamic networks[END_REF], and the focus on graceful degradation in algorihtms for -set agreement of Biely et al. [START_REF] Biely | Gracefully degrading consensus and k-set agreement in directed dynamic networks[END_REF].

Combinatorial Topology. Combinatorial Topology was first applied to the problem of -set agreement in wait-free shared memory by Herlihy and Shavit [START_REF] Herlihy | The topological structure of asynchronous computability[END_REF], Saks and Zaharoglou [START_REF] Saks | Wait-free k-set agreement is impossible: The topology of public knowledge[END_REF] and Borowsky and Gafni [START_REF] Borowsky | Generalized FLP impossibility result for Tresilient asynchronous computations[END_REF].

Beyond these first forays, many other results got proved through combinatorial topology. Among others, we can cite the lower bounds for renaming by Castañeda and Rajsbaum [START_REF] Castañeda | New combinatorial topology bounds for renaming: the lower bound[END_REF] and the derivation of lower-bounds for message-passing by Herlihy and Rasjbaum [START_REF] Herlihy | The topology of distributed adversaries[END_REF]; There is even a result by Alistarh et al. [START_REF] Alistarh | Why extension-based proofs fail[END_REF] showing that traditional proof techniques (dubed extension-based proofs) cannot prove the impossibility of -set agreement in specific shared-memory models, whereas techniques from combinatorial topology can.

For a full treatment of combinatorial topology applied to distributed computing, see Herlihy et al. [START_REF] Herlihy | Distributed Computing Through Combinatorial Topology[END_REF].

Combination of Topology and Round-based models. Two papers at least applied topology (combinatorial or not) to general roundbased models in order to study agreement problems: Godard and Perdereau [START_REF] Godard | k-set agreement in communication networks with omission faults[END_REF] used combinatorial topology to study -set agreement in models with omission failures; and Nowak et al. [START_REF] Nowak | Topological characterization of consensus under general message adversaries[END_REF] characterization of consensus for general round-based models (not necessarily oblivious) using point-set topology.

DEFINITIONS 2.1 Communication models

One common feature of many models of distributed computation is the notion of rounds, or layers. Formally, rounds are communicationclosed as defined by Elrad and Francez [START_REF] Elrad | Decomposition of distributed programs into communication-closed layers[END_REF]: at each round, a process only takes into account the messages (or information) sent by other processes at this same round.

Traditionally, rounds are thought of as synchronous: synchrony indeed provides a natural way to implement them. But asynchronous rounds also exist, both in message-passing [START_REF] Charron-Bost | The heard-of model: computing in distributed systems with benign faults[END_REF] and in sharedmemory [START_REF] Borowsky | Immediate atomic snapshots and fast renaming[END_REF][START_REF] Herlihy | The asynchronous computability theorem for t-resilient tasks[END_REF].

Here we abstract away all implementation details, and consider a model with rounds, parameterized with the allowed sequences of communication graphs -directed graphs where each node is a process and each arrow correspond to a delivered message to the destination from the source. There are no crashes, just a specification of which message can be received at which round. This abstracts the Heard-Of model [START_REF] Charron-Bost | The heard-of model: computing in distributed systems with benign faults[END_REF], synchronous message adversaries [START_REF] Afek | Asynchrony from synchrony[END_REF][START_REF] Raynal | Synchrony weakened by message adversaries vs asynchrony restricted by failure detectors[END_REF] and dynamic networks [START_REF] Kuhn | Distributed computation in dynamic networks[END_REF][START_REF] Kuhn | Coordinated consensus in dynamic networks[END_REF], as well as all other models relying only on the properties of rounds.

We fix Π = { 1 , ..., } as our set of processes for the rest of the paper.

Definition 2.1 (Communication model). Let ℎ Π be the set of graphs. Then ⊆ ( ℎ Π ) is a communication model.
Any set of infinite sequences of graphs defines a model. In order to make models more manageable, we focus on a restricted form, where the graph for each round is decided independently of the others. The model is thus entirely characterized by the set of allowed graphs. We call these communication models oblivious, following Coulouma et al [START_REF] Coulouma | A characterization of oblivious message adversaries for which consensus is solvable[END_REF].

Definition 2.2 (Oblivious communication models).

Let be a communication model. Then is oblivious ∃ ⊆ ℎ Π : = .

Intuitively, oblivous models capture safety properties: bad things that must not happen. Or equivalently, good things that must happen at every round. Usually, these good properties are related to connectivity, like containing a cycle or a spanning tree. Since such a property tends to be invariant when more messages are sent, we can look at oblivious models defined by a set of subgraphs. Classical examples of closed-above models are the non-empty kernel predicate (only graphs where at least one process broadcasts) and the non-split predicate (only graphs where each pair of processes hears from a common process), used notably by Charron-Bost et al. [START_REF] Charron-Bost | Approximate consensus in highly dynamic networks: The role of averaging algorithms[END_REF] for characterizing the solvability of approximate consensus (the variant of consensus where the decided value should be less than apart, where > 0 is fixed beforehand). Another closed-above model is the one satisfying the tournament property of Afek and Gafni [START_REF] Afek | Asynchrony from synchrony[END_REF], which they show is equivalent to wait-free read-write shared memory.

One example of an oblivous model which is not closed-above is the one generated by all graphs containing a cycle, except the clique. More generally, the closure-above forces us to have all graphs with more edges than our generators.

Nonetheless, closed-above models capture a fundamental intuition behind distributed computing models: specifying what should not happen. They also have a good tradeoff between expressivity and simplicity, since the "combinatorial data" used to build them is contained in a small number of graphs. Finally, the patterns expected by safety properties tend to be independent of which processes play which roles -what matters is the existence of a ring or spanning tree, not who is where on it.

We call such closed-above models symmetric.

Definition 2.4 (Symmetric models). Let be a closed-above model, and be the set of graphs generating it. Then is symmetric = ( ), where

( ) = { ( ) | ∈ ∧ : Π → Π a permutation on Π}.
In the rest of the paper, we will limit ourselves to closed-above models, both symmetric and not.

Oblivious algorithms

Because most applications of combinatorial topology to distributed computing aim towards impossibility results, the traditional algorithms considered err on the side of power: full information protocols, which exchange at each round the view of everything ever heard by the process. For example, after a couple of rounds, views will contain nested sets of views, themselves containing views, recursively until the initial values.

In contrast, we focus on oblivious algorithms. That is, we limit each process to remember only the initial values it knows, not who sent them or when. This amounts to a function from Π to the set of initial values (with a ⊥ when the value is not known). In turn, these algorithms lose the ability to trace the path of the value.

We can view oblivious algorithms as full-information protocol whose decision map (the function from final view to decision value) depends only on the set of known pairs (process,initial value). The full-information protocol might still be used for deciding when to apply the decision map, but this map loses everything except the known pairs. That is, the decision map is constrained to decide similarly in situations where it received the same values, even when they were from different processes.

Definition 2.5 (Oblivious algorithm).

Let A be a full-information protocol, with decision map . Then A is an oblivious algorithm ∀ a view :

( ) = ( ( )), where ( ) = ( , ) ∈ (( , )) and (( , )) = {( , )} if is a singleton from ( ) otherwise

ONE ROUND UPPER BOUNDS: A START WITHOUT TOPOLOGY

Although lower bounds are our targets, they require upper bounds to gauge their strength. We thus start with upper bounds onset agreement [START_REF] Chaudhuri | More choices allow more faults: Set consensus problems in totally asynchronous systems[END_REF] for closed-above models. Another advantage of starting with our upper bounds is that they rely on concrete algorithms, and allow us to introduce generalizations of the classical domination number that will be used for our lower bounds. Lastly, we also start with bounds for the one round case in this section and the next one. Bounds for multiple rounds depend on these one round bounds.

These bounds follow from a very simple algorithm for solvingset agreement. We assume the set of initial values is totally ordered. Then everybody sends its initial value for one round, and decide the minimum it received.

Simple closed-above models: almost too easy

Recall that the domination number of a graph is the size of its smallest dominating set, that is the size of the smallest set of nodes whose set of outgoing neighbors is Π. Note that the outgoing neighbors of a set ⊆ Π contains -that is, we assume self-loop.

Definition 3.1 (Domination number). Let be a graph. Then its domination number

( ) { ∈ [1, ] | ∃ ⊆ Π : | | = ∧ ∈ ( ) = Π}.
Because the simple closed-above model generated by only allows graphs containing , their domination number is at most ( ). This entails a very simple upper bound on -set agreement.

T 3.2 (U ( ))
. Let be a graph. Then ( )-set agreement is solvable in one round on the simple closed-above model generated by .

P

. The algorithm is just slightly different from the one stated at the start of the section: after one round, each process decides the minimum value of the ones of a fixed minimum dominating set of . Since is known, this minimum dominating set can be computed beforehand. And because it is a dominating set, every process receives at least one value from it, so every process can decide.

Finally, since the minimum dominating set has at most ( ) distinct values, at most ( ) values are decided, and thus our algorithm solves ( )-set agreement.

From Castañeda et al. [START_REF] Castañeda | A topological perspective on distributed network algorithms[END_REF]Thm 5.1], we know this bound is tight: the oblivious model with a single graph cannot solve -set agreement in one round for < ( ). Hence the weaker simple closedabove model generated by cannot solve -set agreement in one round for < ( ).

Still, simple closed-above models are somewhat artificial, as can be seen in the proof: we know exactly the subgraph that must be contained in the actual communication graph. A more realistic take requires to spread the uncertainty to the underlying subgraph; we thus look next at general closed-above models.

General closed-above models: tweaking of upper bounds

For general closed-above models, we must deal with a set of possible underlying subgraphs. This makes our previous approach inapplicable: we cannot hardcode a dominating set because we don't know the underlying subgraph for sure. This new issue motivates the definition of a weakening of the domination number: the equal-domination number of a set of graphs. Intuitively, any set of that much process is a dominating set in all the graphs considered. Definition 3.3 (Equal Domination number of a set of graphs). Let be a set of graphs. Then its equal domination number

( ) max ∈ ( ), where ( ) = { ∈ [1, ] | ∀ ⊆ Π : | | = =⇒ ∈ ( ) = Π}. T 3.4 (U ( ) ).
Let be a set of graphs. Then ( )-set agreement is solvable on the closed-above model generated by .

P

. Let be a set of ( ) processes with the smallest initial values. They have thus at most ( ) distinct initial values. By definition of ( ), dominates every graph in , and thus every graph in the closed-above model generated by .

Thus taking the minimum after one round will result in deciding one of those initial values, and thus one of at most ( ) values. We conclude that our algorithm solves ( )-set agreement after one round on the closed-above model generated by .

Since the equal-domination number is independent of which process does what, it is the same for any permutation of the graph. This entails an upper bound on symmetric models as a corollary.

C

3.5. Let be a set of graphs. Then ( )-set agreement is solvable on the closed-above model generated by ( ).

Now, the natural question to ask is whether we can improve this bound. Or equivalently, is it tight?

The answer depends on the graphs. To see it, let us look at another combinatorial number : covering numbers. Given fewer processes than the equal-domination number of the graph, they do not always form a dominating set. Nonetheless, they might still get heard by some minimum number of processes. We call such minimums the covering numbers of the graph: the -th covering number of is, given any set of processes, the minimum number of processes hearing this set in . These numbers capture the ability of a set of processes to disseminate their values in the graph. If we take the processes with the smallest initial values, we can be sure that at least ( ) processes will hear, and thus choose one of these. This then gives a solution to ( + ( -( )))-set agreement in one round.

P

. For a set of processes with the smallest initial values, they will reach at least ( ) processes after the first round. Thus these processes will decide one of the values when taking the smallest value they received.

As for the rest of the processes, we can't say anything about what they will receive, and thus we consider the worst case, where they all decide differently, and not one of the smallest values. Then the number of decided values is at most + ( -( )), and the theorem follows.

The covering numbers are also independent of processes names; we thus get a similar upper bound on symmetric models as a corollary. When is this new bound better than the one using the equaldomination number? When there is some such that -( ) < ( ) -. Let us take the symmetric models generated by the two graphs in Figure 1.

In the first model, -( ) < ( ) -never happens, because every covering number of a star equals 1 (the biggest set of outgoing neighbors different from Π contains only one process), and its equal-domination number equals (because when taking only -1 processes, the center of the star might not be in there). Thus -( ) = -1 ≥ ( ) -= -. On the other hand, this is the case in the second model, because 2 ( ) = 3 and ( ) = 4. We thus we have

-2 ( ) = 4 -3 = 1 < ( ) -= 4 -2 = 2.
Hence the upper bound with covering numbers ensure 3-set agreement solvability while the upper bound with the equal-domination number only ensures 4-set agreement solvability.

Intuitions on upper and lower bounds

Why do our upper bounds hold? Because we can extract from the underlying graphs some minimal connectivity of sets of processes. Hence, we know from these combinatorial numbers how much the minimal values will spread in the worst case, and thus we bound the maximum number of values decided.

On the other hand, our lower bounds will follow from studying how much values can spread in the best case. Why? Because the more values can spread, the more processes can distinguish between initial configurations, and the more they have a chance to decide correctly. Ensuring enough indistinguishability thus entails an impossibility at solving -set agreement. This indistinguishability is linked to higher-dimension connectivity in combinatorial topology [START_REF] Herlihy | Distributed Computing Through Combinatorial Topology[END_REF]Thm. 10.3.1]; we thus turn to the topological approach to distributed computing for our lower bounds.

ELEMENTS OF COMBINATORIAL TOPOLOGY 4.1 Preliminary definitions

First, we need to introduce the mathematical objects that this approach uses. These are simplexes and complexes. A simplex is simply a set of values, and can be represented as a generalization of a triangle in higher dimensions. Simplexes capture configurations in

1 2 3 (a) Graph ( 1 , { 1 , 3 } ( 2 , { 1 , 2 } ( 3 , { 3 } (b)
( ) = { ∈ | ∃ ∈ : ( , ) ∈ }. We also write ( ) for the ∈ such that ( , ) ∈ The dimension of is | | -1.
Although we define Views to be any set for readability, the traditional view is of sets of pairs, the first element being a process name, and the second being either another view or an initial value. For more details, refer to [START_REF] Herlihy | Distributed Computing Through Combinatorial Topology[END_REF].

Then a complex is a set of simplexes that is closed under inclusion. It captures all considered configurations. Definition 4.2 (Complex). Let and be sets. Then ∈ P ( × ) is a simplicial complex on and (or colored simplicial complex)

• ∀( , ) ∈ × : {( , )} ∈ . • ∀ , simplexes on × : ∈ ∧ ⊆ =⇒ ∈ . The facets of { ∈ | ∀ ∈ : ⊆ =⇒ = }.
The dimension of is the maximum dimension of its facets. is called pure if all its facets have the same dimension.

How can we go from our round-based models, which are generated by graphs, to simplexes and complexes?

Starting with a single graph, we define the uninterpreted simplex induced by this graph. This simplex captures the configuration after a round using graph , simply in terms of who hears from whom. It disregards input values, which makes it uninterpreted. Given a set of graphs representing the possible graphs, we generalize the previous definition to give the uninterpreted complex of . Definition 4.4 (Uninterpreted complex of an oblivious model). Let be an oblivous model defined by a set of graphs . Then the uninterpreted complex of is the complex whose facets are exactly the { | ∈ }. (a) Bipartite graph

( 1 , 1 ) ( 1 , 2 ) ( 2 , 1 ) ( 2 , 2 ) ( 3 , ) (b) Pseudosphere

Uninterpreted complexes of closed-above models

It so happens that closed-above models give rise to uninterpreted complex that are easy to define and study. Indeed, they are unions of pseudospheres, where pseudospheres are colored complexes topologically equivalent to -spheres. These pseudospheres have already been used in the literature to study multiple models of computation [START_REF] Herlihy | Distributed Computing Through Combinatorial Topology[END_REF]Chap. 13].

Definition 4.5 (Pseudospheres [15, Def 13.3.1]). Let 1 , 2 , ..., be sets. Then the pseudosphere complex

(Π; | ∈ [1, ]) • ∀ , ∀ ∈ : ( , ) is a vertex of . • ∀ ⊆ [1, ] : {( , ) | ∈ , ∈ } is a simplex of iff all are distinct.
We can think of these complexes as a generalization of complete bipartite graphs in dimensions. Recall that a complete bipartite graph is a graph that can be split into two sets of nodes, the nodes of each set not linked to each other and each node of one set linked to all nodes of the other set. For example, Figure 3a is a bipartite graph. Now a pseudosphere is the same, except that nodes can be partitioned into sets, no simplex contains more than one element of each set as a vertex, and all the simplexes built from one element of each set are in the complex. Figure 3b is an example of a pseudosphere built from processes 1 , 2 , 3 , and the three sets

1 = { 1 , 2 }, 2 = { 1 , 2 } and 3 = { }.
Among other things, pseudospheres are closed under intersection, and are ( -2) connected.

(Π; | ∈ [1, ]) ∩ (Π; | ∈ [1, ]) = (Π; ∩ | ∈ [1, ]).
One advantage of pseudosphere is that they have high connectivity [START_REF] Herlihy | Distributed Computing Through Combinatorial Topology[END_REF]Def. 3.5.6]. Intuitively, connectivity concerns the (non-)existence of high-dimensional generalisation of holes in the complexes. Since pseudospheres are topologically equivalent to spheres [15, Sect. 13.3], they only have these holes in the highest dimensions.

(Π; | ∈ [1, ]) is ( -2)-connected, where |{ ∈ [1, ] | ≠ ∅}|.
The connectivity of the uninterpreted complex for a simple closed-above model follows, because such a complex is a pseudosphere. Intuitively, for any process , its possible views are exactly the upward closure of its view in the defining graph . Then the -simplexes of the uninterpreted complex are exactly the simplex you can build with one such view for each process.

P

.

• (⊆). Let be a -simplex of . By definition of , it is the uninterpreted simplex of a graph ∈↑ . This in turn means that ∀ ∈ Π :

( ) = ( ) ⊇ ( ). Thus = {( , ( )) | ∈ [1, ]} ⊆ (Π; { | ( ) ⊆ ⊆ Π} | ∈ [1, ]). • (⊇). Let be a n-simplex of (Π; { | ( ) ⊆ ⊆ Π} | ∈ [1, ]). Then ∀ ∈ Π : ( ) ⊇ ( ).
Thus is the uninterpreted simplex of a graph such that ∀ ∈ Π :

( ) ⊇ ( ). We conclude that ∈↑ and thus that ∈ It follows instantly that the uninterpreted complexes of simple closed-above models are (|Π| -2)-connected. 4.9 (C ). Let be a simple closed-above model, and be the graph from which it is built. Then is (|Π| -2)connected.

From this corollary and the closure of pseudospheres by intersection, we now deduce a similar characterization of the connectivity for general closed-above models.

But to do so, we need to first introduce the main tool in our toolbox for studying connectivity of simplicial complexes: the nerve lemma. This result uses a cover of a complex: a set of subcomplexes such that their union gives the initial complex.

Intuitively, the nerve lemma says that if you provide a cover of a complex that is "nice enough", then the connectivity of the initial complex can be deduced from the way that the cover elements intersects. This is usually easier to determine than computing the connectivity directly. [20, T 15.24]). Let be a simplicial complex, ( ) ∈ a cover of and ≥ 0. Then ∀ ⊆ :

( ∈ ) ≥ ( -| | + 1) ∈ = ∅ =⇒ ( is -connected ⇐⇒ N ( | ) is -connected).
Now we can prove the connectivity of uninterpreted complexes for general closed-above models. T 4.12 (C ). Let be a closed-above model, and be the set of graphs from which it is built. Then is (|Π| -2)-connected.

P

. From the proof of Theorem 4.8, we know that is a union of pseudospheres: = ∈ . We want to apply the nerve lemma to this cover. First, by Theorem 4.9, is ( -2)-connected. As for the intersection of any set of , we have two properties. First, it cannot be empty, since all must contains the uninterpreted simplex of the complete graph on Π, by definition of ↑ . This gives us that the nerve complex is a simplex, and thus ∞-connected.

And second, the intersection is also a pseudosphere, by application of Lemma 4.6. Indeed, these are intersections of pseudospheres with the same processes which have an non-empty intersection for each color : the view of this process in the complete graph.

We can thus conclude by application of the nerve lemma and Theorem 4.9.

Interpretation of uninterpreted complexes

We can only go so far with uninterpreted complexes; at some point, we need to consider initial values. Then the same intuition can be applied to a full uninterpreted complex.

Definition 4.14 (Interpretation of uninterpreted complex). Let A be an uninterpreted complex on Π and I be a pure ( -1) complex colored by Π. Then the interpretation of A on I, A (I)

a facet of I a facet of A ( )
These interpretations give us protocol complexes, on which known result on computability are applicable.

A Powerful Tool

On the combinatorial topology front, our results leverage two main tools: the impossibility result on -set agreement based on connectivity [15, Thm. 10.3.1], and a way to compute the connectivity of a complex from the way it is built. This section develops the second idea.

Let A be a pure complex of dimension . We say that A is shellable if there is an ordering 1 , . . . , of its facets such that for every

1 ≤ ≤ -1, =1 ∩ +1
is a pure subcomplex of dimension -1 of the boundary complex of +1 , i.e., of skel -1 +1 . The intuition here is that the complex is the union of simplexes of dimension , and there is an order in which to add simplexes, so that the new simplex is connected to the rest by -1 simplexes, some of its own facets. In the concrete case of 2-simplexes (triangles), they must be connected to the rest by 1-simplexes (edges).

Here, unions and intersections apply to the complexes induced by the facet and all its faces. Such a sequence of facets is a shelling order of A.

For example, the complex in Figure 4a is shellable, but the one in Figure 4b is not.

Given a shelling order 1 , . . . , of a complex A, ( =1 )∩ +1 is the union of the complexes induced by some ( -1)-faces 1 , . . . , of +1 , by definition of shellability. Each is a face of a facet of ∪ =1 , hence +1 and share a ( -1)-face. Then,

=1 ∩ +1 = =1 ( +1 ∩ ).
We also use the following technical result.

2]).

Let A be a pure ( -1)-dimensional sub-complex of the boundary complex of a simplex of dimension . Then A is shellable, and any sequence of its facets is a shelling order for A.

Finally, we rely on the straightforward corollary of the nerve lemma for a cover with two elements.

C 4.16 (T ).

Let and be -connected complexes. If ∩ is ( -1)-connected, then ∪ is -connected.

We now state the main technical result of the section. It extends the result from Castañeda et al. [START_REF] Castañeda | A topological perspective on distributed network algorithms[END_REF] and adapts it to the interpretation of complexes we need here. While Castañeda et al. studied the complex given by the interpretation of a single graph (to capture models like LOCAL and CON GEST [START_REF] Peleg | Distributed Computing: A Locality-sensitive Approach[END_REF]), we care about the complex resulting of the interpretation of a set of graphs.

We thus send each input simplex into a complex, and show that if both the output complexes and the mapping are "nice", the interpreted complex is highly connected. L 4.17. Let A be a pure shellable complex of dimension , B a complex, (B ) ∈ a cover of B, and ℓ ≥ 0 an integer. Suppose that there is a bijection between the facets of A and the elements of (B ) ∈ such that:

(1) For every facet ′ of A and every pure -subcomplex =1 ⊆ A satisfying that =1 ∩ ′ = =1 ( ′ ∩ ) for some of A's facets 1 , . . . , , with each and ′ sharing a ( -1)-face, it holds that

=1 ( ) ∩ ( ′ ) = =1 ( ( ′ ) ∩ ( )).
(2) For every ≥ 0 and every collection 0 , 1 , . . . , of + 1 facets of A with each and 0 sharing a ( -1)-face, it holds that =0 ( ) is least (ℓ -)-connected.

Then, B is ℓ-connected.

ONE ROUND LOWER BOUNDS FOR ONE ROUND: A TOUCH OF TOPOLOGY

As before, we start with the simple closed-above case, where the model is the closure of a single graph. In this case the tight lower bound follows from Castañeda et al. [START_REF] Castañeda | A topological perspective on distributed network algorithms[END_REF]Thm 5.1], as mentionned above.

T 5.1 (L ). Let a simple closed-above model generated by the graph . Let ≤ ( ). Then -set agreement is not solvable on in a single round.

We thus focus on general closed-above models. Here we have to leverage the underlying structure of the protocol complex. We do so through two tools : the main theorem from Section 4.4, as well as two graph parameters: the equal-domination number over a set of graphs, and the max-covering numbers of a set of graphs.

MULTIPLE ROUNDS

Given that we focus on oblivious algorithms, a natural approach to extending our lower bounds to the multiple rounds case is to look at the product of our graphs. By product, we mean the graph of the paths with one edge per graph. Thus the products of graphs capture who will hear who after corresponding communication rounds. Definition 6.1 (Graph path product). Let and be graphs with auto-loops (∀ ∈ Π : ( , ) ∈ ( ) ∧ ( , ) ∈ ( )). Then their graph path product the graph (Π, ) such that ∀ , ∈ Π : ( , ) ∈ Π =⇒ ∃ ∈ Π : ( , ) ∈ ( ) ∧ ( , ) ∈ ( ).

Since we have a graph as the result, we can apply our lower bounds for one round. At least, if the resulting graph still satisfy the hypotheses of our lower bounds. It does, although product doesn't maintain closure-above. This subtlety is explained in the next subsection.

Closure-above is not invariant by product, but its still works

What is the pitfall mentionned above? Quite simply, that the product of two closed-above models does not necessarily gives a closedabove model. This follows from the fact that the closure-above of a product of graphs doesn't always equal the product of the closure-above of the graphs. Let's take an example: the product of a cycle with itself. Then we cannot build the following graph by extending the cycles and taking the product: Why? Simply put, adding the new edge to either of the two cycles necessarily creates other edges in the product. Adding an edge from 2 to any other node than 3 and 4 also creates new edges; so does adding an edge to 4 and then an edge from 4 to 6 , or an edge from 3 to 6 in the second graph.

Hence the product of the closure above of this cycle with itself is not the closure-above of the squared cycle. To put it differently, closure-above is not invariant by the product operation.

Nonetheless, the bell does not toll for our hopes of extending our properties. What is used in the lower bound proofs above is not closure-above itself, but its consequences: being a union of pseudospheres containing the full simplex, such that for each pseudosphere, all graphs contain the smallest graph.

All three properties are present in a specific subset of the product of two simple closed-above models: all products where edges might be added to the last graph in the product but not to the other. Each added edge only alters the view of its destination, since it is in the second graph, and multiple added edges don't interfere because they are all added to the same graph. Hence we can change the views of processes one at a time, and thus we get a pseudosphere. Since adding no edge gives the original product and adding all missing edges gives the clique, we get the other two properties. Then taking this subset of the product of two general closed-above models result in a union of pseudosphres, one for each product of the underlying graphs.

Therefore we can extract relevant subcomplexes from the product of closed-above models, and then the lower bounds only depends on the properties of the underlying product of graphs.

Upper bounds for multiple rounds

Even if we just explained how to deal with lower bounds for multiple bounds, we still start by giving upper bounds for multiple rounds. This is for the same reason as in the one round case: the upperbounds require no combinatorial topology, and they allow us to introduce concepts needed for the lower bounds.

First, we need to prove a little result that is enough for our upper bounds: that the product of closed-above models is included in the closure-above of the product. L 6.2 (P ). Let and be two graphs. Then ↑ ↑ ⊆↑ ( ).

P . Let ∈↑ ↑ . Thus ∃ ′ ∈↑ , ∃ ′ ∈↑ : = ′ ′ . Let , ∈ Π such that ( , ) ∈
. We show that ( , ) ∈ ; this will entail that ∈↑ (

). Because ( , ) ∈ , ∃ ∈ Π : ( , ) ∈ ( ) ∧ ( , ) ∈ ( ). But ′ ∈↑ and ′ ∈↑ , therefore ( , ) ∈ ( ′ )∧( , ) ∈ ( ′ ). We conclude that ( , ) ∈ ′ ′ = .

What this means is that taking the closure-above of the products of our graphs over-approximate the actual model after rounds. And thus, algorithms working on these approximations work on the actual model. Now, let us start with simple closed-above models. Just like for the one round case, they are completely characterized by the domination number of their underlying graph. T 6.3 (U ( ) ). Let be a simple closed-above model defined by the graph . Let > 0. Then ( )-set agreement is solvable in rounds in .

P

. We have that ( )-set agreement is solvable on ↑ by Theorem 3.2. This then implies by Lemma 6.2 that it is solvable on (↑ ) , that is on . But for general closed-above models, one cannot use the domination number itself, because one cannot know which of the underlying graphs will be there. As in the one round case, we use the equal-domination number and covering numbers. ). Let be a general closed-above model generated by the set of graphs . Let > 0. Then ( )-set agreement is solvable in rounds on .

P

. We have that ( )-set agreement is solvable on ↑ , that is on .

One issue with these bounds is that they require the computation of possibly many products, as well as the computation of the combinatorial numbers for a lot of graphs. One alternative is to forsake the best bound we can get for one that can be computed using only the numbers for the initial graphs.

This hinges on covering number sequences. Recall that the -th covering number of a graph is the minimum number of processes hearing a set of processes that do not broadcast. In a sense, it gives the guaranty of propragation of information by a set of processes.

That's the whole story for one round. But what happens when you do multiple rounds? Then, if the -th covering number of the graph is greater than , this means that in the next rounds, the minimum number of people who will hear the value of the initial processes is the -th covering number. And if this number is greater than , this repeats. Covering number sequences capture this process. One can also see them as the sequences of covering numbers for powers of the graph. Definition 6.6 (Covering number sequences). Let be a graph. Then the -th covering numbers sequence of ( ) ∈N * such that 1 = ( ) and ∀ ≥ 1 :

+1 = |Π| if ≥ ( ) ( ) if < ( )
Armed with these sequences, we get an upper bound directly from . T 6.7 (U ). Let be a simple closed-above model defined by the graph on Π. Then if the -th covering sequence of reaches at some point, -set agreement is solvable on the model .

We can adapt this bound for general closed-above models by generalizing the covering numbers sequences to a set of graphs. Definition 6.8 (Covering numbers sequences for sets of graphs). Let be set of graphs. Then the -th covering numbers sequence of ( ) ∈N * such that 1 = min ∈ ( ) and

∀ ≥ 1 : +1 = if ≥ max ∈ - ( ) min ∈ ( ) if < max ∈ - ( )
T 6.9 (U ). Let be a set of graph on Π. Then if the -th covering sequence of reaches at some point, -set agreement is solvable on the oblivious closed-above model generated by .

P

. If the -th covering number sequence of reaches after step , this means that every set of processes is heard by everyone after rounds. In particular, the processes with the smallest initial values will be heard by everyone.

Hence sending all the values heard for now for rounds, and then deciding the smallest value received, ensures that one of the -th smallest values will be chosen, and thus solves -set agreement. As a concrete applications of these bounds, we consider a classical family of subgraphs: stars. Definition 6.12 (Star graphs). Let be a graph. Then is a star graph ∃ ⊆ Π : = ( , × Π). T 6.13 (L ). Let be the set of graphs which are unions of stars with different centers. Then -set agreement is not solvable in the closed-above model generated by .

Lower bounds for multiple rounds

CONCLUSION

We provided upper and lower bounds on -set agreement for closedabove models, the subset of round-based models defined by subgraphs that must be present in the communication graph at each round. These models encompass many message-passing models of distributed computing focused on safety properties. Regarding the bounds themselves, although their proofs leverage combinatorial topology, all our bounds are expressed in terms of combinatorial numbers of the graphs. That is, these bounds can be used without any knowledge of combinatorial topology. Yet combinatorial topology was instrumental in showing such sweeping results.

Definition 2 . 3 (

 23 Closed-above communication models). Let be an oblivious communication model. Then is closed-above ∃ We call the graphs in the generators of . If is a singleton, then is simple closed-above.

Definition 3 . 6 (

 36 Covering numbers of a set of graphs). Let be a set of graphs. Then ∀ < ( ), its -th covering number (

T 3 .

 3 7 (U ). Let be a set of graphs. Then ∀ ∈ [1, ( ) [: ( + ( -( )))-set agreement is solvable on the oblivious closed-above model generated by .

Figure 1 : 3 . 8 .

 138 Figure 1: Two examples communication graphs

Definition 4 . 3 (

 43 Uninterpreted simplex of a graph). Let be a graph. Then the uninterpreted simplex of is the colored simplex {( , ( ) | ∈ Π }.

Figure 3 :

 3 Figure 3: A bipartite graph and a pseudosphere

  Let be a simple closed-above model, and be the graph from which it is built. Then= (Π; { | ( ) ⊆ ⊆ Π} | ∈ [1, ]).

C

Definition 4 . 10 (

 410 Nerve complex). Let be a simplicial complex, ( ) ∈ a cover of . Then the nerve complex of this cover, N ( | ) the complex generated by• the vertices are the ; • and the simplexes are the sets { | ∈ } for ⊆ such that

Definition 4 . 13 (

 413 Interpretation of uninterpreted simplex). Let be an uninterpreted simplex on Π and be a ( -1)-simplex colored by Π. Then the interpretation of on , ( ) {( , ) | ∈ Π ∧ ( ∈ =⇒ (∃ ∈ ( ) : = ( )))}

Figure 4 :

 4 Figure 4: Examples of a complex that is shellable and one that is not

  Let > 0 and let a simple closed-above model generated by the graph . Then ( ( ) -1)-set agreement is not solvable on in rounds by an oblivious algorithm.

  Let > 0 and let be a closed-above model generated by the set of graphs . Let = min( ( )-2, min{ + ( )-2 | ∈ [1, ( )-1]}) Then ( + 1)-set agreement is not solvable on in rounds by an oblivious algorithm.

ACKNOWLEDGMENTS

Adam Shimi was supported by the Agence Nationale de Recherche under Grant No.: PARDI ANR-16-CE25-0006 and Armando Castañeda was supported by project PAPIIT IN108720.

Definition 5.2 (Distributed domination number of a set of graphs). Let be a set of graphs. Then the distributed domination number , ( )

The difference between ( ) and ( ) is that a set of ( ) processes dominates each graph of separately, whereas a set of ( ) processes might not dominate any graph of , but it dominates every subset of graphs of together. Thus ( ) ≤ ( ). Fitting, considering the former is used in lower bounds and the latter in upper bounds.

Next, the max-covering numbers are quite subtle. For < ( ), the -th max-covering number of is the maximum number of processes hearing a set of procs, summed over graphs in .

That is, the max-covering numbers capture how much values can be disseminated in the best case. They serve in lower bounds by giving a best case scenario on which we can focus to prove impossibility. Definition 5.3 (Max covering numbers of a set of graphs). Let be a set of graphs and < ( ). Then the -th max-covering number of , max-cov ( ) max

We also define the -th max-covering coefficients on , ( ) The term depending on ( ) in the lower bound serves when the max-covering numbers are not sufficient to distinguish adversaries with different properties. Consider for example the symmetric models of all unions of stars, with ≤ . Then for those graphs, for < ( ), we have -( ) = t, and thus ( ) = -. Hence the minimum over the + ( ) -2 is -2.

But this would mean that ( -1)-set agreement is impossible for < , whereas we can clearly solve 2-set agreement for = -1, for example. What depends on is ( ) itself. More precisely, ( ) = -+ 1, because given , we can consider only the graph where the centers of stars are in Π \ , up until the point where | | > -.

Hence our lower bound shows that for the symmetric union of stars, ( -)-set agreement is impossible in one round. Given that our upper bounds above tell that ( -+ 1)-set agreement is possible in one round for this model, the bound is tight.

Finally, the bound can be specialized for symmetric models.

+ --1 ( -

}) = Then ( + 1)-set agreement is not solvable on (↑ ) in a single round.

Notice that all these lower bounds are valid for general algorithms, not only oblivious ones. The reason is that a one round full information protocol is an oblivious algorithm.


Xavier Thirioux
email: xavier.thirioux@enseeiht.fr

Alexis Maffart
email: alexis.maffart@gmail.com

Taylor Series Revisited

Keywords: Taylor expansion, Certification, PDE _

We propose a renovated approach around the use of Taylor expansions to provide polynomial approximations. We introduce a coinductive type scheme and finely-tuned operations that altogether constitute an algebra, where our multivariate Taylor expansions are firstclass objects. As for applications, beyond providing classical expansions of integro-differential and algebraic expressions mixed with elementary functions, we demonstrate that solving ODE and PDE in a direct way, without external solvers, is also possible. We also discuss the possibility of computing certified errors within our scheme.

Motivations

Taylor Expansions

Our principal motivation is to provide an automatic way of approximating arbitrary multivariate numerical expressions, involving elementary functions, integrations, partial derivations and arithmetical operations. In terms of features, we propose an approach where Taylor expansions are first-class objects of our programming language, computed lazily on demand at any order. Finally, we also wish to obtain certified errors, which will by the end include errors of approximation and numerical errors, expressed in any suitable user-provided error domain, such as zero-centered intervals, intervals, zonotopes, etc. From a user's perspective, a typical workflow is first to compute a certified approximation at some order of some expression, second to evaluate the maximum error for the given domains of variables, and maybe third to compute a finer approximation at some higher order (without recomputing previous values) if the error is too coarse, and so on, until the approximation meets the user's expectations in terms of precision. We postulate that the expressions at hand are indeed analytical and possess a valid Taylor expansion around a given point and within variables' domains. If it is not the case, then the error computed at every increasing order won't show any sign of diminishing and could even diverge. Last but not least, our approach yields a direct means to express solutions to ODEs and PDEs and thus solve them, without complex numerical methods based on domains discretization.

Furthermore, we aim at bringing as much robustness and correction as possible to our library through a correct-by-construction approach. The type system is in charge of the correction as it ensures, at compile time, that dimensions of various tensors, functions, convolutions and power series conform to their specifications. This is of a particular importance in a complex and error-prone context involving a vast number of numerical computations such as ODEs and PDEs resolution. The type system which validates all dimension related issues greatly helps in reducing the focus on purely numerical concerns: correctness of approximation, precision, convergence. Moreover, correction could be proved more formally with a proof assistant such as COQ. This idea could be addressed in the future even if this work is likely to be laborious.

As a disclaimer, the current state of our contribution doesn't allow yet the computation of certified errors in the presence of differential equations, so we mainly focus here on infinite Taylor expansions without remainders. Still, as one of our prominent future goals, certified errors were taken into account in the design stage of our framework and we discuss them along this paper.

Applications

Among many possible applications, we more specifically aim at formally verifying systems dealing with complex numerical properties, such as controllers for embedded systems. Moreover, through certified integration of ODE, we may also consider hybrid systems, such as a continuous plant coupled to a discrete controller.

Outline

We start by recalling some related works around formalization and mechanization of Taylor expansions in Sect. 2. Then, we state a mathematical formulation of our on-demand multivariate Taylor expansions with errors in Sect. 3 before introducing our implementation of data structures and operations that form an algebra in Sect. 4. We separately discuss the more complex case of composition in Sect. 5. In Sect. 6, we present some experiments done on solving differential equations in a direct way. Finally, we open up some perspectives, notably about errors, then conclude, respectively in Sects. 7 and 8.

Related Works

Taylor Series

Although Taylor expansions are well known and form a very rich and interesting algebra, their realizations as software items are not widespread. From a mathematical perspective, some weaknesses may explain this lack of success: they only support analytical functions, a rather limited class of functions; they don't possess good convergence properties, uniform convergence is hardly guaranteed for instance; typical applications for polynomial approximations are usually not concerned with certified errors, mean error or integrated square error (through various norms) are more important and don't easily fit into Taylor expansion schemes. Finally, from a programming perspective, Taylor expansions are: hard to implement as they require many different operations to be implemented, from low-level pure numbers to high-level abstract Taylor expansions seen as first-class citizens; error-prone with lots of complex floating-point computations on non-trivial data structures; heavily resource demanding in our multidimensional setting because data structures rapidly grow as the precision order increases.

Here are a few works dealing with Taylor expansions. In [START_REF] Karczmarczuk | Functional differentiation of computer programs[END_REF], the author presents an early application of laziness to cleanly obtain Taylor polynomial approximations. Laziness allows to augment the degree of the resulting polynomial on demand. Yet, the setting is much simpler as it is strictly one-dimensional and certified errors are not in scope. With these restrictions, the author obtains nice formulations of automatic differentiation and polynomial approximations of classical phenomena in physics. Speaking about implementation, related works come in many flavors and date back to the now well established folklore of automatic differentiation (forward or backward modes). As for symmetric tensor algebra, which forms a well-suited representation basis for partial derivatives, a huge menagerie of (mostly C++) libraries exists, for tensors of arbitrary orders and dimensions (but some libraries put a very low upper-bound on these values). These implementations are clearly not oriented towards reliability and proof of correctness, but towards mere efficiency. This also comes at the expense of some user-friendliness, as memory management and user interface are more complex and error-prone than in our own library. Still, we may consider interfacing our code base with a trusted and stable tensor library, for much better performance.

One of the most prominent implementation of Taylor expansions is the COSY tool, cf. [START_REF] Makino | Rigorous integration of flows and ODEs using Tylor models[END_REF]8]. This tool has been used in industrial-scale engineering and scientific contexts, to modelize and predict the complex dynamics of particles in accelerators for instance. This tool supports 1D Taylor expansions with intervalbased certified errors. Polynomial degree is not refinable on demand and Taylor expansions are not handled per se (i.e. not first-class citizens). The authors managed anyway to implement an error refinement scheme for solved form ordinary differential equations, that allows solving them with tight certified errors. Experiments show that this tool compares favorably to other traditional approximations and bounding techniques, such as branch-and-bound approaches and interval arithmetics, in terms of speed and precision. We also aim at implementing differential equation solving in our multi-dimensional setting.

At the other end of the spectrum, [START_REF] Martin-Dorel | Certified, efficient and sharp univariate Taylor models in COQ[END_REF] proposes correct-by-construction univariate Taylor expansions with certified errors, which appears as a huge step. Integration of floating-point errors into this scheme is also a concern addressed in [START_REF] Martin-Dorel | Formally verified certificate checkers for hardest-to-round computation[END_REF]. Still, apart from its limitation to the 1D case, this approach suffers from weaknesses: expansion degree is fixed and differential equations cannot be handled. The underlying algorithm won't be so easily turned into a co-inductive (lazy) equivalent version.

And in the middle of the spectrum comes [START_REF] Pearlmutter | Lazy multivariate higher-order forward-mode AD[END_REF], where the author defines a way to handle multivariate Taylor series and presents its implementation featuring on demand computation thanks to Scheme laziness. The few points he did not implement and that we will try to cope with in our library are: errors certification which is not handled and efficiency which is not optimal. For instance, the author's method to multiply multivariate power series is to define a generic composition between a bivariate function and a power series and to instantiate it with the multiplication. This method is simply built upon the chain rule but has some drawbacks. First, the generic equation given can usually be drastically simplified for instance in the case of multiplication and second, such a generic scheme implies that some parts of the resulting coefficients will be computed several times differently. Conversely, in our solution, the pervasive multiplication operation is implemented with a strong concern on optimality.

Our work and specifically our data-structure is based on the dissertation [9, Part 2], with the nuance that a single unbounded tree will be used instead of an infinite sequence of finite trees, each such tree representing a symmetric tensor of a given order. This choice notably enables the resolution of partial differential equations, which was impossible in the setting of [START_REF] Thirioux | Verifying embedded systems[END_REF].

Differential Equations

Iterative methods are pervasive in integrating differential equations because they often provide an efficient way to find an approximation of an ODE solution. Some of them own validation aspects, such as [START_REF] Sandretto | Validated explicit and implicit Runge-Kutta methods[END_REF] which relies on Runge-Kutta method to integrate ODE with a numerical validation. The main difference between these methods and our work as a direct method is that we don't need these next level iterations. We are able to yield a result in the equivalent of the first iteration.

Formalization

We recall the canonical presentation of a multivariate Taylor expansion at order R in dimension N . This expansion converges to f (x) when R → +∞ for an analytical function f only in a chosen neighbourhood of point 0.

f (x) = |α|<R D α f (0) • x α α! + |α|=R D α f (λ * x) • x α α!
In the above formulation, x = (x 0 , . . . , x N -1 ) ∈ R N , α = (α 0 , . . . , α N -1 ) ∈ N N indexes the derivation order of f in the symmetric tensor of partial derivatives D α f and λ ∈ [0, 1] is an unknown coefficient that characterizes the exact Taylor remainder. We have to compute derivatives both at point 0 for the polynomial part and at point λ * x for the error part. We choose to use a single co-inductive data-structure that encodes all possible derivatives, indexed by some α. As for the elements of this structure, we handle value, error pairs. Our framework is error-agnostic as the value-error domain is user-defined and only requires arithmetical operations. Several solutions are available in the literature: zero-centered intervals, intervals, zonotopes, etc. In the remainder, we only assume that elements of our structures form an algebra (including addition, multiplication and some elementary functions), disregarding whether they are pure values or values with errors.

This co-inductive structure, that we coin a "cotensor", enables to compute finer approximations on demand and also to lazily represent expansions of solutions to ODEs and PDEs, when they are expressed in solved form, i.e. not implicit (as it would for instance be the case if the solution were specified as a zero of a polynomial form in a functional space).

4 An Algebra of Taylor Series

Data Structure

Coefficients are present in each node of a unique tree structure and are written as s o 0 ,...,o N -1 where every o i is the number of occurrences of the variable x i in the path that leads to the considered coefficient

s o 0 ,...,o N -1 .
The principle is quite simple: at each node, we choose either to keep the same variable accounting for the final Taylor series, or we drop it and repeat the same process for lower dimension variables. This is pictured in tree branches of the following example as x i for the first case and x i for the second case. The variable at the root of the tree is X n if the dimension is n + 1. This tree is developed below and represents a symmetric cotensor s of dimension 4: s 0,0,0,0 s 0,0,0,0 s 0,0,0,0 s 0,0,0,0 s 1,0,0,0 ...

s 0,1,0,0 s 0,1,0,0 ... s 0,2,0,0 ... ... s 0,0,1,0 s 0,0,1,0 s 0,0,1,0 ... s 0,1,1,0 ... ... s 0,2,0,0 s 0,2,0,0 ... ... ... s 0,0,0,1 s 0,0,0,1 s 0,0,0,1 s 0,0,0,1 ... s 0,1,0,1 ... ... s 0,0,1,1 s 0,0,1,1 ... ... ... s 0,0,0,2 s 0,0,0,2

Structural Decomposition

We will introduce for this co-inductive structure a few notations inspired from the computation of the quotient and the remainder with respect to variable X n . We will call a left cotensor a cotensor which is the left branch of another cotensor and we will denote L n+1 the set of left symmetric cotensors and R n+1 the set of right symmetric cotensors in dimension n + 1. If V is the set of labels at the root of the tree, we have the following definitions:

L n+1 L n + X n .R n+1 R n+1 L n + X n .R n+1 + V Hence : R n+1 = L n+1 + V
We note that the only difference between left and right cotensors is the constant part v ∈ V and from now, we are going to consider that right case is the general one and that left case is the specification of the right case with constant part equal to 0. This will prevent us from writing similar redundant equations for all algebraic operations we will describe later. A cotensor is, then, considered a right cotensor by default, even if it has no parent because it contains a significant value v ∈ V which is the constant part of the Taylor series. It comes then that a tree is interpreted as a Taylor series by adding together the term for the left tree, X n times the term for the right tree and the label value of the root.

Implementation

Finally, in terms of OCaml implementation, this decomposition scheme naturally translates into the slightly relaxed following type definition, where L n and R n have been conflated in a single type: Here the type of symmetric cotensors tree has two type parameters: the type of elements 'a and the dimension type 'n. The last parameter not being constant through recursion, it appears as _ in the type declaration. Then, the two cases for the dimension N : N = 0 and N = 0, are respectively handled with Nil and Leaf/Node constructors. Leaf is only a special case of Node where all the coefficients are zeros. Handling this particular case with a different constructor aims at saving some computations, for instance all polynomial forms will be represented by finite trees, not by unbounded ones with trailing zeros. And Leaf constructor is used to mark the end of a branch when the dimension has decreased to 0, namely all the variables has been consumed. Type parameters of constructors' arguments behave accordingly to the decomposition of R n+1 .

The Nat.zero and Nat.succ type constructors encode the dimensions of manipulated cotensors, as we use GADT 1 allowed by OCaml. We use a standard type-level encoding of Peano numbers and operations that we don't detail here. We hereby enforce a correct-by-construction use of our data-structures.

Component-Wise Operations

From this section onward, we assume cotensor elements form a field, with arithmetical operations on it. It may be in practice a field of coefficients or/and errors. These elements are denoted by V A and V B . Functions "λ. " and " + " straightforwardly witness the vector space structure of cotensors. The Hadamard product " ⊙ " is the component-wise product of two cotensors of same dimension. Hence with the notation

A n+1 A L n + X n .A R n+1 + V A : A n+1 + B n+1 = (A L n + B L n ) + X n .(A R n+1 + B R n+1 ) + (V A + V B ) λ.A n+1 = λ.A L n + λ.X n .A R n+1 + λ.V A A n+1 ⊙ B n+1 = (A L n ⊙ B L n ) + X n .(A R n+1 ⊙ B R n+1 ) + (V A * V B )

Multiplication

Let us define a new notation for cotensors in order to specify the multiplication.

We are now going to consider that the error term is no longer separated in a precise term of the equation but is distributed in all the terms of the equation. Which gives:

S(X 0 , ..., X N ) = (S 0 + S 1 ⊙ X + S 2 ⊙ X 2 + ... + S m ⊙ X m + ...) shortened in = (S 0 + S 1 X + S 2 X 2 + ... + S m X m + ...)
where X = (X 0 , ..., X N )

This notation is inspired by derivation order; even if we do not consider order of cotensors; because it will be of a great help when defining the multiplication and introducing the convolution product. Product of Taylor expansions is really pervasive and appears in many operations (derivation formulas, composition of Taylor series, etc.). It is naturally defined with an explicit convolution. Concretely:

S(X 0 , ..., X N ) × T (X 0 , ..., X N ) = (S 0 + S 1 X + ... + S p X p + ...) × (T 0 + T 1 X + ... + T q X q + ...) = R 0 + R 1 X + R 2 X 2 + ... + R k X k + ... where ∀k ∈ N, R k = k i=0 S i T k-i
1 Generalized Algebraic Data Types.

To compute the coefficients at order k, we need to consider every product that will produce an order k, i.e. every coefficient of order i by every coefficient of order ki, i ranging from 0 to k.

In our setting, we maintain a typed convolution structure to express computation of the term k i=0 S i T k-i . This structure, while geared towards static guarantees and proof of correctness, still allows for some efficient implementation. Informally, we may specify our structure as an array containing couples of cotensors of a specific dimension and that will represent absolute paths. The same structure is used to represent relative paths. We introduce a path notation, illustrated by the following examples in dimension n:

-() is the considered tree -(n) is the tree we get when we take the n-th variable (X n ) once in the considered tree -(n.n.n-1) when we take the X n variable twice and then X n-1 once

The so called "considered tree" is the original tree given in parameter if considering absolute paths or a specific tree (descendant of the original one) if considering relative paths. Through the relative paths (left part of the semi-colon), we will store the number of times we went down a right branch since the last left branch, namely relative paths are about the current variable and absolute paths are about all previous variables with respect to the order. Initially, the structure contains a couple of the two original trees given in parameter for both absolute paths (and the part for relative paths is empty):

() ; () () ()
Then, at each step of the algorithm:

-If the current node is a right branch, we will update the relative paths by adding the current node (k.k here) and shifting the lines as follows:

() (k) ; ... becomes () (k) (k.k) ; ... (k) () ... (k.k) (k) () ...
-if the current node is a left branch, we will combine the relative paths with the absolute ones, store the result as the new absolute paths and empty the new relative paths:

() (n-1) ; () (n) becomes () ; () (n) (n-1) (n-1.n) (n-1) () (n) () () (n-1.n) (n-1) (n) ()
Folding this structure to compute a term of a product simply consists in combining relative paths with absolute paths, multiplying cotensors roots column-wise and then summing these intermediate results altogether. Associating a relative path to an absolute one means concatenating them. Speaking in terms of trees, it means that the relative path begins where the absolute one ends in the tree.

Differential Operations

Cotensors of dimension N may not only be structurally decomposed on X N -1 but also on any other X k , which we would call a non-structural decomposition. For that purpose, the " [ ]" function specializes a cotensor, i.e. drops some index by specializing it to a specific dimension k, and therefore represents the division by a monomial X k . Conversely, the " ↑ " function represents the multiplication by a monomial X k . For a cotensor of dimension N , they are defined in terms of polynomials as:

(S[k])(X 0 , . . . , X N -1 ) S(X 0 ,...,X N -1 )-S(X 0 ,...,X k-1 ,0,X k+1 ,...,X N -1 ) X k (S ↑ k)(X 0 , . . . , X N -1 ) X k .S(X 0 , . . . , X N -1 )
Using the same notations as for component-wise operations, we show how these operators simply fit the structural decomposition:

S[k] = (S L + X N -1 .S R + V S )[k] =      S R , for k = N -1 S L +X N -1 .S R -S L |X k ←0 -X N -1 .S R |X k ←0 +V S -V S X k = S L [k] + X N -1 .S R [k], for k < N S ↑ k =    0 + X N -1 .S, for k = N -1 (S L + X N -1 .S R + V S ).X k = S L ↑ k + X N -1 .(S R ↑ k) + V S .X k , for k < N
Differential operations introduce partial differentiation and integration in the cotensor algebra. These differentiation and integration operators respectively refer to S[ ] and S ↑ . They also use the cotensor of integration/derivation factors "∆ k ", where the o i are the variable occurrence number, such that:

(∆ k ) (o 0 ,...,o N -1 ) 1 + o k , for i o i = R dS(X 0 ,...,X N -1 ) dX k S[k] ⊙ ∆ k X k 0 S(X 0 , . . . , x k , . . . , X N -1 )dx k (S ⊙ ∆ -1 k ) ↑ k 5
The Composition Operator

Differential Method

Principle. The Taylor series algebra with the previous operations still remains basic, and that is why we are now interested in composing Taylor series with elementary functions. To do so, we only need to apply elementary functions to arbitrary arguments, i.e. to compose univariate Taylor series with multivariate ones. A general composition scheme of Taylor series is also possible in our setting but out of the scope of our current concerns. This method lies on a differential decomposition, namely a function is the sum of the integrals of its derivatives with respect to all its variables, plus a constant term:

H : R N → R, H = H(0) + i<N X i ∂H ∂X i X k =0 k>i dX i
Example. We need to partially evaluate the derivatives at 0 to avoid counting several times the parts shared by different variables, as illustrates the following concrete example: let

F : R 3 → R, F (x, y, z) = x 3 + 2x 2 y + xz + 5y 2 + 3yz 2        ∂f ∂x = 3x 2 + 4xy + z x 0 ∂f ∂x dx = x 3 + 2x 2 y + xz ∂f ∂y = 2x 2 + 10y + 3z 2 y 0 ∂f ∂y dy = 2x 2 y + 5y 2 + 3yz 2 ∂f ∂z = x + 6yz z 0 ∂f ∂z dz = xz + 3yz 2
The blue terms are redundant and that is why we have:

F (x, y, z) = F (0, 0, 0) + x 0 ∂f ∂x dx + y 0 ∂f ∂y x=0 dy + z 0 ∂f ∂z x=0 y=0 dz
Composition. As we are in the specific case of composition, we will use the classic chain rule:

∂(f • g) ∂X i i<N = ( ∂g ∂X i ) i<N × (f ′ • g)
Hence :

f • g = f • g(0) + i<N X i ( ∂g ∂X i × f ′ • g) X k =0 k>i dX i
The computation of the partial derivatives ∂(f •g) ∂X i i<N is done case by case with respect to the elementary function f at use, each such function having a well-known derivative f ′ . The cases where f = exp, sin, cos, log, atan, x a , . . . are easily handled. So, according to the above equation, we only need to partially evaluate these derivatives, to integrate them then and to finally sum the results. This method will bring us satisfying results as detailed below, but one must bear in mind that despite the method is very short in terms of code and then easily implemented, it is not optimal in terms of computation. This differential method for the composition is not canonical in that it does not compute the minimum number of operations to produce the coefficients of the result. As a witness of non canonicity in the definition of composition, the ∆ k coefficients will be used for multiplication and division consecutively, which could be avoided. Besides, as long as we do not handle certified errors, the method does not need an additive decomposition of f but it will be the case as soon as we handle the errors and we will have to deal with this constraint.

Elementary Functions

Elementary functions, limited to one argument functions, are specified as univariate Taylor series. Therefore, as only one branch of the cotensor will be meaningful, such series are treated separately. This is only a matter of efficiency and obviously not mandatory. To obtain a Taylor expansion of an elementary function, we need to be able to compute any n-th derivative. Taylor series for elementary functions are well known, so the first way to produce such a series is to compute the coefficients iteratively and lazily with respect to the known formulas, such as the following ones:

exp(x) = i∈N x i i! log(1 + x) = i∈N -(-x) i i (1 + x) p = i∈N p i x i sin(x) = i∈N (-1) i (2i+1)! x 2i+1 cos(x) = i∈N (-1) i (2i)! x 2i
Similar formulations are available for elementary functions not presented here.

Experimentation

Now that the main operations are available in our algebra, we can start using it. Differential equations are pervasive in dynamical systems and our point is to propose a direct (i.e. non-iterative) way to solve them. By direct method, we mean that coefficients are computed once and for all and therefore there is no need to iterate over their values until a specific precision is reached. Precision in our case is seen differently: coefficients are computed only once and if the user wants a finer precision, the user will increase the order of derivation which means that new and deeper coefficients will be computed.

Airy Equation

To illustrate this direct approach for solving ODEs and PDEs, we will use the first dimension Airy equation which stands as follows:

f ′′ -xf = 0
As the equation contains a second derivative, we split it for convenience in two first order equations introducing f dot as f derivative:

       f dot = f dot 0 + x xf f = f 0 + x f dot
Then, thanks to OCaml laziness, we express and solve this mutually recursive system directly, with the following principle:

-According to the second equation, computing the first coefficient of f , the constant part, means summing the constant part of f 0 with the constant part of x f dot. We know that the constant part of an integral will be 0, whatever the integrand is. -the first coefficient of f dot, or equivalently the second coefficient of f , is computed the same way (no need to evaluate the argument of the integral). -then the mutual recursion works and the third coefficient of f , or the second one of f dot, is simply the result of integrating the constant part of xf , actually 0. The other coefficients are also computed in finite time.

So the trick is to stay a step ahead by computing a first coefficient of a recursive Taylor series without having to evaluate itself, thanks to the integral operator, and then to keep this advance all along the computation so that the recursion will always end. Indeed, if the computation scheme respects the causality, for example in one dimension: computing a coefficient requires only strictly lower order coefficients, then we can ensure the recursion will end.

Once we get the solution up to a specific order, we evaluate it as a polynomial function so that we can draw its graph (Figs. 1 and2): We can observe that the approximation is reliable on a specific interval and diverge outside of it. We can have this conclusion because we know the theoretical result in this case, but we won't know it in most cases. This is what will motivate the necessary handling of certified errors. Intervals of errors, which are only an example of error representation, will give the user information about how far the theoretical function could be from the returned approximation.

Heat Equation

In order to explain the principle of causality more precisely and to show a more general case, we are going to present the 2-dimensional heat equation example:

∂u ∂t = α ∂ 2 u ∂x 2
There are 2 different ways of integrating this equation and we chose to integrate it with respect to variable t so that initial conditions are a function of variable x at initial time t = 0. Here is the new form of the equation:

u(x, t) = u 0 (x) + α × t ∂ 2 u(x, t) ∂x 2
where u 0 (x) will be a data we have. The causality is respected if computing any derivative ∂ i+j u ∂x i ∂t j boils down to compute elements of initial condition u 0 (x). And in the case of the heat equation, we can ensure it will be possible thanks to Schwarz 's theorem about switching partial derivatives:

∂ i+j u ∂x i ∂t j = ∂ i+j-1 u ∂x i ∂t j-1 ∂u ∂t = ∂ i+j-1 u ∂x i ∂t j-1 ∂ 2 u ∂x 2 = ∂ i+j+1 u ∂x i+2 ∂t j-1 = ... = ∂ i+2j u ∂x i+2j
This graph illustrates the dependencies between the partial derivatives and we see that all arrows will end up on the vertical axis which represents the derivatives with respect to x only, namely the different parts of u 0 (x). The causality being respected ensures that the recursion will end. This example in 2 dimensions shows how the principle of causality is more flexible than it was presented with the Airy equation. Indeed, we said that coefficients of specific order should require strictly lower order coefficients, which is graphically represented by arrows crossing the blue line from the top right-hand corner down to bottom left-hand corner. But we state now that it is not a necessary condition as we can see with the heat equation where higher order coefficients are required but with respect to other variables. So arrows are allowed to cross the blue line in the opposite direction as long as they end on the vertical axis.

Figure 3 shows our heat equation solution developed at order 25. The vertical axis is the temperature. We set the initial conditions to a sinus, which concretely means we impose the temperature on one axis to be an alternation of warm and cold at initial time. The graph converges to a uniform average value along the time which is consistent with the physical interpretation.

What we call order here and denote by R is only the unrolling depth of the infinite tree we build. The graph in Fig. 4 shows the computation times (in seconds, on a common laptop computer) of the heat equation solution according to order and the graph in Fig. 5 shows this computation time divided by the number of coefficients of the solution, which lies in θ(R N ) with N the dimension, according to [START_REF] Thirioux | Verifying embedded systems[END_REF]. By dividing the computation time by the number of computed coefficients (normalized to 1 for R = 0), we aimed at evaluating the amount of additional computation done per useful coefficient, i.e. the "administrative" overhead induced by the resolution of the equation, due to auxiliary data structures, memory allocations, etc. We observe only a linear overhead and despite the relative simplicity of the heat equation, it comforts us in the decisions taken so far for implementing our framework.

Perspectives

Canonical Method for Composition

As defined in Sect. 5, composition involves the resolution of a partial differential equation. This hinders the computation of error bounds. Indeed, as far as we know, there is no established general method to solve such equations with certified errors, beyond ad-hoc situations such as elliptic, parabolic, etc., equations with specific initial conditions. In order to devise a direct more tractable and non recursive way to compose Taylor series, following schemes such as Faà di Bruno's formula, we first need to handle errors. As in formal power series, composition (f • g) may be achieved only when g has no constant part. To factorize out the constant part of g (so that we fall back to evaluation at point 0), we depend on an additive decomposition of f , when available.

Again, we sum up some decompositions of standard elementary functions. For every A n+1 ∈ R n+1 , we have the following equations, where we remark that their right-hand sides are built from a constant part (V A ) and another term without a constant part:

exp A L n + X n .A R n+1 + V A = exp(V A ) exp A L n + X n .A R n+1 log A L n + X n .A R n+1 + V A = log(V A ) + log 1 + A L n +X n .A R n+1 V A sin A L n + X n .A R n+1 + V A = sin(V A ) cos A L n + X n .A R n+1 + cos(V A ) sin A L n + X n .A R n+1 arctan A L n + X n .A R n+1 + V A = arctan(V A ) + arctan A L n +X n .A R n+1 1+V A .(A L n +X n .A R n+1 )
We are currently developing a canonical composition operator f • g following decomposition schemes that are all well known to strongly involve combinatorial reasoning. Our preliminary results already show that the administrative content of such heavy combinatorial computations, such as iterating over partitions, combinations, permutations and so on, have a great cost and are not yet on a par with the differential approach in terms of efficiency, at least for the tested instances. More investigation is required in that respect. We still expect to obtain an efficient canonical solution, with a simpler error propagation scheme and furthermore less computations to reduce such propagation.

Certified Errors

Taylor Models. Differential equations put aside, we are already able to compute certified errors in our framework. It merely requires the introduction of an arithmetical domain for errors. We introduce below a very simple error domain based upon symmetric zero-centered monotonic error functions.

Let us assume K stands for the value domain. Error functions are then elements of the following domain E, assuming we work in dimension N :

E {f ∈ (K + ) N → K + | f (0) = 0, f monotonous}
The error model is then the product K×E. The semantics of an element of this model represents a function from variable bounds to sets of possible values:

v, ǫ X ∈ (K + ) N → {k ∈ K | |k -v| ≤ ǫ(X)}
The error model has N + 1 constructors: (k, 0) for k ∈ K, denoted "k" and the i ∈ [0, N -1] indexed family (0, X → X i ), denoted "X i ". It is endowed with a K-algebra structure and is further turned into an full-fledged domain using suitable definitions of elementary functions on K × E, as illustrated below. Similar definitions may be devised for other elementary functions:

v 1 , ǫ 1 + v 2 , ǫ 2 v 1 + v 2 , ǫ 1 + ǫ 2 α × v, ǫ α × v, |α| × ǫ v 1 , ǫ 1 × v 2 , ǫ 2 v 1 × v 2 , |v 1 | × ǫ 2 + |v 2 | × ǫ 1 + ǫ 1 × ǫ 2 e v,ǫ e v , e v × (e ǫ -1) log v, ǫ log v, log 1 + ǫ v (v = 0)
Taylor models are then built from cotensors of value, error terms. We consider a function f ∈ R N → R, assumed analytical at point 0 and note respectively f α and ǫ α as the value and error at derivation multi-index α.

A Taylor model predicate T M(f, R, δ) at order R in a δ-neighbourhood of point 0 (where δ ∈ R + N ) is defined as the following:

T M(f, R, δ) ∀x ∈ R N .|x| ≤ δ =⇒ |f (x) - |α|≤R α=0 f α x α | ≤ |α|=R ǫ α (δ)|x| α
A Taylor model for parameters R and δ is then the set of functions f such that T M(f, R, δ) holds true.

Issues with Recursive Definitions. We recall that the above definitions must be amended in order to account for errors in (recursive) differential equations. Indeed, in that case, dependencies between errors at different derivation orders do not respect the causality relation fulfilled by pure values. So we need to compute another fixed point, different from the one for pure values. We illustrate this discrepancy between values and errors, considering the following partial development of a Taylor series with errors for a bivariate function f :

f (X, Y ) f 0 , ǫ 0 + X. f X , ǫ X + Y. f Y , ǫ Y + . . .
Then, integrating f along X, accounting for errors, yields the following series:

X f = 0, |X|.(|f 0 | + ǫ 0 ) + X. f 0 , ǫ 0 + Y. 0, |X|.(|f Y | + ǫ Y ) + . . .
Unfortunately, we remark that the error term |X|.(|f 0 | + ǫ 0 ) at order 0, while still a zero-centered monotonic error function, directly depends on ǫ 0 , the error function of f at order 0. The same problem occurs at order Y . On the contrary, the value part of the integrand is always 0, so is independent of f . As we wish to define f recursively through such an integrand, setting for instance f = X f , we face the necessity to find a different computation scheme for errors than for values. This is left for future work, but we feel that it would probably imply to transpose in our multivariate setting the kind of argumentation found in the Picard-Lindelôf theorem (that determines existence and unicity of solutions to ODEs in solved forms).

Going Further. Many other sensible choices for computing errors are also possible such as arbitrary intervals, zonotopes, etc., but we haven't experimented with these solutions yet. We chose to stick to the lightweight zero-centered error domain, giving up some precision to save computation time, mostly because it is much simpler to implement and also because we rely on on-demand cotensor exploration to increase precision, by computing deeper coefficients of Taylor expansions. We nevertheless plan to address the problem of finding a well-suited error domain, in terms of precision with respect to computation time.

Accounting for numerical errors is also on our roadmap. As a first approach, we postulate that we would only have to represent every real number with an interval of lower and upper approximations given as two floating-point numbers, lifting every computation from an algebra of real numbers to an algebra of floating-point intervals. The main question will be to test whether accumulating numerical errors along a huge number of computations could significantly degrade precision, as the derivation order increases, jeopardizing the core feature of our framework.

Another method, closely related to our own functional language framework exploiting laziness, would be to consider using a setup for exact real number algebra, as illustrated for instance in [START_REF] Geuvers | Constructive analysis, types and exact real numbers[END_REF]. Besides its lack of efficiency wrt. floatingpoint numbers, it would not suffer from a potential untamable accumulation of errors and would also open the way for a complete formal verification (including tensorial structure and numerical aspects). This is left for future work.

Conclusion

With a renovated view on Taylor series, we provide an implementation of a genuine full-fledged algebra of such series, in the multivariate case. Even if the work is far from being completed, it has been proven useful already as we are able to deal smoothly with partial differential equations in solved form, without any input from domain expert. To the best of our knowledge, implementing such an algebra of Taylor series with a concern on efficiency through carefully crafted algorithmics but also on correctness through strong typing has not been tried before. Indeed, although not presented here, our implementation puts an emphasis on strong typing, through extensive use of advanced OCaml GADT features. This proved really helpful in designing correct-by-construction code, at least with respect to dimensions and derivation orders, while implementing complex and error-prone numerical computations.

The next big challenges to take up are: first, the introduction of a better composition scheme; second, error domains and computation schemes compatible with every construction of our algebra. This would pave the way for applying our library in the paradigm of guaranteed integration for instance, notwithstanding other pervasive usages of Taylor series in various scientific fields.

  type ('a, _) st = | Nil: ('a, Nat.zero) st | Leaf: ('a, 'n Nat.succ) st | Node: ('a, 'n) st Lazy.t * 'a * ('a, 'n Nat.succ) st Lazy.t -> ('a, 'n Nat.succ) st and ('a, 'n) tree = ('a, 'n) st Lazy.t

Fig. 1 .

 1 Fig. 1. Our function (at order 150) Fig. 2. Theoretical result

Fig. 3 .Fig. 4 . 5 .

 345 Fig. 3. Heat equation solution

s 0,0,0,2 ... ... ... ...x 3 x 2 x 1 x 0 x 0 x 0 x 0 x 0 x 0 x 1 x 1 x 0 x 0 x 0 x 0 x 1 x 1 x 0 x 0 x 1 x 2 x 2 x 1 x 0 x 0 x 0 x 0 x 1 x 1 x 0 x 0 x 1 x 2 x 2 x 1 x 0 x 0 x 1 x 2 x 3 x 3 x 2 x 1 x 0 x 0 x 0 x 0 x 1 x 1 x 0 x 0 x 1 x 2 x 2 x 1 x 0 x 0 x 1 x 2 x 3 x 3 x 2 x 1x 0 x 0 x 1x 2x 3