Peter Csaba Ölveczky

Formalizing and Validating the P-Store Replicated Data Store in Maude

P-Store is a well-known partially replicated transactional data store that combines wide-area replication, data partition, some fault tolerance, serializability, and limited use of atomic multicast. In addition, a number of recent data store designs can be seen as extensions of P-Store. This paper describes the formalization and formal analysis of P-Store using the rewriting logic framework Maude. As part of this work, this paper specifies group communication commitment and defines an abstract Maude model of atomic multicast, both of which are key building blocks in many data store designs. Maude model checking analysis uncovered a non-trivial error in P-Store; this paper also formalizes a correction of P-Store whose analysis did not uncover any flaw.

Introduction

Large cloud applications-such as Google search, Gmail, Facebook, Dropbox, eBay, online banking, and card payment processing-are expected to be available continuously, even under peak load, congestion in parts of the network, server failures, and during scheduled hardware or software upgrades. Such applications also typically manage huge amounts of (potentially important user) data. To achieve the desired availability, the data must be replicated across geographically distributed sites, and to achieve the desired scalability and elasticity, the data store may have to be partitioned across multiple partitions.

Designing and validating cloud storage systems are hard, as the design must take into account wide-area asynchronous communication, concurrency, and fault tolerance. The use of formal methods during the design and validation of cloud storage systems has therefore been advocated recently [START_REF] Newcombe | How Amazon Web Services uses formal methods[END_REF][START_REF] Ölveczky | Design and validation of cloud computing data stores using formal methods[END_REF]. In [START_REF] Newcombe | How Amazon Web Services uses formal methods[END_REF], engineers at the world's largest cloud computing provider, Amazon Web Services, describe the use of TLA+ during the development of key parts of Amazon's cloud infrastructure, and conclude that the use of formal methods at Amazon has been a success. They report, for example, that: (i) "formal methods find bugs in system designs that cannot be found though any other technique we know of"; (ii) "formal methods [...] give good return on investment"; (iii) "formal methods are routinely applied to the design of complex real-world software, including public cloud services"; (iv) formal methods can analyze "extremely rare" combination of events, which the engineer cannot do, as "there are too many scenarios to imagine"; and (v) formal methods allowed Amazon to "devise aggressive optimizations to complex algorithms without sacrificing quality."

This paper describes the application of the rewriting-logic-based Maude language and tool [START_REF] Clavel | All About Maude[END_REF] to formally specify and analyze the P-Store data store [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF]. P-Store is a well-known partially replicated transactional data store that provides both serializability and some fault tolerance (e.g., transactions can be validated even when some nodes participating in the validation are down).

Members of the University of Illinois Center for Assured Cloud Computing have used Maude to formally specify and analyze complex industrial cloud storage systems such as Google's Megastore and Apache Cassandra [START_REF] Grov | Formal modeling and analysis of Google's Megastore in Real-Time Maude[END_REF][START_REF] Liu | Formal modeling and analysis of Cassandra in Maude[END_REF]. Why is formalizing and analyzing P-Store interesting? First, P-Store is a well-known data store design in its own right with many good properties that combines widearea replication, data partition, some fault tolerance, serializability, and limited use of atomic multicast. Second, a number of recent data store designs can be seen as extensions and variations of P-Store [START_REF] Sovran | Transactional storage for georeplicated systems[END_REF][START_REF] Ardekani | Non-monotonic snapshot isolation: Scalable and strong consistency for geo-replicated transactional systems[END_REF][START_REF] Ardekani | G-DUR: a middleware for assembling, analyzing, and improving transactional protocols[END_REF]. Third, it uses atomic multicast to order concurrent transactions. Fourth, it uses "group communication" for atomic commit. The point is that both atomic multicast and group communication commit are key building blocks in cloud storage systems (see, e.g., [START_REF] Ardekani | G-DUR: a middleware for assembling, analyzing, and improving transactional protocols[END_REF]) that have not been formalized in previous work. Indeed, one of the main contributions of this paper is an abstract Maude model of atomic multicast that allows any possible ordering of message reception consistent with atomic multicast.

I have modeled (both versions of) P-Store, and performed model checking analysis on small system configurations. Maude analysis uncovered some significant errors in the supposedly-verified P-Store algorithm, like read-only transactions never getting validated in certain cases. An author of the original P-Store paper [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF] confirmed that I had indeed found a nontrivial mistake in their algorithm and suggested a way of correcting the mistake. Maude analysis of the corrected algorithm did not find any error. I also found that a key assumption was missing from the paper, and that an important definition was very easy to misunderstand because of how it was phrased in English. All this emphasizes the need for a formal specification and formal analysis in addition to the standard prose-and-pseudo-code descriptions and informal correctness proofs.

The rest of the paper is organized as follows. Section 2 gives a background on Maude. Section 3 defines an abstract Maude model of the atomic multicast "communication primitive." Section 4 gives an overview of P-Store. Sections 5 and 6 present the Maude model and the Maude analysis, respectively, of P-Store, and Section 7 describes a corrected version of P-Store. Section 8 discusses some related work, and Section 9 gives some concluding remarks.

Due to space limitations, only parts of the specifications and analyses are given. I refer to the longer report [START_REF] Ölveczky | Formalizing and validating the P-Store replicated data store in Maude[END_REF] for more details. Furthermore, the executable Maude specifications of P-Store, together with analysis commands, are available at http://folk.uio.no/peterol/WADT16.

Preliminaries: Maude

Maude [START_REF] Clavel | All About Maude[END_REF] is a rewriting-logic-based formal language and simulation and model checking tool. A Maude module specifies a rewrite theory (Σ, E ∪ A, R), where:

-Σ is an algebraic signature; that is, a set of declarations of sorts, subsorts, and function symbols. -(Σ, E ∪ A) is a membership equational logic theory, with E a set of possibly conditional equations and membership axioms, and A a set of equational axioms such as associativity, commutativity, and identity. The theory (Σ, E∪ A) specifies the system's state space as an algebraic data type. -R is a set of labeled conditional rewrite rules3 l : t -→ t if m j=1 u j = v j specifying the system's local transitions. The rules are universally quantified by the variables in the terms, and are applied modulo the equations E ∪ A. 4I briefly summarize the syntax of Maude and refer to [START_REF] Clavel | All About Maude[END_REF] for more details. Operators are introduced with the op keyword: op f : s 1 . . . s n -> s. They can have user-definable syntax, with underbars '_' marking the argument positions, and equational attributes, such as assoc, comm, and id, stating, for example, that the operator is associative and commutative and has a certain identity element. Equations and rewrite rules are introduced with, respectively, keywords eq, or ceq for conditional equations, and rl and crl. The mathematical variables in such statements are declared with the keywords var and vars, or can be introduced on the fly having the form var:sort. An equation f (t 1 , . . . , t n ) = t with the owise ("otherwise") attribute can be applied to a term f (. . .) only if no other equation with left-hand side f (u 1 , . . . , u n ) can be applied. A class declaration class C | att 1 : s 1 , ... , att n : s n .

declares a class C with attributes att 1 to att n of sorts s 1 to s n . An object of class C is represented as a term < O : C | att 1 : val 1 , ..., att n : val n > of sort Object, where O, of sort Oid, is the object's identifier, and where val 1 to val n are the current values of the attributes att 1 to att n . A message is a term of sort Msg.

The state is a term of the sort Configuration, and is a multiset made up of objects and messages. Multiset union for configurations is denoted by a juxtaposition operator (empty syntax) that is declared associative and commutative, so that rewriting is multiset rewriting supported directly in Maude.

The dynamic behavior of concurrent object systems is axiomatized by specifying each of its transition patterns by a rewrite rule. For example, the rule A subclass inherits all the attributes and rules of its superclasses.

Formal Analysis in Maude. A Maude module is executable under some conditions, such as the equations being confluent and terminating, modulo the structural axioms, and the theory being coherent [START_REF] Clavel | All About Maude[END_REF]. Maude provides a range of analysis methods, including simulation for prototyping, search for reachability analysis, and LTL model checking. This paper uses Maude's search command

(search [[n]] t0 =>* pattern [such that cond ] .)
which uses a breadth-first strategy to search for at most n states that are reachable from the initial state t 0 , match the pattern pattern (a term with variables), and satisfy the (optional) condition cond . If '[n]' is omitted, then Maude searches for all solutions. If the arrow '=>!' is used instead of '=>*', then Maude searches for final states; i.e., states that cannot be further rewritten.

Atomic Multicast in Maude

Messages that are atomically multicast from (possibly) different nodes in a distributed system must be read in (pairwise) the same order: if nodes n 3 and n 4 both receive the atomically multicast messages m 1 and m 2 , they must receive (more precisely: "be served") m 1 and m 2 in the same order. Note that m 2 may be read before m 1 even if m 2 is atomically multicast after m 1 . Atomic multicast is typically used to order events in a distributed system. In distributed data stores like P-Store, atomic multicast is used to order (possibly conflicting) concurrent transactions: When a node has finished its local execution of a transaction, it atomically multicasts a validation request to other nodes (to check whether the transaction can commit). The validation requests therefore impose an order on concurrent transactions.

Atomic multicast does not necessarily provide a global order of all events. If each of the messages m 1 , m 2 , and m 3 is atomically multicast to two of the receivers A, B, and C, then A can read m 1 before m 2 , B can read m 2 before m 3 , and C can read m 3 before m 1 . These reads satisfy the pairwise total order requirement of atomic multicast, since there is no conflict between any pair of receivers. Nevertheless, atomic multicast has failed to globally order the messages m 1 , m 2 , and m 3 . If atomic multicast is used to impose something resembling a global order (e.g., on transactions), it should also satisfy the following uniform acyclic order property: the relation < on (atomic-multicast) messages is acyclic, where m < m holds if there exists a node that reads m before m . Atomic multicast is an important concept in distributed systems, and there are a number of well-known algorithms for achieving atomic multicast [START_REF] Guerraoui | Genuine atomic multicast in asynchronous distributed systems[END_REF]. To model P-Store, which uses atomic multicast, I could of course formalize a specific algorithm for atomic multicast and include it in a model of P-Store. Such a solution would, however, have a number of disadvantages, including:

1. Messy non-modular specifications. Atomic multicast algorithms involve some complexity, including maintaining Lamport clocks during system execution, keeping buffers of received messages that cannot be served, and so on. This solution could also easily yield a messy non-modular specification that fails to separate the specification of P-Store from that of atomic multicast. 2. Increased state space. Running a distributed algorithm concurrently with P-Store would also lead to much larger state spaces during model checking analyses, since also the states generated by the rewrites involving the atomic multicast algorithm would contribute to new states. 3. Lack of generality. Implementing a particular atomic multicast algorithm might exclude behaviors possible with other algorithms. That would mean that model checking analysis might not cover all possible behaviors of P-Store, but only those possible with the selected atomic multicast algorithm.

I therefore instead define, for each of the two "versions" of atomic multicast, a general atomic multicast primitive, which allows all possible ways of reading messages that are consistent with the selected version of atomic multicast. In particular, such a solution will not add states during model checking analysis.

