Enrico Formenti
email: enrico.formenti@univ-cotedazur.fr

Jean-Charles Régin
email: jean-charles.regin@univ-cotedazur.fr

Sara Riva
email: sara.riva@univ-cotedazur.fr

MDDs boost equation solving on discrete dynamical systems

Keywords: Multi-valued decision diagrams, Discrete Dynamical Systems, Graphs semiring

. The boost in performance of our algorithms comes essentially from a clever usage of Multi-valued decision diagrams. These results are an important step forward in the analysis of complex dynamics graphs as those appearing, for instance, in biological regulatory networks or in systems biology.

Introduction

Multi-valued Decision Diagrams (MDD) are a generalization of Binary Decision Diagrams (BDD) [START_REF] Akers | Binary decision diagrams[END_REF][START_REF] Bergman | Decision diagrams for optimization[END_REF] used to obtain efficient representations of functions (with finite domains) or (finite) sets of tuples. An MDD is a Directed Acyclic Graph (DAG) created from a finite set of variables with specific (finite) domains. Associating each variable with a level of the structure, the MDD represents a set of feasible assignments as a path from the root to the final node. A crucial aspect of MDDs is the exponential compression power of the reduction operation and the fact that many classical operations (intersection, union, etc.) can be performed without decompression. In the last years, MDDs have been applied in many disparate research domains proving the potential of this structure. MDDs are used, for instance, to improve random forest algorithms replacing the classic binary decision trees [START_REF] Nakahara | A random forest using a multivalued decision diagram on an fpga[END_REF], to represent and analyze automotive product data of valid/invalid product configurations [START_REF] Berndt | Multi-valued decision diagrams for the verification of consistency in automotive product data[END_REF], and to perform trust analysis in social networks [START_REF] Zhang | Multivalued decision diagrams-based trust level analysis for social networks[END_REF]. There are also applications related to mathematical models like Multi-State Systems (MSS) [START_REF] Zaitseva | A multi-valued decision diagram for estimation of multi-state system[END_REF]. In this case, MDDs represent the MSS structure in terms of Multi-Valued Logic to compute some measures. MDDs find applications also in the analysis of the discrete dynamical systems. Indeed, in [START_REF] Naldi | Decision diagrams for the representation and analysis of logical models of genetic networks[END_REF], MDD represents logic rules to analyze some properties and dynamics aspects in the case of regulatory networks (for example, to perform a stable states identification).

In this work, we propose to apply MDDs to equations over Discrete Dynamical Systems (DDS). DDS are a model to represent complex phenomena which evolve in (discrete) time coming from different domains such as genetic regulatory networks, boolean automata networks, population dynamics, and many others. DDS consist in a finite set of states and a next-state function. In the finite case, DDS correspond to a particular class of graphs, called dynamic graphs (i.e. graphs with outgoing degree one). Therefore, DDS are simple structures that can be applied to any phenomena that evolve according to a function. Often, when representing phenomena with DDS, one needs to validate "macro" dynamics coming from experimental results, or to identify the set of "micro" dynamics which generate a given "macro" behavior observed.

In [START_REF] Dennunzio | Polynomial equations over finite, discrete-time dynamical systems[END_REF], operations of sum and multiplication have been introduced to explain how DDS can be combined to generate new dynamical behaviors. Equipping the set of DDS with these operations provides a commutative semiring, and this naturally leads to write polynomial equations over DDS. This is interesting because these equations (with a constant right-hand term) can model hypotheses over the dynamics, and if one can prove properties for DDS then they can be automatically lifted to the application models. The idea is to solve these equations to validate the corresponding hypotheses. This paper introduces a pipeline to solve equations over the cycles of dynamics graphs (i.e. equations/hypotheses over the long term behavior of a phenomenon), formally introduced in [START_REF] Dennunzio | Solving equations on discrete dynamical systems[END_REF]. The first part of this paper proposes a new algorithm which outperforms the stateof-art technique (Colored-Tree method) using MDDs to enumerate the solutions set of simple equations. The second part introduces a pipeline to solve general equations using MDDs to solve basic cases and to limit the exploration of the solutions space. The overall purpose is to introduce the first complete pipeline, based on MDD, to solve the equations introduced in [START_REF] Dennunzio | Polynomial equations over finite, discrete-time dynamical systems[END_REF][START_REF] Dennunzio | Solving equations on discrete dynamical systems[END_REF] to study the hypothesis over the long term behavior of a phenomenon modeled by DDS.

Preliminaries

Multi-valued decision diagrams (MDDs)

Multi-valued decision diagrams are the extension of Binary decision diagrams (BDD) in which the diagram consists of a rooted acyclic graph able to represent a multi-valued function f : {0...d -1} r → {true, f alse}. Considering a generic MDD, each level represents a variable and a final layer contains the true terminal node (tt). Therefore, the data-structure contains r + 1 layers and a path from the root to the tt node represents a valid set of assignments. Each node is characterised by the variable represented in the layer and at most d outgoing edges. In general, the edges are directed from an upper layer to a lower one. An edge corresponds to an assignment of the variable to a certain value specified over the edge. This data-structure presents only one root and one leaf. The false node, and the corresponding paths to reach this node, are omitted.

