# Reconfigurable Hardware Implementation of a Multivariate Polynomial Interpolation Algorithm

Authors: Rafael A. Arce-Nazario; Edusmildo Orozco; Dorothy Bollman
Journal: International Journal of Reconfigurable Computing
 (2010)

---

## Abstract

Multivariate polynomial interpolation is a key computation in many areas of science and engineering and, in our case, is crucial for the
solution of the reverse engineering of genetic networks modeled by finite fields. Faster implementations of such algorithms are needed to
cope with the increasing quantity and complexity of genetic data. We present a new algorithm based on Lagrange interpolation for multivariate
polynomials that not only identifies redundant variables in the data and generates polynomials containing only nonredundant variables, but
also computes exclusively on a reduced data set. Implementation of this algorithm to FPGA led us to identify a systolic array-based architecture
useful for performing three interpolation subtasks: Boolean cover, distinctness, and polynomial addition. We present a generalization of these tasks
that simplifies their mapping to the systolic array, and control and storage considerations to guarantee correct results for input sequences
longer than the array. The subtasks were modeled and implemented to FPGA using the proposed architecture, then used as building blocks to
implement the rest of the algorithm. Speedups up to172× and 67× were obtained for the subtasks and complete application, respectively, when compared to a software implementation, while
achieving moderate resource utilization.


---

## Body

## 1. Introduction

Recent years have seen a significant increase in methods and tools to collect genetic data from which important information can be extracted using a number of techniques [1]. For instance, microarray data collected at various steps in organism development can help geneticists understand its developmental process and response to environmental stimuli [2]. Several models such as ordinary differential equation models [3], continuous models [4], stochastic models [5], and discrete models, of which Boolean models [6] are special cases, have been proposed. Essentially, each node represents the expression level of a gene (or genes) of interest at a time point. A directed edge from node a to node b symbolizes the influence of gene(s) in a to the gene(s) in b.Our research group focuses on multivariate finite field gene network (MFFGN) models, in which multiple genes are monitored at each time step and their expression levels are discretized to a predefined set of values{0,1,2,…,p-1}, where p is a prime number, that is, the expression level of each gene is an element of the prime field ℤp [7, 8]. One way to deduce gene interaction from these models is to solve the reverse engineering problem, that is, find functions that describe the network's state transitions. An important step in the solution of this problem is determining an interpolating polynomial for the points by using such methods as a multivariate version of Lagrange's interpolation formula [8]. Polynomial interpolation over finite fields also has applications in error correcting codes and is a major building block of many numerical methods [9, 10].Fast algorithms and implementations are needed to sustain rapidly increasing bioinformatics computational demands [11]. Reconfigurable logic represents an attractive platform for the high performance implementation of finite field algorithms, with many designs achieving significant speedups over their software equivalents. However, despite the availability of tools to ease the hardware design flow, there is still a notable gap between the bioinformatics and hardware experts [12]. Our group developed a reconfigurable logic implementation of a multivariate polynomial interpolation algorithm suitable for reverse engineering in genetic networks. Rather than striving for a specialized punctual design, we used this opportunity to identify and develop common structures that might be of use to other researchers and applications. Well-documented and parameterizable components that perform functionalities that are common in many algorithms will alleviate the effort spent on new implementations.In this paper, we discuss the interpolation algorithm and our proposed architecture. We emphasize several computational substructures that appear repeatedly throughout the design and which can be implemented effectively using a hardware systolic array-based architecture. These tasks are of the type where we perform a certain reduction or rearrangement of a sequence of elements from multivariate polynomials or boolean expressions. The systolic array concurrently manages data receipt and parallel processing, making it an amenable structure when dealing with streamed data. The simplicity of the array cells, storage, and control unit, allows the instantiation of multiple cells while maintaining competitive clock frequencies, thus achieving high performance. Several tasks critical to interpolation were modeled and implemented to FPGA using the proposed architecture, obtaining speedups up to 172× when compared to a software implementation, while achieving low resource utilization. These implementations were used as components to develop a complete, high-performance multivariate polynomial interpolation methodology on hardware.Paper Outline
Section2 discusses the interpolation methodology, while Section 3 details data representation for the implementation. Section 4 describes the hardware implementation in general. Sections 5 and 6 describe an algorithmic generalization for the reduction tasks and present the proposed hardware architecture. Section 7 explains the rest of the implementation blocks. Section 8 reports the results of FPGA implementations, both of reduction tasks and the complete implementation. Section 9 provides our conclusions.


## 2. Interpolation Methodology

Discrete models, in particular finite field models, have been proposed for regulatory processes such as genetic networks [7] and biochemical networks [13]. In the setting of a finite field model a genetic network can be seen as a finite system whose states are governed by the iterations of a tuple of polynomials (P1,P2,…,Pn), where n is the number of genes.Given a sequence ofk  n-tuples of values from a set F, representing the states of n genes at times t0 through tk-1, the reverse engineering problem seeks to find a function f:Fn→Fn for which f(s0)=s1,    f(s1)=s2,…,f(sk-2)=sk-1.WhenF is a finite field, function f can be represented by a tuple of polynomials (P1,P2,…,Pn) where Pi:Fn→F for each i=1,2,…,n. As pointed out in [8] such polynomials can be found by the following variation of Lagrange's interpolation formula:(1)Pi(x)=∑m=0k-1fi(am)∏j=0,j≠mk-1xljm-ajljmamljm-ajljm,
where each am denotes an n-tuple and aml denotes the lth component of am. Given aj and am, where aj≠am, ljm is the first component in which aj and am differ.In the context of genetic networks modeled by polynomials over a finite field, the resulting polynomials give a sense to the biologist about the interactions between genes. Nevertheless, results from (1) may contain redundant variables even after algebraic simplifications.Example 1.
Let us consider the functionf:ℤ33→ℤ3 which is incompletely specified by Table 1.
Using (1), an interpolating polynomial for f is given by P=x12x2x3+x12x22+x12x3+2x12x2+x1x2+2x2+2x12+2x1+1. This polynomial depends on x1,x2,   and x3 and cannot be further simplified by any algebraic means. However, the polynomial P′=2x2+2x2x3 also interpolates f at the given points but depends only on variables x2 and x3. In this context, x1 is a redundant variable, since we found another interpolating polynomial for f that does not include x1 in its expression.Table 1
Table forf:ℤ33→ℤ3.
x1x2x3f01201120120102121102Redundant variables are undesirable as they introduce complexity into the polynomials without adding information valuable to the biologist. Furthermore, empirical data suggest that genetic networks are sparsely connected [14]. Redundant variables can be identified using Sasao's algorithm for the detection of dependent variables in incompletely specified multiple valued functions [15].


### 2.1. Essential Variables and Bases

For any setS and any positive integer n, we denote each s∈Sn by (s1,s2,…,sn). Let U={x1,x2,…,xn} be a set of n variables. For any set of variables X={xi1,xi2,…,xim}⊂U and any s∈Sn, we define projXs={si1,…,sim}. For instance, if U={x1,x2,x3,x4,x5,x6}, X={x2,x3,x6}, and s=(1,0,3,2,0,1), the projection of s on X,projXs, is {0,3,1}.LetA={0,1,…,a-1} and B={0,1,…,b-1} where a and b are integers,  a,b≥2 and let f:D→B, where D⊂An. We say that a variable xi is essential in f or depends onf if there exist c,d∈D such that ci≠di,cj=dj for all j≠i, and f(c)≠f(d). A set of variables {xi1,xi2,…,xim} is a basis for f if (1) for every c,d∈D, projXc=projXd if and only if f(c)=f(d) and (2) there is no other set of variables that properly contains X with this property.Lemma 1.
Ifx is an essential variable and X is a basis, then x∈X.Proof.
Letx be an essential variable. Then there exists c,d∈D and i such that ci≠di, projU-{xi}c=projU-{xi}d and f(c)≠f(d). Let X be any basis. Then X⊄U-{xi}, but x⊂U. Hence xi∈X.For any functionf:D→B and any basis X,f can be expressed solely in terms of the variables of X. For any basis X, any variable x∉X is redundant (with respect to X.)Our goal is to determine interpolating polynomials over finite fields in terms of the variables of any of its bases. To this end, let us first review an algorithm of Sasao [15] which determines all bases for any multiple valued function, not just those over finite fields.Letf:D→B, where D⊂An. For each i=0,1,…,b-1, let Sf,i={s∈D∣f(s)=i}.The set of variables appearing in each disjunct ofR constitutes a basis.Thus, the algorithm consists of two stages. First, determineR, which requires O(b2) steps and second, convert R to disjunctive normal form, which requires O(2n) steps.The following lemma, whose proof is immediate, is useful in simplifying expressions such asx1∧x1x2=x1 and x1∨x1x2=x1 in the computation of R.Lemma 2.
LetS={xi1,xi2,…,xid} be a subset of the variables {x1,x2,…,xn}, let E be the disjunction of the variables in S and let F be a disjunction of some subset of S. Also let G be the conjunction of some subset of S. Then  (a)
E∧F=F∧E=F,(b)
E∧G=G∧E=G.In what follows we  abbreviateE∧F by EF.Example 2.
Letf:{0,1,2,3}5→{1,2,3} be a function whose values are included  in Table 2.
Applying Algorithm1 and the reduction rules of Lemma 2, we obtain
(2)R=(x2∨x5)(x3)(x1∨x5)(x1∨x4)=x1x2x3∨x1x3x5∨x3x4x5.
Thus, the bases are X1={x1,x2,x3}, X2={x1,x3,x5}, and X3={x3,x4,x5}.Table 2x1x2x3x4x5f021231021311211011011202212102212012212212012113222223Algorithm 1:Algorithm to determine the collection of bases.
for  i:=0 to b-1  doR=1;for  j:=0 to b-1,j≠i  dor(i,j)=⋁(u,v)∈Sf,i×Sf,j{xk|uk≠vk}endR=R∧r(i,j);ExpressR in disjunctive normal form (DNF)end


### 2.2. Interpolating Polynomials with No Redundant Variables