Atomic Multicast in Maude: "User Interface"

To define an atomic multicast primitive, the system maintains a "table" of read and sent-but-unread atomic-multicast messages for each node. This table must be consulted before reading an atomic-multicast message, to check whether it can be read/served already, and must be updated when the message is read.

The "user interface" of my atomic multicast "primitive" is as follows:

-Atomically multicasting a message. A node n that wants to atomically multicast a message m to a set of nodes {n 1 , . . . , n k } just "sends" the "message" -The user must add the term [emptyAME] (denoting the "empty" atomic multicast table) to the initial state.

Maude Specification of Atomic Multicast

To keep track of atomic-multicast messages sent and received, the table The "wrapper" used for atomic multicast takes as arguments the message (content), the sender's identifier, and the (identifiers of the) set of receivers: "distributes" such an atomic-multicast msg from o to o 1 ... o n message by: (1) "dissolving" the above multicast message into a set of messages The update function, which updates the atomic-multicast table when O reads a message MC, just moves MC from the set of unread messages to the end of the list of read messages in O's record in the table.

The expression okToRead(mc, o, amTable) is used to check whether the object o can read the atomic-multicast message mc with the given global atomicmulticast table amTable. The function okToRead is defined differently depending on whether atomic multicast must satisfy the uniform acyclic order requirement.

okToRead for Pairwise Total Order Atomic Multicast. The following equations define okToRead by first characterizing the cases when the message cannot be read; the last equation uses Maude's owise construct to specify that the message can be read in all other cases: In the first equation, O wants to read MC, and its AM-entry shows that O has not read message MC2. However, another object O2 has already read MC2 before MC, which implies that O cannot read MC. In the second equation some object O2 has read MC2 and has MC in its sets of unread atomic-multicast messages, which implies that O cannot read MC yet (it must read MC2 first).

okToRead for Uniform Acyclic Order Atomic Multicast. To define atomic multicast which satisfies the uniform acyclic order requirement, the above definition must be generalized to consider the induced relation < instead of pairwise reads. The above definition checks whether a node o can read a message m 1 by checking whether it has some other unread message m 2 pending such reading m 1 before m 2 would conflict with the m 1 /m 2 -reading order of another node. This happens if another node has read m 2 before reading m 1 , or if it has read m 2 and has m 1 pending (which implies that eventually, that object would read m 2 before m 1 ). In the more complex uniform acyclic order setting, that solution must be generalized to check whether reading m 1 before any other pending message m 2 would violate the current or the (necessary) future "global order." That is, is there some m 1 elsewhere that has been read or must eventually be read after m 2 somewhere? If so, node o obviously cannot read m 1 at the moment.

The function receivedAfter takes a set of messages and the global AMtable as arguments, and computes the < * -closure of the original set of messages; i.e., the messages that cannot be read before the original set of messages: In the above equation, there is a message MC in the current set of messages in the closure. Furthermore, the global atomic-multicast table shows that some node O2 has read MC2 right after reading MC, and MC2 is not yet in the closure. Therefore, MC2 is added to the closure.

op receivedAfter : MsgContSet AM-Table -> MsgContSet .
In the following equation, there is a message MC in the closure; furthermore, some object O2 has already read MC. This implies that all unread messages MCS2 of O2 must eventually be read after MC, and hence they are added to the closure: Finally, the current set is returned when it cannot be extended:

eq receivedAfter(MCS, AM-TABLE) = MCS [owise] .
The function okToRead can then be defined as expected: O can read the pending message MC if MC is not (forced to be) read after any other pending message (in the set MCS): I have model-checked both specifications of atomic multicast on a number of scenarios and found no deadlocks or inconsistent multicast read orders.

P-Store

P-Store [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF] is a partially replicated data store for wide-area networks developed by Schiper, Sutra, and Pedone that provides transactions with serializability. P-Store executes transactions optimistically: the execution of a transaction T at site s (which may involve remote reads of data items not replicated at s) proceeds without worrying about conflicting concurrent transactions at other sites. After the transaction T has finished executing, a certification process is executed to check whether or not the transaction T was in conflict with a concurrent transaction elsewhere, in which case T might have to be aborted. More precisely, in the certification phase the site s atomically multicasts a request to certify T to all sites storing data accessed by T . These sites then perform a voting procedure to decide whether T can commit or has to be aborted.

P-Store has a number of attractive features: (i) it is a genuine protocol: only the sites replicating data items accessed by a transaction T are involved in the certification of T ; and (ii) P-Store uses atomic multicast at most once per transaction. Another issue in the certification phase: in principle, the sites certify the transactions in the order in which the certification requests are read. However, if for some reason the certification of the first transaction in a site's certification queue takes a long time (maybe because other sites involved in the voting are still certifying other transactions), then the certification of the next transaction in line will be delayed accordingly, leading to the dreaded convoy effect. P-Store has an "advanced" version that tries to mitigate this problem by allowing a site to start the certification also of other transactions in its certification queue, as long as they are not in a possible conflict with "older" transactions in that queue.

The authors of [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF] claim that they have proved the P-Store algorithm correct.

P-Store in Detail

This section summarizes the description of P-Store in [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF].

System Model and Assumptions. A database is a set of triples (k, v, ts), where k is a key, v its value, and ts its time stamp. Each site holds a partial copy of the database, with Items(s) denoting the keys replicated at site s. I do not consider failures in this paper (as failure treatment is not described in the algorithms in [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF]). A transaction T is a sequence of read and write operations, and is executed locally at site proxy(T ). Items(T ) is the set of keys read or written by T ; WReplicas(T ) and Replicas(T ) denote the sites replicating a key written, respectively read or written, by T . A transaction T "is local iff for any site s in Replicas(T ), Items(T ) ⊆ Items(s); otherwise, T is global." Each site ensures order-preserving serializability of its local executions of transactions. As already mentioned, P-Store assumes access to an atomic multicast service that guarantees uniform acyclic order.

Executing a Transaction. While a transaction T is executing (at site proxy(T )), a read on key k is executed at some site that stores k; k and the item time stamp ts read are stored as a pair (k, ts) in T 's read set T.rs. Every write of value v to key k is stored as a pair (k, v) in T 's set of updates T.up. If T reads a key that was previously updated by T , the corresponding value in T.up is returned.

When T has finished executing, it can be committed immediately if T is read-only and local. Otherwise, we need to run the certification protocol, which also propagates T 's updates to the other (write-) replicating sites.

If the certification process, described next, decides that T can commit, all sites in WReplicas(T ) apply T 's updates. In any case, proxy(T ) is notified about the outcome (commit or abort) of the certification. Certification Phase. When T is submitted for certification, T is atomically multicast to all sites storing keys read (to check for stale reads) or written (to propagate the updates) by T . When a site s reads such a request, it checks whether the values read by T are up-to-date by comparing their versions against those currently stored in the database. If they are the same, T passes the certification test; otherwise T fails at s. The site s may not replicate all keys read by T and therefore may not be able to certify T . In this case there is a voting phase where each site s replicating keys read by T sends the result of its local certification test to all sites s w replicating a key written by T . A site s w can decide on T 's outcome when it has received (positive) votes from a voting quorum for T , i.e., a set of sites that together replicate all keys read by T . If some site votes "no," the transaction must be aborted. The pseudo-code description of this certification algorithm in [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF] is shown in Fig. 1.

As already mentioned, a site does not start the certification of another transaction until it is done certifying the first transaction in its certification queue. To avoid the convoy effect that this can lead to, the paper [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF] also describes a version of P-Store where different transactions in a site's certification queue can be certified concurrently as long as they do not read-write conflict.

Formalizing P-Store in Maude

I have formalized both versions of P-Store (i.e., with and without sites initiating multiple concurrent certifications) in Maude, and present parts of the formalization of the simpler version. The executable specifications of both versions, with analysis commands, are available at http://folk.uio.no/peterol/ WADT16, and the longer report [START_REF] Ölveczky | Formalizing and validating the P-Store replicated data store in Maude[END_REF] provides more detail.

Class Declarations

Transactions. Although the actual values of keys in the databases are sometimes ignored during analysis of distributed data stores, I choose for purposes of illus-tration to represent the concrete values of keys (or data items). This should not add new states that would slow down the model checking analysis.

A transaction (sometimes also called a transaction request) is modeled as an object of the following class Transaction: The operations attribute denotes the list of read and write operations that remain to be executed. Such an operation is either a read operation x := read k, where x is a "local variable" that stores the value of the (data item with) key k read by the operation, or a write operation write(k, expr ), where expr in our case is a simple arithmetic expression involving the transaction's local variables. waitRemote(k, x) is an "internal operation" denoting that the transaction execution is awaiting the value of a key k (to be assigned to the local variable x) which is not replicated by the transaction's proxy. An operation list is a list of such operations, with list concatenation denoted by juxtaposition. destination denotes the (identity of the) proxy of the transaction; that is, the site that should execute the transaction. The readSet attribute denotes the ','-separated set of pairs versionRead(k, version), each such pair denoting that the transaction has read version version of the key k. The writeSet attribute denotes the write set of the transaction as a map (k 1 |-> val 1 ), ..., (k n |-> val n ). The status attribute denotes the commit state of the transaction, which is either commit, abort, or undecided. Finally, localVars is a map from the transaction's local variables to their current values.

Replicas. A replicating site (or site or replica) stores parts of the database, executes the transactions for which it is the proxy, and takes part in the certification of other transactions. A replica is formalized as an object instance of the following subclass Replica: The datastore attribute represents the replica's local database as a set < key 1 , val 1, ver 1 > , . . . , < key l , val l , ver l > of triples < key i , val i, ver i > denoting a version of the data item with key key i , value val i , and version number ver i . 5 The attributes executing, submitted, committed, and aborted denote the transactions executed by the replica and which are/have been, respectively, currently executing, submitted for certification, committed, and aborted. The queue holds the certification queue of transactions to be certified by the replica (in collaboration with other replicas). transToCertify contains data used for the certification of the first element in the certification queue (in the simpler algorithm), and decidedTranses show the status (aborted/committed) of the transactions that have previously been (partly) certified by the replica.

Clients. Finally, I add an "interface/application layer" to the P-Store specification in the form of clients that send transactions to be executed by P-Store:

class Client | txns : ObjectList, pendingTrans : TransIdSet .

txns denotes the list of transaction (objects) that the client wants P-Store to execute, and pendingTrans is either the empty set or (the identity of) the transaction the client has submitted to P-Store but whose execution is not yet finished.

Initial State.

The following shows an initial state init4 (with some parts replaced by '...') used in the analysis of P-Store. This system has: two clients, c1 and c2, that want P-Store to execute the two transactions t1 and t2; three replicating sites, r1, r2, and r3; and three data items/keys x, y, and z. Transaction t1 wants to execute the operations (xl :=read x) (yl :=read y) at replica r1, while transaction t2 wants to execute write(y, 5) write(x, 8) at replica r2. The initial state also contains the empty atomic multicast table and the table which assigns to each key the sites replicating this key. Initially, the value of each key is [START_REF] Ardekani | G-DUR: a middleware for assembling, analyzing, and improving transactional protocols[END_REF] and its version is 1. Site r2 replicates both x and y.