According to [START_REF] Amilhastre | Compiling csps: A complexity map of (non-deterministic) multivalued decision diagrams[END_REF][START_REF] Darwiche | A knowledge compilation map[END_REF] , an MDD is deterministic if all the nodes have pairwise distinct labels on outgoing edges. Moreover, an MDD is ordered if given two nodes A and B such that A is the parent of B, we have that the variable represented into the node A is the smallest (w.r.t. a given total order over the set of variables). Fig. 1: The MDD corresponding to the valid assignments {(1, 2, 3), (2, 2, 3), [START_REF] Amilhastre | Compiling csps: A complexity map of (non-deterministic) multivalued decision diagrams[END_REF][START_REF] Andersen | An introduction to binary decision diagrams[END_REF][START_REF] Andersen | An introduction to binary decision diagrams[END_REF], (1, 2, 2)} for the variables x 0 , x 1 and x 2 . The structure presents 4 layers. The root and its edges represent the variable x 0 and its possible assignments.

One of the most interesting aspects of MDDs is their reduction. According to this procedure, they can gain an exponential factor in representation space. The reduction of MDDs aims at merging equivalent nodes i.e. that have equivalent outgoing paths. In particular, two nodes are equivalent if they have the same label and destination for each edge. The idea is to identify and merge equivalent nodes from the bottom to the top of the MDD, see Figure 2 for an example. Several reduction algorithms for MDDs have been introduced [START_REF] Andersen | An introduction to binary decision diagrams[END_REF][START_REF] Cheng | An mdd-based generalized arc consistency algorithm for positive and negative table constraints and some global constraints[END_REF]. In our software we use the pReduce algorithm [START_REF] Perez | Efficient operations on mdds for building constraint programming models[END_REF] since the time complexity per node is bounded by its number of outgoing arcs and the complexity is linear on the size of the MDD.

Another big advantage of using MDDs is that classical set operations such as Cartesian product, complement, intersection, union, difference, and many others can be performed without decompressing the structure. For instance, the cartesian product over MDDs can be performed just by transforming the root of an MDD into the tt node of another one (see Figure 3 for an example). Given two MDDs, the intersection algorithm creates a new MDD to represent the solutions contained in both the original ones (see Figure 4). To achieve this goal, the process starts with the creation of a new root, and only the outgoing edges that are common between the two structures are recreated. In this way, each new node corresponds to two original nodes and the process is iterated. It is important to remember that in the end, it is necessary to verify if there are nodes without children; if any, they must be deleted. The drawback of this algorithm is that the resulting MDD can be larger than the original ones. For more details, we refer the reader to [START_REF] Perez | Efficient operations on mdds for building constraint programming models[END_REF] and [START_REF] Bergman | Mdd propagation for sequence constraints[END_REF]. Fig. 4: Two MDDs (left and center) and their intersection (right).

Discrete Dynamical Systems and dynamics graphs

A Discrete Dynamical System (DDS) is a structure χ, f where χ is a finite set of states and f : χ → χ is a function called the next state map. When modelling a phenomenon by a DDS χ, f , χ is its set of states and f is the law which brings from state α ∈ χ at time t to the state f (α) at time t + 1. When χ is finite, any DDS χ, f can be identified with its dynamics graph G ≡ V, E where V = χ and E = {(α, β) ∈ V × V, f (α) = β}. Therefore, all the properties of the DDS can be deduced from the properties of its dynamics graph. As a first property, one can remark that they are graphs with outgoing degree one and hence each strongly connected components of such graphs is made by a single cycle (or loop). From now on, we will turn all the discussion about DDS in terms of dynamics graphs. Call DG the set of all dynamics graphs up to (graph) isomorphism. One can define on DG two operations: sum and product as follows. Given two graphs,

G 1 = V 1 , E 1 ∈ DG and G 2 = V 2 , E 2 ∈ DG, the sum G 1 + G 2 is the graph G = V 1 V 2 , E 1 E 2 ∈ DG, where is the disjoint union operator. The product G 1 • G 2 is the graph V , E ∈ DG with V = V 1 × V 2 and E = {((α 1 , α 2 ), (β 1 , β 2 )) ∈ V × V , (α 1 , β 1 ) ∈ E 1 and (α 2 , β 2 ) ∈ E 2 }.
The product defined above consists in the parallel synchronous execution of the two dynamics graphs. The sum is the mutually exclusive alternative between the two behaviours.

R := DG, +, • is a commutative semiring in which ∅, ∅ is the neutral element w.r.t. + and {α} , {(α, α)} is the neutral element w.r.t. multiplication. Now, consider the semiring R[x 1 , x 2 , . . . , x s ] of polynomials over R in the variables x i , naturally induced by R. Polynomial equations of the form (1) model hypotheses about a certain dynamics deduced from experimental data.