Given any partially defined functionf over a finite field and any basis X for f, we would like to determine interpolating polynomials involving only those variables of X, that is, polynomials containing no redundant variables.Definition 1.
Letf:D→B where D⊂An and let X be a basis for f. For all c,d∈An define the relation c≡Xd if and only if ci=di, for each xi∈X, f(c), and f(d) are defined, and f(c)=f(d).In other words,c≡Xd if and only if (projXc,f(c))=(projXd,f(d)). It is straightforward to verify the following.Lemma 3.
≡X is an equivalence relation.LetD/≡X be a set of representatives of each equivalence class under ≡X.Theorem 1.
LetF be a field, let (a0,b0), (a1,b1), …, (ak-1,bk-1) be a set of points in Fn×F, and let X be a basis for the function f defined by f(ai)=bi, i=0,1,…,k-1. Then an interpolating polynomial for f whose only variables are those of X is
(3)P(X)=∑m=0,am′∈D/≡Xk-1bm∏j≠m(xljm-ajljmamljm′-ajljm),
where ℓjm is the smallest index such that xℓjm∈X and aj′ and am′ differ.Proof.
For eacham, m=0,1,…,k-1, let am′ be the representative of the equivalence class Cm containing am. Then
(4)xljm-ajljmamljm′-ajljm={0ifx=aj1ifx=am
and so P(am)=bm, where bm is the value of each tuple of Cm.Note that for each of the factors(5)cjm(x)=xljm-ajljmamljm′-ajljm
in (3), cjm(aj)=0 and cjm(am)=1 and so for all i=0,1,…,k-1 and all r=1,2,…(6)cjmr(ai)=cjm(ai).Equation (3) gives us a new algorithm for polynomial interpolation that involves only those variables in the chosen basis X, and furthermore, it also avoids redundant tuples by computing only on the representatives of the equivalence classes given by the relation ≡X.Using (3), O(k2) operations are needed to interpolate k points for each function f:Fn→F. Given k points of a function f:Fn→Fn, the bases for the component functions fi need not be the same. In this case, we can apply (3) n times, thus giving an algorithm with complexity O(nk2).Example 3.
Let us assume that the0,1,2,3 of Example 2 represent the corresponding elements of ℤ5. We found that the bases are X1={x1,x2,x3}, X2={x1,x3,x5}, and X3={x3,x4,x5}. Now let us use the new algorithm given by (3) to find an interpolating polynomial over ℤ5 in terms of the basis X1. The equivalence classes under X1 are
(7)C1={(0,2,1,2,3),(0,2,1,3,1)},C2={(2,1,1,0,1)},C3={(0,1,1,2,0)},C4={(2,1,2,1,0),(2,1,2,0,1),(2,1,2,2,1)},C5={(0,1,2,1,1)},C6={(2,2,2,2,2)}.
and the set of class representatives is
(8)D≡X1={(0,2,1,2,3),(2,1,1,0,1),(0,1,1,2,0),(2,1,2,1,0),(0,1,2,1,1),(2,2,2,2,2)}.
Applying (3) and making use of (6), we have
(9)P(X1)=1*x1-20-2x2-12-1+1*x1-02-0x3-21-2x2-21-2+2*x2-21-2x1-20-2x3-21-2+2*x1-02-0x3-12-1x2-21-2+3*x2-21-2x1-20-2x3-12-1+3*x1-02-0x2-12-1=2(x1-2)(x2-1)+3x1(x2-2)(x3-2)+4(x1-2)(x2-2)(x3-2)+4x1(x3-1)(x2-2)+4(x1-2)(x2-2)(x3-1)+4x1(x2-1).
Thus, a polynomial which depends only on the variablesx1, x2, and x3 and interpolates the tuples given by Example 2 is
(10)P(X1)=4x1x2+4x2x3+3x1+2x3+1.


## 2.1. Essential Variables and Bases

For any setS and any positive integer n, we denote each s∈Sn by (s1,s2,…,sn). Let U={x1,x2,…,xn} be a set of n variables. For any set of variables X={xi1,xi2,…,xim}⊂U and any s∈Sn, we define projXs={si1,…,sim}. For instance, if U={x1,x2,x3,x4,x5,x6}, X={x2,x3,x6}, and s=(1,0,3,2,0,1), the projection of s on X,projXs, is {0,3,1}.LetA={0,1,…,a-1} and B={0,1,…,b-1} where a and b are integers,  a,b≥2 and let f:D→B, where D⊂An. We say that a variable xi is essential in f or depends onf if there exist c,d∈D such that ci≠di,cj=dj for all j≠i, and f(c)≠f(d). A set of variables {xi1,xi2,…,xim} is a basis for f if (1) for every c,d∈D, projXc=projXd if and only if f(c)=f(d) and (2) there is no other set of variables that properly contains X with this property.Lemma 1.
Ifx is an essential variable and X is a basis, then x∈X.Proof.
Letx be an essential variable. Then there exists c,d∈D and i such that ci≠di, projU-{xi}c=projU-{xi}d and f(c)≠f(d). Let X be any basis. Then X⊄U-{xi}, but x⊂U. Hence xi∈X.For any functionf:D→B and any basis X,f can be expressed solely in terms of the variables of X. For any basis X, any variable x∉X is redundant (with respect to X.)Our goal is to determine interpolating polynomials over finite fields in terms of the variables of any of its bases. To this end, let us first review an algorithm of Sasao [15] which determines all bases for any multiple valued function, not just those over finite fields.Letf:D→B, where D⊂An. For each i=0,1,…,b-1, let Sf,i={s∈D∣f(s)=i}.The set of variables appearing in each disjunct ofR constitutes a basis.Thus, the algorithm consists of two stages. First, determineR, which requires O(b2) steps and second, convert R to disjunctive normal form, which requires O(2n) steps.The following lemma, whose proof is immediate, is useful in simplifying expressions such asx1∧x1x2=x1 and x1∨x1x2=x1 in the computation of R.Lemma 2.
LetS={xi1,xi2,…,xid} be a subset of the variables {x1,x2,…,xn}, let E be the disjunction of the variables in S and let F be a disjunction of some subset of S. Also let G be the conjunction of some subset of S. Then  (a)
E∧F=F∧E=F,(b)
E∧G=G∧E=G.In what follows we  abbreviateE∧F by EF.Example 2.
Letf:{0,1,2,3}5→{1,2,3} be a function whose values are included  in Table 2.
Applying Algorithm1 and the reduction rules of Lemma 2, we obtain
(2)R=(x2∨x5)(x3)(x1∨x5)(x1∨x4)=x1x2x3∨x1x3x5∨x3x4x5.
Thus, the bases are X1={x1,x2,x3}, X2={x1,x3,x5}, and X3={x3,x4,x5}.Table 2x1x2x3x4x5f021231021311211011011202212102212012212212012113222223Algorithm 1:Algorithm to determine the collection of bases.
for  i:=0 to b-1  doR=1;for  j:=0 to b-1,j≠i  dor(i,j)=⋁(u,v)∈Sf,i×Sf,j{xk|uk≠vk}endR=R∧r(i,j);ExpressR in disjunctive normal form (DNF)end


## 2.2. Interpolating Polynomials with No Redundant Variables

Given any partially defined functionf over a finite field and any basis X for f, we would like to determine interpolating polynomials involving only those variables of X, that is, polynomials containing no redundant variables.Definition 1.
Letf:D→B where D⊂An and let X be a basis for f. For all c,d∈An define the relation c≡Xd if and only if ci=di, for each xi∈X, f(c), and f(d) are defined, and f(c)=f(d).In other words,c≡Xd if and only if (projXc,f(c))=(projXd,f(d)). It is straightforward to verify the following.Lemma 3.
≡X is an equivalence relation.LetD/≡X be a set of representatives of each equivalence class under ≡X.Theorem 1.
LetF be a field, let (a0,b0), (a1,b1), …, (ak-1,bk-1) be a set of points in Fn×F, and let X be a basis for the function f defined by f(ai)=bi, i=0,1,…,k-1. Then an interpolating polynomial for f whose only variables are those of X is
(3)P(X)=∑m=0,am′∈D/≡Xk-1bm∏j≠m(xljm-ajljmamljm′-ajljm),
where ℓjm is the smallest index such that xℓjm∈X and aj′ and am′ differ.Proof.
For eacham, m=0,1,…,k-1, let am′ be the representative of the equivalence class Cm containing am. Then
(4)xljm-ajljmamljm′-ajljm={0ifx=aj1ifx=am
and so P(am)=bm, where bm is the value of each tuple of Cm.Note that for each of the factors(5)cjm(x)=xljm-ajljmamljm′-ajljm
in (3), cjm(aj)=0 and cjm(am)=1 and so for all i=0,1,…,k-1 and all r=1,2,…(6)cjmr(ai)=cjm(ai).Equation (3) gives us a new algorithm for polynomial interpolation that involves only those variables in the chosen basis X, and furthermore, it also avoids redundant tuples by computing only on the representatives of the equivalence classes given by the relation ≡X.Using (3), O(k2) operations are needed to interpolate k points for each function f:Fn→F. Given k points of a function f:Fn→Fn, the bases for the component functions fi need not be the same. In this case, we can apply (3) n times, thus giving an algorithm with complexity O(nk2).Example 3.
Let us assume that the0,1,2,3 of Example 2 represent the corresponding elements of ℤ5. We found that the bases are X1={x1,x2,x3}, X2={x1,x3,x5}, and X3={x3,x4,x5}. Now let us use the new algorithm given by (3) to find an interpolating polynomial over ℤ5 in terms of the basis X1. The equivalence classes under X1 are
(7)C1={(0,2,1,2,3),(0,2,1,3,1)},C2={(2,1,1,0,1)},C3={(0,1,1,2,0)},C4={(2,1,2,1,0),(2,1,2,0,1),(2,1,2,2,1)},C5={(0,1,2,1,1)},C6={(2,2,2,2,2)}.
and the set of class representatives is
(8)D≡X1={(0,2,1,2,3),(2,1,1,0,1),(0,1,1,2,0),(2,1,2,1,0),(0,1,2,1,1),(2,2,2,2,2)}.
Applying (3) and making use of (6), we have
(9)P(X1)=1*x1-20-2x2-12-1+1*x1-02-0x3-21-2x2-21-2+2*x2-21-2x1-20-2x3-21-2+2*x1-02-0x3-12-1x2-21-2+3*x2-21-2x1-20-2x3-12-1+3*x1-02-0x2-12-1=2(x1-2)(x2-1)+3x1(x2-2)(x3-2)+4(x1-2)(x2-2)(x3-2)+4x1(x3-1)(x2-2)+4(x1-2)(x2-2)(x3-1)+4x1(x2-1).
Thus, a polynomial which depends only on the variablesx1, x2, and x3 and interpolates the tuples given by Example 2 is
(10)P(X1)=4x1x2+4x2x3+3x1+2x3+1.


## 3. Data Representation

As can be deduced from the previous section, an implementation of the multivariate polynomial interpolation algorithm must be able to represent and compute on disjunctive normal form (DNF), conjunctive normal form (CNF), and multivariate polynomial expressions. The following subsections establish our method of representation, which ultimately determines the techniques used for implementing the algorithms.


### 3.1. DNF and CNF Expressions

A Boolean expressionD in DNF using only uncomplemented variables is a disjunction of conjunctive clauses D=⋁i=0q-1ci, where ci=x0ti0∧x1ti1∧⋯∧xm-1tim-1, tij∈{0,1}. We shall refer to each ci term as a cube. For example, the following is a Boolean expression in DNF: x1x2∨x1x3∨x0. In our scheme, a DNF expression is represented as a set of cubes {c0,c1,…,cq-1} where each ci is represented as a sequence 〈ti0,ti1,…,tim-1〉. Observe that complemented variables are not used in the algorithm for finding essential variables and bases, hence their representation is not necessary in this context.Example 4.
The DNF expressionD(x0,x1,x2,x3)=x1x2∨x1x3∨x0 is represented as {〈0,1,1,0〉,〈0,1,0,1〉,〈1,0,0,0〉}A Boolean expressionC in CNF using only uncomplemented variables is represented as a set of disjunctions of literals (DOL) {d0,d1,…,dq-1} where each DOL di is represented as a sequence 〈ti0,ti1,…,tim-1〉. Within our implementation there is no need to distinguish the CNF from the DNF representation since no individual component needs to operate on both simultaneously.Example 5.
The CNF expressionC(x0,x1,x2)=(x0∨x2)∧(x0∨x1)∧(x1) is represented as {〈1,0,1〉,〈1,1,0〉,〈0,1,0〉}


### 3.2. Polynomial Representation

A multivariate polynomialP over ℤp={0,1,…,p-1}, where p is prime is defined as P(x0,x1,…,xn-1)=∑i=0q-1mi, where mi=αix0ti0x1ti1⋯xn-1tin-1 and αi,tij∈ℤp. For example, one multivariate polynomial over ℤ5 is P(x0,x1,x2)=4+x03x1+3x2+x24. In our scheme, a polynomial is represented as a set of monomials {m0,m1,…,mq-1} where each mi is represented as a tuple (αi,〈ti0,ti1,…,tin-1〉). Binary representation requires ⌈log2p⌉ bits for each αi and tij.Example 6.
The multivariate polynomialP(x0,x1,x2)=4+x03x1+3x2+x24 over ℤ5 is represented as {(4,〈0,0,0〉),(1,〈3,1,0〉),(3,〈0,0,1〉),(1,〈0,0,4〉)}. For the three variable polynomial over ℤ5 each monomial requires 4×⌈log25⌉=12 bits for binary representation.


## 3.1. DNF and CNF Expressions

A Boolean expressionD in DNF using only uncomplemented variables is a disjunction of conjunctive clauses D=⋁i=0q-1ci, where ci=x0ti0∧x1ti1∧⋯∧xm-1tim-1, tij∈{0,1}. We shall refer to each ci term as a cube. For example, the following is a Boolean expression in DNF: x1x2∨x1x3∨x0. In our scheme, a DNF expression is represented as a set of cubes {c0,c1,…,cq-1} where each ci is represented as a sequence 〈ti0,ti1,…,tim-1〉. Observe that complemented variables are not used in the algorithm for finding essential variables and bases, hence their representation is not necessary in this context.Example 4.
The DNF expressionD(x0,x1,x2,x3)=x1x2∨x1x3∨x0 is represented as {〈0,1,1,0〉,〈0,1,0,1〉,〈1,0,0,0〉}A Boolean expressionC in CNF using only uncomplemented variables is represented as a set of disjunctions of literals (DOL) {d0,d1,…,dq-1} where each DOL di is represented as a sequence 〈ti0,ti1,…,tim-1〉. Within our implementation there is no need to distinguish the CNF from the DNF representation since no individual component needs to operate on both simultaneously.Example 5.
The CNF expressionC(x0,x1,x2)=(x0∨x2)∧(x0∨x1)∧(x1) is represented as {〈1,0,1〉,〈1,1,0〉,〈0,1,0〉}


## 3.2. Polynomial Representation

A multivariate polynomialP over ℤp={0,1,…,p-1}, where p is prime is defined as P(x0,x1,…,xn-1)=∑i=0q-1mi, where mi=αix0ti0x1ti1⋯xn-1tin-1 and αi,tij∈ℤp. For example, one multivariate polynomial over ℤ5 is P(x0,x1,x2)=4+x03x1+3x2+x24. In our scheme, a polynomial is represented as a set of monomials {m0,m1,…,mq-1} where each mi is represented as a tuple (αi,〈ti0,ti1,…,tin-1〉). Binary representation requires ⌈log2p⌉ bits for each αi and tij.Example 6.
The multivariate polynomialP(x0,x1,x2)=4+x03x1+3x2+x24 over ℤ5 is represented as {(4,〈0,0,0〉),(1,〈3,1,0〉),(3,〈0,0,1〉),(1,〈0,0,4〉)}. For the three variable polynomial over ℤ5 each monomial requires 4×⌈log25⌉=12 bits for binary representation.


## 4. General Description of Implementation

Our algorithm for determining interpolating polynomials that do not contain any redundant variables consists of two stages. The first stage identifies the dependent variables using Sasao's algorithm, while the second uses this information along with (3) to determine Pi. Thus, a key part of the reconfigurable logic implementation is choosing efficient hardware structures to perform the critical parts of both stages. For the first stage, the most time-consuming task is the conversion of a Boolean expression in conjunctive normal form (CNF) to disjunctive normal form (DNF). The second step relies heavily on the identification of distinct binomials as they are generated from (3) and multivariate polynomial multiplication/addition.To take advantage of the FPGA's fine-grained parallelism and considering their limitations in I/O bandwidth, these computations were implemented in a pipelined manner. In other words, computational blocks were designed to sustain stream processing as allowed by the FPGA's resources, rather than accumulate all needed data and then perform block processing.Figure1 shows a block diagram of the implementation. The CNF Generator receives the n-tuples and performs Algorithm 1 to generate the disjuntions of literals, that is, the r(i,j) terms. These are streamed to the Cover DOLs block which eliminates redundant disjuntions of literals by using Lemma 2. The nonredundant DOLs are then fed to block CNF2DNF which computes the conjunctions of literals and uses Lemma 2 to eliminate redundant conjunctions. The output of CNF2DNF is the set of bases. A priori biological knowledge about the organism under study is required to choose the most adequate base. The n-tuples are also fed to the second stage which uses the chosen base to determine representatives of the equivalence classes from the n-tuples. The representatives are then streamed to the modified Lagrange interpolation block which implements (3) by generating the binomials for each pair of class representatives (Binomial Generator), eliminates duplicate binomials (Filter Duplicates), then multiplies and adds the product results.Figure 1
Block diagram of the implementation, highlighting the subtasks that are implemented using the systolic array architecture.The highlighted blocks in Figure1, that is, Cover DOLs, Cover COLs, Determine Class Representatives, and Polynomial Addition, have two characteristics that make them suitable candidates for systolic array processing: (1) each performs an operation on a sequence of elements that are streamed from the preceeding block, (2) the operation is such that it performs reduction of a sequence of elements from multivariate polynomials or Boolean expressions. For example, the Cover Cubes receives preliminary results from the CNF to DNF conversion and eliminates redundant cubes using Lemma 2. This is an operation that can be described as given a sequence C of cubes (C=c1,c2,…,cm) determine a subsequence C′⊂C such that it contains only cubes ci∈C where there exists no other cj≠ci∈C such that ci∧cj=cj.The following sections describe a generalized algorithm and systolic array-based architecture for the type of reduction operations common throughout our interpolation methodology. This is followed by a description of how they are utilized as part of the architectural data path.


## 5. Generalization of Reduction Tasks

The subtasks that deal with reduction in our interpolation algorithm can be generally described as follows. Given the sequenceX=x1,x2,…,xm of elements from Σ, compute the sequence Y=y1,y2,…,ym where yi∈Σɛ=Σ∪ɛ by eliminating or combining elements based on binary comparisons between them. The empty element (ɛ) is introduced for the representation of operations where groups of elements can be reduced to a single element. Algorithm 2 shows a general form of these problems, where W and S are binary functions of the form Σɛ×Σɛ→Σɛ, and the symbol ∥ denotes parallel computation.Algorithm 2:Nested loop algorithm for reduction tasks.
for  i:=1 to n-1  dofor  j:=i+1 to n  doxi:=W(xi,xj)∥xj:=S(xi,xj);endendThe tasks of distinctness, polynomial addition, and redundant Boolean term elimination (hereon referred to ascover) can be implemented by specifying W and S, as explained in the following subsections.


### 5.1. Distinctness

The task of identifying the distinct elements of a sequenceX=x1,x2,…,xm is needed within our interpolation algorithm to determine class representatives and filter out binomial duplicates. It can be implemented by Algorithm 2 with(11)W(xi,xj):={xjifxi=ɛxi,otherwise,S(xi,xj):={ɛifxi=xjorxi=ɛxj,otherwise.Example 7.
AssumeX=5,8,8,5,8,1. After computation with Algorithm 2 using (11), the result is X′=5,8,1,ɛ,ɛ,ɛ,ɛ.We provide a correctness proof for this operation. Proofs of the other reduction operations given in this section are similar.Lemma 4.
For an input sequencex=x1,x2,…,xm, where xi∈Σ, the output of Algorithm 2 with functions W and S defined as in (11) outputs xi1,xi2,…,xik,ɛ,ɛ,…,ɛ where xi1,xi2,…,xik is the subsequence of distinct elements of the input sequence.Proof.
We show by induction that afterj≤k iterations of the outer loop,
(12)xi,x2,…,xj,xj+1,…,xm=xi1,xi2,…,xij,yj+1,…,ym,
where for each p≥j+1,yp∈{xij+1,…,xik}∪{ɛ}.
Afterj=1 iterations, x1=xi1. Next suppose that after j iterations,
(13)x1,x2,…,xj,xj+1,…,xm=xi1,…,xij,yj+1,…,ym.
for some j<k where yp∈{xij+1,…,xik}∪{ɛ} for each p≥j+1. Then either (1) yj+1=xj+1 or (2) yj+1=ɛ. In the first case, yj+1=xij+1 since xj+1 is not equal to any xp,p≤j and is not changed by the j+1st iteration. In the second case, yj+1=ɛ will be replaced by xij+1 during the j+1st iteration. Hence in either case,
(14)x1,x2,…,xj,xj+1,…,xm=xi1,xi2,…,xij+1,yj+2,…,ym,
where for each p≥j+2,yp is either ɛ or an element xip, where p≥j+2.
Hence afteri=k iterations,
(15)x1,x2,…,xm=xi1,xi2,…,xik,ɛ,…,ɛ
and remains the same after i>k iterations since W(ɛ,ɛ)=S(ɛ,ɛ)=ɛ.


### 5.2. Polynomial Addition

Assume that each of the elements inX=x1,x2,…,xm is a monomial. A polynomial addition operation can be implemented by Algorithm 2 with(16)W(xi,xj):={xi+xjifmon(xi)=mon(xj),orxi=ɛ,xi,otherwise,S(xi,xj):={ɛifmon(xi)=mon(xj),orxi=ɛ,xj,otherwise,
where mon(xi)=〈ti0,ti1,…,tin-1〉, that is, the monomial without the coefficient.Example 8.
AssumeX=x2,xy2,3x2,xy,2xy. After computation with Algorithm 2 using (16), the result is X′=4x2,xy2,3xy,ɛ,ɛ.


### 5.3. Cover

Assume that each of the elements inX is a cube. A cube ci covers another cube cj (represented as ci⊃cj) if all Boolean variables present in ci are present in cj (e.g., x2⊃x1x2x3, but x2x4⊅x1x2x3). The cover operation can be implemented by Algorithm 2 with(17)W(xi,xj):={xjifxj⊃xi,xi,otherwise,S(xi,xj):={ɛifxi⊃xj,xj,otherwise,
where qi⊃ɛ for any qi.Example 9.
AssumeX=x1x2,x2x4,x1,x1x2x4. After computation with Algorithm 2 using (17), the result is X′=x1,x2x4,x1,ɛ. Notice that there may be duplicate terms in the result but these can be easily eliminated using a distinctness operation.Algorithm2 is also amenable to other problems that can be expressed as the result of binary comparisons between every two elements of a sequence. For example, the problem of sorting a sequence S can be implemented with(18)W(xi,xj):=min⁡⁡(xi,xj),S(xi,xj):=max⁡⁡(xi,xj).


## 5.1. Distinctness

The task of identifying the distinct elements of a sequenceX=x1,x2,…,xm is needed within our interpolation algorithm to determine class representatives and filter out binomial duplicates. It can be implemented by Algorithm 2 with(11)W(xi,xj):={xjifxi=ɛxi,otherwise,S(xi,xj):={ɛifxi=xjorxi=ɛxj,otherwise.Example 7.
AssumeX=5,8,8,5,8,1. After computation with Algorithm 2 using (11), the result is X′=5,8,1,ɛ,ɛ,ɛ,ɛ.We provide a correctness proof for this operation. Proofs of the other reduction operations given in this section are similar.Lemma 4.
For an input sequencex=x1,x2,…,xm, where xi∈Σ, the output of Algorithm 2 with functions W and S defined as in (11) outputs xi1,xi2,…,xik,ɛ,ɛ,…,ɛ where xi1,xi2,…,xik is the subsequence of distinct elements of the input sequence.Proof.
We show by induction that afterj≤k iterations of the outer loop,
(12)xi,x2,…,xj,xj+1,…,xm=xi1,xi2,…,xij,yj+1,…,ym,
where for each p≥j+1,yp∈{xij+1,…,xik}∪{ɛ}.
Afterj=1 iterations, x1=xi1. Next suppose that after j iterations,
(13)x1,x2,…,xj,xj+1,…,xm=xi1,…,xij,yj+1,…,ym.
for some j<k where yp∈{xij+1,…,xik}∪{ɛ} for each p≥j+1. Then either (1) yj+1=xj+1 or (2) yj+1=ɛ. In the first case, yj+1=xij+1 since xj+1 is not equal to any xp,p≤j and is not changed by the j+1st iteration. In the second case, yj+1=ɛ will be replaced by xij+1 during the j+1st iteration. Hence in either case,
(14)x1,x2,…,xj,xj+1,…,xm=xi1,xi2,…,xij+1,yj+2,…,ym,
where for each p≥j+2,yp is either ɛ or an element xip, where p≥j+2.
Hence afteri=k iterations,
(15)x1,x2,…,xm=xi1,xi2,…,xik,ɛ,…,ɛ
and remains the same after i>k iterations since W(ɛ,ɛ)=S(ɛ,ɛ)=ɛ.


## 5.2. Polynomial Addition

Assume that each of the elements inX=x1,x2,…,xm is a monomial. A polynomial addition operation can be implemented by Algorithm 2 with(16)W(xi,xj):={xi+xjifmon(xi)=mon(xj),orxi=ɛ,xi,otherwise,S(xi,xj):={ɛifmon(xi)=mon(xj),orxi=ɛ,xj,otherwise,
where mon(xi)=〈ti0,ti1,…,tin-1〉, that is, the monomial without the coefficient.Example 8.
AssumeX=x2,xy2,3x2,xy,2xy. After computation with Algorithm 2 using (16), the result is X′=4x2,xy2,3xy,ɛ,ɛ.


## 5.3. Cover

Assume that each of the elements inX is a cube. A cube ci covers another cube cj (represented as ci⊃cj) if all Boolean variables present in ci are present in cj (e.g., x2⊃x1x2x3, but x2x4⊅x1x2x3). The cover operation can be implemented by Algorithm 2 with(17)W(xi,xj):={xjifxj⊃xi,xi,otherwise,S(xi,xj):={ɛifxi⊃xj,xj,otherwise,
where qi⊃ɛ for any qi.Example 9.
AssumeX=x1x2,x2x4,x1,x1x2x4. After computation with Algorithm 2 using (17), the result is X′=x1,x2x4,x1,ɛ. Notice that there may be duplicate terms in the result but these can be easily eliminated using a distinctness operation.Algorithm2 is also amenable to other problems that can be expressed as the result of binary comparisons between every two elements of a sequence. For example, the problem of sorting a sequence S can be implemented with(18)W(xi,xj):=min⁡⁡(xi,xj),S(xi,xj):=max⁡⁡(xi,xj).


## 6. Hardware Structure for Reduction Tasks

This section explains our proposed hardware structure and outlines the mapping for the reduction tasks. We first discuss the systolic array and cells with the assumption that the array is deep enough to process the input sequence without overflowing. We proceed by discussing the additional components that must be added to guarantee correct processing even if the overflow condition does not hold.


### 6.1. Systolic Array and Cells

The proposed structure is a linear systolic array ofn computational cells, each performing an operation deduced from W and S of the above discussion. Figure 2 illustrates the array and contents of the basic cell. Each cell contains two registers F and M, where each has enough resources to hold any single element from the alphabet Σɛ. Each cell also contains the logic to perform functions W(F,M) and S(F,M). All cells are initialized to the empty symbol. At each clock step a new element from the sequence X is input to the first cell and every cell performs a comparison between its stored value and the value from the previous cell and assigns (possibly) new values to the registers according to W and S. After all elements have been processed the M registers contain the resulting sequence and can be shifted out of the array using the connections between the M registers.Figure 2
Block diagram of systolic array.More formally, our linear systolic array computation can be modeled as a simplified version of the generic systolic array presented in [16], as shown in Algorithm 3. The notation Fi[r] symbolizes the content of register F in cell r at iteration i.Algorithm 3:Model of the linear systolic array computation.
i:=0∥r=1n{M0[r]:=ɛF0[r]:=ɛwhile stop_condition = false doi:=i+1if  i<m  then  Fi-1[0]:=xielse  Fi-1[0]:=ɛend  if∥r=1n{Fi[r]:=S(Fi-1[r-1],Mi-1[r])Mi[r]:=W(Fi-1[r-1],Mi-1[r])endFigure3 illustrates the operation of the array using cells that implement the distinctness task by defining W and S as in (11). Observe that the operations implied by W and S and the connections between neighboring cells accomplish the comparisons established in Algorithm 2, albeit in a different order and in a parallel manner. In fact, the computation implemented by the systolic array can be interpreted as a reindexing of Algorithm 2 as presented in Algorithm 4. Table 3 shows the comparison indexes generated by Algorithm 4 for the example in Figure 3. Notice that even though Algorithm 4 has reordered the indexes, the same data dependence is maintained as in Algorithm 2. In other words, if f(xi,xj) is performed before a f(xi,xk) (where j≠k) in Algorithm 2 the same order is preserved in Algorithm 4. The same condition holds for the order of f(xi,xj) and f(xk,xj).Table 3
Illustration of index generation by Algorithm4 for the example shown in Figure 3.
j=1j=2j=3j=4j=5j=6j=7j=8j=9x0, x1x0, x2x0, x3x0, x4x0, x5x1, x5x2, x5x3, x5x4, x5x1, x2x1, x3x1, x4x2, x4x3, x4x2, x3Algorithm 4:Reindexed algorithm for reduction tasks.
for  j:=1 to 2n-3  doparallel  for  i:= iif (j≤n,1,j-n) to |(j)/2|  doxi=W(xi,xj-i+1)∥xj-i+1=S(xi,xj-i+1)endEndFigure 3
Systolic array implementing the distinctness operation forX=0,6,2,2,7,6. Each row represents a cell. Within each cell, the left and right columns represent the S and M registers, respectively. Elements are input to the top row and flow downwards.


### 6.2. Processing Sequences Longer than the Array Depth

The functionality of the basic cells mandates their implementation on user logic inside the FPGA (rather than on embedded units, such as block RAMS). Depending on the application, the characteristics of the data sequence, and FPGA model, the FPGA might not have enough resources to instantiate a systolic array whose depthdsa is greater than or equal to m (the sequence length). In such cases, there is no guarantee that the array by itself will be able to compare each element against each other to obtain the correct solution. To illustrate this anomaly, assume that in the example presented in Figure 3 the input is a sequence of distinct elements X=x0,x1,…,x5. The array would be filled once x3 is registered in the last cell, hence none of the cells would have operated on the combination of x4 and x5. A control and temporary storage mechanism must be provided to reiterate the data through the pipeline if needed and guarantee that all elements are compared against each other. This can be accomplished by a scheme such as the one shown in Figure 4 through the addition of a FIFO, control unit and several datapath control elements (i.e., multiplexers).Figure 4
Pipeline with added overflow FIFO and control unit to support sequences longer thann.TheOverflow FIFO (OF) stores elements that make it through the array without having been eliminated or registered. These elements will have to be reiterated through the array to test their relation to other elements in the OF. For example, suppose we have an array of depth dsa=4 that implements the distinctness operation. The sequence 〈3,5,1,3,2,6,7,6〉 will fill the pipeline with 〈3,5,1,2〉 so that the final 〈6,7,6〉 overflows the pipeline. These 3 elements are stored in the OF for a later pass through the array. The FIFO depth df≥m-dsa where m is the number of elements in the sequence and dsa is the array depth. Although both systolic array cells and FIFO spend mostly on register resources, a reason why one may be able to increase FIFO depth and not necessarily array depth is that in FPGAs FIFOs are easily mapped to the embedded block RAM, whereas processing cells map to user logic.Thecontrol unit controls the data flow between the array and the OF and determines how to reiterate the overflowed data through the array. For the reduction tasks presented in Section 5 we can deduce three basic modes of operation for the control unit. Figures 5 and 6 illustrate the need for these modes. Figure 5 shows the end of a first pass of a sequence through systolic arrays that implement (a) a sorting algorithm, and (b) sum of polynomials.Contents of the systolic array and OF after a first pass of the sequences through (a) sorting and (b) polynomial addition.
(a)(b)Figure 6
Contents of systolic array and overflow FIFO for the shown sequenceX (a) after first pass, (b) after second pass (OF to array refeed), (c) systolic array is emptied for new OF refeed, (d) after final refeed from OF. Dotted lines illustrate the paths of data refeed.(i)
In the case of a nonreducing operation like sorting (Figure5(a)), the array will contain dsa sorted elements. Thus, elements in the array can be shifted out while elements in the OF are reinputted into the array. This must be repeated until the OF is empty after a pass, that is, at most ⌈m/dsa⌉ times.(ii)
When adding polynomials (Figure5(b)), the array will contain at most dsa monomials with their final coefficients. Thus, elements in the array can be shifted out while elements in the OF are reinputted into the array. This must be repeated until the OF is empty after a pass, that is, at most ⌈m/dsa⌉ times.Figure6 illustrates the various passes needed for the cover operation. At the end of a first pass the array will contain at most dsa cubes. However, in this case a first pass through the array does not guarantee that all cubes in the OF have been compared to all the cubes registered in the array. For example, in Figure 6(a) cube x1x6 passed through the array without being registered. Meanwhile, the last cube x1 covered x1x3 and x1x2 in two cells of the array. Thus, we have at least one cube (x1x6) in the OF which was not compared to a cube in the array (x1). The elements in the OF must be refed through the (unflushed) array to allow every ci stored in the FIFO to be compared against every cj registered in the array. The refeed is repeated until a complete pass of the OF cubes through the array does not modify the cubes registered in the array. After this (Figure 6(b)), if any ci in the OF is going to be part of the result, for example, x4, it must only be still compared to other elements in the OF. The contents of the array may be shifted out to a memory that stores the results, the array is reset and the OF cubes are refed (Figure 6(c)). The three-step process continues until the OF is empty after an OF refeed.


## 6.1. Systolic Array and Cells

The proposed structure is a linear systolic array ofn computational cells, each performing an operation deduced from W and S of the above discussion. Figure 2 illustrates the array and contents of the basic cell. Each cell contains two registers F and M, where each has enough resources to hold any single element from the alphabet Σɛ. Each cell also contains the logic to perform functions W(F,M) and S(F,M). All cells are initialized to the empty symbol. At each clock step a new element from the sequence X is input to the first cell and every cell performs a comparison between its stored value and the value from the previous cell and assigns (possibly) new values to the registers according to W and S. After all elements have been processed the M registers contain the resulting sequence and can be shifted out of the array using the connections between the M registers.Figure 2
Block diagram of systolic array.More formally, our linear systolic array computation can be modeled as a simplified version of the generic systolic array presented in [16], as shown in Algorithm 3. The notation Fi[r] symbolizes the content of register F in cell r at iteration i.Algorithm 3:Model of the linear systolic array computation.
i:=0∥r=1n{M0[r]:=ɛF0[r]:=ɛwhile stop_condition = false doi:=i+1if  i<m  then  Fi-1[0]:=xielse  Fi-1[0]:=ɛend  if∥r=1n{Fi[r]:=S(Fi-1[r-1],Mi-1[r])Mi[r]:=W(Fi-1[r-1],Mi-1[r])endFigure3 illustrates the operation of the array using cells that implement the distinctness task by defining W and S as in (11). Observe that the operations implied by W and S and the connections between neighboring cells accomplish the comparisons established in Algorithm 2, albeit in a different order and in a parallel manner. In fact, the computation implemented by the systolic array can be interpreted as a reindexing of Algorithm 2 as presented in Algorithm 4. Table 3 shows the comparison indexes generated by Algorithm 4 for the example in Figure 3. Notice that even though Algorithm 4 has reordered the indexes, the same data dependence is maintained as in Algorithm 2. In other words, if f(xi,xj) is performed before a f(xi,xk) (where j≠k) in Algorithm 2 the same order is preserved in Algorithm 4. The same condition holds for the order of f(xi,xj) and f(xk,xj).Table 3
Illustration of index generation by Algorithm4 for the example shown in Figure 3.
j=1j=2j=3j=4j=5j=6j=7j=8j=9x0, x1x0, x2x0, x3x0, x4x0, x5x1, x5x2, x5x3, x5x4, x5x1, x2x1, x3x1, x4x2, x4x3, x4x2, x3Algorithm 4:Reindexed algorithm for reduction tasks.
for  j:=1 to 2n-3  doparallel  for  i:= iif (j≤n,1,j-n) to |(j)/2|  doxi=W(xi,xj-i+1)∥xj-i+1=S(xi,xj-i+1)endEndFigure 3
Systolic array implementing the distinctness operation forX=0,6,2,2,7,6. Each row represents a cell. Within each cell, the left and right columns represent the S and M registers, respectively. Elements are input to the top row and flow downwards.


## 6.2. Processing Sequences Longer than the Array Depth

The functionality of the basic cells mandates their implementation on user logic inside the FPGA (rather than on embedded units, such as block RAMS). Depending on the application, the characteristics of the data sequence, and FPGA model, the FPGA might not have enough resources to instantiate a systolic array whose depthdsa is greater than or equal to m (the sequence length). In such cases, there is no guarantee that the array by itself will be able to compare each element against each other to obtain the correct solution. To illustrate this anomaly, assume that in the example presented in Figure 3 the input is a sequence of distinct elements X=x0,x1,…,x5. The array would be filled once x3 is registered in the last cell, hence none of the cells would have operated on the combination of x4 and x5. A control and temporary storage mechanism must be provided to reiterate the data through the pipeline if needed and guarantee that all elements are compared against each other. This can be accomplished by a scheme such as the one shown in Figure 4 through the addition of a FIFO, control unit and several datapath control elements (i.e., multiplexers).Figure 4
Pipeline with added overflow FIFO and control unit to support sequences longer thann.TheOverflow FIFO (OF) stores elements that make it through the array without having been eliminated or registered. These elements will have to be reiterated through the array to test their relation to other elements in the OF. For example, suppose we have an array of depth dsa=4 that implements the distinctness operation. The sequence 〈3,5,1,3,2,6,7,6〉 will fill the pipeline with 〈3,5,1,2〉 so that the final 〈6,7,6〉 overflows the pipeline. These 3 elements are stored in the OF for a later pass through the array. The FIFO depth df≥m-dsa where m is the number of elements in the sequence and dsa is the array depth. Although both systolic array cells and FIFO spend mostly on register resources, a reason why one may be able to increase FIFO depth and not necessarily array depth is that in FPGAs FIFOs are easily mapped to the embedded block RAM, whereas processing cells map to user logic.Thecontrol unit controls the data flow between the array and the OF and determines how to reiterate the overflowed data through the array. For the reduction tasks presented in Section 5 we can deduce three basic modes of operation for the control unit. Figures 5 and 6 illustrate the need for these modes. Figure 5 shows the end of a first pass of a sequence through systolic arrays that implement (a) a sorting algorithm, and (b) sum of polynomials.Contents of the systolic array and OF after a first pass of the sequences through (a) sorting and (b) polynomial addition.
(a)(b)Figure 6
Contents of systolic array and overflow FIFO for the shown sequenceX (a) after first pass, (b) after second pass (OF to array refeed), (c) systolic array is emptied for new OF refeed, (d) after final refeed from OF. Dotted lines illustrate the paths of data refeed.(i)
In the case of a nonreducing operation like sorting (Figure5(a)), the array will contain dsa sorted elements. Thus, elements in the array can be shifted out while elements in the OF are reinputted into the array. This must be repeated until the OF is empty after a pass, that is, at most ⌈m/dsa⌉ times.(ii)
When adding polynomials (Figure5(b)), the array will contain at most dsa monomials with their final coefficients. Thus, elements in the array can be shifted out while elements in the OF are reinputted into the array. This must be repeated until the OF is empty after a pass, that is, at most ⌈m/dsa⌉ times.Figure6 illustrates the various passes needed for the cover operation. At the end of a first pass the array will contain at most dsa cubes. However, in this case a first pass through the array does not guarantee that all cubes in the OF have been compared to all the cubes registered in the array. For example, in Figure 6(a) cube x1x6 passed through the array without being registered. Meanwhile, the last cube x1 covered x1x3 and x1x2 in two cells of the array. Thus, we have at least one cube (x1x6) in the OF which was not compared to a cube in the array (x1). The elements in the OF must be refed through the (unflushed) array to allow every ci stored in the FIFO to be compared against every cj registered in the array. The refeed is repeated until a complete pass of the OF cubes through the array does not modify the cubes registered in the array. After this (Figure 6(b)), if any ci in the OF is going to be part of the result, for example, x4, it must only be still compared to other elements in the OF. The contents of the array may be shifted out to a memory that stores the results, the array is reset and the OF cubes are refed (Figure 6(c)). The three-step process continues until the OF is empty after an OF refeed.


## 7. Further Implementation Considerations

Aside from the reduction tasks, the computational blocks in our implementation can be classified as performing one of two functionalities: element pair generation or multiple term multiplication. Element pair generators, such as theCNF Generator and the Binomial Generator in Figure 1, input a set of tuples and output a distinct pair of tuples at each computational step. This is accomplished by using the architecture illustrated in Figure 7. As tuples are input (each accompanied by a load signal) both RAMs store the tuples and an Address Control Generation unit counts the number of tuples. When a status signal indicates that all tuples have been received, the ACG begins generating all possible address combinations (a0,a1) for 0≤a0<n, a0+1≤a1≤n. The distinct pairs of n-tuples produced by the CNF Generator are fed to n  ⌈log2p⌉-comparators to compute the r(i,j) terms in Algorithm 1, as shown Figure 8. In the case of the Binomial Generator, the corresponding terms within the distinct pairs of n-tuples are subtracted from each other and the results are used to determine the terms of the binomial for (3), as shown in Figure 9.Figure 7
Element pair generator.Figure 8
Computation of the disjunction of literalsr(i,j) from a pair of n-tuples.Figure 9
Computation of the terms of a binomial for (3) from a pair of n-tuples.Multiple-term multiplication/conjunction, such as needed in blocksCNF2DNF and Lagrange, is performed using the architectures shown in Figures 10 and 11, respectively. The architecture depicted in Figure 10 receives DOLs from the CNF Generator, the DOL2Literals block generates a literal expression for each found in the DOL. The literals are queued in a circular FIFO and are conjuncted to each of the results in the Preliminary Results FIFO. After each DOL is conjuncted, the resulting cubes are passed through the Cover Cubes block which eliminates redundant cubes. A control unit monitors and generates signals to coordinate the data path activities.Figure 10
Illustration of a step in the conversion of CNF(x1∨x2∨x3)∧(x1∨x4)∧(x3∨x5)∧(x1∨x2) to DNF. The first two conjunctions have been computed and their redundant terms eliminated and stored in the Preliminary Results FIFO.Figure 11
Illustration of a step in the multiplication of binomials(x2+1)(2x3+3)(3x1+2) in ℤ5. The multiplication of the first two binomials has completed and the results are stored in the Preliminary Results FIFO.The architecture depicted in Figure11 performs the computations of (3). Each binomial that is received from the Binomial Generator is registered. Then each of its two monomials is multiplied by each of the monomials that has resulted from previous binomial multiplications. The monomials that result from the multiplication are input to the PrelimPolyAdd block which adds similar monomials to reduce the size of the preliminary product. PrelimPolyAdd outputs its result to the PrelimRes FIFO. The iterative process continues until the last binomial of a product is received, in which case each monomial of the preliminary product is output to the final Polynomial Addition block (see Figure 1). Once the last binomial of the last product has been multiplied, the (final) Polynomial Addition block is signaled and begins to output the monomials of the final result.


## 8. Results and Discussion

This section presents and discusses results for the individual reduction tasks as well as the complete interpolation algorithm implementation.


### 8.1. Reduction Tasks

The systolic array cells for the subtasks of distinctness, polynomial addition, and Boolean cover were modeled using Verilog HDL, based on their function definitions presented in Section5. The array and control units were also modeled using Verilog based on the descriptions provided in Section 6. Cell, systolic array and control unit models were designed as parameterizable modules in terms of width of elements, array depth and control mode, respectively. Models were behaviorally simulated using ModelSim SE 6.0 to validate their results against software versions of the algorithms. Xilinx ISE Design Suite 11.1 was used to synthesize the architectures for a Virtex-4 XC4VLX200ff1513-11 in the DRC Dev Systems DS2000. Default synthesis parameters were used, such as mapping FIFOs to block RAMs and speed as the optimization goal. All the reported results for the individual tasks are before Place and Route.For each of the three subtasks, several systolic array depths (dsa) were implemented, then randomly generated sequences of various bit widths and lengths were input. Table 2 shows timing and resource results for (1) the cover operation on sequences of length 4096 and width of 32 bits, (2) the polynomial addition operation on sequences of length 4096, where each element represents a monomial of 8 variables in ℤ5, and (3) the distinctness operation on the same sequences as (3). Results are compared against a software implementation compiled using the GNU project c++ compiler with optimization level 3 and running on a 3.40 GHz Intel Pentium D CPU with 2 MB cache and 3 GB RAM. CPU execution times were measured by accessing the processor’s cycle counter. In all cases, we report the smallest measured CPU execution time. The freq column indicates the FPGA clock frequency, tF and tC are the FPGA and CPU run times, and speedup is computed as tC/tF. The numbers shown in parenthesis next to each resource quantity are the resource percentages for the target device.Several important observations can be highlighted from the results in Table4. First, the maximum operation frequency is independent of the array depth, allowing performance to scale adequately with depth. In fact, given the simplicity of the required control unit and storage, the longest path was mostly determined by the implementation of functions W and S inside the basic cells. Second, the parallelism achieved through the use of the systolic array greatly compensates for the refeeds that must be performed when the result is longer than the array depth. For instance, the result for the cover operation was of length 3386 on average, yet there is a speedup of more than 40× for dsa=256. Finally, the resource utilization to achieve competitive speedups is minimal and easy to estimate when scaling since it is almost completely attributable to the systolic array cells and overflow FIFO. Hence when using this structure as part of a bigger design, the designer could readily determine which parameters are best suited for the area/performance constraints and objectives.Table 4
Experimental results for subtask blocks.
dsaSlices*Slice F/F*4 input LUTs*FIFO16/RAMB16*freq (Mhz)tF (s)tC (s)speedup3*Cover644141 (4.65%)4348 (2.44%)7708 (4.33%)9 (2.68%)1761.12E-031.29E-0211.51×1288279 (9.29%)8700 (4.88%)15054 (8.45%)9 (2.68%)1765.90E-041.29E-0221.93×25616554 (18.58%)17404 (9.77%)30796 (17.28%)9 (2.68%)1763.21E-041.29E-0240.33×3*P-add642967 (3.33%)3516 (1.97%)4119 (2.31%)7 (2.08%)2006.84E-043.36E-0249.17×1285885 (6.61%)7455 (4.18%)7156 (4.02%)7 (2.08%)2003.58E-043.36E-0293.99×25611741 (13.18%)14880 (8.35%)14228 (7.99%)7 (2.08%)2001.95E-043.36E-02172.28×3*Distinct642680 (3.01%)3741 (2.10%)3173 (1.78%)7 (2.08%)1907.20E-043.37E-0246.87×1285342 (6.00%)7465 (4.19%)6261 (3.51%)7 (2.08%)1903.77E-043.37E-0289.60×25610659 (11.96%)14903 (8.36%)12437 (6.98%)7 (2.08%)1902.05E-043.37E-02164.25×R*esource percentages are for a Virtex-4 XC4VLX200ff1513-11.Part of our purpose here is to argue that the reduction tasks can be implemented efficiently using the systolic array design. With Table4, we intend to demonstrate the effect of each particular operation for users that might want to use any of them as part of a bigger design. Clearly, the ultimate speedup may be dictated by other components of the final design, as shown by Tables 5 and 6.Table 5
Experimental results for stage 1.
ndSASlicesSlice FFs*4-input LUTS*FIFO16/ RAMB16*tF (ms)tC (ms)Speedup8321651 (1.85%)1908 (1.07%)3047 (1.71%)25 (7.44%)0.041.2729.59×9643011 (3.38%)3596 (2.02%)5515 (3.1%)27 (8.04%)0.192.6413.87×101285860 (6.58%)7201 (4.04%)11186 (6.28%)33 (9.82%)0.497.2214.78×111286191 (6.95%)7728 (4.34%)11755 (6.60%)43 (12.80%)1.0122.9722.74×1225611880 (13.34%)14907 (8.37%)22573 (12.67%)55 (16.37%)2.5482.2832.46×1325614101 (15.83%)16981 (9.53%)25676 (14.41%)95 (28.27%)8.89596.967.13×R*esource percentages are for a Virtex-4 XC4VLX200ff1513-11.Table 6
Experimental results for stage 2.
nBasis VarsClass RepsSlices*Slice FFs*4-input LUTs*FIFO16/ RAMB16*tF (ms)tC (ms)Speedup82237745 (8.69%)8335 (4.68%)13846 (7.77%)28 (8.33%)0.101.5916.67×83637950 (8.92%)8352 (4.69%)14205 (7.97%)28 (8.33%)0.447.7617.68×84907955 (8.93%)8359 (4.69%)14213 (7.98%)30 (8.93%)0.8416.6019.71×92258520 (9.56%)9220 (5.17%)15205 (8.53%)28 (8.33%)0.101.8517.64×93778539 (9.58%)9226 (5.18%)15239 (8.55%)28 (8.33%)0.6911.2016.13×94948582 (9.63%)9260 (5.2%)15340 (8.61%)34 (10.12%)1.1420.2617.81×102249188 (10.31%)10098 (5.67%)16710 (9.38%)31 (9.23%)0.101.7817.95×103659192 (10.32%)10106 (5.67%)16714 (9.38%)31 (9.23%)0.478.3017.69×104929201 (10.33%)10130 (5.69%)16742 (9.40%)33 (9.82%)0.8618.4321.36×1122510130 (11.37%)10978 (6.16%)18024 (10.12%)34 (10.12%)0.112.0118.80×1137510156 (11.4%)11001 (6.17%)18067 (10.14%)36 (10.71%)0.5812.8022.22×1149210115 (11.35%)11008 (6.18%)18025 (10.12%)37 (11.01%)1.0523.0021.93×1222510861 (12.19%)11861 (6.66%)19077 (10.71%)36 (10.71%)0.102.0519.93×1237010806 (12.13%)11881 (6.67%)19012 (10.67%)38 (11.31%)0.5210.6820.41×1249310877 (12.21%)11895 (6.68%)19121 (10.73%)38 (11.31%)1.0319.9019.26×1322511405 (12.8%)12747 (7.15%)20240 (11.36%)39 (11.61%)0.102.1320.74×1337111469 (12.87%)12760 (7.16%)20358 (11.43%)40 (11.90%)0.5312.5023.68×1349111416 (12.81%)12771 (7.17%)20274 (11.38%)41 (12.20%)0.8521.3825.21×R*esource percentages are for a Virtex-4 XC4VLX200ff1513-11.


### 8.2. Interpolation Algorithm

We modeled the proposed interpolation algorithm in Verilog HDL using the reduction tasks and the architectures described in previous sections as building blocks. The behavioral simulation and synthesis tools as well as the software compilation parameters and platform were the same as in the reduction task experiments.The results are shown separately for each of the two interpolation stages since, once the bases have been computed by the first stage, a scientist's intervention would be necessary to choose the most appropriate before proceeding to the second stage. For both stages, the synthesis tool determined maximum operating frequencies above 150 MHz, thus we chose 150 MHz for the experiments. The Xilinx ISE 11.1 place and route tools (with default effort level) were able to meet the timing requirements reported for Stages 1 and 2. We attribute this to the fact that the great majority of connections in our design are local and grow only linearly withn.


#### 8.2.1. Stage 1

Table5 reports timing and resource results for the first stage. Randomly generated sets of 100 n-tuples (8≤n≤13) were input to implementations with diverse systolic array depths. The table shows results for the dSA case with the shortest execution time for each n. The smaller values of n generated fewer intermediate terms during computation, hence smaller values of dSA produced better results since they incurred in less unnecessary overhead for filling and emptying. As the number of variables increases, the number of terms produced by Algorithm 1 increases exponentially, which is clearly evidenced by the run times, that is, columns tF and tC. The FPGA implementation achieves speedups of 13× to 67×, which mostly increase with n. The increase in speedup at the higher values of n can be attributed to a higher percentage of time of full pipeline utilization as well as the increased benefit of parallelism in the comparisons between the n-tuple variables, that is, the computation of the r(i,j) terms in Algorithm 1. As expected, the use of necessarily serial processing tasks such as the multiplication of the DOLs to obtain the DNF (see Example 2) limit the speedup achieved by the FPGA implementation when compared to the results obtained using the cover subtask by itself.Although the ultimate goal of reverse engineering might be understanding complete complex biological systems, recent reported models limit their exploration over specific subsystems or punctual mechanisms, for example, the study cell-cycle regulatory network of fission yeast and apoptosis (programmed cell death) [17, 18], using a moderate number of genes between 8 and 20. Our results lead us to believe that the advantages of using the proposed architecture in reconfigurable logic will be sustained even at bigger sizes. As the sizes increase our architecture will have to rely increasingly on iterations rather than parallel computations. Nevertheless, due to the custom nature (finite field) of the data, even at these sizes the parallelism over each iteration will be enough to offer significant speedups over the alternative, fully serial implementation in general purpose processors.Speedup is maintained at a cost of increasing resource utilization. For the considered cases, BRAMs, which are used to implement the various FIFOs and RAMs of stage 1, is the fastest growing resource. This could pose a challenge to maintain performance at even highern values. However, this issue can potentially be resolved by complementing the BRAM capacity with a dedicated external memory bank, as is commonly included in modern reconfigurable computer platforms [19, 20]. A study of the impact of external memory on the performance of the proposed architecture is future work. An advantage of our proposed architecture is that all the required memories behave as FIFOs. Thus, their accesses patterns are predictable and (as part of future work) we can use this information to device a buffering scheme that hides the latency of accessing the external memories.The increase in user logic (Slices, LUTs) in highern can be circumvented by using shallower systolic arrays (smaller dSA), which still maintain competitive speedups. For example, although not shown in the table, for n=13 speedups of 30× and 45× were obtained by using dSA values of 64 and 128.


#### 8.2.2. Stage 2

As explained in Section2.2, our algorithm uses a base Xi computed in Stage 1 to determine a set of equivalence classes D/≡Xi from the original set of n-tuples. The equivalence classes are then used to compute an interpolating polynomial using (3). A greater number of variables in the chosen Xi and a high variability among the n-tuples translates to larger cardinalities of D/≡Xi. The greater the cardinality of D/≡Xi the higher the number of binomials computed by (3), which constitutes the bulk of computation for Stage 2.The timing results shown in Table6 support the previous statements. The sets of 100 n-tuples from Stage 1 along with bases of 2 to 4 variables were input to Stage 2. Observe that since the data sets were randomly generated there was high variability in the sets of n-tuples, which resulted in similar D/≡Xi cardinalities (column Class Reps) and run times (columns tF and tC) for the various sizes of n (for a given number of variables in the chosen basis). The FPGA implementation achieves speedups of 16× to 25×, which mostly increase with n. The moderate increment in speedup as n increases is attributed to the increased benefit of parallelism for the generation of binomials, and the comparison and addition of monomials.Resource utilization distribution in Stage 2 is roughly uniform between the user logic, that is, Slices, and BRAMs. Furthermore, the increment in resource utilization for increasing values ofn and numbers of variables in the basis is quite moderate, for example, the percentage of slices goes from 8.69% for n=8 to 12.87% for n=13.


## 8.1. Reduction Tasks

The systolic array cells for the subtasks of distinctness, polynomial addition, and Boolean cover were modeled using Verilog HDL, based on their function definitions presented in Section5. The array and control units were also modeled using Verilog based on the descriptions provided in Section 6. Cell, systolic array and control unit models were designed as parameterizable modules in terms of width of elements, array depth and control mode, respectively. Models were behaviorally simulated using ModelSim SE 6.0 to validate their results against software versions of the algorithms. Xilinx ISE Design Suite 11.1 was used to synthesize the architectures for a Virtex-4 XC4VLX200ff1513-11 in the DRC Dev Systems DS2000. Default synthesis parameters were used, such as mapping FIFOs to block RAMs and speed as the optimization goal. All the reported results for the individual tasks are before Place and Route.For each of the three subtasks, several systolic array depths (dsa) were implemented, then randomly generated sequences of various bit widths and lengths were input. Table 2 shows timing and resource results for (1) the cover operation on sequences of length 4096 and width of 32 bits, (2) the polynomial addition operation on sequences of length 4096, where each element represents a monomial of 8 variables in ℤ5, and (3) the distinctness operation on the same sequences as (3). Results are compared against a software implementation compiled using the GNU project c++ compiler with optimization level 3 and running on a 3.40 GHz Intel Pentium D CPU with 2 MB cache and 3 GB RAM. CPU execution times were measured by accessing the processor’s cycle counter. In all cases, we report the smallest measured CPU execution time. The freq column indicates the FPGA clock frequency, tF and tC are the FPGA and CPU run times, and speedup is computed as tC/tF. The numbers shown in parenthesis next to each resource quantity are the resource percentages for the target device.Several important observations can be highlighted from the results in Table4. First, the maximum operation frequency is independent of the array depth, allowing performance to scale adequately with depth. In fact, given the simplicity of the required control unit and storage, the longest path was mostly determined by the implementation of functions W and S inside the basic cells. Second, the parallelism achieved through the use of the systolic array greatly compensates for the refeeds that must be performed when the result is longer than the array depth. For instance, the result for the cover operation was of length 3386 on average, yet there is a speedup of more than 40× for dsa=256. Finally, the resource utilization to achieve competitive speedups is minimal and easy to estimate when scaling since it is almost completely attributable to the systolic array cells and overflow FIFO. Hence when using this structure as part of a bigger design, the designer could readily determine which parameters are best suited for the area/performance constraints and objectives.Table 4
Experimental results for subtask blocks.
dsaSlices*Slice F/F*4 input LUTs*FIFO16/RAMB16*freq (Mhz)tF (s)tC (s)speedup3*Cover644141 (4.65%)4348 (2.44%)7708 (4.33%)9 (2.68%)1761.12E-031.29E-0211.51×1288279 (9.29%)8700 (4.88%)15054 (8.45%)9 (2.68%)1765.90E-041.29E-0221.93×25616554 (18.58%)17404 (9.77%)30796 (17.28%)9 (2.68%)1763.21E-041.29E-0240.33×3*P-add642967 (3.33%)3516 (1.97%)4119 (2.31%)7 (2.08%)2006.84E-043.36E-0249.17×1285885 (6.61%)7455 (4.18%)7156 (4.02%)7 (2.08%)2003.58E-043.36E-0293.99×25611741 (13.18%)14880 (8.35%)14228 (7.99%)7 (2.08%)2001.95E-043.36E-02172.28×3*Distinct642680 (3.01%)3741 (2.10%)3173 (1.78%)7 (2.08%)1907.20E-043.37E-0246.87×1285342 (6.00%)7465 (4.19%)6261 (3.51%)7 (2.08%)1903.77E-043.37E-0289.60×25610659 (11.96%)14903 (8.36%)12437 (6.98%)7 (2.08%)1902.05E-043.37E-02164.25×R*esource percentages are for a Virtex-4 XC4VLX200ff1513-11.Part of our purpose here is to argue that the reduction tasks can be implemented efficiently using the systolic array design. With Table4, we intend to demonstrate the effect of each particular operation for users that might want to use any of them as part of a bigger design. Clearly, the ultimate speedup may be dictated by other components of the final design, as shown by Tables 5 and 6.Table 5
Experimental results for stage 1.
ndSASlicesSlice FFs*4-input LUTS*FIFO16/ RAMB16*tF (ms)tC (ms)Speedup8321651 (1.85%)1908 (1.07%)3047 (1.71%)25 (7.44%)0.041.2729.59×9643011 (3.38%)3596 (2.02%)5515 (3.1%)27 (8.04%)0.192.6413.87×101285860 (6.58%)7201 (4.04%)11186 (6.28%)33 (9.82%)0.497.2214.78×111286191 (6.95%)7728 (4.34%)11755 (6.60%)43 (12.80%)1.0122.9722.74×1225611880 (13.34%)14907 (8.37%)22573 (12.67%)55 (16.37%)2.5482.2832.46×1325614101 (15.83%)16981 (9.53%)25676 (14.41%)95 (28.27%)8.89596.967.13×R*esource percentages are for a Virtex-4 XC4VLX200ff1513-11.Table 6
Experimental results for stage 2.
nBasis VarsClass RepsSlices*Slice FFs*4-input LUTs*FIFO16/ RAMB16*tF (ms)tC (ms)Speedup82237745 (8.69%)8335 (4.68%)13846 (7.77%)28 (8.33%)0.101.5916.67×83637950 (8.92%)8352 (4.69%)14205 (7.97%)28 (8.33%)0.447.7617.68×84907955 (8.93%)8359 (4.69%)14213 (7.98%)30 (8.93%)0.8416.6019.71×92258520 (9.56%)9220 (5.17%)15205 (8.53%)28 (8.33%)0.101.8517.64×93778539 (9.58%)9226 (5.18%)15239 (8.55%)28 (8.33%)0.6911.2016.13×94948582 (9.63%)9260 (5.2%)15340 (8.61%)34 (10.12%)1.1420.2617.81×102249188 (10.31%)10098 (5.67%)16710 (9.38%)31 (9.23%)0.101.7817.95×103659192 (10.32%)10106 (5.67%)16714 (9.38%)31 (9.23%)0.478.3017.69×104929201 (10.33%)10130 (5.69%)16742 (9.40%)33 (9.82%)0.8618.4321.36×1122510130 (11.37%)10978 (6.16%)18024 (10.12%)34 (10.12%)0.112.0118.80×1137510156 (11.4%)11001 (6.17%)18067 (10.14%)36 (10.71%)0.5812.8022.22×1149210115 (11.35%)11008 (6.18%)18025 (10.12%)37 (11.01%)1.0523.0021.93×1222510861 (12.19%)11861 (6.66%)19077 (10.71%)36 (10.71%)0.102.0519.93×1237010806 (12.13%)11881 (6.67%)19012 (10.67%)38 (11.31%)0.5210.6820.41×1249310877 (12.21%)11895 (6.68%)19121 (10.73%)38 (11.31%)1.0319.9019.26×1322511405 (12.8%)12747 (7.15%)20240 (11.36%)39 (11.61%)0.102.1320.74×1337111469 (12.87%)12760 (7.16%)20358 (11.43%)40 (11.90%)0.5312.5023.68×1349111416 (12.81%)12771 (7.17%)20274 (11.38%)41 (12.20%)0.8521.3825.21×R*esource percentages are for a Virtex-4 XC4VLX200ff1513-11.


## 8.2. Interpolation Algorithm

We modeled the proposed interpolation algorithm in Verilog HDL using the reduction tasks and the architectures described in previous sections as building blocks. The behavioral simulation and synthesis tools as well as the software compilation parameters and platform were the same as in the reduction task experiments.The results are shown separately for each of the two interpolation stages since, once the bases have been computed by the first stage, a scientist's intervention would be necessary to choose the most appropriate before proceeding to the second stage. For both stages, the synthesis tool determined maximum operating frequencies above 150 MHz, thus we chose 150 MHz for the experiments. The Xilinx ISE 11.1 place and route tools (with default effort level) were able to meet the timing requirements reported for Stages 1 and 2. We attribute this to the fact that the great majority of connections in our design are local and grow only linearly withn.


### 8.2.1. Stage 1

Table5 reports timing and resource results for the first stage. Randomly generated sets of 100 n-tuples (8≤n≤13) were input to implementations with diverse systolic array depths. The table shows results for the dSA case with the shortest execution time for each n. The smaller values of n generated fewer intermediate terms during computation, hence smaller values of dSA produced better results since they incurred in less unnecessary overhead for filling and emptying. As the number of variables increases, the number of terms produced by Algorithm 1 increases exponentially, which is clearly evidenced by the run times, that is, columns tF and tC. The FPGA implementation achieves speedups of 13× to 67×, which mostly increase with n. The increase in speedup at the higher values of n can be attributed to a higher percentage of time of full pipeline utilization as well as the increased benefit of parallelism in the comparisons between the n-tuple variables, that is, the computation of the r(i,j) terms in Algorithm 1. As expected, the use of necessarily serial processing tasks such as the multiplication of the DOLs to obtain the DNF (see Example 2) limit the speedup achieved by the FPGA implementation when compared to the results obtained using the cover subtask by itself.Although the ultimate goal of reverse engineering might be understanding complete complex biological systems, recent reported models limit their exploration over specific subsystems or punctual mechanisms, for example, the study cell-cycle regulatory network of fission yeast and apoptosis (programmed cell death) [17, 18], using a moderate number of genes between 8 and 20. Our results lead us to believe that the advantages of using the proposed architecture in reconfigurable logic will be sustained even at bigger sizes. As the sizes increase our architecture will have to rely increasingly on iterations rather than parallel computations. Nevertheless, due to the custom nature (finite field) of the data, even at these sizes the parallelism over each iteration will be enough to offer significant speedups over the alternative, fully serial implementation in general purpose processors.Speedup is maintained at a cost of increasing resource utilization. For the considered cases, BRAMs, which are used to implement the various FIFOs and RAMs of stage 1, is the fastest growing resource. This could pose a challenge to maintain performance at even highern values. However, this issue can potentially be resolved by complementing the BRAM capacity with a dedicated external memory bank, as is commonly included in modern reconfigurable computer platforms [19, 20]. A study of the impact of external memory on the performance of the proposed architecture is future work. An advantage of our proposed architecture is that all the required memories behave as FIFOs. Thus, their accesses patterns are predictable and (as part of future work) we can use this information to device a buffering scheme that hides the latency of accessing the external memories.The increase in user logic (Slices, LUTs) in highern can be circumvented by using shallower systolic arrays (smaller dSA), which still maintain competitive speedups. For example, although not shown in the table, for n=13 speedups of 30× and 45× were obtained by using dSA values of 64 and 128.


### 8.2.2. Stage 2

As explained in Section2.2, our algorithm uses a base Xi computed in Stage 1 to determine a set of equivalence classes D/≡Xi from the original set of n-tuples. The equivalence classes are then used to compute an interpolating polynomial using (3). A greater number of variables in the chosen Xi and a high variability among the n-tuples translates to larger cardinalities of D/≡Xi. The greater the cardinality of D/≡Xi the higher the number of binomials computed by (3), which constitutes the bulk of computation for Stage 2.The timing results shown in Table6 support the previous statements. The sets of 100 n-tuples from Stage 1 along with bases of 2 to 4 variables were input to Stage 2. Observe that since the data sets were randomly generated there was high variability in the sets of n-tuples, which resulted in similar D/≡Xi cardinalities (column Class Reps) and run times (columns tF and tC) for the various sizes of n (for a given number of variables in the chosen basis). The FPGA implementation achieves speedups of 16× to 25×, which mostly increase with n. The moderate increment in speedup as n increases is attributed to the increased benefit of parallelism for the generation of binomials, and the comparison and addition of monomials.Resource utilization distribution in Stage 2 is roughly uniform between the user logic, that is, Slices, and BRAMs. Furthermore, the increment in resource utilization for increasing values ofn and numbers of variables in the basis is quite moderate, for example, the percentage of slices goes from 8.69% for n=8 to 12.87% for n=13.


## 8.2.1. Stage 1

Table5 reports timing and resource results for the first stage. Randomly generated sets of 100 n-tuples (8≤n≤13) were input to implementations with diverse systolic array depths. The table shows results for the dSA case with the shortest execution time for each n. The smaller values of n generated fewer intermediate terms during computation, hence smaller values of dSA produced better results since they incurred in less unnecessary overhead for filling and emptying. As the number of variables increases, the number of terms produced by Algorithm 1 increases exponentially, which is clearly evidenced by the run times, that is, columns tF and tC. The FPGA implementation achieves speedups of 13× to 67×, which mostly increase with n. The increase in speedup at the higher values of n can be attributed to a higher percentage of time of full pipeline utilization as well as the increased benefit of parallelism in the comparisons between the n-tuple variables, that is, the computation of the r(i,j) terms in Algorithm 1. As expected, the use of necessarily serial processing tasks such as the multiplication of the DOLs to obtain the DNF (see Example 2) limit the speedup achieved by the FPGA implementation when compared to the results obtained using the cover subtask by itself.Although the ultimate goal of reverse engineering might be understanding complete complex biological systems, recent reported models limit their exploration over specific subsystems or punctual mechanisms, for example, the study cell-cycle regulatory network of fission yeast and apoptosis (programmed cell death) [17, 18], using a moderate number of genes between 8 and 20. Our results lead us to believe that the advantages of using the proposed architecture in reconfigurable logic will be sustained even at bigger sizes. As the sizes increase our architecture will have to rely increasingly on iterations rather than parallel computations. Nevertheless, due to the custom nature (finite field) of the data, even at these sizes the parallelism over each iteration will be enough to offer significant speedups over the alternative, fully serial implementation in general purpose processors.Speedup is maintained at a cost of increasing resource utilization. For the considered cases, BRAMs, which are used to implement the various FIFOs and RAMs of stage 1, is the fastest growing resource. This could pose a challenge to maintain performance at even highern values. However, this issue can potentially be resolved by complementing the BRAM capacity with a dedicated external memory bank, as is commonly included in modern reconfigurable computer platforms [19, 20]. A study of the impact of external memory on the performance of the proposed architecture is future work. An advantage of our proposed architecture is that all the required memories behave as FIFOs. Thus, their accesses patterns are predictable and (as part of future work) we can use this information to device a buffering scheme that hides the latency of accessing the external memories.The increase in user logic (Slices, LUTs) in highern can be circumvented by using shallower systolic arrays (smaller dSA), which still maintain competitive speedups. For example, although not shown in the table, for n=13 speedups of 30× and 45× were obtained by using dSA values of 64 and 128.


## 8.2.2. Stage 2

As explained in Section2.2, our algorithm uses a base Xi computed in Stage 1 to determine a set of equivalence classes D/≡Xi from the original set of n-tuples. The equivalence classes are then used to compute an interpolating polynomial using (3). A greater number of variables in the chosen Xi and a high variability among the n-tuples translates to larger cardinalities of D/≡Xi. The greater the cardinality of D/≡Xi the higher the number of binomials computed by (3), which constitutes the bulk of computation for Stage 2.The timing results shown in Table6 support the previous statements. The sets of 100 n-tuples from Stage 1 along with bases of 2 to 4 variables were input to Stage 2. Observe that since the data sets were randomly generated there was high variability in the sets of n-tuples, which resulted in similar D/≡Xi cardinalities (column Class Reps) and run times (columns tF and tC) for the various sizes of n (for a given number of variables in the chosen basis). The FPGA implementation achieves speedups of 16× to 25×, which mostly increase with n. The moderate increment in speedup as n increases is attributed to the increased benefit of parallelism for the generation of binomials, and the comparison and addition of monomials.Resource utilization distribution in Stage 2 is roughly uniform between the user logic, that is, Slices, and BRAMs. Furthermore, the increment in resource utilization for increasing values ofn and numbers of variables in the basis is quite moderate, for example, the percentage of slices goes from 8.69% for n=8 to 12.87% for n=13.


## 9. Conclusion

This paper presents a new methodology based on Lagrange interpolation with two important properties: (1) it identifies redundant variables and generates polynomials containing only nonredundant variables, and (2) it computes exclusively on a reduced data set. The analysis of the methodology for its hardware implementation led us to the identification of several reduction tasks which were generalized to a simple algorithm. The generalized algorithm can be efficiently mapped to a systolic array in which each processing cell implements a pair of binary operations between an incoming and a stored value. The tasks of Boolean cover, distinctness, and multivariate polynomial addition were implemented and served as building blocks to the rest of the application. The FPGA implementation of the reduction operations and the complete application achieved speedups of up to 172× and 67×, respectively, as compared to software implementations run on a contemporary CPU, with moderate resource utilization.



---

*Source: 313479-2010-12-29.xml*