Local Execution of a Transaction

The execution of a transaction has two phases. In the first phase, the transaction is executed locally by its proxy: the transaction performs its reads and writes, but the database is not updated; instead, the reads are recorded in the transaction's read set, and its updates are stored in the writeSet attribute.

The second phase is the certification (or validation) phase, when all appropriate nodes together decide whether or not the transaction can be committed or must be aborted. If it can be committed, the replicas update their databases.

This section specifies the first phase, which starts when a client without pending transactions sends its next transaction to its proxy. I do not show the variable declarations (see [START_REF] Ölveczky | Formalizing and validating the P-Store replicated data store in Maude[END_REF]), but follow the convention that variables are written with (all) capital letters. P-Store assumes that the local executions of multiple transactions on a site are equivalent to some serialized executions. I model this assumption by executing the transactions one-by-one. Therefore, a replica can only receive a transaction request if its set of currently executing transactions is empty (none): There are three cases to consider when executing a read operation X :=read K: (i) the transaction has already written to key K; (ii) the transaction has not written K and the proxy replicates K; or (iii) the key K has not been read and the proxy does not replicate K. I only show the specification for case (i). I do not know what version number should be associated to the read, and I choose not to add the item to the read set. (The paper [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF] does not describe what to do in this case; the problem disappears if we make the common assumption that a transaction always reads a key before updating it.) As an effect, the local variable X gets the value V: Write operations are easy: evaluate the expression EXPR to write and add the update to the transaction's writeSet:

Certification Phase

When all the transaction's operations have been executed by the proxy, the proxy's next step is to try to commit the transaction. If the transaction is readonly and local, it can be committed directly; otherwise, it must be submitted to the certification protocol. Some colleagues and I all found the definition of local in [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF] (and quoted in Section 4) to be quite ambiguous. We thought that "for any site s in Replicas(T ), Items(T ) ⊆ Items(s)" means either "for each site s . . . " or that proxy(T ) replicates all items in T . The first author of [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF], Nicolas Schiper, told me that it actually means "for some s . . . ." In hindsight, we see that this is also a valid interpretation of the definition of local . To avoid misunderstanding, it is probably good to avoid the phrase "for any" and use either "for each" or "for some."

If the transaction T cannot be committed immediately, it is submitted for certification by atomically multicasting a certification request-with the transaction's identity TID, read set RS, and write set WS-to all replicas storing keys read or updated by T (lines 9-10 in Fig. 1 According to lines 7-8 in Fig. 1, a replica's local certification succeeds if, for each key in the transaction's read set that is replicated by the replica in question, the transaction read the same version stored by the replica: op certificationOk : ReadSet DataStore -> Bool . eq certificationOk((versionRead(K, VERSION) , READSET), (< K, V, VERSION2 > , DS)) = (VERSION == VERSION2) and certificationOk(READSET, (< K, V, VERSION2 > , DS)) . eq certificationOk(RS, DS) = true [owise] .

If the transaction to certify is not local, the certifying sites must together decide whether or not the transaction can be committed. Each certifying site therefore checks whether the transaction passes the local certification test, and sends the outcome of this test to the other certifying sites (lines 13 and 19-22): If the local certification fails, the site sends an abort vote to the other write replicas and also notifies the proxy of the outcome. Otherwise, the site sends a commit vote to all other site replicating an item written by the transaction. The voting phase ends when there is a voting quorum; that is, when the voting sites together replicate all keys read by the transaction. This means that a certifying site must keep track of the votes received during the certification of a transaction. The set of sites from which the site has received a (positive) vote is the fourth parameter of the certify record it maintains for each transaction. If a site receives a positive vote, it stores the sender of the vote (lines [START_REF] Ölveczky | Design and validation of cloud computing data stores using formal methods[END_REF][START_REF] Ölveczky | Formal modeling, performance estimation, and model checking of wireless sensor network algorithms in Real-Time Maude[END_REF]. If a site receives a negative vote, it decides the fate of the transaction and notifies the proxy if it replicates an item written by the transaction (lines 28-29).

If a write replica has received positive votes from a voting quorum (lines 23-27 and 29), the transaction can be committed, and the write replica applies the updates and notifies the proxy. The following rule models the behavior when a site has received votes from a voting quorum RIDS for transaction TID: Finally, the proxy of transaction TID receives the outcome from one or more sites in TID's certification set (the abort case is similar): ---notify client

6 Formal Analysis of P-Store

In the absence of failures, P-Store is supposed to guarantee serializability of the committed transactions, and that a decision (commit/abort) is made on all transactions.

To analyze P-Store, I search for all final states-i.e., states that cannot be further rewritten-reachable from a given initial state, and inspect the result. This analysis therefore also discovers undesired deadlocks. In the future, I should instead automatically check serializability, possibly using the techniques in [START_REF] Grov | Formal modeling and analysis of Google's Megastore in Real-Time Maude[END_REF], which adds to the state a "serialization graph" that is updated whenever a transaction commits, and then checks whether the graph has cycles.

The search for final states reachable from state init4 in Section 5.1 yields a state which shows that t1's proxy is not notified about the outcome of the certification (see [START_REF] Ölveczky | Formalizing and validating the P-Store replicated data store in Maude[END_REF] for details). The problem seems to be line 29 in the algorithm in Fig. 1: only sites replicating items written by transaction T (WReplicas(T )) send the outcome of the certification to T 's proxy. It is therefore not surprising that the outcome of the read-only transaction t1 does not reach t1's proxy.

The transactions in init4 are local. What about non-local transactions? The initial state init5 is the same as init4 in Section 5.1, except that item y is only replicated at site r3, which means that t1 and t2 become non-local transactions.

Searching for final states reachable from init5 shows a result where the certification process cannot reach a decision on the outcome of transaction t1: The fate of t1 is not decided: both r2 and r3 are stuck in their certification process. The problem seems to be lines 22 and 23 in the P-Store certification algorithm: why are only write replicas involved in sending and receiving votes during the certification? Shouldn't both read and write replicas vote? Otherwise, it is impossible to certify non-local read-only transactions, such as t1 in init5. are read from the same site and some additional concurrency control is used to ensure serializability), but admitted that this is indeed not mentioned anywhere in their paper. My specifications consider the algorithm as given in [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF], without taking the unstated assumptions into account, and also subjects the local read-only transactions to certification.

Analysis of the Updated Specification. I have analyzed the corrected specification on five small initial configurations (3 sites, 3 data items, 2 transactions, 4 operations). All the final states were correct: the committed transactions were indeed serializable.

The Advanced Algorithm. I have also specified and successfully analyzed the (corrected) version of P-Store where multiple transactions can be certified concurrently. It is beyond the scope of this paper to describe that specification.

Related Work

Different communication forms/primitives have been defined in Maude, including wireless broadcast that takes into account the geographic location of nodes and the transmission strength/radius [START_REF] Ölveczky | Formal modeling, performance estimation, and model checking of wireless sensor network algorithms in Real-Time Maude[END_REF], as well as wireless broadcast in mobile systems [START_REF] Liu | Modeling and analyzing mobile ad hoc networks in Real-Time Maude[END_REF]. However, I am not aware of any model of atomic multicast in Maude. Maude has been applied to a number of industrial and academic cloud storage systems, including Google's Megastore [START_REF] Grov | Formal modeling and analysis of Google's Megastore in Real-Time Maude[END_REF], Apache Cassandra [START_REF] Liu | Formal modeling and analysis of Cassandra in Maude[END_REF], and UC Berkeley's RAMP [START_REF] Liu | Formal modeling and analysis of RAMP transaction systems[END_REF]. However, that work did not address issues like atomic multicast and group communication commit.

Lamport's TLA+ has also been used to specify and model check large industrial cloud storage systems like S3 at Amazon [START_REF] Newcombe | How Amazon Web Services uses formal methods[END_REF] and the academic TAPIR transaction protocol targeting large-scale distributed storage systems.

On the validation of P-Store and similar designs, P-Store itself has been proved to be correct using informal "hand proofs" [START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF]. However, such hand proofs do not generate precise specifications of the systems and tend to be error-prone and rely on missing assumptions, as I show in this paper. I have not found any model checking validation of related designs, such as Jessy [START_REF] Ardekani | Non-monotonic snapshot isolation: Scalable and strong consistency for geo-replicated transactional systems[END_REF] and Walter [START_REF] Sovran | Transactional storage for georeplicated systems[END_REF].

Concluding Remarks

Cloud computing relies on partially replicated wide-area data stores to provide the availability and elasticity required by cloud systems. P-Store is a well-known such data store that uses atomic multicast, group communication commitment, concurrent certification of independent transactions, etc. Furthermore, many other partially replicated data stores are extensions and variations of P-Store.

I have formally specified and analyzed P-Store in Maude. Maude reachability analysis uncovered a numbers of errors in P-Store that were confirmed by one of the P-Store developers: both read and write replicas need to participate in the certification of transactions; write replicas are not enough. I have specified the proposed fix of P-Store, whose Maude analysis did not uncover any error.

Another main contribution of this paper is a general and abstract Maude "primitive" for both variations of atomic multicast.

One important advantage claimed by proponents of formal methods is that even precise-looking informal descriptions tend to be ambiguous and contain missing assumptions. In this paper I have pointed at a concrete case of ambiguity in a precise-looking definition, and at a crucial missing assumption in P-Store.

This work took place in the context of the University of Illinois Center for Assured Cloud Computing, within which we want to identify key building blocks of cloud storage systems, so that they can be built and verified in a modular way by combining such building blocks in different ways. Some of those building blocks are group communication commitment certification and atomic multicast. In more near term, this work should simplify the analysis of other state-of-the-art data stores, such as Walter and Jessy, that can be seen as extensions of P-Store.

The analysis performed was performed using reachability analysis; in the future one should also be able to specify the desired consistency property "directly."

  rl [l] : m(O,w) < O : C | a1 : x, a2 : O', a3 : z > => < O : C | a1 : x + w, a2 : O', a3 : z > m'(O',x) .

  op atomic-multicast_from_to_ : MsgCont Oid OidSet -> Configuration . The equation eq (atomic-multicast MC from O to OS) [AM-ENTRIES] = (distribute MC from O to OS) [insert(MC, OS, AM-ENTRIES)] .