a 1 • x w1 1 + a 2 • x w2 2 + . . . + a k • x ws s = C (1)
Equation ( 1) can be interpreted as follows. The constant term C on the righthand side is the dynamical system deduced from experimental data. The point is that C comes just from experimental data and hence it might be the "macro" result of many cooperating hidden variables at a "micro" level. On the left hand side of (1), we have a hypothesis on the "micro" structure based on partial information (the coefficients) and unknown information (the variables). In other words, the coefficients a i are hypothetical sub-dynamical systems that should cooperate to produce the observed dynamics C. Finding valid values for the unknown variables provides a finer structure for C which can bring further knowledge about the observed phenomenon. More generally, one can interpret Equation (1) as a question over dynamics graphs (i.e. directed graphs with outgoing degree 1). The constant right-hand side represents the current graph and the left-hand side is a question (hypothesis) about a possible decomposition (according to the semi-ring operations).

In [START_REF] Dennunzio | Polynomial equations over finite, discrete-time dynamical systems[END_REF], it has been proved that finding solutions to generic polynomial equations (i.e. in which both left and right-hand side of the equation are made by polynomials) over DDS is undecidable, while the problem is decidable when the right-hand side of the equation is constant. However, even in the decidable case, the complexity of the problem is beyond NP, except for very particular cases. Some abstractions are introduced to progressively filter the solutions space according to features of the real solutions. Therefore, the general solutions are found in the intersection of the solutions of these abstractions. For this reason, the solutions enumeration of each abstraction is fundamental.

At least three abstractions can be devised on dynamics graphs, namely, abstraction on cycles, on cardinality (of the vertex set) and on paths. Studying each abstraction separately allows to study different aspects of a dynamics. In particular, solving equations over cycles leads to the validation of hypotheses over the long term behavior of a phenomenon.

The abstraction on cycles

In this paper we analyze equations over the cyclic part of the dynamics. Considering a generic Equation (1), we introduce a pipeline to solve the abstraction over the long term behavior but before we need to recall some notation and some concepts that will be useful in the sequel.

For a dynamics graph G, let G be the subgraph of G which contains only the cycles and the loops. Denote R the restriction of R such that for any G ∈ R, G ∈ R. It is not difficult to see that R, +, • is a sub-semiring of R, +, • . Thus, solving Equation (1) over R is a necessary step for solving it over R. This also implies that we have to enumerate all solutions in R first and then filter them out using the other abstractions to find the general solutions. Notation 1 A cycle {α 1 , α 2 , . . . , α p } of length p can be conveniently denoted by

C 1 p . A subgraph, with K different lengths of cycles is denoted K i=1 C ni pi .
Fig. 5: A dynamics graph G, in which the subgraph G is the red part and it is denoted

C 2 2 ⊕ C 1 3 ⊕ C 1 5
according to our notation. The operations of sum and product over R can be conveniently applied to the new notation. Graphically, we can consider the result of a sum as a new system composed by all the cycles of the input systems, and the product one as a new system generated with the Cartesian product of the cycles of the input systems.

Definition 1 (Sum). Consider two dynamics graphs Å

≡ K A i=1 C n Ai p Ai and B ≡ K B j=1 C n Bj p Bj , Å⊕ B is K A i=1 C n Ai p Ai ⊕ K B j=1 C n Bj p Bj = C n A1 p A1 ⊕. . .⊕C n AK A p AK A ⊕C n B1 p B1 ⊕. . .⊕C n BK B p BK B .
Considering two sets of cycles C n Ai p Ai and C n Bj p Bj , if they have the same cycles length (p Ai = p Bj ), then they can be rewritten like C

n Ai +n Bj p Ai . Definition 2 (Product). Consider Å ≡ K A i=1 C n Ai p Ai and B ≡ K B j=1 C n Bj p Bj , Å B is K A i=1 C n Ai p Ai K B j=1 C n Bj p Bj = K A i=1 K B j=1 C n Ai p Ai C n Bj p Bj = K A i=1 K B j=1 C n Ai •n Bj •gcd(p Ai ,p Bj ) lcm(p Ai ,p Bj )
where gcd is the greatest common divisor and lcm is the least common multiple.

With the help of the previous notation, Equation [START_REF] Akers | Binary decision diagrams[END_REF] can be rewritten as

( K1 i=1 C n1i p1i x1 w1 ) ⊕ ( K2 i=1 C n2i p2i x2 w2 ) ⊕ . . . ⊕ ( Ks i=1 C nsi psi xs ws ) = m j=1 C nj pj (2)
where K z is the number of distinct cycles size in the system a z with z ∈ {1, . . . , s} (cf. Equation ( 1)) and n zi is the number of cycles of length p zi of a z . In the right term C, there are m different periods, where for the j th different period there are n j cycles of period p j . However, Equation ( 2) is still hard to solve in the present form. We can simplify it further by performing a contraction step which consists in cutting Equation (2) into two simpler equations

             C n11 p11 x1 w1 = m j=1 C uj pj C 1 1 ẙ = m j=1 C vj pj (3a) (3b) with ẙ = ( K1 i=2 C n1i p1i x1 w1 ) ⊕ ( K2 i=1 C n2i p2i x2 w2 ) ⊕ . . . ⊕ ( Ks i=1
C nsi psi xs ws ) and n j = u j + v j for j ∈ {1, . . . , m}. By recursively applying contraction steps, and for all possible u j , v j values in Equations (3a) and (3b), solving Equation ( 2) boils down to solve multiple times simple equations i.e. equations of the form

