Abstract-In this paper, we consider the problem of assigning a set of clients with demands to a set of servers with capacities and degree constraints. The goal is to find an allocation such that the number of clients assigned to a server is smaller than the server's degree and their overall demand is smaller than the server's capacity, while maximizing the overall throughput. This problem has several natural applications in the context of independent tasks scheduling or virtual machines allocation. We consider both the offline (when clients are known beforehand) and the online (when clients can join and leave the system at any time) versions of the problem. We first show that the degree constraint on the maximal number of clients that a server can handle is realistic in many contexts. Then, our main contribution is to prove that even if it makes the allocation problem more difficult (NP-Complete), a very small additive resource augmentation on the servers degree is enough to find in polynomial time a solution that achieves at least the optimal throughput. After a set of theoretical results on the complexity of the offline and online versions of the problem, we propose several other greedy heuristics to solve the online problem and we compare the performance (in terms of throughput) and the cost (in terms of disconnections and reconnections) of proposed algorithms through a set of extensive simulation results.

I. INTRODUCTION

In a client-server computing platform, where servers have capacity and degree constraints and clients have demands, we consider the problem of finding an allocation of clients to servers such that each server's degree and capacity constraints are satisfied while fulfilled demand is maximized. For instance, this models the problem of scheduling a very large number of identical tasks on a server-client platform [START_REF] Hong | Distributed adaptive task allocation in heterogeneous computing environments to maximize throughput[END_REF]. Initially, several servers hold or generate tasks that are transferred and processed by clients. The goal is to maximize the overall throughput achieved using this platform, i.e., the (fractional) number of tasks that can be processed within one time unit. Since QoS mechanisms for bandwidth control have to be used in order to cope with the heterogeneity of the clients [START_REF] Brown | Traffic Control HOWTO. Chapter 6. Classless Queuing Disciplines[END_REF], [START_REF] Hubert | Linux Advanced Routing & Traffic Control[END_REF], the degree constraint is related to the maximal number of TCP connections that a server can handle using QoS and the capacity of the server is defined as its overall outgoing bandwidth. This resource allocation problem also has applications in the context of Cloud Computing [START_REF] Zhang | Cloud computing: stateof-the-art and research challenges[END_REF], [START_REF] Van | SLA-aware virtual resource management for cloud infrastructures[END_REF], [START_REF] Calheiros | A heuristic for mapping virtual machines and links in emulation testbeds[END_REF]. In this case, servers represent physical machines and clients represent services, which can be deployed on the servers by using one or more virtual machines (VMs). Each service comes with its demand and a physical machine can host at most a given number of virtual machines (see Section II). In this context, the resource allocation problem can be used to find the allocation that allows the maximal fraction (the same for all services) that can be processed on a set of physical machines.

In the general setting, each server S j is characterized by its capacity b j (i.e., the quantity of data that it can send, or the number of flops that it can process during one time-unit, depending on the context) and its degree d j (i.e., the maximal number of open TCP connections, or the number of virtual machines that it can handle simultaneously). On the other hand, each client C i is characterized by its demand w i (i.e., the number of tasks that it can process during one timeunit, or its computational demand per time unit). Our goal is to build a bipartite graph between servers and clients, so that capacity, degree and demand constraints are satisfied.

Formally, let us denote by w j i the capacity allocated by server S j to client C i . Then, a valid allocation must satisfy the following conditions

∀j i w j i ≤ b j (1)
∀j Card{i : w j i > 0} ≤ d j (2) ∀i j w j i ≤ w i (3)
where Equation (1) refers to the capacity constraint at server S j , Equation (2) refers to the degree constraint at server S j and Equation (3) refers to the demand constraint at client C i .

Therefore, as introduced in [START_REF] Beaumont | Allocation of clients to multiple servers on large scale heterogeneous platforms[END_REF], Maximize-Throughput-Bounded-Degree (MTBD) problem is defined as follows:

Maximize j i
w j i under constraints (1), ( 2) and [START_REF] Armbrust | Above the clouds: A berkeley view of cloud computing[END_REF].

Due to the dynamic nature of the clients participating into a large scale volunteer computation or the virtual machines running on a Cloud, it is both interesting to study MTBD when the set of clients is known in advance or when clients join and leave the system at any moment, the offline [START_REF] Beaumont | Allocation of clients to multiple servers on large scale heterogeneous platforms[END_REF] and online scenarios [START_REF] Beaumont | On-line allocation of clients to multiple servers on large scale heterogeneous systems[END_REF] respectively. In the online context, it makes sense to compare the algorithms according to their cost, the number of changes in the allocation induced by a client arrival or departure, and their performance, the achieved throughput, as discussed in Section III.

The rest of the paper is organized as follows. In Section II, we present the applications of MTBD to the scheduling of independent tasks in the context of large scale volunteer computing platforms and to the allocation of services to physical machines in the context of Cloud computing. In Section III, we justify the used model, we formalize the allocation problem and we visit the results of this paper. In Section IV, we prove that MTBD is NP-Complete in the strong sense but that a small additive resource augmentation (of 1) on the servers degrees is enough to find in polynomial time a solution that achieves at least the optimal throughput. Then, we consider in Section V the more realistic setting where the set of clients is not known in advance but clients rather join and leave the system at any time, i.e., the online version of MTBD. We prove that no fully online algorithm (where only one change is allowed for each event) can achieve a constant approximation ratio, whatever the resource augmentation on servers degrees. Then, we prove that it is possible to maintain the optimal solution at the cost of at most 4 changes per server each time a new node joins or leaves the system. At last, we propose in Section VI several other greedy heuristics to solve the online problem and we compare performance in terms of throughput, and cost in terms of disconnections and reconnections, of proposed algorithms through a set of extensive simulation results based on realistic datasets. Concluding remarks are given in Section VII.

II. APPLICATIONS AND RELATED WORK

A. Independent Tasks Scheduling on Large Scale Platforms Scheduling computational tasks on a given set of processors is a key issue for high-performance computing, especially in the context of large scale computing platforms such as BOINC [START_REF] Anderson | BOINC: A System for Public-Resource Computing and Storage[END_REF] or Folding@home [START_REF] Larson | Folding@ Home and Genome@ Home: Using distributed computing to tackle previously intractable problems in computational biology[END_REF]. These platforms are characterized by their large scale, their heterogeneity and the performance variations of the participating resources. These characteristics strongly influence the set of applications that can be executed on these platforms. First, the running time of the application has to be large enough to benefit from the platform scale, and to minimize the influence of start-up times due to sophisticated middleware. Second, the applications should consist of many small independent tasks in order to minimize the influence of variations in resource performances and to limit the impact of resource failures. From a scheduling point of view, the set of applications that can be efficiently executed is therefore restricted, and we can concentrate on "embarrassingly parallel" applications consisting in many independent tasks.