(

  msg msg from o to o1) ... (msg msg from o to on), one for each receiver o k in the set {o 1 , . . . , o n }; and (2) by adding, for each receiver o k , the message (content) msg to the set unread k of unread atomicmulticast messages in the atomic-multicast table.

  vars MC MC2 : MsgContent . vars MCS MCS2 : MsgContSet . vars MCL MCL2 MCL3 MCL4 : MsgContList . eq okToRead(MC, O, [am-entry(O, MCL, MCS MC MC2) am-entry(O2, MCL2 :: MC2 :: MCL3 :: MC :: MCL4, MCS2) AM-ENTRIES]) = false . eq okToRead(MC, O, [am-entry(O, MCL, MCS MC MC2) am-entry(O2, MCL2 :: MC2 :: MCL4, MCS2 MC) AM-ENTRIES]) = false . eq okToRead(MC, O, [AM-ENTRIES]) = true [owise] .

  ceq receivedAfter(MC MCS, [am-entry(O2, MCL :: MC :: MC2 :: MCL2, MCS2) AM-ENTRIES]) = receivedAfter(MC MCS MC2, [am-entry(O2, MCL :: MC :: MC2 :: MCL2, MCS2) AM-ENTRIES]) if not (MC2 in MCS) .

  ceq receivedAfter(MC MCS, [am-entry(O2, MCL2 :: MC :: MCL4, MCS2) AM-ENTRIES]) = receivedAfter(MC MCS MCS2, [am-entry(O2, MCL2 :: MC :: MCL4, emptyMsgContSet) AM-ENTRIES]) if MCS2 =/= emptyMsgContSet .

  eq okToRead(MC, O, [am-entry(O, MCL, MCS MC) AM-ENTRIES]) = not (MC in receivedAfter(MCS, [am-entry(O, MCL, MCS) AM-ENTRIES])) .

Fig. 1 .

 1 Fig. 1. The P-Store certification algorithm in [14].

  x, r2) ;; replicatingSites(y, (r2 , r3)) ;; replicatingSites(z, r1)] < c1 : Client | txns : < t1 : Transaction | operations : ((xl :=read x) (yl :=read y)), destination : r1, readSet : emptyReadSet, status : undecided, writeSet : emptyWriteSet, localVars : (xl |-> [0] , yl |-> [0]) >, pendingTrans : empty > < c2 : Client | txns : < t2 : Transaction | operations : (write(y, 5) write(x, 8)), destination : r2, ... > pendingTrans : empty > < r1 : Replica | datastore : (< z, [2], 1 >), committed : none, aborted : none, executing : none, submitted : none, queue : emptyTransList, transToCertify : noTrans, decidedTranses : noTS > < r2 : Replica | datastore : ((< x, [2], 1 >) , (< y, [2], 1 >)), ... > < r3 : Replica | datastore : (< y, [2], 1 >), ... > .

  rl [sendTxn] : < C : Client | pendingTrans : empty, txns : < TID : Transaction | destination : RID > ; TXNS > => < C : Client | pendingTrans : TID, txns : TXNS > (msg executeTrans(< TID : Transaction | >) from C to RID) .

  rl [receiveTxn] : (msg executeTrans(< TID : Transaction | >) from C to RID) < RID : Replica | executing : none > => < RID : Replica | executing : < TID : Transaction | > > .

  rl [executeRead1] : < RID : Replica | executing : < TID : Transaction | operations : (X :=read K) OPLIST, writeSet : (K |-> V), WS, localVars : VARS > > => < RID : Replica | executing : < TID : Transaction | operations : OPLIST, localVars : insert(X, V, VARS) > > .

  rl [executeWrite] : < RID : Replica | executing : < TID : Transaction | operations : write(K, EXPR) OPLIST, localVars : VARS, writeSet : WS > > => < RID : Replica | executing : < TID : Transaction | operations : OPLIST, writeSet : insert(K, eval(EXPR, VARS), WS) > > .

  rl [readCommit] : (msg commit(TID) from RID2 to RID) < RID : Replica | submitted : < TID : Transaction | >, committed : TRANSES > => < RID : Replica | submitted : none, committed : (TRANSES < TID : Transaction | >) > done(TID) .

  Maude> (search init5 =>! C:Configuration .) ... Solution 4 ... < r1 : Replica | submitted : < t1 : Transaction | localVars :(xl |->[8], yl |->[5]), operations : nil, readSet : versionRead(x,2), versionRead(y,2), ... > , transToCertify : noTrans > < r2 : Replica | committed : < t2 : Transaction | writeSet : (x |-> [8], y |-> [5]), ... >, datastore : < x,[8],2 >, decidedTranses : transStatus(t2,commit), transToCertify : certify(t1,r1,(versionRead(x,2),versionRead(y,2)), emptyWriteSet,r2) , ... > < r3 : Replica | aborted : none, committed : none, datastore : < y,[5],2 >, decidedTranses : transStatus(t2,commit), transToCertify : certify(t1, r1, ..., emptyWriteSet, r3) , ... >

  defines a family of transitions in which a message m, with parameters O and w, is read and consumed by an object O of class C, the attribute a1 of the object

	O is changed to x + w, and a new message m'(O',x) is generated. Attributes
	whose values do not change and do not affect the next state of other attributes
	or messages, such as a3, need not be mentioned in a rule. Likewise, attributes
	that are unchanged, such as a2, can be omitted from right-hand sides of rules.

  ):

	crl [commit/submit2] :
	< RID : Replica | executing :
	< TID : Transaction | operations : nil, readSet : RS, writeSet : WS >,
	submitted : TRANSES >
	REPLICA-TABLE
	=>
	< RID : Replica | executing : none, submitted : TRANSES < TID : Transaction | > >
	REPLICA-TABLE
	(atomic-multicast certify(TID, RS, WS) from RID
	to replicas((keys(RS) , keys(WS)), REPLICA-TABLE))
	if WS =/= emptyWriteSet or not localTrans(keys(RS), REPLICA-TABLE) .

  Nicolas Schiper confirmed that the errors pointed out in Section 6 are indeed errors in P-Store. He also suggested the fix alluded to in Section 6: replace WReplicas(T ) with Replicas(T ) in lines 22, 23, and 29. The Maude specification of the proposed correction is given in http://folk.uio.no/peterol/WADT16/.Missing Assumptions. One issue seems to remain: why can read-only local transactions be committed without certification? Couldn't such transactions have read stale values? Nicolas Schiper kindly explained that local read-only transactions are handled in a special way (all values

	7 Fixing P-Store

An equational condition ui = wi can also be a matching equation, written ui:= wi, which instantiates the variables in ui to the values that make ui = wi hold, if any.

Operationally, a term is reduced to its E-normal form modulo A before a rewrite rule is applied.

The paper[START_REF] Schiper | P-Store: Genuine partial replication in wide area networks[END_REF] does not specify whether a replica stores multiple versions of a key.

Acknowledgments. I would like to thank Nicolas Schiper for quick and friendly replies to my questions about P-Store, the anonymous reviewers for helpful comments, and Si Liu and José Meseguer for valuable discussions about P-Store and atomic multicast.

This work was partially supported by AFOSR/AFRL Grant FA8750-11-2-0084 and NSF Grant CNS 14-09416.


Tom Lyche

Jean-Louis Merrien
email: jmerrien@insa-rennes.fr

J.-L Merrien

Simplex-Splines on the Clough-Tocher Element

Keywords: Triangle Mesh, Piecewise polynomials, Interpolation, Simplex Splines, Marsden-like Identity

We propose a simplex spline basis for a space of C 1 -cubics on the Clough-Tocher split on a triangle. The 12 elements of the basis give a nonnegative partition of unity. We derive two Marsden-like identities, three quasi-interpolants with optimal approximation order and prove L ∞ stability of the basis. The conditions for C 1 -junction to neighboring triangles are simple and similar to the C 1 conditions for the cubic Bernstein polynomials on a triangulation. The simplex spline basis can also be linked to the Hermite basis to solve the classical interpolation problem on the Clough-Tocher split.

Introduction

Piecewise polynomials over triangles have applications in several branches of the sciences ranging from finite element analysis, surfaces in computer aided design and other engineering problems. For many of these applications, piecewise linear C 0 surfaces do not suffice. In some cases, we need smoother surfaces for modeling, or higher degrees to increase the approximation order. To obtain C 1 smoothness on an arbitrary triangulation, one needs piecewise quintic polynomials, [START_REF] Lai | Spline Functions on Triangulations[END_REF]. We can use lower degrees if we are willing to split each triangle into a number of subtriangles. Examples are the Clough-Tocher split (CT), [START_REF] Clough | Finite element stiffness matrices for analysis of plate bending[END_REF] and the Powell-Sabin 6 and 12-splits (PS6, PS12), [START_REF] Powell | Piecewise quadratic approximation on triangles[END_REF]. The 1 number of subtriangles is 3, 6 and 12 for CT, PS6 and PS12, respectively. A B-spline like basis both for C 1 cubics and C 2 quintics has been constructed for PS6, [START_REF] Grošelj | Construction and analysis of cubic Powell-Sabin B-splines[END_REF][START_REF] Speleers | A normalized basis for quintic Powell Sabin splines[END_REF] and references therein. Recently a B-spline like basis has also been proposed for a 9 dimensional subspace of C 1 cubics on CT, [START_REF] Grošelj | Construction and analysis of cubic Powell-Sabin B-splines[END_REF]. The PS12-split can be defined as the complete graph obtained by connecting vertices and edge midpoints of each triangle. A B-spline basis for PS12 and the full space of C 1 -cubics on CT seem difficult. An alternative to the B-spline basis is the Hermite basis. Since it uses both values and derivatives it is not as stable as the B-spline basis and it does not form a nonnegative partition of unity.

Here we construct a B-spline basis for one triangle in the coarse triangulation and connect to neighboring triangles using Bernstein-Bézier techniques. This was done for PS12 using C 1 quadratics, [START_REF] Cohen | A B-spline-like basis for the Powell-Sabin 12-split based on simplex splines[END_REF], and C 2 and C 3 quintics, [START_REF] Lyche | A Hermite interpolatory subdivision scheme for C 2 -quintics on the Powell-Sabin 12-split[END_REF][START_REF] Lyche | Stable Simplex Spline Bases for C 3 Quintics on the Powell-Sabin 12-Split[END_REF]. These bases, consisting of simplex splines (see for example [START_REF] Micchelli | On a numerically efficient method for computing multivariate B-splines[END_REF] for a general introduction), all share attractive properties of univariate B-splines such as In this paper we consider the full 12 dimensional space of C 1 cubics on the CT-split. We will define a simplex spline basis for this split and show that it has all the B-spline and Bernstein-Bézier properties mentioned above.

The CT-split is interesting for many reasons. To obtain a space of C 1 piecewise polynomials of degree at most 3 on an arbitrary triangulation, we only need to divide each triangle into 3 subtriangles, while 6 and 12 subtriangles are needed for PS6 and PS12. Moreover, the approximation order of the space S 3 of piecewise C 1 cubics on CT is 4 and this is at least as good as for the spaces S 6 and S 12 of piecewise cubics on PS6 and piecewise quadratics on PS12. The degrees of freedom for S 6 are values and gradients of the vertices of the coarse triangulation while for S 3 and S 12 we need in addition cross boundary derivatives at the midpoint of the edges, see Figure 1 (left). For further comparisons of these three spaces see Section 6.6 in [START_REF] Lai | Spline Functions on Triangulations[END_REF].