C 1 p X = C n q (4)
where p ∈ {p 11 , p 12 , . . . , p sKs }, q ∈ {p 1 , p 2 , . . . , p m }, and n is smaller than the number of cycles of length q in the right part. Assume that X is xz wz such that z ∈ {1, . . . , s}. It is necessary to find xz from X i.e. we need to compute the w z root of X. Given 2M integers p i , k i ∈ N, with p i > 0 for all i ∈ {1, . . . , M }, let l(p 1 , p 2 , ..., p t , k 1 , k 2 , ..., k t ) be the lcm between the p i for which

k i = 0 and t ≤ M (with l(p 1 , p 2 , ..., p t , k 1 , k 2 , ..., k t ) = 1 iff ∀1 ≤ i ≤ t, k i = 0). Assume Å ≡ C 1 p1 ⊕ C 1 p2 ⊕ ... ⊕ C 1 p M . Then, Å w ≡ M i=1 C p w-1 i pi ⊕ k1+k2+...+kM =w 0≤k1,k2,...,kM <w w k1,k2,...,kM C M t=1 kt =0 p k t -1 t • M t=2 kt =0
gcd(l(p1,...,pt-1,k1,...,kt-1),pt) l(p1,p2,...,pM ,k1,k2,...,kM )

.

(5)

The state-of-art method

In [START_REF] Dennunzio | Solving equations on discrete dynamical systems[END_REF], two computational problems are devised concerning simple equations.

The SOBFID (SOlve equation on BIjective Finite DDS ) problem is a decision problem which takes in input p, n, q ∈ N \ {0} and returns true iff C 1 p X = C n q admits a solution. EnumSOBFID is the problem which takes the same input as SOBFID and outputs the list of all solutions of C 1 p X = C n q . To the best of our knowledge, the Colored-Tree Method (CTM) proposed in [START_REF] Dennunzio | Solving equations on discrete dynamical systems[END_REF] is the best current technique to solve this problem. The method exploits a connection between EnumSOBFID and the well-known Change-making problem [START_REF] Martello | Knapsack Problems: Algorithms and Computer Implementations[END_REF] coupled with a completeness-check (running in exponential time) to explore the feasible solutions space. Essentially, the right term n is decomposed in every possible way and these possibilities are represented in a tree. The method comprises two main phases: tree building and solutions aggregation. In the first one, the algorithm uses a tree to decompose the right part of the equation in a certain number of subgraphs, searching for each node of the tree (each subgraph) the minimum number of product operations (minimum number of cycles in the variable or minimum number of children) which are necessary to produce it. The subgraphs found are then divided into subsets of children. Iterating this idea on each subset produced, the method arrives to enumerate all the possible ways to generate the cycles involved in the right part. During the second phase, the method computes (bottom-up) the real solutions of the equation represented in the tree.

Boosting everything up with MDDs

The enumeration of solutions for Equation ( 2) is one of the main objective of this paper. In order to achieve this, we solve the EnumSOBFID problem and show that MDDs boost up the enumeration. In other words, we are interested to find all the possible ways to generate C n q cycles starting with one cycle of length p. This last problem is similar to the well-known Change-making problem (in its enumeration version) in which one aims at finding all the possible ways to change a total amount with a given set of coins. In our case, the total amounts are the C n q cycles but to complete the similarity we need to find the coins which can be part of a solution. Remark that a cycle C 1 u generates C r q cycles in the right part iff r divides q, u = q p • r, gcd(p, q p • r) = r and lcm(p, q p • r) = q. A cycle C 1 u with the previous properties is called feasible and r is a feasible divisor of q.

Let D p,q = {d 1 , . . . , d l } be the set of feasible divisors (w.r.t. Equation (4) of course). There is a solution to (4) iff there exist x 1 , . . . , x l such that l i=1 d i x i = n (i.e. there is a solution to the Change-making problem for a total amount n and a coins system D p,q ). To solve EnumSOBFID we need to enumerate all solutions of the previous Change-making problem. At this point MDDs come into play. We are going to use them to have a compact and handful representation of the set of all possible solution to Equation (4).

SB-MDD

The MDD M p,q,n containing all solutions to Equation ( 4) is a labelled digraph V, E, where

V = Z i=1 V i with Z = n min Dp,q + 1 and V 1 = {root}, V i is a multiset of {1, .
. . , n -1}, and, finally, V Z = {tt}. For any node α ∈ V , let val(α) = α if α = root and α = tt, val(root) = 0 and val(tt) = n. For any i ∈ {1, . . . , Z -1} and for any α ∈ V i and β ∈ V i+1 ∪ {tt}, (α, β) ∈ E iff val(β) -val(α) ∈ D p,q and val(β) ≤ val(tt). The labelling function : E → D p,q is such that for any (α, β) ∈ E, ((α, β)) = val(β) -val(α) ∈ D p,q .