Even in the context of independent tasks on heterogeneous resources [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-completeness[END_REF], makespan minimization, i.e., minimizing the time to process a given number of tasks, is intractable. An idea to circumvent the difficulty of makespan minimization is to lower the ambition of the scheduling objective. Instead of aiming at the absolute minimization of the execution time, it is generally more efficient to consider asymptotic optimality only (when the number of tasks is large). The goal is then to optimize the throughput. i.e., the fractional number of tasks that can be processed in one time-unit once steady-state has been reached. This approach has been pioneered by Bertsimas and Gamarnik [START_REF] Bertsimas | Asymptotically optimal algorithm for job shop scheduling and packet routing[END_REF] and has been extended to task scheduling [START_REF] Banino | Scheduling Strategies for Master-Slave Tasking on Heterogeneous Processor Platforms[END_REF] and collective communications [START_REF] Beaumont | Pipelining Broadcasts on Heterogeneous Platforms[END_REF]. Steady-state scheduling allows to relax the scheduling problem in many ways, and aims at characterizing the activity of each resource during each time-unit by deciding which (rational) fraction of time is spent sending and receiving tasks and to which client tasks are delegated, that is to focus on resource allocation rather than scheduling.

Independent task scheduling on large scale computing platforms can be modeled using MTBD problem. Following MTBD notation, each server S j is characterized by its capacity b j , the number of tasks it can send during one time-unit, and its maximal degree d j , the number of open connections that it can handle simultaneously. On the other hand, each volunteer is considered as a client C i and it's characterized by its demand w i , the number of tasks it can handle during one time-unit. In this case, client's capacity w i encompasses both its processing and communication capacities. More specifically, if comp i denotes the number of tasks C i can process during one time-unit, and comm i denotes the number of tasks it can receive during one time-unit, then we set w i = min(comp i , comm i ).

To model contentions, we rely on the bounded multiport model, that has already been advocated by Hong et al. [START_REF] Hong | Distributed adaptive task allocation in heterogeneous computing environments to maximize throughput[END_REF] for independent tasks distribution on heterogeneous platforms. In this model, a server S j can serve any number of clients simultaneously, each using a bandwidth w i ≤ w i provided that its outgoing bandwidth is not exceeded, i.e., i w i ≤ b j . This corresponds to modern network infrastructure, where each communication is associated to a TCP connection. This model strongly differs from the traditional oneport model used in scheduling literature, where connections are made in exclusive mode: the server can communicate with a single client at any time-step. Previous results obtained in steady-state scheduling of independent tasks [START_REF] Banino | Scheduling Strategies for Master-Slave Tasking on Heterogeneous Processor Platforms[END_REF] have been obtained under this model, which is easier to implement. For instance, Saif and Parashar [START_REF] Saif | Understanding the Behavior and Performance of Non-blocking Communications in MPI[END_REF] report experimental evidence that achieving the performances of bounded multi-port model may be difficult, since asynchronous sends become serialized as soon as message sizes exceed a few megabytes. Their results hold for two popular implementations of MPI, the message-passing standard: MPICH on Linux clusters and IBM MPI on the SP2. Nevertheless, in the context of large scale platforms, the networking heterogeneity ratio may be high, and it is unrealistic to assume that a 100MB/s server may be kept busy for 10 seconds while communicating a 1MB data file to a 100kB/s DSL node. Therefore, in our context, all connections must directly be handled at TCP level, without using high level communication libraries.

It is worth noting that at TCP level, several QoS mechanisms enable a prescribed sharing of the bandwidth [START_REF] Brown | Traffic Control HOWTO. Chapter 6. Classless Queuing Disciplines[END_REF], [START_REF] Hubert | Linux Advanced Routing & Traffic Control[END_REF]. In particular, it is possible to handle simultaneously several connections and to fix the bandwidth allocated to each connection. In our context, these mechanisms are particularly useful since w i encompasses both processing and communication capabilities of C i and therefore, the bandwidth allocated to the connection between S j and C i may be lower than both b j and w i . Nevertheless, handling a large number of connections at server S j with prescribed bandwidths consumes a lot of kernel resources, and it may therefore be difficult to reach b j by aggregating a large number of connections. In order to circumvent this problem, we introduce another parameter d j in the bounded multiport model, that represents the maximal number of connections that can be simultaneously opened at server S j .

Therefore, the model we propose encompasses the benefits of both the bounded multi-port model (by setting ∀i, d i = +∞) and the one-port model (by setting ∀i, d i = 1). It enables several communications to take place simultaneously, what is compulsory in the context of large scale distributed platforms, and practical implementation is achieved by using TCP QoS mechanisms and by bounding the maximal number of connections.

B. Virtualization in Cloud Computing Platforms

Cloud Computing [START_REF] Zhang | Cloud computing: stateof-the-art and research challenges[END_REF], [START_REF] Armbrust | Above the clouds: A berkeley view of cloud computing[END_REF] has recently emerged as a new paradigm for service providing over the Internet. Among the challenges associated to Cloud Computing is the efficient use of virtualization technologies such as Xen [30], KVM [START_REF]Kernel based virtual machine[END_REF] and VMware [29] and the migration of Virtual Machines (VMs) onto Physical Machines (PMs). Using virtualization, it is possible to run several Virtual Machines on top of a given Physical Machine. Since each VM hosts its complete software stack (Operating System, Middleware, Application), it is possible to migrate VMs from a PM to another. The ability to move virtual machines is crucial in order to achieve good load balancing [START_REF] Van | SLA-aware virtual resource management for cloud infrastructures[END_REF], [START_REF] Calheiros | A heuristic for mapping virtual machines and links in emulation testbeds[END_REF] in a dynamic context where VMs are added and removed from the system. It is also crucial for energy minimization [START_REF] Berl | Energy-efficient cloud computing[END_REF], [START_REF] Beloglazov | Energy efficient allocation of virtual machines in cloud data centers[END_REF] in order to determine if some PM can be switched off.

The mapping problem of services having heterogeneous computing demands onto PM having heterogeneous capacities can be modeled using MTBD. In this context, each physical machine S j is characterized by its computing capacity b j (i.e., the number of flops it can process during one time-unit) and its maximal degree d j (i.e., the number of different VMs that it can handle simultaneously, given that each VM comes with its complete software stack). On the other hand, each service C i is characterized by its demand w i (i.e., its overall processing demand during one time-unit). Then, a valid solution of MTBD provides a valid mapping of services onto PMs. The online version of MTBD corresponds to the case where services are added to or removed from the Cloud, or to the case when their demands change over time. In this case, the property that we prove in Section V stating that the online algorithm we propose bounds the number of changes on any PM is crucial as it enables to bound the number of VM migrations.

C. Related Works

A closely related problem is Bin Packing with Splittable Items and Cardinality Constraints, where the goal is to pack a given set of items in as few bins as possible. The items may be split, but each bin may contain at most k items or pieces of items. This is very close to the problem we consider, with two main differences: in our case the number of servers (corresponding to bins) is fixed in advance, and the goal is to maximize the total used capacity of the servers (corresponding to the total packed size), whereas the goal in Bin Packing is to minimize the number of bins used to pack all the items (corresponding to the number of used servers). Furthermore, we consider heterogeneous servers (what would correspond to bins with heterogeneous capacities and heterogeneous cardinality constraints).

Bin Packing with splittable items and cardinality constraints was introduced in the context of memory allocation in parallel processors by Chung et al. [START_REF] Chung | Parallelism versus Memory Allocation in Pipelined Router Forwarding Engines[END_REF], who considered the special case when k = 2. They showed that even in this simple case, this problem is NP-Complete, and they proposed a 3/2-approximation algorithm. Epstein and van Stee [START_REF] Epstein | Improved results for a memory allocation problem[END_REF] showed that Bin Packing with splittable items and cardinality constraints is NP-Hard for any fixed value of k, and that the simple NEXT-FIT algorithm achieves an approximation ratio of 2 -1/k. They also design a PTAS and a dual PTAS [START_REF] Epstein | Approximation Schemes for Packing Splittable Items with Cardinality Constraints[END_REF] for the general case where k is a constant.

Other related problems were introduced by Shachnai et al. [START_REF] Shachnai | Approximation schemes for packing with item fragmentation[END_REF]. They propose to model the size of an item as increasing when it is split and to ask for a global bound on the number of fragmentations. The authors prove that this problem does not admit a PTAS, and provide a dual PTAS and an asymptotic PTAS. In a multiprocessor scheduling context, another related problem is scheduling with allotment and parallelism constraints [START_REF] Shachnai | Multiprocessor Scheduling with Machine Allotment and Parallelism Constraints[END_REF]. The goal is to schedule a certain number of tasks, where each task comes with a bound on the number of machines that can process it simultaneously and a bound on the overall number of machines that can participate in its execution. This problem can also be seen as a splittable packing problem, but this time with a bound k i on the number of times an item can be split. In [START_REF] Shachnai | Multiprocessor Scheduling with Machine Allotment and Parallelism Constraints[END_REF], an approximation algorithm of ratio max

i (1 + 1/k i ) is presented.
In a related context, resource augmentation techniques have already been successfully applied to online scheduling problems [START_REF] Kalyanasundaram | Speed is as powerful as clairvoyance[END_REF], [START_REF] Phillips | Optimal Time-Critical Scheduling via Resource Augmentation[END_REF], [START_REF] Chekuri | Multi-processor scheduling to minimize flow time with resource augmentation[END_REF], [START_REF] Funk | On-line scheduling on uniform multiprocessors[END_REF] in order to prove optimality or good approximation ratio. More precisely, it has been established that several well-known online algorithms, that have poor performance from an absolute worst-case perspective, are optimal for these problems when allowed moderately more resources [START_REF] Phillips | Optimal Time-Critical Scheduling via Resource Augmentation[END_REF]. In this paper, we consider a slightly different context, since the off-line solution already requires resource augmentation on the servers degrees. We prove that it is possible in the on-line context to maintain at relatively low cost a solution that achieves the optimal throughput with the same resource augmentation as in the off-line context.

III. MODEL AND SUMMARY OF RESULTS

A. Platform Model and Maintenance Costs

Let us denote by b j the capacity of server S j and by d j the maximal number of clients that it can handle simultaneously (its degree). The capacity of client C i is denoted by w i . All capacities are normalized and expressed in terms of (fractional) number per timeunit. Moreover, let us denote by w j i the allocated value by server S j to client C i . Then, we have noticed in the introduction that MTBD can be expressed as a maximization problem under constraints (1), ( 2) and ( 3).

In the online version of MTBD, we introduce the notion of (virtual) rounds. A new round starts when a client joins or leaves the system, so that no duration is associated to a round. We denote by LC t the set of clients present at round t (with their respective capacities). Client C joins (resp. leaves) the system at round t if C ∈ LC t \LC t-1 (resp. C ∈ LC t-1 \LC t ). The arrival or departure of a client can therefore only take place at the beginning of a round and ∀t, |LC

t \LC t-1 | + |LC t-1 \LC t | ≤ 1.
Let us denote by LS the set of servers (with their respective capacity and degree constraints).

Solving the online version of MTBD comes into two flavors. First, one may want to maintain the optimal throughput at a minimal cost in terms of changes in existing connections between clients and servers. Second, one may want to achieve a minimal number of changes in existing connections at each server and to obtain the best possible throughput. In order to compare online solutions, we need to define precisely the cost of changing the existing allocation of clients to servers due to the arrival or departure of a new client.

Let us denote by w j i (t) the allocated value by server S j to client C i at round t. We say that client C i is connected to server S j at round t if w j i (t) > 0. We say that the connection between server S j and client C i changes at round t if w j i (t -1) = w j i (t), and we denote by N t j = |{i, w j i (t-1) = w j i (t)}| the number of changes occurring at server S j at round t.

This notion of change covers three different situations. If w j i (t -1) = 0 and w j i (t) > 0, then this change corresponds to a new connection to the server. Symmetrically, if w j i (t -1) > 0 and w j i (t) = 0, then client C i was disconnected from server S j . Finally, if both w j i (t -1) and w j i (t) are positive, this corresponds to a change in the allocation between client C i and server S j .

In the context of independent tasks scheduling, since we rely on complex QoS mechanisms to achieve the prescribed bandwidth sharing between clients and servers, any change in bandwidth allocation induces some cost. If a new client connects to a server, a new TCP connection needs to be opened, what also induces some cost. On the other hand, all modifications in bandwidth connections made by the different server nodes can take place in parallel. Similarly, in the context of Virtualization, adding or removing a VM from a PM induces some cost, due to migration. On the other hand, the different migration operations can be done in parallel.

Therefore, we introduce the following definition to measure and compare algorithms that solve online MTBD.

Definition 3.1: Let A be an algorithm solving the online version of MTBD. A induces l changes in connections per round if

max t max Sj ∈LS N t j = l.

B. Main Results

In the offline context, we first prove that MTBD is NP-Complete due to the degree constraint at the server nodes. On the other hand, we propose a sophisticated polynomial time algorithm, based on a slight resource augmentation to solve MTBD. More specifically, we prove that, if d j denotes the degree constraint at node S j , then the throughput achieved using this algorithm and degree d j + 1 is at least the same as the optimal one with degree d j (Theorem 4.4 in Section IV).

In the online context, the first result we present is that no online algorithm with cost less than 2 can achieve a constant approximation ratio, whatever the resource augmentation on the degree (Theorem 5.1 in Section V).

The second result presented in the online context shows that there exists a polynomial time online algorithm whose cost is at most 4 (see Theorem 5.5), with a resource augmentation of 1 (Lemma 5.2), and that maintains the optimal throughput at any round. Indeed, we know that Algorithm SEQ (Algorithm 1 in Section IV) provides at least the optimal throughput allowing the smallest possible additive resource augmentation α = 1. Hence, we transform algorithm SEQ into an online algorithm, and we use it to solve the online version of MTBD. The online version of Algorithm SEQ is called OSEQ (see Algorithm 2 in section V).

Therefore, in our context, maintaining the optimal throughput (with resource augmentation) is not more expensive, in terms of online cost and up to a constant ratio smaller than 2, than maintaining a constant approximation ratio of the optimal throughput.

IV. OFFLINE CASE ANALYSIS

We start the study of MTBD with the analysis of its complexity. Let us consider the corresponding decision problem, Throughput-Bounded-Degree-Dec (TBD-DEC), where the goal is to decide whether a throughput K can be achieved given a set of servers and a set of clients.

Lemma 4.1: TBD-DEC is NP-Complete in the strong sense.

Proof: To prove this result, we use a reduction to the 3-Partition problem [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-completeness[END_REF]. Indeed, let us consider an instance of 3-Partition consisting of 3m items a i such that a i = mB and ∀i, B 4 < a i < B 2 and let us set ∀j, d j = 3, b j = B, n = 3m, ∀i, w i = a i and K = mB. Since the overall out degree of the servers is at most 3m and since all 3m clients must be used in order to reach throughput mB, each server must be connected to exactly 3 clients and no client should be connected to more than one server. Since the overall capacity of the servers is m × B, each server must be connected to 3 clients whose aggregated capacity is exactly B, what achieves the NP-Completeness proof.

A. A Resource Augmentation Based Algorithm

Let us now present Algorithm SEQ, that relies on resource augmentation to provide a solution to MTBD problem. Due to the mentioned resource augmentation, SEQ outputs a non-valid solution in the sense that the number of clients allocated to a server S j may be d j + 1 instead of d j as stated in constraint (2). SEQ is described precisely in Algorithm 1.

In the following, we will consider lists of clients sorted by increasing capacities, and if LC = {C i } denotes such a list, we will denote by LC(l, k) = k i=l w i the sum of the capacities of the clients between C l and C k , both of them included.

Throughout a whole computation, Algorithm SEQ maintains an ordered list of remaining clients. At each step, it picks up a server S j arbitrarily and goes through the list to find a suitable set of clients for this server. A suitable set of clients is a set of d j +1 consecutive clients in the ordered list, called an interval of length d j +1, with total capacity at least b j , and such that the sum of the capacities of the first d j clients is less than the capacity b j of the server. These constraints ensure that the whole capacity and the maximum out-degree of the server are used. If such an interval [l, l + d j ] exists (there may be several, but any of them does the trick), Algorithm SEQ selects the rightmost one, i.e., the interval [l, l + d j ] such that LC(l, l + d j -1) < b j and LC(l

+ 1, l + d j ) ≥ b j .
This choice ensures that clients C l , C l+1 , . . . , C l+dj -1 are served completely by server S j by setting w j i = w i for all i ∈ {l, l + 1, . . . , l + d j -1}. If the total capacity of the interval exceeds b j , the last client can only be partially served. In that case, client C l+dj is served with capacity w j l+dj = b j -LC(l, l+d j -1) and then reinserted in the list of remaining clients with new capacity w l+dj equal to LC(l, l + d j ) -b j . In that case, client C l+dj will be connected to more than one server in the final solution. The list of clients is then updated and the algorithm goes on with the next server.

With respect to the ordering of the updated list of clients, let us point out that the choice of the rightmost interval ensures an ordering property. That is: the position of the modified client C l+dj in the sorted list remains the same. Indeed, C l+dj 's new capacity is equal to w l+dj = LC(l, l+d j )-b j = w l +LC(l+1, l+d j )-b j , and then the constraint LC(l + 1, l + d j ) ≥ b j ensures that w l+dj ≥ w l . Hence, w l-1 ≤ w l+dj ≤ w l+dj +1 , and thus the updated list of clients is already ordered. This property will be crucial in Section V. Indeed, among all possible valid intervals that can be allocated to S j , only the rightmost one produces an allocation that does not require many changes when a client joins or leaves the system (see Section VI).

It may happen that there exists no suitable interval for two reasons. The first one is that any set of d j +1 clients has not enough capacity to use all the bandwidth b j (i.e., the overall capacity of the d j + 1 largest clients is not big enough). In this case, SEQ allocates to server S j the d j largest clients (the last d j clients in the ordered list). Note that in that case, SEQ would be allowed to allocate one more client to server S j . But no valid solution could allocate more bandwidth to this server, and the extra connection may actually be useful later on.

On the other hand, if any set of d j clients has overall capacity larger than b j (i.e., the overall capacity of the d j smallest clients is already too large), then the algorithm simply allocates the k smallest clients, where k is the smallest index such that LC(1, k) ≥ b j . In this case also, the last client may be split, and its remaining capacity will be LC(1, k) -b j (clearly, the new client is the smallest one in the list and is reinserted at the same place in this case also).

B. Approximation Results

Let us now prove that the throughput allocated by Algorithm SEQ is at least as much as the throughput Split C l+dj in C l+dj and C l+dj with w l+dj = w l+dj +w l+dj and w l+dj = b j -LC(l, l+d j -1)

7: Set A j = {C l , C l+1 , . . . , C l+dj -1 , C l+dj } 8: Remove C l , C l+1 , . . . , C l+dj and insert C l+dj in LC 9: end if 10: if LC(1, d j ) ≥ b j then 11: Search for the smallest k such that LC(1, k) ≥ b j 12: Split C k in C k and C k with w k = w k + w k and w k = b j -LC(1, k -1) 13: Set A j = {C 1 , C 2 , . . . , C k-1 , C k } 14: Remove C 1 , C 2 , . . . , C k and insert C k in LC 15: end if 16: if LC(n -d j , n) < b j then 17: Set A j = {C n-dj +1 , C n-dj +2 , . . . , C n } 18: Remove C n-dj +1 , C n-dj +2 , . . . , C n from LC 19:
end if 20: end for 21: RETURN A = {A j } m j=1 provided by any valid solution. For the sake of simplicity, we consider that the length of the list of clients remains n during the execution of the algorithm. Without loss of generality, we assume that removed clients will thus be considered as 0-capacity clients and reinserted at the beginning of the list. To prove the result, we need to introduce an order between two lists of clients. Intuitively, if two lists of clients LC and R satisfy LC R, then, whatever the remaining servers, list LC will be easier to allocate than list R. Definition 4.2: Let LC and R be two lists of clients with the same length n and ordered by increasing capacities. We say that LC is easier than R (denoted by LC R), if

∀k ≤ n, LC(1, k) ≤ R(1, k)
Let us now consider a given step of the algorithm SEQ in which the considered server has capacity b and degree d. Let LC and R be two lists of clients. The application of this step of algorithm SEQ to the list LC leads to a remaining list LC . Similarly, a valid allocation 1 of this server to the list R yields a list of remaining clients R . The following lemma states that this atomic operation preserves the order .

Lemma 4.3: If LC R, and LC and R are obtained from LC and R as described above, then LC R .

Proof: We begin by proving two lower bounds for R (1, k). Since R is obtained from R by a valid allocation, there exists a set C ⊆ [1, n] of chosen clients, and assigned values v i for i ∈ C such that: Card(C) ≤ d, ∀i, v i ≤ R i (where R i denotes the capacity of the i th client in R), and v i ≤ b. There also exists a sorting permutation

σ such that R σ(i) = R i if i / ∈ C, and R σ(i) = R i -v i if i ∈ C. We can then write R (1, k) in two different ways, R (1, k) = i:i / ∈C∧σ(i)≤k R i + i:i∈C∧σ(i)≤k R i -v i = i:σ(i)≤k R i - i:i∈C∧σ(i)≤k v i For k > d, since there are at least k -d indexes i such that i / ∈ c ∧ σ(i) ≤ k, and since R(1, k -d) is the sum of the k -d smallest R i values, then i:i / ∈c∧σ(i)≤k R i ≥ R(1, k -d). Together with R i -v i ≥ 0, we obtain the first upper bound R (1, k) ≥ R(1, k -d) ∀k > d. (4)
Similarly, since there are k indexes i such that σ(i) ≤ k, then i:σ(i)≤k R i ≥ R(1, k). Together with i∈C v i ≤ b, we obtain the second upper bound

R (1, k) ≥ R(1, k) -b. (5)
To complete the proof, we need to evaluate LC (1, k). Since we identified three main situations when adding a server, we evaluate LC (1, k) in each possible situation.

Case 1 ∃l such that LC(l, l + d -1) < b and LC(l, l + d) ≥ b: In this case (see lines 4 to 8 in Algorithm 1) the algorithm allocates completely clients C l , C l+1 , . . . , C l+d-1 to S and only partially client C l+d , whose remaining capacity is w l+d . The first d clients of the list LC will thus have zero capacity, and C l+d will be reinserted at the same position as pointed out earlier. Then, the updated list LC is equal to

{C l = 0, . . . , C l+d-1 = 0, C 1 , . . . , C l-1 , C l+d = LC(l, l + d) -b, C l+d+1 , . . . , C n }.
1 Remember that the number of clients allocated to the server may be as high as d + 1 with SEQ, whereas it is limited to d in the valid solution.

Then, for k ≤ d, LC (1, k) is a sum over the completely allocated reinserted clients, and thus LC (1, k) = 0. For the second interval d < k ≤ l -1 + d, LC (1, k) is a sum of the first k -d capacities in LC, since they were shifted by d positions (due to the insertion of d clients at the beginning of the list), and so LC (1, k) = LC(1, k -d). If the interval includes one more client, i.e., d < k ≤ l +d, the sum is the same than in the previous interval, but the last element in the sum is replaced by the size of the split client that has been inserted, LC (1, k) = LC(1, k -d -1) + w l+d . Finally when l + d < k, the sum is equal to the sum in the original list, decreased by the total capacity allocated to

S, LC (1, k) = LC(1, k) -b.
Now, using Equations ( 4) and ( 5), and the fact that LC R, we have:

LC (1, k) = 0 ≤ R (1, k) for k ≤ d LC (1, k) = LC(1, k -d -1) + w l+d ≤ LC(1, k -d) ≤ R(1, k -d) ≤ R (1, k) for d < k ≤ l + d LC (1, k) = LC(1, k) -b ≤ R(1, k) -b ≤ R (1, k) for l + d < k. Case 2 A(1, d) ≥ b:
In this case (see lines 9 to 14 in Algorithm 1), since SEQ uses the first l ≤ d clients, there is no reordering of the list. The new list LC can therefore be written as {C 1 , . . . , C l-1 , C l , C l+1 , . . . , C n }, where C i has zero capacity for i < l. Moreover, since the overall allocated capacity is equal to b, then

LC (1, k) = 0 when k ≤ l -1 and LC (1, k) = LC(1, k) -b for k > l -1. Hence, Equation (5) combined to LC R leads to LC (1, k) ≤ R (1, k).
Case 3 A(n -d, n) < b: In this case (see lines 15 to 18 in Algorithm 1), SEQ allocates completely the d last clients to S, and therefore all reinserted clients C i will have zero capacity and will be reinserted at the beginning of the list. The new list LC can therefore be written as

{C n-d+1 , . . . , C n , C 1 , . . . , C n-d }. Therefore, LC (1, k) = 0 when k ≤ d and LC (1, k) = LC(1, k - (d + 1)) for k > d. Once again, Equation (4) combined with LC R leads to LC (1, k) ≤ R (1, k).
We can now state and prove the main result of this section.

Theorem 4.4: Let A be any valid solution of an instance I, and SEQ(I) be the solution given by algorithm SEQ. Then the throughput of SEQ(I) is at least as much as the throughput of A.

Proof: Using the ordering and Lemma 4.3, the proof of Theorem 4.4 becomes straightforward. Indeed, let us start with the initial list of clients LC 0 = LR 0 = L and let us denote by LC j (resp., LR j ) the list of remaining clients after the j-th first steps of Algorithm SEQ (resp., not fully allocated to servers S 1 , . . . , S j in the valid allocation A).

Then, a trivial induction, based of successive applications of Lemma 4.

C. Approximation algorithms

SEQ can easily be turned into a valid approximation algorithm with ratio ρ = dmin dmin+1 , where d min is the smallest degree of all servers. At the end of algorithm SEQ, we can disconnect one client from each server whose out-degree has been exceeded. Removing the smallest connected client cannot decrease the average quantity of resource allocated per connection. Thus, if we denote by w j the average quantity of resource allocated per connection of server S j at the end of SEQ, and by w j the average quantity of resource allocated per connection after the modification, we have w j dj ≥ w j dj +1 . Hence w j ≥ dj dj +1 w j ≥ ρw j . Since the overall throughput T is equal to the sum of all w j (and therefore is larger than the optimal throughput T * ), we obtain T ≥ ρT * .

This resource augmentation result can also be seen as an approximation result for the problem MDGT (Minimize Degree for a Given Throughput). Indeed, if we are given a bound T ≤ min( j b j , i w i ) on the throughput, a simple dichotomic search finds the minimum value α SEQ of α such that the throughput of SEQ(I(α)) is at least T on the modified instance I(α) in which server S j has degree d j + α. Theorem 4.4 states that if there exists a solution A of throughput T for instance I(α -1), then SEQ(I(α -1)) provides a valid solution for instance I(α) whose throughput is at least T .

Therefore, α SEQ ≤ α * + 1, where α * is the optimal (integer) value of the problem MDGT for instance I. Since MDGT is NP-complete, this is the best possible approximation result.

V. ONLINE CASE ANALYSIS

In this section, we consider more specifically the online case, where the set of clients is not known in advance, but clients can join and leave the system at any time.

Let us start the analysis of the online case by proving that no online algorithm whose cost is less than 2 (see Definition 3.1) can achieve a constant approximation ratio for the online MTBD problem. This result holds true even if we allow any constant resource augmentation ratio on the degree of the servers, what strongly differs from the offline setting, where a constant additive resource augmentation of 1 is enough to achieve optimal throughput. The proof is by counter-example.

An algorithm A α uses α ≥ 1 resource augmentation ratio when the maximal degree used by a server S j is d j + α, while its original degree is d j . Moreover, let us denote by OP T (I) the optimal throughput on instance I, and by A α (I) the throughput provided by Algorithm A α on instance I.

Theorem 5.1: Given a resource augmentation ratio α and a constant k, there exists an instance I of online MTBD, such that for any algorithm A α with cost less than 2,

A α (I) < 1 k OP T (I).
Proof: The proof is by exhibiting an instance I on which any online algorithm with cost less than 2 will fail to achieve the required approximation ratio.

This platform consists in only one server S with bandwidth b = (2k) α+1 and degree constraint d = 1. On the other hand, let us consider a set of clients C 0 , C 1 , . . . , C α+1 whose capacities are 1, 2k, (2k) 2 , . . . , (2k) α+1 . In the online instance I, clients arrive one after the other, by increasing capacities. More precisely, at round j, for 0 ≤ j ≤ α + 1, client C j with capacity (2k) j is added. Clearly, since the degree of the server is 1, only 1 client can be attached to the server and, since clients arrive by increasing capacity, the optimal solution consists in attaching C j to the server at round j. Note that maintaining this optimal solution at any time step has cost 2, since at each round, client C j is connected to the server and client C j-1 is disconnected.

In fact, any online algorithm that achieves an approximation ratio of at most k must attach C j to the server at round j. Indeed, the capacity of C j is larger than 3 2 k times the overall capacity of all previous clients, since

j-1 0 (2k) i < ( 3 2 k)(2k) j .
Therefore, any online algorithm whose approximation ratio is at most k needs to connect a new client at each round. Therefore, if its cost is strictly less than 2, it cannot disconnect clients, so that after round α + 1, the degree of the server would be α + 2, thus violating the maximal resource augmentation on the degree of the server node.

A. OSEQ Algorithm

Let us now present OSEQ Algorithm, the online version of Algorithm SEQ. OSEQ Algorithm retains the performance guarantee of SEQ by achieving the optimal throughput with only one extra connection per server. Moreover, OSEQ guarantees that each time a client joins or leaves the platform it produces at most 4 changes at each server, i.e., the cost of OSEQ Algorithm is 4.

OSEQ Algorithm can at first be seen as a pseudoonline algorithm in the sense that it produces the same solution as if SEQ was computed from the start at each round. In fact, even if it is easier to present and analyze OSEQ in this way, we will show in Section V-C how to re-use some of the computations to lower the complexity. A global view of the naive version of OSEQ is presented in Algorithm 2.

Algorithm 2 Algorithm OSEQ (naive version)

UPON a new round starts; SET LS the list of servers; SET LC = sort(LC) the ordered available clients at the current round;

APPLY algorithm SEQ to the instance (LS, LC); RETURN SEQ(LS, LC), the allocation at the current round;

Corollary 5.2: The throughput provided by algorithm OSEQ at every round is at least as much as the optimal throughput when the degree constraint is satisfied.

Above corollary follows directly from Theorem 4.4

B. Guarantee on the number of changes

We proceed now by proving that the solution provided by OSEQ Algorithm at every round (i.e., when a client joins or leaves the platform) produces at most 4 changes per server. To this end, we will keep track of the differences between the lists of remaining clients throughout the execution of OSEQ.

Definition 5.3: Let LC and R be two ordered lists of clients. We will say that R is an augmented version of C if it is obtained from LC by the insertion of a new client and possibly the increase of the capacity of the next client. Formally, LC is augmented to R if there exists an integer p ≤ n, a new client X and a value y ≥ 0 such that R = {C 1 , . . . , C p-1 , X , C p , C p+1 , . . . , C n }, where the capacity of X is smaller or equal to w p (the new capacity of client C p ) and w p = w p + y ≤ w p+1 .

The following lemma shows that a list of clients and any augmented version of it, when allocated to the same server, produces almost the same allocation. Let I be an instance consisting on one server with capacity b and degree d. Let LC be the updated list of clients after OSEQ is applied to instance I with list of clients LC, and let us denote with A the output of OSEQ. Similarly, let R denote the updated list of clients and B be the output of OSEQ when applied to instance I with a different list of clients R.

Lemma 5.4: If R is an augmented version of LC, then R is an augmented version of LC , and the allocations A and B differ by at most 4 changes.

Proof:

Let R consist of {C 1 , . . . , C p-1 , X , C p , C p+1 , . . . , C n },
where w p = w p +y and the capacity of X (w X ) is equal to x. The first step of the proof consists in computing the partial sums R(u, v) for any u ≤ v ≤ n. A quick case study shows that

R(u, v) =                LC(u -1, v -1) if p < u -1, LC(u -1, v -1) + y if p = u -1, LC(u, v -1) + x + y if u ≤ p < v, LC(u, v -1) + x if p = v, LC(u, v) if p > v.
In particular, since by hypothesis x ≤ w p and x + y ≤ w p+1 , then in all cases, R(u, v) ≤ LC(u, v). Furthermore, since x ≥ w p-1 , R(u, v) ≥ LC(u -1, v -1) also holds in all cases.

Let us now consider the application of OSEQ(d, b) to LC. Without loss of generality, we consider that a suitable interval [l, l + d] has been found, i.e., that LC(l, l+d-1) < b and LC(l+1, l+d) ≥ b. In that case, allocation A is (C l , . . . , C l+d-1 , C If the change from LC to R lies outside of the interval [l, l + d] (i.e., p < l or p > l + d), then this change does not affect to the execution of Algorithm OSEQ, and allocations A and B are the same. In that case, thus, the result holds.

Otherwise, the resulting allocation B depends on the value of R(l + 1, l + d). Indeed, our previous bounds for R(u, v) shows that R(l, l + d -1) < b, and R(l

+ 2, l + d + 1) ≥ b. Thus, either [l, l + d] or [l + 1, l + d + 1]
is the suitable interval for the application of OSEQ to R.

Let us suppose that R(l + 1, l + d) ≥ b: In this case, the resulting allocation B is (C l , . . . , X , C p , . . . , C If p = l + d, then X is the split client. This case actually results in only two changes in the allocations: the addition of (one part of) X , and the removal of C If p = l, then X is actually not included in B. Thus, this case results in two changes only: the modification of the capacity of C p , and the fact that C l+d is split differently. The updated list R is also an augmented version of LC , with X inserted and the capacity of C (b) l+d increased.

Theorem 5.5: The cost of Algorithm OSEQ' is at most 4.

Proof: Let us prove that if two lists of clients LC and R differ by the addition of a new client, then the resulting allocations to each server computed by OSEQ differ by at most 4 changes.

Denote by LC j the current list of clients after the first j rounds of OSEQ starting from LC 0 = LC, and similarly for R. It is clear that R is an augmented version of LC, Lemma 5.4 shows that if R j is an augmented version of LC j , then R j+1 is an augmented version of LC j+1 . Then, a trivial induction based on the application of Lemma 5.4 proves that resulting allocations differ by at most 4 changes.

In the case of the removal of a client, we can simply swap the role of LC and R in the previous statements.

C. Efficient Implementation Issues

We have first presented the naive version of OSEQ, as an algorithm that recomputes from scratch the whole solution at each round. However, the proof of Lemma 5.4 shows that it is possible to compute only the changes between the previous allocation and the new one. The implementation is made more complex by the analysis of many cases depending on the values of p, l, l + d and so on. However, it is also more efficient, since for each server, we only have to decide whether the suitable interval is [l, l + d j ] or [l + 1, l + d j + 1], instead of going through the whole list of clients. It is thus possible to do it in constant time, what leads to a global complexity of θ(m) for Algorithm OSEQ.

VI. EXPERIMENTAL EVALUATION

A. Heuristics for comparison

As already mentioned in Section II, related work has mostly been done in the context of Bin Packing, where there is an infinite amount of identical bins, and the goal is to pack all items in as few bins as possible 2 . Interestingly, in this setting, the NEXT-FIT algorithm has a worst-case approximation ratio of 2 -1/k [START_REF] Epstein | Improved results for a memory allocation problem[END_REF], but it can easily be observed that it does not exhibit a constant approximation ratio for the total packed size when the number of bins is fixed. Moreover, most of existing algorithms in this context are approximation schemes, with prohibitive running times. To provide a basis of comparison, we thus introduce online versions of the natural greedy heuristics that performed best in the offline setting.

• LCLS (Largest Client Largest Server) At each step, the client with the largest w i is associated with the server with the largest available capacity b j = b ji w j i . The client is split if necessary, in which case the remaining w i = w i -b j is inserted in the ordered list.

• LCBC (Largest Client Best Connection) In this heuristic, we also consider the largest client first, but servers are ordered according to their remaining capacity per connection, which is defined as the ratio between the remaining capacity b j and the remaining available degree d j . The server with the largest capacity per connection is selected. Here also, the client is split if necessary. We analyze here a worst-case instance on which LCBC achieves a throughput significantly lower than SEQ, for the same resource augmentation on the degree. For fixed even m and B, let us consider an instance with m servers, each of capacity B and degree d = B 2 , mB 2 small clients of size 1, with an additional big client of size mB 2 . On this instance, the solution of SEQ is to assign to each server d small clients and a part of size B 2 of the big client. This solution achieves a throughput of mB.

On the other hand, LCBC assigns first the big client to as few servers as possible ( m 2 of them), with parts of size B. Once this is done, m 2 servers remain unused, and each can only accommodate d+1 small clients because of the degree constraint. The total throughput of this solution is thus B m 2 + (d + 1) m 2 = mB 1 2 + 1 4 + 1 2B . Hence, when B grows, the ratio between the throughput of SEQ and LCBC tends to 3 4 .

B. Random Instance Generation

We generate instances randomly, trying to focus at the same time on realistic scenarios and difficult instances. Instances are more difficult to solve when the sum of server capacities is roughly equal to the sum of client capacities. Indeed, the minimum of both is a trivial upper bound on the total achievable throughput, and a large difference between them provides a lot of freedom on the largest component to reach this upper bound. Based on the same idea, we generate instances where the sum of the server degrees j d j is roughly equal to the number n of clients.

In order to get a realistic distribution of server and client capacities, we have used information available from the volunteer computing project GIMPS [START_REF]The great internet mersenne prime search (gimps)[END_REF] that provides the average computing power of all its participants. A simple statistical study shows that the computational power (based on the 7,000 largest participants) follows a power-law distribution with exponent α ≈ 2.09. We have thus used this distribution and this exponent to generate the capacities of both clients and servers. The resulting values are then scaled so that their sums ( i w i and j b j ) are roughly equal. Furthermore, the degree d j of server S j is chosen proportional to its capacity b j (it seems reasonable to assume that a server with a larger capacity can accommodate more clients), with a Gaussian multiplicative factor of mean 1 and variance 0.1. We generate instances with m servers and n = pm clients, where p is chosen as 10 or 50, and m varies between 10 and 160.

To generate online instances, we start from a complete instance. Two kinds of random events are then generated: departure of a client (picked uniformly at random), or arrival of a newly generated client. We generate 300 such events, each kind having probability 1/2. The time intervals between two successive events are generated as a Poisson process.

C. Results

We ran simulations for different instance sizes by varying the number m of servers. For each value of m, 250 instances were generated, and we plot on the figures the average, median, and the first and last decile over these 250 instances. For each algorithm, the line connects the average values, the upper error bar shows the last decile (which means that on 10% of the instances, the value was higher), the lower error bar shows the first decile (the value was lower on 10% of the instances), and the lonely mark in between is the median (half of the instances had lower values).

1) Offline simulations: In the first set of experiments, we have measured the throughput of the solutions proposed by each algorithm. All values are normalized against the previously mentioned upper bound min( j b j , i w i ). Figure 1 shows the average results on 250 instances when the number of servers varies from 10 to 160. We can already make some remarks:

• For these instances, algorithm SEQ performs consistently better than the others. In fact, it almost always reaches the upper bound. In a second set of experiments, we have computed for each algorithm A the minimum value α * that needs to be added to the degree of each server so that algorithm A reaches the upper bound B = min( j b j , i w i ). Note that the results of Section IV do not imply that α * ≤ 1 for algorithm SEQ, since it may well be the case that the upper bound cannot be reached with the original degree sequence. Average results for all algorithms and for varying m are depicted in Figure 2.

We can see that, as expected, algorithm SEQ makes very good use of the additional degree, and can almost always reach the upper bound with an increase of 1 or 2. As expected also, the ranking of algorithms observed for the total throughput is still the same when considering α * . We see that with LCBC, one needs about 5 more connections to reach the bound for p = 10, and between 10 and 15 when p = 50 (notice that since the sum of the server degrees j d j is roughly equal to the number n of clients, p represents the average degree of the servers).

A more precise look at the results for m = 160 is shown on Figure 3, where the value of α * for each instance is plotted against the dispersion of the clients capacities, measured by the relative mean difference of these values 3 . We can see that most of the values for LCBC are between 2 and 5 for p = 10, and between 6 and 12 for p = 50. However, it can be as high as 2p for instances with very large dispersion in client capacities, and these high values tend to increase the average. The results for algorithm LCLS exhibit the same kind of behavior, with larger values of α * for the most heterogeneous instances, and this explains larger average values. Therefore, for these difficult heterogeneous instances, we can see the benefit of the guarantee proved in Section IV for algorithm SEQ. Indeed, in these simulations the mean value of α is p, so that a value of α * of order more than 5 is expected to degrade significantly the networking performances of the servers. Thus, greedy algorithms fail to use the whole capacity of the platform in strongly heterogeneous cases, whereas 1 or 2 extra connections are enough using SEQ.

2) Online simulations: In the online simulations, we compare Online SEQ with OBC, an online version of LCBC. LCBC consistently outperforms LCLS in offline simulations, and this ranking still holds about their online version, and thus the online version of LCLS is not analyzed here. However, we also consider another version of SEQ, named SEQLEFT, which selects the leftmost suitable set of clients (instead of the rightmost one for SEQ). SEQ andSEQLEFT have very similar performance in the offline case, but their cost in online situations is quite different.

On Figure 4, we plot the total number of computed tasks throughout the instance, which is simply the integral over time of the instantaneous throughput (assuming for simplicity that changes from one solution to the other take no time). The value obtained is then normalized against the upper bound min( j b j , i w i ) (so that an average over 250 instances make sense). We can see that the offline results can be observed in this situation as well: the performance of OBC is about 5% worse than that of SEQ, which is always very close to the upper bound. Furthermore, higher values of p lower the variability of the results.

Figure 5 shows the cost of the algorithms. We can see that the cost of SEQ is always 4, while the cost of OBC is between 10 and 20 on average. However, it is once again very variable and reaches 35 on roughly 10% of the instances. Remember that the average outdegree of the servers is p (instances with m servers contain pm clients), so this result means that it is quite likely that, using OBC, at some point in the execution, one server has to change more than half of the clients it is connected to. This figure also shows the importance of the locality obtained by selecting the rightmost suitable interval in SEQ: the cost of SEQLEFT is not bounded by 4, and can get as high as the cost of OBC.

On the other hand, Figure 6 shows the average of the costs over the 300 events. We can see that the average cost of an event for SEQ is between 3 and 4, while it is around 1.5 (varying between 1.3 and 2) for OBC. This shows that events that incur many changes for OBC are relatively rare, and are compensated by many events that generate no or very few changes. However, we feel that this cost for maintaining the guarantees are justified by the higher performance and use of computing resources, and by the stability of SEQ.

VII. CONCLUSIONS

In this work, we have considered a resource allocation problem that models both independent tasks scheduling and virtual machines allocation problems. With respect to existing literature, our main contribution is to introduce a degree constraint, that is crucial for realism in both contexts. We prove that even if this additional constraint makes the resource allocation problem NP-Complete, only a very small resource augmentation on the degree is sufficient to achieve optimality. We also analyze the online setting, where the resources can change during the execution, as expected in mentioned applications. In the online context, we prove that maintaining optimality is not more expensive (up to a ratio of 2) than achieving a constant approximation ratio. Finally, we provide an extensive set of simulation results based on realistic data.

Algorithm 1 5 :

 15 Algorithm SEQ 1: Set S = {S j } m j=1 and LC = sort({C i } n i=1 ); 2: Set A = {A j = {∅}} m j=1 and j = 1; 3: for j = 1 to m do 4: if ∃l such that LC(l, l +d j -1) < b j and LC(l, l + d j ) ≥ b j then Pick l s.t. LC(l, l + d -1) < b and LC(l + 1, l + d) ≥ b 6:

  3 proves that LC m LR m . This means that ∀k ≤ n, LC m (1, k) ≤ LR m (1, k), and in particular LC m (1, n) ≤ LR m (1, n), where LC m (1, n) and LR m (1, n) respectively denote the overall unused capacity of the clients in the solution computed respectively by SEQ and A. Hence, the throughput obtained using Algorithm SEQ is larger than the throughput obtained in solution A, what achieves the proof of Theorem 4.4.

  l+d ), and the updated list LC is (C 1 , . . . , C l-1 , C (b) l+d , . . . , C n ), where C (a) l+d and C (b) l+d are the two parts of the split client C l+d .

  It differs from A by 4 changes: the addition of X , the removal of C (a) l+d , and the modification of both C p and C l+d-1 . The updated list of clients R is (C 1 , . . . , C l-1 , C (b) l+d-1 , C l+d , . . . , C n ). It is thus an augmented version of LC , since C (b) l+d-1 is inserted between C l-1 and the capacity of C (b) l+d is increased to w l+d . Indeed, by definition of the splitting process, w(C (b) l+d-1 ) = R(l, l + d) -b and w(C (b) l+d ) = LC(l, l + d) -b, which implies w(C (b) l+d-1 ) ≤ w(C (b) l+d ).

  updated list R is also an augmented version of LC , with the remaining part of X inserted and the capacity of C (b) l+d increased to w l+d + y. Let us now suppose that R(l + 1, l + d) < b: In this case, the suitable interval is [l + 1, l + d + 1] and thus the resulting allocation B is (C l+1 , . . . , X , C p , . . . , C (a ) l+d ). Once again, it differs from A by 4 changes: the addition of X , the removal of C l , and the modification of both C p and C (a) l+d . The updated list of clients R is (C 1 , . . . , C l , C (b ) l+d , . . . , C n ). It is therefore an augmented version of LC , since C l is inserted in right after C l-1 and the capacity of C (b) l+d is increased to w(C (b ) l+d ). In this case, w l ≤ w(C (b) l+d ) comes from the ordering property of SEQ, see Section V-A.

Figure 1 .Figure 2 .

 12 Figure 1. Offline simulations: Average normalized throughput for p = 10 and p = 50.

Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 3456 Figure 3. Offline simulations: α * against client dispersion for m = 160, p = 10 and p = 50.

•

  We also define an online version of this heuristic: leaves, OBC tries to use the newly available bandwidth to reduce the indegree of other clients. Assume that X was connected to server S, and that client Y is connected to both S and S . When X leaves, S can reallocate the corresponding bandwidth to client Y, what can be of interest if this allows to disconnect Y from S , since this lowers the outdegree of S . OBC selects as many such incident connections as possible, starting from the smallest ones. If there are some unconnected clients remaining, OBC then acts as if they had just arrived and tries to connect them with the procedure described earlier.

	newly arrived client with a smaller one, selecting
	the server which yields the largest gain in total
	throughput.
	When a client X
	Online Best Connection (OBC). Servers are still
	ordered by their remaining capacity per connection.
	When a new client arrives, it is connected to the
	server whose capacity per connection is closest to
	the client's capacity. If no server is available, OBC
	goes through all servers which have some band-
	width remaining but no degree left, and swaps the

The work of Olivier Beaumont, Lionel Eyraud-Dubois and Hejer Rejeb was partially supported by French ANR project Alpage. The work of Christopher Thraves was supported in part by Spanish MICINN grant Juan de la Cierva.

In our context, servers are bins and clients are items

The mean difference of values {y i } is the average absolute difference of all couples of values. The relative mean difference is the mean difference divided by the arithmetic mean.


Jean-David Génevaux

Julien Narboux

Pascal Schreck

Formalization of Wu's simple method in Coq

Keywords: formalization, automation, geometry, Wu's method, proof assistant, Coq

We present in this paper the integration within the Coq proof assistant, of a method for automatic theorem proving in geometry. We use an approach based on the validation of a certificate. The certificate is generated by an implementation in Ocaml of a simple version of Wu's method.

Introduction

Bringing new automation techniques to interactive theorem provers while preserving safety is an important goal in order to spread the use of proof assistants.

In this paper, we focus on one of the most successful methods for automated theorem proving in geometry, namely Wu's method. We integrate this method into the Coq proof assistant [START_REF]The. The Coq proof assistant reference manual[END_REF]. Indeed, Wu's method can prove hundreds of geometry theorems [START_REF] Wu | On the Decision Problem and the Mechanization of Theorem Proving in Elementary Geometry[END_REF] and some conjectures were first proved using this method [START_REF] Chou | Mechanical Geometry Theorem Proving[END_REF][START_REF] Wang | A new theorem discovered by computer prover[END_REF].

One could formally prove correct the implementation of this decision procedure within a proof assistant such as Coq or Isabelle, but this would require a tedious work. Another approach is to modify the decision procedure to make it produce a witness of its correctness (the certificate) which can be checked by an external tool called a validator. The certificate is such that the validator is simpler to prove than the original decision procedure. Using such an approach, the completeness of the method is not guaranteed as the certifying algorithm may generate an erroneous certificate or even fail to generate one. But if we prove formally the correctness of the validator, when the validator confirms a result we get the same level of confidence as if we had proven the correctness of the decision procedure. This approach has several advantages: first the algorithm which generates the certificate is independent of the validator and hence can be written by different people using different languages, second the decision procedure does not need to be proved formally, hence the implementation can be optimized more easily. The difficulty is to generate a certificate which can be checked in a reasonable amount of time.

Gröbner basis is another well known tool for automated deduction in geometry [START_REF] Kapur | Geometry Theorem Proving using Hilbert's Nullstellensatz[END_REF]. Gröbner basis have already been integrated into Coq by Benjamin Grégoire, Loïc Pottier and Laurent Théry [START_REF] Pottier | Connecting Gröbner Bases Programs with Coq to do Proofs in Algebra, Geometry and Arithmetics[END_REF][START_REF] Grégoire | Proof Certificates for Algebra and their Application to Automatic Geometry Theorem Proving[END_REF]. Formalizing Wu's method brings us a tool to decide ideal membership as the Gröbner basis method does. But, in the context of automatic deduction in geometry, Wu's method and its variants developed by Shang-Ching Chou, Xiao-Shan Gao and Dongming Wang have been shown to be generally more efficient [START_REF] Shang | Ritt-Wu's Decomposition Algorithm and Geometry Theorem Proving[END_REF][START_REF] Wang | Elimination Method[END_REF][START_REF] Wang | Elimination Practice[END_REF]. Moreover, in geometry, the statements as given by a user are usually false in some special cases, such as when three points are collinear for example. These conditions under which a statement is true are called non degeneracy conditions (ndgs). An advantage of Wu's method is that it can generate the ndgs. This is crucial for the applications we aim at in the context of education. In the future, we would like to integrate this method in a tool to prove geometry theorems within a dynamic geometry software following the work presented in [START_REF] Tuan | An Additional Tool About the Orientation for Theorem Proving in the Coq Proof Assitant[END_REF][START_REF] Tuan | A Coq-based Library for Interactive and Automated Theorem Proving in Plane Geometry[END_REF].

Related Work

Wu's method has been implemented by Chou, he used his implementation to prove hundreds of geometric problems [START_REF] Chou | Mechanical Geometry Theorem Proving[END_REF]. Xiao-Shan Gao created Geometry Expert (GEX) [START_REF] Gao | MMP/Geometer -a software package for automated geometry reasoning[END_REF][START_REF] Gao | Geometry expert, software package[END_REF] and Zheng Ye produced Java Geometry Expert (Java GEX) [START_REF] Ye | An Introduction to Java Geometry Expert[END_REF]. Dongming Wang developed further elimination algorithms based on pseudo-division and implemented them as a Maple library called Geother [START_REF] Wang | GEOTHER 1.1: Handling and Proving Geometric Theorems Automatically[END_REF]. Judit Robu implemented Wu's method in Theorema [START_REF] Robu | Geometry Theorem Proving in the Frame of the Theorema Project[END_REF]. Predrag Janičić implemented the method within GCLC [START_REF] Janičić | System Description: GCLCprover + GeoThms[END_REF].

We are not aware of any formalization of Wu's method inside a proof assistant. But other methods for automatic deduction in geometry have been integrated into proof assistants: Gröbner basis can be used in Coq, HOL-Light and Isabelle [GPT11,Har07,CW07], the area method for non oriented euclidean geometry has been formalized in Coq as a tactic [CGZ94,Nar04,JNQ10] and geometric algebras are also available in Coq and can be used to prove automatically theorems in projective geometry [START_REF] Fuchs | A Formalisation of Grassmann-Cayley Algebra in Coq[END_REF].

The closest work to ours is the implementation of a Gröbner basis method within Coq by Benjamin Grégoire, Loïc Pottier and Laurent Théry [START_REF] Grégoire | Proof Certificates for Algebra and their Application to Automatic Geometry Theorem Proving[END_REF]. In this paper, we present an extension of their work to deal with Wu's method.

This paper is organized as follows. We first give a quick overview of Wu's method and we highlight the facts we need to prove to show the correctness of the method. Then we describe our formalization within the Coq proof assistant.

Overview of Wu's Method

Cartesian geometry

We are interested in formally proving theorems in geometry, but there are many geometries! Wu's method focuses on geometry with coordinates in a field. To be more precise, we consider here the plane F 2 where F is a field where all predicates on points, lines, circles, etc. can be expressed under the form P (X) = 0, X being a vector of variables corresponding to the coordinates of the involved objects.

Let us give a simple (and very classical) example:

Example 1.

In a parallelogram the diagonals intersect in their midpoints:

A B C D I Fig. 1. Classical parallelogram example
In the usual Euclidean geometry, the fact that lines (AB) and (CD) are parallel can be expressed by the polynomial equation1

(x B -x A )(y C -y D ) -(y B -y A )(x C -x D ) = 0.
The fact that point I is the midpoint of [AC] and the midpoint of [BD] can be expressed by the polynomial equations (

x B + x D ) -(x A + x C ) = 0 and (y B + y D ) -(y A + y C ) = 0.
The theorem stating that in a parallelogram the diagonals intersect in their midpoints is then directly stated by

(x B -x A )(y C -y D ) -(y B -y A )(x C -x D ) = 0 (x C -x B )(y D -y A ) -(y C -y B )(x D -x A ) = 0 ⇒ (x B + x D ) -(x A + x C ) = 0 ∧ (y B + y D ) -(y A + y C ) = 0
where all the variables are implicitly universally quantified.

The goal of Wu's method is precisely to prove geometric theorems which can be put under this algebraic form.

Rings and ideals

The general question is the following:

(1) Given k polynomials H 1 , . . . , H k and G in F(x 1 , . . . x m ), is it true that ∀x 1 , . . . x m ∈ F :

k i=1 H i (x 1 , . . . x m ) = 0 ⇒ G(x 1 , . . . , x m ) = 0 ?
This question gets obviously a positive answer when G belongs to the radical of the ideal I =< H 1 , . . . , H k > generated by polynomials H 1 , . . . , H k , that is, when there exists an integer r and

k polynomials Q 1 , . . . , Q k such that G r = Σ i Q i H i .
The famous Hilbert's Nullstellensatz theorem states that if F is algebraically closed, then the converse is also true. That is, we can always find such polynomials.

Then, in this framework, proving a geometric theorem consists in proving that the polynomial belongs to an ideal for which a set of generator is known. Note that the Nullstellensatz has to be used only to prove that a polynomial does not belong to an ideal: if we want to prove a known theorem, producing a linear combination of polynomials H i equal to G is enough.

It would be easy to test ideal membership by using Euclidean division. But, unfortunately, the ring of multivariate polynomials over a field is not Euclidean when more than one variable is involved. Nonetheless, the two most famous methods use a kind of Euclidean division to perform such a test. Buchberger's method consists in transforming the generating set into a so-called Gröbner basis in which a division algorithm can be used to give the good result while in the Wu's method a pseudo-division is used which closely mimics the Euclidean division.

Pseudo-division and pseudo-remainder

The idea of pseudo-division consists in multiplying the polynomial to be divided, say P , by the leading coefficient of the dividing polynomial, say Q, to a certain power in order that all the coefficients of this product are divisible by the leading coefficient of Q. More formally, let P and Q be polynomials of F[x 1 , . . . , x m ], and I ∈ F[x 1 , . . . , x m-1 ] be the leading coefficient of Q in the variable x m , the pseudo-division of P by Q in the variable x m yields two polynomials T and R such that

I r P = T Q + R
where r is the number of non zero coefficients of P and

R ∈ F[x 1 , . . . , x m ] with deg(x k , R) < deg(x k , Q).
Establishing the existence of the polynomials T and R and the correction of the pseudo-division algorithm is very similar to the proof of the correction of the Euclidean division algorithm. R is called the pseudoremainder of P by Q into the variable x k and it obviously belongs to the ideal < P, Q >.

Pseudo-division and pseudo-remainder can be used to perform a triangulation of the system i H i (x 1 , . . . , x m ) = 0. The triangulation produces an algebraic system i T i (x 1 , . . . , x m ) = 0 where each

T i belongs to F[x 1 , . . . x i ] ∩ < H 1 , . . . , H k >. This way, the implication i H i (x 1 , . . . , x m ) = 0 ⇒ i T i (x 1 , . . . , x m ) = 0
holds, but the converse is true if and only if the leading coefficient considered in the used pseudo-divisions do not vanish. But for the formalization, we do not the need the converse.

Note that, if the statement is given as a ruler-and-compass construction then the triangulation is tractable. Moreover, after choosing a reference, it may even become trivial. For instance, by fixing point A at (0, 0) and point B on Ox axis and observing that y C = y D , the statement of example 1 becomes :

(x C -x B -x D ).y C = 0 ⇒ x B + x D -x C = 0
where x B , x C and y C are parameters of the figure and x D is a dependent variable. Note that, by making this choice, we implicitly assume that A = B. The hypotheses are trivially under a triangular form. Next Section explains in details how to fix a reference, why this is correct, and how to construct the figure corresponding to the hypotheses.

Eventually, using successive pseudo-divisions of G by T k , T k-1 , . . . , T 1 , the method tests if IG belongs to the ideal generated by polynomials H i where I is a product of the leading coefficients of the triangulated polynomials. Summarizing the calculi, Wu's method, in his simplest form, allows to compute some polynomials I, S i and R such that:

IG = i (H i S i ) + R
It is clear that if the polynomial R is null, then the theorem is proved under the assumption that I(x 1 , . . . , x m ) = 0. When the theorem is stated as a construction this assumption corresponds to the non degeneracy conditions [START_REF] Shang | A class of geometry statements of constructive type and geometry theorem proving[END_REF].

In our example, the degenerated case occurs when y C = 0. This matches the case where the four points A, B, C and D are collinear : in this case, the theorem does not hold. The converse is not true: R = 0 does not mean that the theorem is false. The simple method of Wu is not complete but according to Chou [START_REF] Chou | Mechanical Geometry Theorem Proving[END_REF] it is powerful enough to prove hundred of classical theorems. A complete method would need to use ascending chains which are considered in the Ritt-Wu principle [START_REF] Chou | Mechanical Geometry Theorem Proving[END_REF].

The main steps of the Wu's method for geometry theorem proving can be summarized as follows:

1. Transforming the statements into an algebraic form. 2. Choosing the origin and a direction for the system of coordinates. 3. Showing that to prove the statement in general it is sufficient to show that it holds in the given system of coordinates. 4. Simplifying the polynomials thanks to the choice of the reference. 5. Triangulating the list of hypotheses using pseudo-division.

6. Pseudo-dividing successively the goal by the triangulated hypotheses. If the final remainder is the null polynomial, then the statement holds under the condition that some polynomials do not vanish (we call those polynomials the non degeneracy conditions) 7. Re-interpreting non degeneracy conditions as geometric predicates.

Formalization in Coq

In this section we describe our formalization of Wu's method in Coq. To formalize a decision procedure within Coq, there are several solutions: 1. We could write a tactic in the implementation language of Coq (Ocaml) which generates a proof term. 2. We could write a tactic in Ltac, the tactic language of Coq: a domain specific language which allows pattern matching on the context of the proof and backtracking. 3. We could implement the decision procedure within Coq itself, prove its correctness and use it within Coq using reflection. 4. As explained in the introduction we can implement the decision procedure as an external tool which generates a certificate and then check the certificate by reflection using a procedure written in Coq itself.

The first and third solutions can be seen as special (extreme) cases of the fourth one. For the first solution, the certificate can be considered as the proof itself and the validator is trivial, and for the third solution the certificate contains no information and the validator is the whole decision procedure.

As explained in the introduction for the core of the method we choose an approach based on the validation of a certificate, and we have been able to reuse the machinery developed by Benjamin Grégoire, Loïc Pottier and Laurent Théry for the Gröbner basis. But for some steps of the method we use the tactic language of Coq (Ltac). We use Ltac for the first steps, to put the problem in algebraic form and to choose a system of coordinates. Figure 2 gives an overview of the different parts of the development and their implementation language. The simplification phase consists only perfoming some small simplifications using substitution in Ltac before starting the actual procedure when equations are trivial. The geometrization phase is not yet implemented. The simple procedure of Wu that we use may produce ndgs which are not necessary. We plan to implement geometrization only for the predicates corresponding to ndgs conditions of ruler and compass constructions (collinear, parallel, point equality).

In the next sections, we describe first the algebraization process and then the generation of the certificate. We do not describe fully the validator here because we could reuse the validator of Benjamin Grégoire, Loïc Pottier and Laurent Théry.

Algebraization

Putting the statement in algebraic form, choosing the right system of coordinates, and simplifying the generated coordinates, is the first step of the method, and it is crucial. In this section we describe the formalization of this step using the tactic language of Coq.

Stating the conjecture As most automatic theorem provers in geometry, to state a theorem, we will assume that the user provides the assumptions as a list of geometric predicates taking only points as parameters. This has the advantage to simplify the formalization and as shown in [START_REF] Narboux | A graphical user interface for formal proofs in geometry[END_REF], it is possible to transform a statement containing points and lines into a statement containing only points. Figure 3 gives the Coq definitions of some common geometric predicates. X A and Y A denotes respectively the x and y coordinates of point A. Figure 3 lists some of the Coq definitions for the common geometric predicates. These definitions have the advantage that degenerated cases are not excluded. For instance, parallel A B C D holds when A = B. This leads to statements which are as general as possible.

Definition collinear A B C := (X A -X B) * (Y B -Y C) = (Y A -Y B) * (X B -X C). Definition parallel A B C D := (X A -X B) * (Y C -Y D) = (Y A -Y B) * (X C

Fig. 3. Definition of some geometric predicates

Then the statement corresponding to our example (without the ndg) is the following:

Lemma parallelogram : forall A B C D E F:Point, B <> A -> parallel A B C D -> parallel A D B C -> is_midpoint E A C -> is_midpoint F B D -> equal E F.
To put the statement in algebraic form we just need to unfold the definitions of the geometric predicates. But in practice, even for simple examples, if we do not fix a specific coordinate system the computations will take much more time or may even fail due to lack of memory. The classical solution consists in adding additional assumptions to fix the coordinate system, assuming that one point is the origin and another belongs to the x axis:

Lemma parallelogram : forall A B C D E F:Point, X A = 0 -> Y A = 0 -> Y B = 0 -> B <> A -> parallel A B C D -> parallel A D B C -> is_midpoint E A C -> is_midpoint F B -> equal E F.
This solution is sufficient to perform benchmarks in a context where the user is an expert. As we aim to apply this method in an educational context, we want to provide a procedure which works on the original statement.

We need to show that, without loss of generality, we can assume that A is of coordinates (0, 0) and B of coordinates (x b , 0). Following the idea of John Harrison [START_REF] Harrison | Without loss of generality[END_REF], this requires to show that the predicates we use are invariant under translation and rotation. For that purpose, we define a translation function taking as arguments a point and a vector (trans) and a rotation function (rot) taking as arguments a point and the sine and cosine of the angle of rotation. Then for each predicate (collinear in this example), we prove lemmas of the following form:

Lemma collinear_inv_translation: forall A B C V, collinear A B C <-> collinear (trans A V) (trans B V) (trans C V). Lemma collinear_inv_rotation: forall A B C cos sin, cos*cos + sin*sin = 1 -> (collinear A B C <-> collinear (rot A cos sin) (rot B cos sin) (rot C cos sin)).

Finally, we have a tactic Algebraization O I H which takes as input a point O to be put at the origin, a point I to be put on the x-axis and a proof that O = I, and which makes use of the above lemmas to perform the required simplifications. The algebraization tactic only work for goals which are stated using function and predicates which take as arguments only points. The following predicates/functions are available: collinear, parallel, orthogonal, midpoint, intersection of lines, square of length, equality of points, angles or lengths. The tactic can not deal with user defined predicates automatically. Adding a new predicate requires to add the lemmas for invariance under translation and rotation and to update the tactic.

We leave for future work, the automatic choice of a reference. Our experiments shows that this choice is crucial : choosing different references can lead to very different computation times. The heuristic proposed by Chou is to choose as an axis a line which contains many points and also a perpendicular line. We noticed that choosing as the origin one of the points of the goal is also often a good choice, because it simplifies the starting polynomial of the successive pseudo-division.

After application of the tactic Algebraization A B H, we get the following goal:

H : 0 = -X P2 * P3 + X P0 * P3 H0 : -P3 * X P0 + Y P2 * X P0 = X P2 * Y P0 H4 : 2 * X P = X P0 H5 : 2 * Y P = P3 + Y P0 H2 : 2 * X P1 = X P2 H6 : 2 * Y P1 = Y P2 ______________________________________(1/1) X P = X P0

Certified implementation

In the section we describe our implementation of a simple version of Wu's method and how it generates certificates.

Implementation of the method To integrate Wu's method in Coq, we need to modify the implementation to generate a certificate. Unfortunately, the implementations of Wu's method which are already available are either not open source or rely on proprietary Computer Algebra Systems. Moreover, we hope that in the future our tactic could be distributed with Coq, hence we decided to write our own version of the method in Ocaml, the implementation language of Coq. We base our implementation of Wu's method on the Ocaml libraries for dealing with multivariate polynomials developed by Loïc Pottier for the integration of Gröbner basis in Coq. He used two different data structures for polynomials : first a list of monomials, each monomial is represented by a coefficient and an array containing the degree of each variable and the total degree of the monomial, second a recursive data structure considering polynomials with several variables as a polynomial with one variable but whose coefficients are polynomials in the other variables. We reused his data-structures, and slightly optimized the one based on list of monomials by changing the representation of monomials: in the arrays of degrees we store only the degrees of the variables until the highest variable which appears in the monomial. This variant reduces both the memory footprint and the number of needed comparisons of degree. This is crucial, because 30% of the computing time is spent in the monomial comparison function which is the core of the multiplication function. We implement the simple algorithm presented by Chou in [START_REF] Chou | Mechanical Geometry Theorem Proving[END_REF], hence it is incomplete because we do not check polynomials for irreducibility. But in practice many theorems can be proved using this simple method. Our certifying implementation consists in 3000 lines of Ocaml half of which are tests and examples.

Certificate generation Certificates are pieces of information, which allow to check if the result of a computation is correct. In our case, the certificate we generate is based on the Nullstellensatz, hence we need to produce a proof that the polynomial which represents the goal (multiplied by a coefficient representing the ndgs) is in the ideal generated by the hypotheses. We need to show that there are some integer r, polynomials I and S i such that:

I r G = i (S i H i )
As explained in Sec. 2, the core of Wu's method can be decomposed into two main steps: the triangulation and the successive pseudo-divisions.

To generate the certificate we generate three kinds of intermediary certificates: one for the pseudo-divisions, one for the triangulation, and another one for the successive pseudo-division procedure.

Certificates for pseudo-division

The certificates for the pseudo-division will be used both by the triangulation and successive pseudo-division procedure.

Recall that the pseudo-division of A and B return a remainder R2 such that there exists an integer d and a polynomial Q such that:

I d A = QB + R
where I is the leading coefficient of B.

So to certify that R is in the ideal generated by A and B, instead of just returning the remainder R we also return I, d and Q which verify that R =

I d × A + (-Q) × B.
To ease implementation we number our polynomials. Then the implementation of the pseudo-division is encapsulated into a function which stores in the list of certificates a certificate for each polynomial identified by its number. The certificate for each polynomial is in the form of a list of pairs of a polynomial and the number of a polynomial which we already know to be in the ideal : Certificates for the triangulation phase For the triangulation, we do not need to show that the result is triangulated, we just need to show that the polynomials in the triangulation are in the ideal generated by the hypotheses. We notice that the triangulation phase of Wu's method and its variants based on pseudo-division rely on an invariant: the method maintains a list of polynomials which belong to the ideal generated by the hypotheses. At the beginning of the triangulation procedure, the list consists of the original hypotheses themselves, hence they are in the ideal. At each step of the triangulation, the list of polynomials l is replaced by a new list l such that l and l differ only in one polynomial p ∈ l (let's call it p ∈ l ). This polynomial p is in the ideal generated by l.

The fact that p is in the ideal generated by l follows from the fact that p is computed using a pseudo-division p = prem(p, h) for some h in l.

Hence, to generate the certificate of the triangulation, we keep track of every intermediate polynomial thanks to its number, and we combine the certificates of the pseudo-divisions.

Certificates for the successive pseudo-division For the last phase (the successive pseudo-division of the goal by the triangulated polynomials), we just save the S i , and I as:

IG = i (S i T i )
where But, if during the combination of the different auxiliary certificates we evaluate the S i , then we get polynomials with thousands of terms. Polynomials of this size cannot be managed by Coq. Grégoire, Pottier and Théry proposed to use certificates involving not simply the list of the S i but a straight line program to evaluate the S i . This provides the possibility to have some let ... in instructions inside the certificates to allow to decompose and to share polynomial expressions. The validator is a function which computes IG and i (S i T i ) and which test the equality of these two polynomials.

S i = q i i-
Using their certificate structure we can define auxiliary polynomials which can then be used to define several other polynomials, etc.

To reduce the time required to check certificates, we make use of those let ... in instructions as most as possible. But as shown in the Sec. 4 the verification time of the certificates is still significant.

Toward more optimized certificates for ideal membership tests In this section we present some ideas for improvement (which are not yet implemented) of the structure for certificates.

One limitation of the certificates we use relies on the fact that every polynomial which is defined by a let ... in needs to be in the ideal generated by the previous polynomials. This bring no restriction for the polynomials generated by the Gröbner basis method, but in the context of Wu's method, it would be useful to be able to define polynomials using let ... in which are not in the ideal generated by the hypotheses. This would allow some sharing between the definition of I and the definition of the S i . Another idea to reduce drastically the size of the certificates would be to allow to have unevaluated pseudo-divisions in the certificates. This would require to prove formally the pseudo-division. But the pseudo-division is a time consuming task of Wu's method, hence this would imply to do a lot of computations during the validation of the certificate.

Benchmark

Table 1 presents our results compared to those of L. Pottier et al. based on Gröbner basis. The results show no clear winner, Wu's method is faster in some cases and Gröbner basis is faster in some other cases. We were surprised by the time needed to check the certificate in Coq. To study this, we reimplemented the validator in Ocaml and used it with our implementation of the method. The first column provides the computation time to produce the certificate and to check it when using an Ocaml validator : those timing are encouraging. The percentage of time required to check the certificates using the Ocaml validator varies from 1% to 80% but it is about 50% on average.

We did not expect such a difference of computation time between Ocaml and Coq. The data structure used by the Coq validator is a a modified Hörner form provided by the ring tactic [START_REF] Grégoire | Proving Equalities in a Commutative Ring Done Right in Coq[END_REF]. We tried another implementation based on lists of monomials, but the computation time was similar. The results could be improved in the future by the use of Native-Coq a Coq version using the native Ocaml compiler to perform strong reduction [START_REF] Boespflug | Full reduction at full throttle[END_REF]. In this paper, we described the first formalization of Wu's simple method in an interactive proof assistant. We implemented a version of Wu's method which generates certificates to ensure the correctness of the results generated by the method. We reused and extended the framework for certificate validation introduced by Benjamin Grégoire, Loïc Pottier and Laurent Théry, and we present some ideas for improvement. Thanks to the power of Wu's method, the new tactic we obtain does not only prove geometry theorems but also generates non degeneracy conditions necessary to prove a conjecture. Compared to other implementations of the method, as we formalize the invariance by translation and rotation of the statements, we prove the real geometry theorems and not a special case of the algebraic version of the statement given in a specific coordinate system.

In the future, we plan to extend our Ocaml implementation with a more optimized version of the pseudo-division (using GCD computations) along with the complete method of Wu-Ritt [START_REF] Wu | On the Decision Problem and the Mechanization of Theorem Proving in Elementary Geometry[END_REF][START_REF] Chou | Mechanical Geometry Theorem Proving[END_REF][START_REF] Shang | Ritt-Wu's Decomposition Algorithm and Geometry Theorem Proving[END_REF] or its variants developed by Dongming Wang [START_REF] Wang | Elimination Method[END_REF][START_REF] Wang | Elimination Practice[END_REF]. For educational applications our implementation still needs to be polished for example by extending the input language and by choosing automatically an origin and a direction for the system of coordinates.

The formalization of Wu's method in Coq enriches our framework for interactive and automatic theorem proving in geometry. We have now four methods for automated deduction in geometry available in Coq. It opens the door to several directions of research. Among them, it would be interesting to study the combination of Wu's method with other methods inside the same framework, for instance with the cylindrical algebraic decomposition method (CAD) whose formalization in Coq is under way by Assia Mahboubi [START_REF] Mahboubi | Contributions à la certification des calculs dans R : théorie, preuves, programmation[END_REF]. We could also study the different applications of Wu's method, such as solving geometric constraint systems [START_REF] Gao | Solving geometric constraint systems. I. A global propagation approach[END_REF][START_REF] Gao | Solving geometric constraint systems. II. A symbolic approach and decision of Rc-constructibility[END_REF] or finding loop invariants in the context of program verification.

Availability. The current version of the plug-in is available here: http://dpt-info.u-strasbg.fr/ ~narboux/nsatzwu.tar.gz

Fig. 2 .

 2 Fig. 2. Overview of the development

  -X D). Definition orthogonal A B C D := (X A -X B) * (X C -X D) + (Y A -Y B) * (Y C -Y D) = 0. Definition is_midpoint I A B := 2 * X I = X A + X B /\ 2 * Y I = Y A + Y B. Definition length_eq A B C D := (X A -X B) * (X A -X B) + (Y A -Y B) * (Y A -Y B) = (X C -X D) * (X C -X D) + (Y C -Y D) * (Y C -Y D). Definition is_in_intersection A B C D E := collinear A B C /\ collinear A D E.

  let pseudo_div_num a b x certif = let (r,c,d,s) = pseudo_div (a.p) (b.p) x in let new_n = new_num () in certif := (new_n, r, [(c^^d, a.n);(p_zero --s, b.n)]) ::(!certif); {p=r ; n= new_n}

  c j are the leading coefficients of the polynomials in the triangulation.