This paper is organized as follows: In the remaining part of the introduction, we review some properties of CT, introduce our notation and recall the main properties of simplex splines. In Section 2, we construct a cubic simplex spline basis for CT, from which, in Section 3, we derive two Marsden identities and then, in Section 4, three quasi-interpolants, and show L ∞ stability of the basis. In Section 5, conditions to ensure C 0 and C 1 continuity through an edge between two triangles are derived. The conversion between the simplex spline basis and the Hermite basis for CT is considered in Section 6. Lagrange and Hermite interpolation on triangulations using C 1 cubics, quartics and higher degrees have also been considered in [START_REF] Davydov | Interpolation by Splines on Triangulations[END_REF]. We end the paper with numerical examples of interpolation on a triangulation.

The Clough-Tocher split

To describe this split, let T := p 1 , p 2 , p 3 be a nondegenerate triangle in R 2 . Using the barycenter p T := (p 1 + p 2 + p 3 )/3 we can split T into three subtriangles T 1 := p T , p 2 , p 3 , T 2 := p T , p 3 , p 1 and T 3 := p T , p 1 , p 2 . On T we consider the space

S 1 3 ( ) := {f ∈ C 1 (T ) : f |T i is a polynomial of degree at most 3, i = 1, 2, 3}.
(1) This is a linear space of dimension 12, [START_REF] Lai | Spline Functions on Triangulations[END_REF]. Indeed, each element in the space can be determined uniquely by specifying values and gradients at the 3 vertices and cross boundary derivatives at the midpoint of the edges, see Figure 1,(right).

We associate the half open edges

p i , p T ) := {(1 -t)p i + tp T : 0 ≤ t < 1}, i = 1, 2, 3,
with subtriangles of T as follows

p 1 , p T ) ∈ T 2 , p 2 , p T ) ∈ T 3 , p 3 , p T ) ∈ T 1 , (2)
and we somewhat arbitrarily assume p T ∈ T 2 .

Notation

We let N be the set of natural numbers and N 0 := N ∪ {0} the set of nonnegative integers. For a given degree d ∈ N 0 , the space of polynomials of total degree at most d will be denoted by P d . The Bernstein polynomials of degree d on T are given by

B d ijk (p) := B d ijk (β 1 , β 2 , β 3 ) := d! i!j!k! β i 1 β j 2 β k 3 , i, j, k ∈ N 0 , i+j +k = d, (3)
where p ∈ R 2 and β 1 , β 2 , β 3 , given by

p = β 1 p 1 + β 2 p 2 + β 3 p 3 , β 1 + β 2 + β 3 = 1, (4)
are the barycentric coordinates of p. The set

B d := {B d ijk : i, j, k ∈ N 0 , i + j + k = d} (5)
is a partition of unity basis for P d . The points