Graphically, M p,q,n can be represented by layers as usual, interpreting each layer as the usage of the i-th coin. Each node represents the sum of coins from the root to the node. Moreover, the longest path in M p,q,n is Z = n min Dp,q + 1. Remark that M p,q,n contains duplicated solutions, indeed, it contains all the permutations of a solution but according to Equation (4) different permutations lead to the same solution. For this reason, we impose a symmetry breaking constraint: for any node α (different from tt), let e be the label of the incoming edge; the only allowed outgoing edges of α are those with label ≤ e. In this way all the paths of the MDD will be ordered and the size of the MDD will be smaller. An SB-MDD is an MDD which satisfies the symmetry breaking constraint. The building of M p,q,n ends with a pReduction which merges equivalent nodes and deletes all nodes (and the corresponding edges) which are not on a path from root to tt. Let us illustrate the construction with an example.

Example 1. Consider the simple equation C 1 2 X = C 6 6 .
The set of divisors of q (smaller or equal to n) is {6, 3, 2, 1}. However, D p,q = {1, 2}. Indeed, we have

r = 6 ∧ u = 18 → gcd(2, 18) = 6 ∧ lcm(2, 18) = 6 r = 3 ∧ u = 9 → gcd(2, 9) = 3 ∧ lcm(2, 9) = 6 r = 2 ∧ u = 6 → gcd(2, 6) = 2 ∧ lcm(2, 6) = 6 r = 1 ∧ u = 3 → gcd(2, 3) = 1 ∧ lcm(2, 3) = 6
Figure 6 shows M 2,6,6 in its classic form (left) and in its SB-MDD form (right). Remark that in Figure 6 (left) many solutions are duplicated. For example, the solution [2, 2, 1, 1] (in red) is represented (more than) twice. Once M 2,6,6 is built and reduced, reading solutions correspond to paths labels from root to tt: The solution [2, 2, 2] says that 6 is changed with 3 coins of value r = 2. Recalling that u = q p • r we can express the solution in term of dynamics graphs as C 3 6 . Operating similarly for all the other solutions we find

{[2, 2, 2], [2, 2, 1, 1], [2, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1]}.
C 3 6 , C 2 6 ⊕ C 2 3 , C 1 6 ⊕ C 4 3 , C 6 3 .

Equations over dynamics graphs

We have shown how MDDs can be used to compute the solutions set of simple equations. In this section, we introduce a pipeline to solve Equations of form (2) to validate hypotheses over discrete systems represented by DDS. The goal is the enumeration of solutions. The pipeline consists in the following steps:

identification and resolution of the necessary equations (each necessary equation corresponds to an SB-MDD); enumeration of the contractions steps by an MDD structure; computation of the solutions of each contraction step with a particular technique to compute the intersection between SB-MDDs. Necessary equations. As we have already seen, to solve a generic equation over dynamics graphs we need to find the solutions of a certain number of contraction steps (Equations (3a) and (3b)). Each of these contraction steps ends up with a system of equations of the form

C nzi pzi X z = m j=1 C vj pj (6)
where z ∈ {1, . . . , s}, i ∈ {1, . . . , K z }, and v j ≤ n j . This means that each monomial is responsible for the generation of a subgraph of the right term. In other words, the solution of ( 6) is the Cartesian product of the solutions of a certain number of simple equations. Remark that a single simple equation might occur several times while searching for a solution of a single contraction step or in multiple contraction steps.

According to the product rules, each equation

C nzi pzi X z = C vj pj is equivalent to C 1 pzi X z = C v j n zi
pj . Remark that if v j /n zi is not an integer, then the equation has no solutions.

Since there are many contraction steps that need to be explored, we aim at limiting the exploration by ignoring those which involve simple equations without solutions. Therefore, the method starts with the computation of the set of simple equations that can be involved in the result of a contractions step. This amounts to compute all the SB-MDDs M pzi,pj , n n zi with p zi ∈ {p 11 , . . . , p sKs }, p j ∈ {p 1 , . . . , p m }, for all n ∈ {1, . . . , n j }. In this way, even if a simple equation is involved in many contraction steps, it is solved only once. A simple equation with at least one solution is called necessary equation. It is important to notice that it is not necessary to explore an SB-MDD to decide if a solution exists. Indeed, by construction, an SB-MDD is generated iff it contains at least a solution.

Contractions steps. Once the set of necessary equations has been computed, we create an MDD CS to enumerate all the contractions steps that must be taken into account to enumerate the solutions of Equation ( 2). This MDD will be a Cartesian product of other MDDs, one for each p j of the right term, i.e. CS = × m j=1 CS j . Consider a cycle length p j in the constant term. The MDD CS j accounts for all the feasible ways (according to the set of necessary equations) to generate n j cycles of length p j using the monomials of the equation. Therefore, CS j is a labelled digraph V j , E j , j where V j = s+1 z=1 Kz i=1 V j,zi with K s+1 = 1 (see Equation ( 2) for the meaning of K i ) and V j,11 = {root}, V j,zi ⊆ {0, ..., n j }, and, finally, V j,(s+1)1 = {tt}. For any node α ∈ V j , let val(α) = α if α = root and α = tt, val(root) = 0 and val(tt) = n j . The set of possible outgoing edges of level V j,zi is D pzi,pj = {g ∈ N | 1 ≤ g ≤ n j and M pzi,pj , g n zi ∈ necessary equations} ∪ {0} . For any α, β ∈ V j , (α, β) ∈ E j iff 1. α ∈ V j,zi and either β ∈ V j,z(i+1) for some i < K z or β ∈ V j,(z+1)1 ; 2. val(β) -val(α) ∈ D pzi,pj and val(β) ≤ val(tt).