Table 1 .

 1 Benchmark : Computations are done using an Intel(R) Core(TM) i5 CPU 750 @ 2.67GHz with 4Gb RAM.

	Theorem Certificate gener-	Certificate gener-	Certificate gener-	Relative	speed
		ated using Wu	ated using Wu	ated using GB	of Wu's method
		and checking us-	and checking us-	and checking us-	over GB
		ing Ocaml	ing Coq	ing Coq
	Pascal 2	0.013	21	-	-
	Pascal 1	0.024	22	1652	×75
	Ptolemy95 0.010	10	30	×3
	Pappus	0.043	3	8	×2.6
	Altitudes 0.002	3	7	×2.3
	Simson	0.002	5	8	×1.6
	Perp-bisect 0.001	2	3	×1.5
	Pythagore 0.001	1	1	×1
	Feuerbach 0.038	15	15	×1
	Isoceles	0.001	1	1	×1
	Euler Line 0.063	9	6	×0.6
	Medians	0.001	3	2	×0.6
	Chords	0.015	4	2	×0.5
	Thales	0.003	6	3	×0.5
	Bissectors 0.001	6	3	×0.5
	Desargues 0.027	99	10	×0.1
	Ceva	0.025	98	6	×0.06

Note that this equation is actually equivalent to A = B or C = D or (AB) and (CD) are parallel. See the discussion on ndgs (Sec.3).

Note that there is a condition on the degree of R but this is not important in this discussion.

Acknowledgments. We wish to thank Loïc Pottier and Laurent Théry for having made their work publicly available and for the discussions we had.

This work is partially supported by the ANR project Galapagos.