p d ijk := ip 1 + jp 2 + kp 3 d , i, j, k ∈ N 0 , i + j + k = d, (6)
are called the domain points of B d relative to T . In this paper, we will order the cubic Bernstein polynomials by going counterclockwise around the boundary, starting at p 1 with B 3 300 and ending with B 3 111 , see Figure 2 {B 1 , B 2 , . . . , B 10 } := {B 3 300 , B

p 3 , 2p 3 + p 1 3 , p 3 + 2p 1 3 , p T . (8)
The partial derivatives of a bivariate function

f = f (x 1 , x 2 ) are denoted ∂ 1,0 f := ∂f ∂x 1 , ∂ 0,1 f := ∂f ∂x 2 , and ∂ u f := (u 1 ∂ 1,0 + u 2 ∂ 0,1
)f is the derivative in the direction u := (u 1 , u 2 ). We denote by ∂ β j f , j = 1, 2, 3 the partial derivatives of f (β 1 , β 2 , β 3 ) with respect to the barycentric coordinates of f . The symbols S and S o are the closed and open convex hull of a set S ∈ R m . For k ≤ m, we let vol k (S) be the k dimensional volume of S and define 1 S : R m → R by

1 S (x) := 1, if x ∈ S, 0, otherwise.
By the association (2), we note that for any x ∈ T

1 T 1 (x) + 1 T 2 (x) + 1 T 3 (x) = 1 T (x). (9)
We write #K for the number of elements in a sequence K.

Bivariate simplex splines

In this section we recall some basic properties of simplex splines.

For n ∈ N, d ∈ N 0 , let m := n + d and k 1 , . . . , k m+1 ∈ R n be a sequence of points called knots. The multiplicity of a knot is the number of times it occurs in the sequence. Let σ = k 1 , . . . , k m+1 with vol m (σ) > 0 be a simplex in R m whose projection π : R m → R n onto the first n coordinates satisfies π(k i ) = k i , for i = 1, . . . , m + 1.

With [K] := [k 1 , . . . , k m+1 ], the unit integral simplex spline M [K] can be defined geometrically by

M [K] : R n → R, M [K](x) := vol m-n σ ∩ π -1 (x) vol m (σ) .
For properties of M [K] and proofs see for example [START_REF] Micchelli | On a numerically efficient method for computing multivariate B-splines[END_REF]. Here, we mention:

• If n = 1 then M [K]
is the univariate B-spline of degree d with knots K, normalized to have integral one.

• In general M [K] is a nonnegative piecewise polynomial of total degree d and support K .

• For d = 0 we have

M [K](x) := 1/vol n ( K ), x ∈ K o , 0, if x / ∈ K . (10)
• The value of M [K] on the boundary of K has to be delt with separately, see below.

• If vol n ( K ) = 0 then M [K] can be defined either as identically zero or as a distribution.

We will deal with the bivariate case n = 2, and for our purpose it is convenient to work with area normalized simplex splines, [START_REF] Lyche | Stable Simplex Spline Bases for C 3 Quintics on the Powell-Sabin 12-Split[END_REF]. They are defined by Q[K](x) = 0 for all x ∈ R 2 if vol 2 ( K ) = 0, and otherwise

Q T [K] = Q[K] := vol 2 (T ) d+2 2 M [K], (11)
where T in general is some subset of R 2 , and in our case will be the triangle T := p 1 , p 2 , p 3 . The knot sequence is [p 1 , p 2 , p 3 , p T ] taken with multiplicities. Using properties of M [K] and [START_REF] Powell | Piecewise quadratic approximation on triangles[END_REF], we obtain the following for

Q[K]:
• It is a piecewise polynomial of degree d = #K -3 with support K

• knot lines are the lines in the complete graph of K

• local smoothness: Across a knot line, Q[K] ∈ C d+1-µ
, where d is the degree and µ is the number of knots on that knot line, including multiplicities

• differentiation formula: ∂ u Q[K] = d d+3 j=1 a j Q[K \ k j ],
for any u ∈ R 2 and any a 1 , . . . , a d+3 such that j a j k j = u, j a j = 0 (A-recurrence)

• recurrence relation: Q[K](x) = d+3 j=1 b j Q[K \ k j ](x), for any x ∈ R 2 and any b 1 , . . . , b d+3 such that j b j k j = x, j b j = 1 (B-recurrence) • knot insertion formula: Q[K] = d+3 j=1 c j Q[K ∪ y \ k j ],
for any y ∈ R 2 and any c 1 , . . . , c d+3 such that j c j k j = y, j c j = 1 (C-recurrence)

• degree zero: From ( 10) and ( 11) we obtain for d = 0

Q[K](x) := vol 2 (T )/vol 2 ( K ), x ∈ K o , 0, if x / ∈ K . ( 12
)
2 A simplex spline basis for the Clough-Tocher split

In this section we determine and study a basis of C 1 cubic simplex splines on the Clough-Tocher split on a triangle. For fixed x ∈ T we use the simplified notation

i j k := Q[p [i] 1 , p [j] 2 , p [k] 3 , p [l] T ](x), i, j, k, l ∈ N 0 , i + j + k + l ≥ 3,
where the notation p

[n] m denotes that p m is repeated n times. When one of the integers i, j, k, l is zero we have Lemma 1 For i, j, k, l ∈ N 0 , i+j +k+l = d ≥ 0 and x ∈ T with barycentric coordinates β 1 , β 2 , β 3 we have

i = 0, j+1 k+1 + = d! j!k!l! (β 2 -β 1 ) j (β 3 -β 1 ) k (3β 1 ) l 1 1 1 , j = 0, i+1 k+1 + = d! i!k!l! (β 1 -β 2 ) i (β 3 -β 2 ) k (3β 2 ) l 1 1 1 , k = 0, i+1 j+1
+ = d! i!j!l! (β 1 -β 3 ) i (β 2 -β 3 ) j (3β 3 ) l 1 1 1 , l = 0, i+1 j+1 k+1 = d! i!j!k! β i 1 β j 2 β k 3 1 1 1 = B d ijk (x), ( 13
)
where the constant simplex splines are given by

1 1 1 = 3 1 T 1 (x), 1 1 1 = 3 1 T 2 (x), 1 1 1 = 3 1 T 3 (x), 1 1 1 = 1 T (x). ( 14
)
Proof: Suppose i = 0. The first equation in (13) holds for d = 0. Suppose it holds for d -1 and let j + k + l = d. Let β 023 j , j = 0, 2, 3 be the barycentric coordinates of x with respect to T 1 = p 0 , p 2 , p 3 , where p 0 := p T . By the B-recurrence

j+1 k+1 + = β 023 2 j k+1 + + β 023 3 j+1 k + + β 023 0 j+1 k+1 .
It is easily shown that

β 023 2 = β 2 -β 1 , β 023 3 = β 3 -β 1 , β 023 0 = 3β 1 .
Therefore, by the induction hypothesis

j+1 k+1 + = (d -1)! j!k!l! (j + k + l)(β 023 2 ) j (β 023 3 ) k (β 023 0 ) l 1 1 1 Since j + k + l = d we obtain the first equation in (13).
The next two equations in (13) follow similarly using

β 031 1 = β 1 -β 2 , β 031 3 = β 3 -β 2 , β 031 0 = 3β 2 , β 012 1 = β 1 -β 3 , β 012 2 = β 2 -β 3 , β 012 0 = 3β 3 .
Using the B-recurrence repeatedly, we obtain the first equality for l = 0. The values of the constant simplex splines are a consequence of [START_REF] Speleers | A normalized basis for quintic Powell Sabin splines[END_REF].

Remark 2 For i = 0 we note that the expression

d! j!k!l! (β 2 -β 1 ) j (β 3 - β 1 ) k (3β 1 ) l in (13) is a Bernstein polynomial on T 1 . Similar remarks hold for j, k = 0. The set C1 := i j k ∈ S 1 3 ( ) : i j k = 0 (15)
of all nonzero simplex splines that can be used in a basis for S 1 3 ( ) contains precisely the following 13 simplex splines.

Lemma 3 We have

C1 = i j k : i, j, k ∈ N, i + j + k = 6 2 2 1 1 , 1 2 2 1 , 2 1 2 1
.

Proof: For l = 0 it follows from Lemma 1 that i j k ∈ S 1 3 ( ) for all i + j + k = 6. Consider next l = 1. By the local smoothness property, C 1 smoothness implies that each of i, j, k can be at most 2. But then

2 2 1 1 , 1 2 2 1 , 2 1 2 1
are the only possibilities. Now if l = 2 then i + j + k = 4 implies that one of i, j, k must be at least 2 and we cannot have C 1 smoothness. Similarly l > 2 is not feasible. Recall that S 1 3 ( ) is a linear space of dimension 12, [START_REF] Clough | Finite element stiffness matrices for analysis of plate bending[END_REF]. Thus, in order to obtain a possible basis for this space, we need to choose 12 of the 13 elements in C1. Since C1 contains the 10 cubic Bernstein polynomials we have to include at least two of

2 2 1 1 , 1 2 2 1 , 2 1 2 1
. We also want a symmetric basis and therefore, we have to include all of them. But then one of the Bernstein polynomials has to be excluded. To see which one to exclude, we insert the

knot p 3 = -p 1 -p 2 + 3p T into 2 2 1 1
and use the C-recurrence to obtain

2 2 1 1 = - 1 2 2 1 - 2 1 2 1 + 3 2 2 2
, or by (13)

2 2 1 1 + 1 2 2 1 + 2 1 2 1 = 3B 3 111 (x). ( 16
)
Thus, in order to have symmetry and hopefully obtain 12 linearly independent functions, we see that B 3 111 is the one that should be excluded. We obtain the following simplex spline basis for S 1 3 ( ).

Theorem 4 (CTS-basis) The 12 simplex splines S 1 , . . . , S 12 , where

S j (x) := B j (x)
, where B j is given by (7) j = 1, . . . , 9,

S 10 (x) := 1 3 2 2 1 1 = (B 3 210 -B 3 300 )1 T 1 + (B 3 120 -B 3 030 )1 T 2 + (B 3 111 -B 3 102 -B 3 012 + 2B 3 003 )1 T 3 S 11 (x) := 1 3 1 2 2 1 = (B 3 111 -B 3 210 -B 3 201 + 2B 3 300 )1 T 1 + (B 3 021 -B 3 030 )1 T 2 + (B 3 012 -B 3 003 )1 T 3 S 12 (x) := 1 3 2 1 2 1 = (B 3 201 -B 3 300 )1 T 1 + (B 3 111 -B 3 120 -B 3 021 + 2B 3 030 )1 T 2 + (B 3 102 -B 3 003 )1 T 3 .
(17) form a partition of unity basis for the space S 1 3 ( ) given by (1). This basis, which we call the CTS-basis, is the only symmetric simplex spline basis for S 1 3 ( ). On the boundary of T the functions S 10 , S 11 , S 12 have the value zero, while the elements of {S 1 , S 2 , . . . , S 9 } reduce to zero, or to univariate Bernstein polynomials.

Proof: By Lemma 1, it follows that the Bernstein polynomials B 1 , . . . , B 9 are cubic simplex splines, and the previous discussion implies that the functions in (17), apart from scaling, are the only candidates for a symmetric simplex spline basis for S 1 3 ( ).

We can find the explicit form of (see definitions at the end of Section 1) . Consider the C-recurrence. Insert-ing p 1 twice and using p 1 = -p 2p 3 + 3p T and (13) we find

2 2 1 1 = - 3 1 1 1 - 3 2 1 + 3 3 2 1 = 4 1 1 + 4 1 1 -3 4 1 1 - 3 2 1 + 3 3 2 1 = (β 1 -β 2 ) 3 1 1 1 + (β 1 -β 3 ) 3 1 1 1 -3β 3 1 1 1 1 -3(β 1 -β 3 ) 2 (β 2 -β 3 ) 1 1 1 + 9β 2 1 β 2 1 1 1 = (β 1 -β 2 ) 3 1 1 1 + [(β 1 -β 3 ) 3 -3(β 1 -β 3 ) 2 (β 2 -β 3 )] 1 1 1 + 3β 2 1 (3β 2 -β 1 ) 1 1 1 .
(18)

Using ( 9) and Lemma 1, we can write 3

1 1 1 = 1 1 1 + 1 1 1 + 1 1 1 , so that 2 2 1 1 = [(β 1 -β 2 ) 3 + β 2 1 (3β 2 -β 1 )] 1 1 1 + β 2 1 (3β 2 -β 1 ) 1 1 1 + [(β 1 -β 3 ) 2 (β 1 -3β 2 + 2β 3 ) + β 2 1 (3β 2 -β 1 )] 1 1 1 = (3β 2 1 β 2 -β 3 1 ) 1 1 1 + (3β 1 β 2 2 -β 3 2 ) 1 1 1 + (6β 1 β 2 β 3 -3β 1 β 2 3 -3β 2 β 2 3 + 2β 3 3 ) 1 1 1 . ( 19
)
By symmetry we obtain

1 2 2 1 = (6β 1 β 2 β 3 -3β 2 1 β 2 -3β 2 1 β 3 + 2β 3 1 ) 1 1 1 + (3β 2 2 β 3 -β 3 2 ) 1 1 1 + (3β 2 β 2 3 -β 3 3 ) 1 1 1 , 2 1 2 1 = (3β 2 1 β 3 -β 3 1 ) 1 1 1 + (3β 1 β 2 3 -β 3 3 ) 1 1 1 + (6β 1 β 2 β 3 -3β 1 β 2 2 -3β 2 2 β 3 + 2β 3 2 ) 1 1 1 . ( 20
)
The formulas for S 10 , S 11 and S 12 in (17) now follows from ( 19) and (20) using ( 3) and ( 14).

By the partition of unity for Bernstein polynomials we find

12 j=1 S j (x) = i+j+k=3 B 3 ijk (x) = 1, x ∈ T .
It is well known that B 3 ijk reduces to univariate Bernstein polynomials or zero on the boundary of T .

Clearly S j ∈ C(R 2 ), j = 10, 11, 12, since no edge contains more than 4 knots. This follows from general properties of simplex splines. By the local support property they must therefore be zero on the boundary. It also follows that S j ∈ C 1 (T ), j = 10, 11, 12, since no interior knot line contains more than 3 knots.

It remains to show that the 12 functions S j , j = 1, . . . , 12 are linearly independent on T . Suppose that 12 j=1 c j S j (x) = 0 for all x ∈ T and let (β 1 , β 2 , β 3 ) be the barycentric coordinates of x. On the edge p 1 , p 2 , where β 3 = 0, the functions S j , j = 5, . . . 12 vanish, and thus

12 j=1 c j S j (x) = c 1 B 3 300 (x) + c 2 B 3 210 (x) + c 3 B 3 120 (x) + c 4 B 3 030 (x) = 0.
On p 1 , p 2 this is a linear combination of linearly independent univariate Bernstein polynomials and we conclude that c 1 = c 2 = c 3 = c 4 = 0. Similarly c j = 0 for j = 5, . . . , 9. It remains to show that S 10 , S 11 and S 12 are linearly independent on T . For x ∈ T o 3 and β 3 = 0 we find

∂S 10 ∂β 3 | β 3 =0 = 6β 1 β 2 = 0, ∂S j ∂β 3 | β 3 =0 = 0, j = 11, 12.
We deduce that c 10 = 0 and similarly c 11 = c 12 = 0 which concludes the proof.

In Figure 3 we show graphs of the functions S 10 , S 11 , S 12 .

Two Marsden identities and representation of polynomials

We give both a barycentric and a Cartesian Marsden-like identity.

Theorem 5 (Barycentric Marsden-like identity) For u := (u 1 , u 2 , u 3 ), The CTS-basis functions S 10 , S 11 , S 12 on the triangle (0, 0), (1, 0), (0, 1) .

β := (β 1 , β 2 , β 3 ) ∈ R 3 with β i ≥ 0, i = 1, 2, 3 and β 1 + β 2 + β 3 = 1 we have (β T u) 3 = u 3 1 S 1 (β) + u 2 1 u 2 S 2 (β) + u 1 u 2 2 S 3 (β) + u 3 2 S 4 (β) + u 2 2 u 3 S 5 (β) + u 2 u 2 3 S 6 (β) + u 3 3 S 7 (β) + u 1 u 2 3 S 8 (β) + u 2 1 u 3 S 9 (β) + u 1 u 2 u 3 S 10 (β) + S 11 (β) + S 12 (β) .
Proof: By the multinomial expansion we obtain

(β 1 u 1 + β 2 u 2 + β 3 u 3 ) 3 = i+j+k=3 3! i!j!k! (β 1 u 1 ) i (β 2 u 2 ) j (β 3 u 3 ) k = i+j+k=3 u i 1 u j 2 u k 3 B 3 ijk (β).
Using B 3 111 = S 10 + S 11 + S 12 and the ordering in Theorem 4 we obtain (21).

Corollary 6 For l, m, n ∈ N 0 with l + m + n ≤ 3 we have an explicit representation for lower degree Bernstein polynomials in terms of the CTSbasis (17).

B l+m+n lmn = 3 l + m + n -1 3 l 0 m 0 n S 1 + 2 l 1 m 0 n S 2 + 1 l 2 m 0 n S 3 + 0 l 3 m 0 n S 4 + 0 l 2 m 1 n S 5 + 0 l 1 m 2 n S 6 + 0 l 0 m 3 n S 7 + 1 l 0 m 2 n S 8 + 2 l 0 m 1 n S 9 + 1 l 1 m 1 n S 10 + S 11 + S 12 , (22)
where 0 0 := 1 and r s := 0 if s > r.

Proof: Differentiating, for any d ∈ N 0 , (β 1 u 1 +β 2 u 2 +β 3 u 3 ) d a total of l, m, n times with respect to u 1 , u 2 , u 3 , respectively, and setting

u 1 = u 2 = u 3 = 1 we find d! (d -l -m -n)! β l 1 β m 2 β n 3 = i+j+k=d i(i -1) . . . (i -l + 1)j . . . (j -m + 1)k . . . (k -n + 1)B d ijk ,
and by a rescaling

B l+m+n lmn = d l + m + n -1 i+j+k=d i l j m k n B d ijk , l + m + n ≤ d.
(23) Using ( 17) with d = 3, we obtain (22).

As an example, we find

B 1 100 = 1 3 3S 1 + 2S 2 + S 3 + S 8 + 2S 9 + S 10 + S 11 + S 12 .
Theorem 7 (Cartesian Marsden-like identity) We have

(1 + x T v) 3 = 12 j=1 ψ j (v)S j (x), x ∈ T , v ∈ R 2 , ( 24
)
where the dual polynomials in Cartesian form are given by

ψ j (v) := 3 l=1 (1 + d T j,l v), j = 1, . . . , 12, v ∈ R 2 . ( 25
)
Here the dual points d j := [d j,1 , d j,2 , d j,3 ], are given as follows.

                     d 1 d 2 d 3 d 4 d 5 d 6 d 7 d 8 d 9 d 10 d 11 d 12                      :=                      p 1 p 1 p 1 p 1 p 1 p 2 p 1 p 2 p 2 p 2 p 2 p 2 p 2 p 2 p 3 p 2 p 3 p 3 p 3 p 3 p 3 p 1 p 3 p 3 p 1 p 1 p 3 p 1 p 2 p 3 p 1 p 2 p 3 p 1 p 2 p 3                      . ( 26
)
The domain points p * j in (8) are the coefficients of x in terms of the CTSbasis

x = 12 j=1 p * j S j (x), ( 27
)
where p * 10 = p * 11 = p * 12 = p T .

Proof: We apply (21) with β 1 , β 2 , β 3 the barycentric coordinates of x and

u i = 1 + p T i v, i = 1, 2, 3.
Then

β 1 u 1 + β 2 u 2 + β 3 u 3 = β 1 + β 2 + β 3 + β 1 p T 1 v + β 2 p T 2 v + β 3 p T 3 v = 1 + x T v.
and ( 24), ( 25), (26) follow from (21). Taking partial derivatives in (24) with respect to v,

∂ v 1 , ∂ v 2 (1 + x T v) 3 = 3x(1 + x T v) 2 = 12 j=1 ∂ v 1 , ∂ v 2 ψ j (v)S j (x),
where

∂ v 1 , ∂ v 2 ψ j (v) := d j,1 (1+d T j,2 v)(1+d T j,3 v)+d j,2 (1+d T j,1 v)(1+d T j,3 v)+ d j,3 (1 + d T j,1 v)(1 + d T j,2 v).
Setting v = 0 we obtain (27). Note that the domain point p T for B 3 111 has become a triple domain point for the CTS-basis.

Following the proof of (27) we can give explicit representations of all the monomials x r y s spanning P 3 . We do not give details here.

Three quasi-interpolants

We consider three quasi-interpolants on S 1 3 ( ). They all use functionals based on point evaluations and the third one will be used to estimate the L ∞ condition number of the CTS-basis.

To start, we consider the following polynomial interpolation problem on T . Find g ∈ P 3 such that g(p * i ) = f i , where f := [f 1 , . . . , f 10 ] T is a vector of given real numbers and the p * i given by ( 8) are the domain points for the cubic Bernstein basis.

Using the ordering (7), we write g in the form 10 j=1 c j B j and obtain the linear system

A = A 1 0 A 2 A 3 , (28)
and if A 1 and A 3 are nonsingular then

A -1 = A -1 1 0 -A -1 3 A 2 A -1 1 A -1 3 = B 1 0 B 2 B 3 . ( 29
)
Using the barycentric form of the domain points in ( 8) we find

A 2 = [1, 3, 3, 1, 3, 3, 1, 3, 3]/27, A 3 = B 3 111 ( 1 3 , 1 3 , 1 3 ) = 2 9 , A 1 := 1 27              
27 0 0 0 0 0 0 0 0 8 12 6 1 0 0 0 0 0 1 6 12 8 0 0 0 0 0 0 0 0 27 0 0 0 0 0 0 0 0 8 12 6 1 0 0 0 0 0 1 6 12 8 0 0 0 0 0 0 0 0 27 0 0 1 0 0 0 0 0 8 12 6 8 0 0 0 0 0 1 6 12

              ∈ R 9×9 (30)
and

B 1 := A -1 1 = 1 6               6 0 0 0 0 0 0 0 0 -5 18 -9 2 0 0 0 0 0 2 -9 18 -5 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 -5 18 -9 2 0 0 0 0 0 2 -9 18 -5 0 0 0 0 0 0 0 0 6 0 0 2 0 0 0 0 0 -5 18 -9 -5 0 0 0 0 0 2 -9 18               , B 3 = [ 9 2 ], B 2 := -B 3 A 2 B 1 = 1 12
λ P i (f )B i , λ P i (f ) := 10 j=1 α i,j f (p * j ), ( 32
)
where the matrix α := A -1 has elements α i,j in row i and column j, i, j = 1, . . . , 10. We have

λ P i (B j ) = 10 k=1 α i,k B j (p * k ) = 10 k=1 α i,k a k,j = δ i,j , i, j = 1, . . . , 10.
It follows that QI P (g) = g for all g ∈ P 3 . Since B j = S j , j = 1, . . . , 9 and B 10 = B 3 111 = S 10 + S 11 + S 12 the quasi-interpolant

QI P : C(T ) → S 1 3 ( ), QI P (f ) := 12 i=1 λ P i (f )S i , λ P 11 = λ P 12 = λ P 10 , ( 33
)
where λ P i (f ) is given by (32), i = 1, . . . , 10, reproduces P 3 . Moreover, since for any f ∈ C(T ) and x ∈ T

|QI P (f )(x)| ≤ max 1≤i≤12 |λ P i (f )| 12 i=1 S i (x) = max 1≤i≤10 |λ P i (f )|,
we obtain

QI P (f ) L∞(T ) ≤ α ∞ f L∞(T ) = 10 f L∞(T ) ,
independently of the geometry of T . Using the construction in [START_REF] Lyche | Stable Simplex Spline Bases for C 3 Quintics on the Powell-Sabin 12-Split[END_REF], we can derive another quasi-interpolant which also reproduces P 3 . It uses more points, but has a slightly smaller norm. Consider the map P : C(T ) → S 1 3 (T ) defined by P (f ) = 12 ℓ=1 M ℓ (f )S ℓ , where

M ℓ (f ) := 1 6 f (d ℓ,1 ) + f (d ℓ,2 ) + f (d ℓ,3 ) + 9 2 f (p * ℓ ) - 4 3 f d ℓ,1 + d ℓ,2 2 + f d ℓ,1 + d ℓ,3 2 + f d ℓ,2 + d ℓ,3 2 .
Here the d ℓ,m are the dual points given by (26) and the p * ℓ are the domain points given by (27). Note that this is an affine combination of function values of f .

We have tested the convergence of the quasi-interpolant, sampling data from the function f (x, y) = e 2x+y + 5x + 7y on the triangle A = [0, 0],

B = h * [1, 0], C = h * [0.2, 1.2]
for h ∈ {0.05, 0.04, 0.03, 0.02, 0.01}. The following array indicates that the error: f -P (f ) L∞(T ) , is O(h 4 ). h 0.05 0.04 0.03 0.02 0.01 error/h 4 0.0550 0.0547 0.0543 0.0540 0.0537 Using a standard argument the following Proposition shows that the error is indeed O(h 4 ) for sufficiently smooth functions.

Proposition 8

The operator P is a quasi-interpolant that reproduces P 3 . For any f ∈ C(T )

P (f ) L∞(T ) ≤ 9 f L∞(T ) , (34)
independently of the geometry of T . Moreover,

f -P (f ) L∞(T ) ≤ 10 inf g∈P 3 f -g L∞(T ) . (35)
Proof: Since d 10 = d 11 = d 12 and B 3 111 = S 10 + S 11 + S 12 , B 3 ijk = S ℓ for (i, j, k) = (1, 1, 1) and some ℓ, we obtain

P (f ) = i+j+k=3 Mijk (f )B 3 ijk
where Mijk = M ℓ for (i, j, k) = (1, 1, 1) and corresponding ℓ and M111 = 3M 10 .

To prove that P reproduces polynomials up to degree 3, i.e., P (B 3 ijk ) = B 3 ijk , whenever i + j + k = 3, it is sufficient to prove the result for B

p 3 , p 2 + p 3 2 , it is easy to compute that M300 (B 3 300 ) = 1, M300 (B 3 ijk ) = 0 for (i, j, k) = (3, 0, 0), M210 (B 3 210 ) = 1, M210 (B 3 ijk ) = 0 for (i, j, k) = (2, 1, 0), M111 (B 3 111 ) = 1, M111 (B 3 ijk ) = 0 for (i, j, k) = (1, 1 , 1).
Therefore, by a standard argument, P is a quasi-interpolant that reproduces P 3 . Since the sum of the absolute values of the coefficients defining M ℓ (f ) is equal to 9, another standard argument shows (34) and (35).

The operators QI P and P do not reproduce the whole spline space S 1 3 ( ). Indeed, since λ P 10 (B 10 ) = M 10 (B 10 ) = 1, we have λ P 10 (S j ) = M 10 (S j ) = 1 3 , j = 10, 11, 12.

To give un upper bound for the condition number of the CTS-basis we need a quasi-interpolant which reproduces the whole spline space. We again use the inverse of the coefficient matrix of an interpolation problem to construct such an operator. We need 12 interpolation points and a natural choice is to use the first 9 cubic Bernstein domain points p * j , j = 1, . . . , 9 and split the barycenter p * 10 = p T into three points. After some experimentation we redefine p * 10 and choose p * 10 := (3, 3, 1)/7, p * 11 := (3, 1, 3)/7 and p * 12 := (1, 3, 3)/7. The problem is to find s = 12 j=1 c j S j such that s(p * i ) = f i , i = 1, . . . , 12. The coefficient matrix for this problem has again the block tridiagonal form (28), where A 1 ∈ R 9×9 and B 1 := A -1 1 are given by ( 30) and (31) as before. Moreover, using the formulas in Theorem 4 we find

A 3 = [S j (p * i )]
α S := A -1 = B 1 0 B 2 B 3 ,
where It follows that the quasi-interpolant QI given by

B 2 = -B 3 A 2 B 1 =    
QI : C(T ) → S 1 3 ( ), QI(f ) := 12 i=1 λ S i (f )S i , λ S i (f ) = 12 j=1 α S i,j f (p * j ), (37)
is a projector onto the spline space S 1 3 ( ). In particular

s := 12 i=1 c i S i =⇒ c i = λ S i (s), i = 1, . . . , 12. (38)
The quasi-interpolant (37) can be used to show the L ∞ stability of the CTS-basis. For this we prove that the condition number is independent of the geometry of the triangle.

We define the ∞-norm condition number of the CTS-basis on T by

κ ∞ (T ) := max c =0 b T c L∞(T ) c ∞ max c =0 c ∞ b T c L∞(T )
,

where b T c := 12 j=1 c j S j ∈ S 1 3 ( ).

) |c i | = |λ S i (b T c)| ≤ α S ∞ b T c L∞(T ) . Therefore, c ∞ b T c L∞(T ) ≤ α S ∞ = 27 - 32 405 ,
and the upper bound κ ∞ < 27 follows.