In this MDD the outgoing edges of a level V j,zi represent the cycles of length p j generated by the monomial C nzi pzi X z of the left part. The labelling function j : E j → s z=1 Kz i=1 D pzi,pj is such that for any (α, β) ∈ E j , j ((α, β)) = val(β)-val(α) ∈ D pzi,pj with α ∈ V j,zi . Starting from each node, a label 0 means that the monomial is not involved in the generation of the cycles of length p j . The sum of the labels of each path from the root to the tt node will be equal to n j , because n j cycles of length p j must be generated.

For each path of CS, it is necessary to perform some some additional steps to understand if it leads to feasible solutions of the equation as explained in the next section.

Solve a system. Each contraction step corresponds to a system of Equations of type [START_REF] Berndt | Multi-valued decision diagrams for the verification of consistency in automotive product data[END_REF]. To compute the solutions set of these equations one needs to compute the Cartesian product between the different solutions of the corresponding simple equations. In its turn, a solution to a simple equation is represented by a SB-MDD. Therefore, their Cartesian product is computed in linear time by placing the SB-MDDs one on top of the other to form a new MDD. We call SB-Cartesian MDD an MDD build in such a way. Remark that an SB-Cartesian MDD is not a SB-MDD.

Recall that X z in Equation ( 6) represents a variable xz wz and if xz is involved in different monomials, then, an intersection operation is required. To compute the solutions set of the variable, we start considering equations involving xz with the same power w z by computing the intersection over the corresponding MDDs. Moreover, remark that each xz wz corresponds to a SB-Cartesian MDD (except for the case in which a monomial is responsible for the generation of only one cycle length). However, notice that the classic algorithm to perform the intersection over MDDs cannot be used if the goal is the intersection between MDDs issued by a Cartesian product (i.e. SB-Cartesian MDDs) because the result depends on the order of the MDDs. In the next section, we propose a new algorithm to perform this task independently from the order.

Once xz wz is assigned with a set of solutions, we need to compute the w z -th root for each value of xz wz and finally the intersection between all the roots found so far.

Finally, we stress that the root procedure is not a trivial step. Indeed, the inverse operations of sum and product are not definable in the commutative semiring of DDS. Therefore, we need an algorithmic technique to compute the result of the w-th root of x. Considering Formula (5), the root can be computed combinatorially or through a finite number of polynomials equations over real numbers of increasing degree. We combine both approaches in order to speedup the computations. Consider a generic system x = C

n 1 p 1 ⊕ C n 2 p 2 ⊕ . . . ⊕ C n l p l such that xw = C n1 p1 ⊕ C n2 p2 ⊕ . . . ⊕ C n h p h
, the number of cycles n 1 , of the minimum length p 1 , involved into the root's solution is computed with the polynomial equation (p 1 ) w-1 • (n 1 ) w = n 1 with p 1 = p 1 . The remaining part of the solution is combinatorially computed knowing that a length of cycle p i may be involved in the root solution iff n i ≥ (p i ) w-1 , and we will have a maximum number of cycles of this length equals to ni (pi) w-1 .

SB-Cartesian Intersection. Consider a set M of MDDs in which some are SB-Cartesian MDDs and others are not. We propose an algorithm (Algorithm 1) which computes the intersection of all the elements in M . Our algorithm needs to be started with an initial set of candidate solutions S (initial guess). If M contains at least a simple SB-MDD (i.e. one which is not a SB-Cartesian MDD), then S is the set of solutions read in the MDD resulting from classical intersection of the SB-MDDs in M ; otherwise S is the set of solutions read in an arbitrarily chosen SB-Cartesian of M . Using S, we will compute the intersection between the remaining SB-Cartesians. The idea is to search the solutions into each SB-Cartesian MDD and update each time the remaining solutions. Each candidate solution is ordered and recursively searched in a SB-Cartesian MDD in M (Algorithms 2 and 3). If it is not found, then it is removed from the set of candidate solutions. Given a SB-Cartesian element M of M , a candidate solution s is validated if it possible to visit each SB-MDD involved in M using a subset of the elements of the solution. Starting from the biggest elements of s, we try to find a path from the root to the first tt node (recall that we are visiting a SB-Cartesian MDD). We proceed in this way because the paths of each SB-MDD are ordered. After having gone through the first SB-MDD, we need to recursively repeat the procedure with the remaining elements of s. Two special cases need our attention:

a generic node in which it is not possible to find a common element between the remaining elements of a solution and the outgoing edges; a node (different from the final tt node) in which there are no more elements of s to compare with the outgoing edges. In both cases, it is necessary to go back in the visited nodes until there is another possible outgoing edge that can be taken into consideration. Once the validation procedure arrives at the last SB-MDD of M , a linear search is performed with the remaining elements. Finally, s ∈ M if there exists a path from the last root to the final tt node following the remaining elements (Algorithm 4). If the linear search fails, then we return to the first node with a different feasible outgoing edge. If no different feasible edges exist, then s is not a solution. The previous procedure is performed over each candidate solution of the initial guess. We stress that in the worst case, for a given candidate solution, our algorithm explores only the subgraph of a SB-Cartesian MDD made of feasible edges. The approach introduced to compute the intersection of SB-Cartesian MDDs is suitable for our needs, its improvement constitutes an interesting future research direction. We need to precise an additional aspect of our application. As explained above, an SB-MDD corresponding to a simple equation has la-bels based on the feasible divisors set (or feasible coins) and each divisor/coin corresponds to a certain cycle C 1 u . However, two different coins of two different simple equations can correspond to the same C 1 u as well as the same coin may correspond to different cycles lengths for different simple equations. Therefore, when searching for a solution in a SB-Cartesian MDD, we must take into account these cases. In our application, two divisors/coins are considered equivalent if they correspond to the same C

Experiments

The experimental evaluation is divided into two parts: one concerns simple equations and the other one is devoted to the complete method which solves generic equations.

Concerning equations of the form C 1 p X = C n q , in our experiments, we set p = q since this grants the existence of at least a solution. Using the MDDs it is possible to outperform the Colored-Tree method (CTM) w.r.t. both memory and time. If we compare the dimension (in terms of nodes) of a colored-tree with the corresponding SB-MDD for a given equation, the second one is smaller. CTM presents some out of memory cases even for equations with n, q, and p smaller than 30 and memory limit of 30GB (see Figure 7 (left)). Using MDDs, we solved equations with n, q, and p up to 100 without any out of memory case and only 6GB RAM limit (see Figure 7 (right)).

Analysing the time to solve equations with parameters up to 30, it turns out that the new technique is faster than the previous one (see Table 1). The reason is that CTM requires a time consuming check procedure to ensure the completeness of the solutions which is not necessary in the MDD case. Due to too high memory and time costs, CTM is unsuitable to solve simple equations coming from contractions steps. The new method fixes these issues allowing to solve generic polynomial equations. Fig. 7: number of nodes for the colored-tree with memory limit of 30GB (left) and for the SB-MDD with memory limit of 4GB (right) in the case of equation for type C 1 q X = C n q . Remark that the black square in the right part of the left diagram are out of memory cases.

Turning to generic equations, we use the proposed pipeline to find the solutions (if any). If roots computation is not performed then everything depends on the number of monomials, the number of distinct sizes of cycles in the right side of the equations. If these quantities grow also the number of contraction steps to consider grows. Solutions spaces are limited by considering only contraction steps that are feasible according to the necessary equations. For example, consider the equation

C 1 5 ẙ ⊕ C 1 4 ẙ2 ⊕ C 1 3 x2 ⊕ C 1 2 x3 = C 3 10 ⊕ C 68 4 ⊕ C 9 3 ⊕ C 9 6 ⊕ C 136 12
, the number of contraction steps is ≈ 2, 42 • 10 16 , but only 6665400 are considered.

The computation of roots is the most time consuming operation. Indeed, consider the equation

C 2 3 x ⊕ C 4 3 ẙ = C 162 3 ⊕ C 20 6 ⊕ C 104 12
, we can find 49329000 solutions considering the 6642 feasible contraction steps in only 58 seconds. If we consider the same equation but with x2 in place of x, then the number of contraction steps explored is the same, but we found 4510 solutions in 5.6 hours (therefore the average computational time per contraction is 3.04 seconds).

Roots computation is expensive and this is reasonable in a sense. Since the division and subtraction operations cannot be performed directly, one needs sophisticated techniques. The pipeline introduced above leads to the first complete technique to solve equations over the asymptotic behavior of DDS. In the end, the pipeline is already proving itself suitable for applications over real experimental data.

Conclusions and perspectives

Equations on DDS are useful to analyze dynamics of phenomena. They allow to model hypotheses on the dynamical behavior and solving them leads to their validation or invalidation. In particular, Equation (2) allow studying the longterm aspects of the dynamics. This paper introduces new algorithms to solve these equations over cycles (i.e. long term behavior). In the case of simple equations, our technique outperforms the CTM. This is an important breakthrough because CTM was used to solve simple equations (generated by contractions steps) and it was practically unusable due to huge memory consumption and time-consuming check procedures. Moreover, this paper proposes a pipeline to solve general equations. The pipeline computes the necessary equations to limit the exploration of the contractions steps and, in the end, it solves each system of equations corresponding to a contractions step. This allows the treatment of much larger dynamics graphs and much more complicated hypotheses.

Future perspectives include the improvement of the pipeline by parallelising the identification of solutions of different feasible contractions steps. Another research direction would try to speedup roots computation by increasing the number of coefficients computed through polynomial equations to reduce the combinatorics. In the end, this work aims to call for further research in MDDs. In fact, studying different ways to perform the intersection between SB-Cartesian MDDs is surely another subject that is worth exploring. This research work leads to a complete and performing pipeline to validate hypotheses over the long term behavior of dynamics graphs adding one more item to the growing list of successful applications of MDDs.