5 C 0 and C 1 -continuity

In the following, we derive conditions to ensure C 0 and C 1 continuity through an edge between two triangles. The conditions are very similar to the classical conditions for continuity of Bernstein polynomials.

Theorem 10

Let s 1 = 12 j=1 c j S j and s 2 = 12 j=1 d j Sj be defined on the triangle T := p 1 , p 2 , p 3 and T := p 1 , p 2 , p3 , respectively, see Figure 4.

The function s

= s 1 on T s 2 on T is continuous on T ∪ T if d 1 = c 1 , d 2 = c 2 , d 3 = c 3 , d 4 = c 4 . (39)
Moreover, s ∈ C 1 (T ∪ T ) if in addition to (39) we have

d 5 = γ 1 c 3 + γ 2 c 4 + γ 3 c 5 , d 9 = γ 1 c 1 + γ 2 c 2 + γ 3 c 9 , d 10 = γ 1 c 2 + γ 2 c 3 + γ 3 c 10 .
(40) where γ 1 , γ 2 , γ 3 are the barycentric coordinates of p3 with respect to T . Suppose next (39) holds and s ∈ C 1 (T ∪ T ). By the continuity property we see that S j , j = 6, 7, 8, 11, 12 are zero and have zero cross boundary derivatives on p 1 , p 2 since they have at most 3 knots on that edge. We take derivatives in the direction u := p3p 1 using the A-recurrence (defined at the end of Section 1) with a := (γ 1 -1, γ 2 , γ 3 , 0) for s 1 and a := (-1, 0, 1, 0)

The Hermite basis

The classical Hermite interpolation problem on the Clough-Tocher split is to interpolate values and gradients at vertices and normal derivatives at the midpoint of edges, see Figure 1. These interpolation conditions can be described by the linear functionals

ρ(f ) = [ρ 1 (f ), . . . , ρ 12 (f )] T := [f (p 1 ), ∂ 1,0 f (p 1 ), ∂ 0,1 f (p 1 ), f (p 2 ), ∂ 1,0 f (p 2 ), ∂ 0,1 f (p 2 ), f (p 3 ), ∂ 1,0 f (p 3 ), ∂ 0,1 f (p 3 ), ∂ n 1 f (p 5 ), ∂ n 2 f (p 6 ), ∂ n 3 f (p 4 )] T ,
where p 4 , p 5 , p 6 , are the midpoints on the edges p 1 , p 2 , p 2 , p 3 , p 3 , p 1 , respectively, and ∂ n j f is the derivative in the direction of the unit normal to that edge in the direction towards p j . We let p j = (x j , y j ) be the coordinates of each point. The coefficient vector c := [c 1 , . . . , c 12 ] T of the interpolant g := 12 j=1 c j S j is solution of the linear system Ac = ρ(f ), where A ∈ R 12×12 with a i,j := ρ i (S j ).

Let H 1 , . . . , H 12 be the Hermite basis for S 1 3 ( ) defined by ρ i (H j ) = δ i,j . The matrix A transforms the Hermite basis to the CTS-basis. Since a basis transformation matrix is always nonsingular, we have

[S 1 , . . . , S 12 ] = [H 1 , . . . , H 12 ]A, [H 1 , . . . , H 12 ] = [S 1 , . . . , S 12 ]A -1 . (44)
To find the elements ρ i (S j ) of A we define for i, j, k = 1, 2, 3

ν ij := p ij 2 , p ij := p i -p j , x ij := x i -x j , y ij := y i -y j , ν ijk := p T i,j p j,k ν ij , for i = j, δ := 1 1 1 x 1 x 2 x 3 y 1 y 2 y 3 . (45)
We note that ν ijk is the length of the projection of p j,k in the direction of p i,j and that δ is twice the signed area of T . By the definition of the unit normals and the chain rule for j = 1, . . . , 12 we find ∂ 1,0 S j = (y 23 ∂ β 1 S j + y 31 ∂ β 2 S j + y 12 ∂ β 3 S j )/δ, ∂ 0,1 S j = (x 32 ∂ β 1 S j + x 13 ∂ β 2 S j + x 21 ∂ β 3 S j )/δ, ∂ n 1 S j = (y 23 ∂ 1,0 S j + x 32 ∂ 0,1 S j )/ν 32 = (ν 32 ∂ β 1 S j + ν 231 ∂ β 2 S j + ν 321 ∂ β 3 S j )/δ, ∂ n 2 S j = (y 31 ∂ 1,0 S j + x 13 ∂ 0,1 S j )/ν 31 = (ν 132 ∂ β 1 S j + ν 31 ∂ β 2 S j + ν 312 ∂ β 3 S j )/δ, ∂ n 3 S j = (y 12 ∂ 1,0 S j + x 21 ∂ 0,1 S j )/ν 21 = (ν 123 ∂ β 1 S j + ν 213 ∂ β 2 S j + ν 21 ∂ β 3 S j )/δ. This leads to

A := A 1 0 A 2 A 3 , with A 1 ∈ R 9×9 , A 2 ∈ R 3×9 , A 3 ∈ R 3×3 ,
where We find

A 1 := 3 δ              
A -1 := B 1 0 B 2 B 3 = [b i,j ] 12 i,j=1 ,
where

B 1 := A -1 1 = 1 3              
3 0 0 0 0 0 0 0 0 3 x 21 y 21 0 0 0 0 0 0 0 0 0 3 x 12 y 12 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 3 x 32 y 32 0 0 0 0 0 0 0 0 0 3 x 23 y 23 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 3 x 13 y 13 3 x 31 y 31 0 0 0 0 0 0

              ∈ R 9×9 ,

•

  a differentiation formula • a stable recurrence relation • a knot insertion formula • they constitute a nonnegative partition of unity • simple explicit dual functionals • L ∞ stability • simple conditions for C 1 and C 2 joins to neighboring triangles • well conditioned collocation matrices for Lagrange and Hermite interpolation using certain sites.

Figure 1 :

 1 Figure 1: The PS12-split (left) and the CT-split (right). The C 1 quadratics on PS-12 and C 1 cubics on CT have the same degrees of freedom as indicated.

Figure 2 :

 2 Figure 2: The cubic Bernstein basis (left) and the CTS-basis (right), where B 3 111 is replaced by S 10 , S 11 , S 12 .

  the B-or C-recurrence

  Figure 3:

10 j=1c

 10 j B j (p * i ) = f i , i = 1, . . . , 10, or in matrix form Ac = f for the unknown coefficient vector c := [c 1 , . . . , c 10 ] T . Since B 10 (p * i ) = B 3 111 (p * i ) = 0 for i = 1, . . . , 9 the coefficient matrix A is block triangular

4 - 9 -9 4 - 9 -9 4 - 9 - 9 .

 4949499

Theorem 9

 9 For any triangle T we have κ ∞ (T ) < 27. Proof: Since the S j form a nonnegative partition of unity it follows that max c =0 b T c L∞(T ) / c ∞ = 1. If s = 12 j=1 c j S j = b T c then by (38

  c 10 ,c 11 ,c 12

Figure 4 :

 4 Figure 4: C 1 -continuity and components

Figure 5 :

 5 Figure 5: C 1 smoothness

Figure 6 :

 6 Figure 6: The Hermite basis functions H 1 , H 2 , H 3 , H 10 on the unit triangle.

  ν 32 , ν 231 , ν 231 -ν 32 , ν 321 -ν 32 , ν 321 , ν 32 , 0 , A 2 (2) := 3 4δ ν 132 , ν 31 , 0, 0, 0, ν 31 , ν 312 , ν 312 -ν 31 ν 132 -ν 31 , A 2 (3) := 3 4δ ν 123 , ν 123 -ν 21 , ν 213 -ν 21 , ν 213 , ν 21 , 0, 0, 0, ν 21 ,

B 3 :

 3 3×3 , and the rows ofB 2 = -B 3 A 2 B 1 ∈ R 3×9 are given by B 2 (1) :=

Figure 7 :

 7 Figure 7: The triangulation and the C 1 surface

Figure 8 : A C 1

 81 Figure 8: A C 1 Hermite interpolating surface on the triangulation

  1 6ν 21 -6ν 123 , x 12 ν 123 + ν 21 x 23 , y 12 ν 123 + ν 21 y 23 , -6ν 213 , x 21 ν 213 + ν 21 x 13 , y 21 ν 213 + ν 21 y 13 , 0, 0, 0 , -6ν 231 , x 23 ν 231 + ν 32 x 31 , y 23 ν 231 + ν 32 y 31 , -6ν 321 , x 23 ν 231 + ν 32 x 21 + ν 32 x 23 , y 23 ν 231 + ν 32 y 21 + ν 32 y 23 , 6ν 132 , x 13 ν 132 + ν 31 x 32 , y 13 ν 132 + ν 31 y 32 , 0, 0, 0, -6ν 312 , x 31 ν 312 + ν 31 x 12 , y 31 ν 312 + ν 31 y 12 .

		P 6
	P 4	P 5
		P 3
	B 2 (2) := 0, 0, 0, B 2 (3) := 1 6ν 32 1 6ν 31 -P 1 P 2

for s 2 . We find with x ∈ p 1 , p 2

The last equality follows from (13) since β 3 = 0 on p 1 , p 2 so that = 3B 2 110 (x). Consider next Sj . By the same argument as for S j , we see that Sj , j = 6, 7, 8, 11, 12 are zero and have zero cross boundary derivatives on p 1 , p 2 . We find for x ∈ p 1 , p 2

We note that on p 1 , p 2 , the polynomials B 2 101 , B2 101 , B 2 011 , B2 011 vanish and

As an example, on the unit triangle (p 1 , p 2 , p 3 ) = ((0, 0), (1, 0), (0, 1)) we find

Some of the Hermite basis functions are shown in Figure 6.

We have also tested the convergence of the Hermite interpolant, sampling again data from the function f (x, y) = e 2x+y + 5x + 7y on the triangle

Examples

Several examples have been considered for scattered data on the CT-split, see for example [START_REF] Farin | A modified Clough-Tocher interpolant[END_REF][START_REF] Mann | Cubic precision Clough-Tocher interpolation[END_REF]. Here, we consider a triangulation with vertices p 1 = (0, 0), p 2 = (1, 0), p 3 = (3/2, 1/2), p 4 = (-1/2, 1), p 5 = (1/4, 3/4), p 6 = (3/2, 3/2), p 7 = (1/2, 2) and triangles T 1 := p 1 , p 2 , p 5 , T 2 := p 2 , p 3 , p 5 , T 3 := p 4 , p 1 , p 5 , T 4 := p 3 , p 6 , p 5 , T 5 := p 6 , p 4 , p 5 , T 6 := p 4 , p 6 , p 7 ,. We divide each of the 6 triangles into 3 subtriangles using the Clough-Tocher split. We then obtain a space of C 1 piecewise polynomials of dimension 3V + E = 3 × 7 + 12 = 33, where V is the number of vertices and E the number of edges in the triangulation. We can represent a function s in this space by either using the Hermite basis or using CTS-splines on each of the triangles and enforcing the C 1 continuity conditions. The function s on T 1 depends on 12 components, while the C 1 -continuity through the edges gives only 5 free components for T 2 ,T 3 and T 4 . Closing the 1-cell at p 5 gives only one free component for T 5 and 5 free components for T 6 , Figure 7 left.

In the following graph, Figure 7, right, once the 12 first components on T 1 were chosen, the other free ones are set to zero. Then, in Figure 8, we have plotted the Hermite interpolant of the function f (x, y) = e 2x+y +5x+7y and gradients using the CTS-splines.