Fig. 2 :

 2 Fig. 2: An MDD before reduction (left) and its reduced version (right).

Fig. 3 :

 3 Fig. 3: The MDD representing the Cartesian product of the MDD in Figure 1 and the MDD in Figure 2 (right). The diagram is drawn horizontally for lack of room.

Fig. 6 : 2 X = C 6 6

 626 Fig. 6: The reduced MDD representing all the solutions of C 1 2 X = C 6 6 in its classic form (left) and in its SB-MDD form (right).

Algorithm 1 :

 1 SB-Cartesian Intersection Input : M , set of MDDs Output: S, solutions of the intersection in M Cartesian← ∅; Traditional← ∅; forall m ∈ M do if m is a SB-Cartesian MDD then Cartesian.add(m) ; else Traditional.add(m) ; S ← ∅; if |Traditional|> 0 then if |Traditional| = 1 then S ←Traditional[0].readSolutions(); else MddIntersected←ClassicIntersection(Traditional[0],Traditional[1]); forall m ∈Traditional \ {T raditional[0], T raditional[1]} do MddIntersected←ClassicIntersection(MddIntersected,m); S ←MddIntersected.readSolutions(); else S ←Cartesian[0].readSolutions(); Cartesian.remove(0); if |S| = 0 then forall m ∈Cartesian do CartesianSearch(S,m); return S

Table 1 :

 1 Computation times (millisec) of CTM (left) and MDD (right) over different input parameters. Symbols '-' represent out of memory cases.

		q,p	n	2	3		4	5	6		7		8	9	10
		2	56|0	56|0		54|0	54|0	53|0	54|0	54|0	54|0	53|0
		3	54|1	54|0		55|1	54|0	55|1	54|0	55|1	54|0	55|1
		4	55|2	54|1		56|2	53|0	57|2	54|0	55|2	55|1	55|2
		5	55|2	56|2		59|2	54|0	61|30	54|0	58|2	55|3	56|2
		6	58|2	56|2		60|2	58|1	63|9	55|0	60|2	56|2	62|2
		7	60|3	58|2		63|19	56|1	78|21	54|0	63|29	60|2	61|2
		8	63|2	59|2		96|10	60|2	107|20	56|1	97|9	61|2	65|2
		9	66|3	60|3		106|21	57|2	153|22	57|1	168|21	60|3	76|20
		10	84|3	62|2		140|11	58|2	185|21	57|2	369|11	70|2	120|17
	q,p	n	11	12		13	14		15	16		17	18	19	20
	11		55|0	824|25		55|0	116|3		74|21	406|22		55|0	1071|22	57|0	1334|23
	12		58|1	17678|26	57|0	132|4		88|21	4105|12		56|0	6022|22	56|0	3672|22
	13		61|2	177894|27	56|0	246|21		92|21	4163|27		55|0	5967|24	56|0	3332|24
	14		59|2	1277979|26	61|1	900|11		116|22	19895|24	56|0	27381|27	56|0	96997|26
	15		60|2	-|28		60|2	3721|22		169|22	19711|25	56|0	637457|26	59|0	419000|25
	16		62|2	-|29		61|2	19900|12	502|23	-|13		57|0	1185947|26	60|0	759365|26
	17		62|2	-|30		62|2	25908|24	554|23	-|26		57|0	-|27	57|0	-|27
	18		64|2	-|46		62|2	164167|13	1102|22	-|26		61|2	-|28	61|0	-|30
	19		66|2	-|32		63|2	226315|25	950|24	-|27		62|2	-|34	57|0	-|39
	20		68|2	-|32		65|2	1707299|25 2749|24	-|16		63|2	-|31	62|2	-|29
	q,p	n	21	22		23	24		25	26		27	28	29	30
	21	2712|23	343542|26	60|0	-|35		95|4	389971|7	2034|12	-|28	58|0	-|37
	22	23399|24	-|14		62|0	-|36		103|4	381929|7	2711|24	-|30	59|0	-|35
	23	27430|24	-|27		59|0	-|38		134|4	-|7		2712|24	-|31	59|0	-|37
	24	149296|25	-|16		64|2	-|39		149|4	-|8	20641|14	-|42	59|0	-|41
	25	162413|25	-|28		65|2	-|40		160|4	-|26	24632|24	-|34	59|0	-|39
	26	212277|25	-|16		66|2	-|42		403|5	-|15	24177|25	-|33	60|0	-|40
	27		-|25	-|27		69|2	-|45		454|4	-|33		-|15	-|34	59|0	-|45
	28		-|26	-|17		70|2	-|47		488|5	-|16		-|25	-|35	60|0	-|44
	29		-|28	-|28		69|2	-|66		506|5	-|28		-|27	-|36	61|1	-|46
	30		-|26	-|27		69|3	-|51		838|6	-|17		-|17	-|38	65|2	-|48

Acknowledgments This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.