Rémi Garcia
email: remi.garcia@etu.univ-nantes.fr

Anastasia Volkova

Alexandre Goldsztejn

Design of optimal multiplierless digital filters

Keywords: digital signal processing, filter design, multiple constant multiplication, ilp models, iir filters

In this work we solve the design of second order infinite impulse response (IIR) filters for the first time. This design is based on integer linear programming (ILP) methods. One of the underlying problem of this design is the multiple constant multiplication (MCM) problem. We propose variants to known ILP models that solve the MCM problem and illustrate our improvements with multiple examples. Our multiplierless filter design is optimal with respect to the number of adders used. Combined with MCM variants this allows addressing this optimal design in multiple contexts, in particular for a field-programmable gate array (FPGA) implementation.

Introduction

In the beginning there was nothing, which exploded.

Terry Pratchett

Digital circuits are replacing analog ones for many years, and analog systems find their digital counterparts. Filters are systems that perform manipulations to reduce or enhance certain characteristics of signals. Digital filters are less dependent of the environment as it is electronic. The design of such filters is of great interest for many reasons, for example, speeding up a filter allows reducing the delay between the input and the output signals. Decreasing the size of the hardware on which digital filters are implemented is of great interest, e. g. for embedded systems like drones. This can be done, in particular, using fixed-point arithmetic: it is a trade-off between precision and efficiency. In fixed-point arithmetic, the size of the circuit can be directly reduced by using smaller wordlengths. However, this worsen the precision and new design algorithms are in place.

Filters are designed to satisfy filter specifications. The traditional filter design process is decomposed into three steps: Find filter coefficients, quantize these coefficients and optimize the hardware implementation of the filter. However, it is hard to infer the best coefficients for hardware implementation as early as the first step. The quantization step makes unlikely the possibility to certify optimality of the hardware implementation according to the filter specification. Given filter specifications, our goal is to directly find filter coefficients in fixed-point arithmetic. Volkova et al. [START_REF] Kumm | Design of Optimal Multiplierless FIR Filters[END_REF] proposed a method that combines into a single procedure this whole design. Dropping the quantization and finding filter coefficients directly in fixed-point arithmetic allows combining the search for filter coefficients with the hardware optimization. This was applied for finite impulse response (FIR) filters and proved to be efficient [START_REF] Kumm | Design of Optimal Multiplierless FIR Filters[END_REF] using integer linear programming (ILP) based models. Obtained solutions are certified optimal to the chosen criterion and our goal is to generalize this method to infinite impulse response (IIR) filters.

Sum of products are involved in filter evaluations. In order to be more efficient, multiplications can be replaced with bit shifts, additions and subtractions. This allows reducing the size of the circuits and speeding up the operations. The multiplierless hardware implementation of a filter through its coefficients relies on the multiple constant multiplication (MCM) problem. Solving this problem allows reducing the hardware cost and/or speed up computations. This problem was tackled in an optimal way with ILP based models by Kumm [START_REF] Kumm | Optimal Constant Multiplication Using Integer Linear Programming[END_REF][START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF] and efforts are being made to improve efficiency of MCM solving. To improve the hardware implementation, MCM builds on bit shifts and addition or subtraction. The cost and space needed for these operations is controlled by the fixed-point representation.

The MCM problem has been tackled in many ways and still has many layers to unfold. We participated to this by analyzing existing ILP models and proposing variants. The optimal design of FIR filters is not close to an end and any improvement is a step in the right direction. We used a similar approach for the design of second order IIR filters and overcame a few issues inherent to ILP modeling of IIR filters. Overall, the contributions of this work are:

• The generalization of a new filter design process to second order IIR filters;

• Contributions to the MCM problem.

Our work on mathematical models has been implemented and experimented using MILP solvers, as CPLEX [START_REF] Cplex | CPLEX User's Manual[END_REF] and Gurobi [START_REF]Gurobi Optimizer Reference Manual[END_REF], and nonlinear programming (NLP) solvers, as SCIP [START_REF] Gleixner | The SCIP Optimization Suite 6.0, software, Optimization Online[END_REF], through JuMP [START_REF] Dunning | JuMP: A Modeling Language for Mathematical Optimization[END_REF], a package for linear and nonlinear modeling in Julia. ScaLP [START_REF] Sittel | ScaLP: A Light-Weighted (MI)LP Library[END_REF] was considered but discarded because of the absence of nonlinear modeling capabilities. For MCM those experiments showed that ILP formulations still have many interesting improvements to be think of and that NLP models are not a promising path to go on. Experimental results for the design of second order IIR filters using ILP models were encouraging and this lead is of great interest.

This work is organized into the following sections: First, the current state of the research about digital signal processing, MCM and FIR filter design is detailed in Section 2; second, in Section 3, we give an analysis of ILP formulations that are used to solve MCM problems and our improvements to these. Then, in Section 4, we extend the FIR filter design process to second order IIR filters by dealing with quadratic terms in the modeling and stability constraints. Finally, in Section 5 we propose our perspectives on leads that we think have great potential.

Notation

For the rest of this work, we will use the following notations.

Background

Wisdom comes from experience. Experience is often a result of lack of wisdom.

Terry Pratchett

Digital Filters

The design of digital filters is an area of research that emerged from the advancements of digital systems. Digital signal processing (DSP), ranging from communications to the analysis of astronomical, seismic signals, has a lot of applicable uses in our modern world. Despite increasing computational power, DSP methods still need to be improved. Actually, in order to obtain a denser time-discretization or faster computation, simply increasing computational power does not suffice and theoretical results are awaited. Furthermore, reducing the circuit area does not only rely on miniaturization advances but on theoretical results too. Systems which do manipulations on spectra of signals are called digital filters 1 . A filter takes an input signal, i. e. a sequence of numbers, and outputs a modified signal. The faster the filter computes its output the more inputs it can take and the smaller is the delay. It allows to work on a closer spectrum of the continuous signal, which is a frequent objective, and to be more reactive to the input. A signal will be seen as a sequence of numbers. Filters that we will focus on are linear time-invariant (LTI) discrete-time systems with a single input and a single output (SISO). A causal filter uses past and present inputs and past outputs to output a single value.

Constant-coefficient difference equations (CCDE) define the filter behavior. A CCDE defines the relation between an input signal u n ∈ (u k ) k∈Z and an output signal y n ∈ (y k ) k∈Z , it states that the output y n can be computed using past and present inputs and past outputs:

N 1 k=0 b k u n-k = N 2 k=0 a k y n-k , a 0 0. ( 1
)
By normalizing, a 0 = 1 can be ensured. Hence CCDE is also the next relation which explicitly shows the output at time n with respect to past data:

y n = N 1 k=0 b k u n-k - N 2 k=1 a k y n-k . (2)
Clearly, max {N 1 , N 2 } points in the past are needed, thus there is a delay which grows with this value which is the filter order. If the filter does not rely on outputs, i. e. a k = 0 ∀k ∈ [[1, N 2 ]], then the filter is called a finite response impulse (FIR) filter. If past output values are needed to compute the present output value, then the form of a such filter is called recursive or infinite impulse response (IIR). For most applications, CCDE is not the method used to compute filter outputs or to find the right design. The impulse signal, denoted by δ, permits to obtain another description of a filter and is defined this way:

δ(k) := 1 if k = 0 0 if k 0 . (3)
The impulse response of a filter is its output to δ [START_REF] Antoniou | Digital Filters: Analysis, Design, and Signal Processing Applications[END_REF]. Given a filter F , applying the z-transform [START_REF] Prandoni | Signal Processing for Communications[END_REF] on its impulse response gives the transfer function. The transfer function denoted by H fully characterizes F and is equal to:

H(z) = b(z) a(z) = N 1 k=0 b k z -k 1 + N 2 k=1 a k z -k , (4)
where

z ∈ C, b N 1 0, a N 2 0, a k ∈ R, ∀k ∈ [[0, N 2 ]] and b k ∈ R, ∀k ∈ [[0, N 1 ]
]. Delays and past inputs or outputs are still visible in the transfer function: z -k corresponds to the k-th delay of the input. In other words z -k is a point in the past, the greater is k the further is this point. The frequency response is the evaluation of its transfer function H on the unit circle: H e iω , ω ∈ [-π, π]. Designing a filter basically means constraining the frequency response for different frequencies. The design relies on filter specifications that are bounds, D and D, on the modulus of the frequency response:

D(ω) ≤ H e iω ≤ D(ω) ∀ω ∈ [0, π] . (5)
A gain, denoted by G ∈ R, can be added to the previous equation:

GD(ω) ≤ H e iω ≤ GD(ω) ∀ω ∈ [0, π] . (6)
This gain might be necessary to find filter coefficients for givens filter type and specifications. It can be fixed to one in some cases.

Example 1. For a lowpass filter specification, low frequencies are preserved while higher frequency are suppressed, filter specifications will resemble to this:

D(ω) = 1 -ε if ω < ω 1 0 if ω > ω 2 , (7)
D(ω) = 1 + ε if ω < ω 1 0 + ε if ω > ω 2 , (8)
where ω 1 , ω 2 ∈ [0, π] and ω 1 < ω 2 . For ω 1 = 0.2π, ω 2 = 0.5π and ε = 0.01 this could lead to the 16th order FIR filter with (b k ) k = (3, 6, 0, -16, -19, 12, 76, 128, 128, 76, 12, -19, -16, 0, 6, 3) and a gain of G = 376.999805 [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF]. This filter is represented by its transfer function in Figure 1. The passband and the stopband of this filter are [0, ω 1 ] and [ω 2 , π] respectively.

Another important aspect of filter design is the stability. Despite that the absence of stability can be wanted for oscillators, stable filters generally are what is needed. A filter is said to be stable if any bounded input results in a bounded output (BIBO). This definition does not lead to a simple way to test for stability in the design process. Zeros and poles are, respectively, the roots of the numerator and the denominator of the transfer function [START_REF] Cplex | CPLEX User's Manual[END_REF]. From poles, a brief characterization of stability can be obtained which is more suitable in practice [10, Section 6.1.4]: a filter F is stable ⇔ poles of F are inside the complex unit circle. [START_REF] Antoniou | Digital Filters: Analysis, Design, and Signal Processing Applications[END_REF] This means that a filter is stable if and only if for all λ, poles of the filter, |λ| < 1. It is unambiguous that FIR filters are inherently stable: since there are no poles, they all are inside the unit circle. On the other hand, IIR filters require a special attention with their design due to this stability requirement. This is usually done by focusing on pole, to ensure the stability, and zero positions as the first design step. The design of filters can be done in compliance with field-programmable gate array (FPGA) requirements. FPGAs are integrated circuits that it is possible to reconfigured after manufacturing. Despite FPGAs have a strong potential with floating-points [START_REF] De Dinechin | When FPGAs Are Better at Floating-Point than Microprocessors[END_REF], fixed-point arithmetic is privileged for filter implementations. This permits to speed up the computation replacing multiplications, b k x n-k 's and a k y n-k 's, with faster operations, Section 2.2 will largely tackle this aspect. It is possible to convert b k 's and a k 's into fixed-point representation with a wordlength that will, in an approximate manner, be included in the filter specifications. Let d be that expected wordlength. To ease the use of optimization methods, it can be ensure that b k 's and a k 's are in the interval [-1, 1]. To achieve that, initial bounds must be found for b k 's and a k 's: b k 's, b k 's, a k 's and a k 's. Then, the transfer function can be "normalized" by multiplying it with:

2 log 2 max k∈[[0,N 1 ]] -b k ,b k 2 log 2 max k∈[[0,N 2 ]] {-a k ,a k } . ( 10
)
The use of log 2 and power of two permits to avoid numerical error since it basically corresponds to a shift which is possible in fixed-point arithmetic. Furthermore, a i 's and b i 's can be multiplied by 2

d such that b k = 2 d b k 's and a k = 2 d a i 's are in [[-2 d , 2 d ]]
. Actually, CCDE can be rewritten into the following:

2 d y n = N 1 k=0 b k x n-k - N 2 k=1 a k y n-k , (11)
where

a k = 2 d a k , ∀k ∈ [[0, N 2 ]] and b k = 2 d b k , ∀k ∈ [[0, N 1 ]].
As a consequence the fixed-point coefficient design can be formulated as an integer coefficient design. FIR filters are the most studied filters for their robustness with respect to fixed-point implementation and, as already mentioned, their unconditional stability [START_REF] Prandoni | Signal Processing for Communications[END_REF]Section 7.1.1]. However FIR filters come with two disadvantages: longer input-output delay and still a high computational cost. On the other hand, IIR filters are faster but stability is not a priori guaranteed thus it has to be guaranteed in the design process. Furthermore, numerical precision is a more frequent issue than with FIR filters and the design is much more complex. Yet, if a way to control those drawbacks is found then the achieved speedup is of great interest for DSP applications. In either case we want to evaluate filters with FPGAs to gain in efficiency and to reduce the memory usage we could optimize the multiplication between the filter coefficients and the inputs and past outputs.

Multiplierless Constant Multiplication

Multiple Constant Multiplication (MCM) is a problem that consists in the speedup of the multiplication of an input number by multiple integer constants. Using generic multipliers in integrated circuits, like FP-GAs, is costly since those multipliers do not have any knowledge of the predefined constants and potential improvements. The initial complexity, which is roughly a cost of n 2 with n the number of bits, can be reduced, especially in fixed-point arithmetic. In fixed-point arithmetic context we mentioned that coefficients a k 's and b k 's can be treated as integers. Next, MCM could be used to optimize the filter evaluation and finally coefficients could be converted back into their fixed-point arithmetic writing. Those integer/fixed-point conversions simply rely on bit shifts which do not degrade the numerical error in fixed-point arithmetic.

Obviously, the research started with a Single Constant Multiplication (SCM). Premisses of SCM problem were given by Bernstein in 1986 [START_REF] Bernstein | Multiplication by integer constants[END_REF]. Since then, SCM has been well studied. It can be tackled using canonical signed digit (CSD) [START_REF] Hewlitt | Canonical signed digit representation for FIR digital filters[END_REF][START_REF] Backenius | Two's Complement Conversion to Minimal Signed Digit Code[END_REF], heuristics, bounds [START_REF] Gustafsson | Lower Bounds for Constant Multiplication Problems[END_REF], optimal approaches [START_REF] Thong | An Optimal and Practical Approach to Single Constant Multiplication[END_REF], etc. Next, MCM problem quickly followed since it is a generalization of SCM problem and it has many applications.

Multiplication has a computational cost that can be reduced using only "shifts", "add" and "subtract" instructions with dedicated hardware. A shift consists in the bit shift of the binary representation of the current number. A left 1-shift corresponds to the multiplication of the initial number by 2. In general a left s-shift of an integer a, written a << s, corresponds to a multiplication of a by 2 s . This operation can be hardwired on FPGAs thus its cost remains negligible. A first straightforward example of a speedup that can be obtained by MCM problem is the replacement of the multiplication by 16 with a 4-shift, multiply by 16 does not have a negligible cost while a shift does. However, 17 or 15 cannot be obtained using only shifts. Add and subtract operations combined with shifts permit to obtain any integer. Yet, those operations have a cost, therefore the overall objective is to replace the initial multiplication with a minimum number of add and subtract operations to reduce the total cost of the multiplications. Those operations, additions and subtractions, will be referred to as adders. In other words, an adder takes two inputs and computes the sum or subtraction of those. To be more visual we will say that an adder has a left input and a right input, this should be clearly visible in figures along the text. Both operations have an equivalent cost, hence SCM and MCM consist in the rewriting of a multiplication using a minimum number of adders and as many shifts as needed. Furthermore, since multiple constants are involved in MCM problem it is possible to share adders between constants to speed up the multiplications even more.

An easy way to rewrite any multiplication is using the binary notation and deducing shifts and adds from it, e. g. for 7x:

7x = 0111 2 x = 2 2 x + 2 1 x + 2 0 x = (x << 2) + (x << 1) + (x << 0) , (12)
so two adders are required to compute 7x this way. This corresponds to the number of 1's in their binary representation minus one. The corresponding circuit is shown in Figure 2a.

There are other ways to represent numbers such as signed digit (SD) or canonical signed digit (CSD) representation [START_REF] Hewlitt | Canonical signed digit representation for FIR digital filters[END_REF]. A ternary representation system in which 1 corresponds to -2 bit position is called SD representation. It is direct that a similar addition/subtraction and shift circuit naturally derives from this representation. A non-adjacent form is the particular case in which no two non-zero values are adjacent. The CSD representation2 is a ternary representation system in which only non-adjacent forms are allowed. For a given integer its CSD representation is unique and can be computed in polynomial time, the CSD representation of 7 is:

1001 CS D . (13)
Figure 2b shows the corresponding circuits to compute 7x, which is optimal. This representation gives a minimal Hamming weight (the number of non-zero values in the representation), and therefore it permits to obtain a good bound on the number of adders in general. Although CSD representations have this important property, it does not ensure a minimal number of adders, e. g. with an example from [3, p. 9-11]:

93x = 0101 1101 2 x = 1010 0101 CS D x, (14)
naturally uses three adders. Yet, 93x = ((x << 1) + x) + (((x << 1) + x) << 5) ,

and despite three addition terms, the term (x << 1) + x is redundant thus can be computed only one time and used twice as shown in Figure 3: two adders suffice. This is why solving SCM problem is needed. However, when multiple constants are involved, the solving of multiple SCM problems does not give an optimal solution as shown in Figure 4.

Usually, MCM problem is presented in its decision form:

Given multiple integers, is there an adder graph that ouputs these integers with N adders,

where an adder graph is an adder circuit for constant multiplication. Adder graphs are directed acyclic graph (DAG) and Figures 2, 3 and4 are examples of such graphs. Each node, except the input node, corresponds to an adder. Nodes have an in-degree of two. Each edge weight, which can go along with an arrow, represents a bit shift. No weight or zero means the absence of a bit shift. Each node has an associated value which corresponds to the multiple computed in the node, this constant is called fundamental. It can be desirable to minimize the number of adders N instead of dealing with a decision problem: this gives an optimization problem. It should be noted that both problems, the decision one and its minimization counterpart, have similar complexities.

x

← 1 (3x) ← 5 - 93x
Figure 3: DAG of 93 in CSD representation needs three adders while minimum adder realization only takes two.

x

← 3 ← 4 (7x) - 23x 7x
Figure 4: One adder is needed for 7x and 23x needs two. Sum of SCMs needs three while an optimal solution for MCM only needs two in total.

We focus our efforts on exact approaches, precisely on ones that use ILP based models to solve this problem [START_REF] Kumm | Optimal Constant Multiplication Using Integer Linear Programming[END_REF][START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF]. In particular, two ILP based models were developed. The first model will be referred to as ILP Formulation 1 and the second is denoted by ILP Formulation 2. ILP Formulation 1 tackles the theoretical MCM problem while ILP Formulation 2 deals with a more practicable version of MCM problem in which a constraint on the adder depth is added. We present here this variant and a few useful terms.

Terminology

The adder depth, denoted by AD(•), is a property of nodes in an adder graph. It is the delay before a node is computed. The adder depth of the input node is 0. Let c be a node, with inputs n 1 and n 2 , its adder depth is defined as follows:

AD(c) = max {AD(n 1 ) , AD(n 2 )} + 1. ( 17
)
The adder depth of a graph is the maximum of its nodes' adder depth. This metric is essential for two reasons: First, note that computation time is proportional to the adder depth; second, a greater adder depth might imply more additions and subtraction to compute a given constant, thus a bigger numerical error inherent to computer arithmetic. For example it is clear that in Figure 2a the second adder depends on the first adder's computation: Second adder has to wait for the first one and the computation error that the first adder outputs is passed on and amplified by the second. For these reasons minimizing the adder depth as a second objective in MCM problem can be considered. This variant is denoted MCM MAD . Furthermore, bounding the adder depth, i. e. adding a constraint on the problem, allows new possible ways to speed up the resolution, e. g. it naturally limits the time-delay. This problem, which is the one that is useful for FPGAs, will be referred to as MCM BAD . This last variant is the one covered by ILP Formulation 2.

Remark. The adder depth of DAGs represented Figures 2a and3 is two while the adder depth of the DAG represented Figure 2b is equal to one.

An odd fundamental graph is an adder graph for which all fundamentals are odd. Dempster and Macleod showed that only odd fundamental graphs can be considered without any loss of generality [START_REF] Dempster | Constant integer multiplication using minimum adders[END_REF]. Additionally, any MCM instance can be transformed into an equivalent MCM problem with only odd positive constants, this transformed MCM problem will be referred to as MCM odd problem. Consequently, solving MCM problem in this restricted scope is possible and reduce the search space. Reducing the search space might speed up the solving and MCM odd can be combined with other variants. For that reason, in general by MCM we will refer to MCM odd . The following provides understanding of how MCM odd problem is derived from an instance of MCM problem. Let T ⊂ Z be the set of target constants of a given MCM problem,

T odd := t ∈ N | t mod 2 = 1 ∧ ∃k ∈ N such that 2 k t ∈ T , (18)
is the set of constants for the equivalent MCM odd instance. This set of positive odd constants T odd corresponds to the elements of T divided by two until they are odd numbers. It can be rewritten in a more compact and readable way using the odd part notion [20, Sequence A000265]. The odd part of a positive integer n is the biggest odd number that divides n. It is written odd(n) and can be defined as follows:

odd(n) := n gcd(2 n , n) . ( 19
)
The odd part of a number n is computed by dividing n by two until the obtained number is odd. With that definition of odd(n) it can be deduced that:

T odd = {odd(|t|) | t ∈ T } . ( 20
)
Example 2. Given a type II low-pass FIR filter with 3, 0, -25, 0, 150, 256, 150, 0, -25, 0, 3 its coefficients, target constants, T = {0, 3, -25, 150, 256} ,

, T odd = {0, 1, 3, 5, 75} . (21) becomes
In both cases, 0 can be exclude from the target constants, in the second case 1 can be exclude too: It corresponds to the absence of shifts and adders. The remaining set to work on is smaller, both in size and in maximum value.

ILP Formulation 1

Kumm gave in [START_REF] Kumm | Optimal Constant Multiplication Using Integer Linear Programming[END_REF] a first model for MCM problem: ILP Formulation 1. This model has no objective function, it corresponds to the decision problem. However it was done to minimize the number of adders, to tackle the optimization problem. Hence, the solving process is embedded in a while loop which generates a model with a fixed number of adders and solves it. While it is infeasible, it increases the number of adders and repeats. An optimal solution for the optimization problem is found when a model has a solution, i. e. when a decision problem has a solution. The initial number of adders could be one or any greater known lower bound of the given instance.

Formulations are defined by a set of data, that partly depends on each instance and solving choices, a set of variables and a set of constraints. Eventually, an objective might be added. The input data and variables

c a ∈ [[0, 2 d ]], ∀a ∈ [[0, N A ]]: constant obtained in adder a; c a,i ∈ [[0, 2 d ]], ∀a ∈ [[1, N A ]], i ∈ {l, r}: constant of adder from input i (left or right) before adder a; c sh a,i ∈ [[0, 2 d ]], ∀a ∈ [[1, N A ]], i ∈ {l, r}: constant of adder from input i before adder a after the shift; c sh,sg a,i ∈ [[-2 d , 2 d ]], ∀a ∈ [[1, N A ]], i ∈ {l, r}: signed constant of adder from input i before adder a after the shift; Φ a,i ∈ {0, 1}, ∀a ∈ [[1, N A ]], i ∈ {l, r}: sign of i input of adder a. 0 for + and 1 for -; c a,i,k ∈ {0, 1}, ∀a ∈ [[1, N A ]], i ∈ {l, r}, k ∈ [[0, N A -1]]: 1 if input i of adder a is adder k, 0 otherwise; φ a,i,s ∈ {0, 1}, ∀a ∈ [[1, N A ]], i ∈ {l, r}, s ∈ [[S min , S max ]]: 1 if shift on input i before adder a is equal to s, 0 otherwise; o a, j ∈ {0, 1}, ∀a ∈ [[1, N A ]], j ∈ [[1, N O ]]: 1 if adder a is equal to the j-th target constant, 0 otherwise.
of ILP Formulation 1 are defined Table 1. Then, constraints follow:

c 0 = 1 (C1.1) c a = c sh,sg a,l + c sh,sg a,r ∀a ∈ [[1, N A ]] (C1.2) c a,i = c k if c a,i,k = 1 ∀a ∈ [[1, N A ]], i ∈ {l, r} , ∀k ∈ [[0, a -1]] (C1.3) a-1 k=0 c a,i,k = 1 ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.4) c sh a,i = 2 s c a,i if φ a,i,s = 1 ∀a ∈ [[1, N A ]], i ∈ {l, r} , s ∈ [[S min , S max ]] (C1.5) S max s=S min φ a,i,s = 1 ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.6) φ a,l,s = 0 ∀s ∈ [[1, S max ]], a ∈ [[1, N A ]] (C1.7) φ a,l,s = φ a,r,s ∀s ∈ [[S min , -1]], a ∈ [[1, N A ]] (C1.8) c sh,sg a,i = -c sh a,i if Φ a,i = 1 ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.9) c sh,sg a,i = c sh a,i if Φ a,i = 0 ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.10) Φ a,l + Φ a,r ≤ 1 ∀a ∈ [[1, N A ]] (C1.11) c a = C j if o a, j = 1 ∀a ∈ [[0, N A ]], j ∈ [[1, N O ]] (C1.12) N A a=0 o a, j = 1 ∀ j ∈ [[1, N O ]] (C1.13)
Constraint (C1.2) states that the value of an adder is equal to the addition of its inputs. Constraints (C1.3) and (C1.5) state that the unsigned i-th input of an adder a is equal to the s shifted constant of the adder k if the input i of adder a is adder k and is shifted by s. Constraints (C1.4) and (C1.6) state that for each input of each adder there is only one previous adder and it has been shifted only once -a 0-shift is allowed. Constraints (C1.7) and (C1.8) ensure that one input is positive shifted while the other shift is zero or that both have an identical negative shift, this is sufficient and reduces the exploration space [START_REF] Dempster | Constant integer multiplication using minimum adders[END_REF]Th. 3]. Constraints (C1.9), (C1.10) and (C1.11) apply the sign for each input and ensure that only one input per adder is negative. This last constraint is redundant since

c a ∈ [[0, 2 d ]] would be contradictory with c a = c sh,sg a,l + c sh,sg a,r if c sh,sg
a,l and c sh,sg a,r are both negative. Finally (C1.12) fixes the value of adder a to the j-th constant and (C1.13) ensures that each constant has exactly one adder that computes it.

Note that this model is not linear: implications are not linear constraints and are present in constraints (C1.3), (C1.5), (C1.9), (C1.10) and (C1.12). However common mixed integer linear programming (MILP) solvers, as CPLEX or Gurobi, handle those constraints very well and this permits to avoid using big M [START_REF] Klotz | Practical guidelines for solving difficult mixed integer linear programs[END_REF]. Yet, the five constraints involving implications could have been replaced by two big M constraints each:

c a,i ≥ c k + c a,i,k -1 × M ∀a ∈ [[1, N A ]], i ∈ {l, r} , ∀k ∈ [[0, a -1]] (C1.3a) c a,i ≤ c k + 1 -c a,i,k × M ∀a ∈ [[1, N A ]], i ∈ {l, r} , ∀k ∈ [[0, a -1]] (C1.3b) c sh a,i ≥ 2 s c a,i + φ a,i,s -1 × M ∀a ∈ [[1, N A ]], i ∈ {l, r} , s ∈ [[S min , S max ]] (C1.5a) c sh a,i ≤ 2 s c a,i + 1 -φ a,i,s × M ∀a ∈ [[1, N A ]], i ∈ {l, r} , s ∈ [[S min , S max ]] (C1.5b) c sh,sg a,i ≥ -c sh a,l + Φ a,i -1 × M ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.9a) c sh,sg a,i ≤ -c sh a,l + 1 -Φ a,i × M ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.9b) c sh,sg a,i ≥ c sh a,l -Φ a,i × M ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.10a) c sh,sg a,i ≤ c sh a,l + Φ a,i × M ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.10b) C j ≥ c a + o a, j -1 × M ∀a ∈ [[0, N A ]], j ∈ [[1, N O ]] (C1.12a) C j ≤ c a + 1 -o a, j × M ∀a ∈ [[0, N A ]], j ∈ [[1, N O ]] (C1.12b)
An equality constraint in an implication constraint has to be decomposed into two inequality constraints with big M. One can replace the equal symbol and add a big M which is activated or not depending of the binary variable used for the implication. Constraint (C1.3) derives into constraints (C1.3a) and (C1.3b). For constraint

(C1.3a), if c a,i,k = 0 then c a,i ≥ c k + c a,i,k -1 × M is inactivated because of the big M that, in a way, implies c a,i ≥ -∞.
However, since M is finite, the real implication of c a,i,k = 0 is c a,i has to be greater than something smaller than its lower bound. The same reasoning holds for constraint (C1.3b). On the other hand, if c a,i,k = 1 then 0 × M is removed from both inequalities which becomes c a,i ≥ c k and c a,i ≤ c k , thus c a,i = c k . Kumm chose to use implications to avoid numerical instability in the solving process as it is suggested by Klotz and Newman [START_REF] Klotz | Practical guidelines for solving difficult mixed integer linear programs[END_REF]Section 3.4.]. Default integrality tolerance on Gurobi and CPLEX is equal to 10 -5 . Hence, only M > 10 5 might lead to numerical instability. Since M < 2 2d with d denoting the wordlength, default integrality tolerance is enough for wordlengths up to 8bit. Lowering integrality tolerance is possible and safe down to 10 -9 . This allows wordlengths up to 14bit but with a downside: an increasing computing time. For any larger wordlengths, big M constraints might not be usable. We give a more detailed analysis in Section 3.1.2 and a solution is proposed to bound big M by 2 d : with the drawback of creating more variables, this allows bigger wordlengths using big M constraints.

ILP Formulation 2

ILP Formulation 2 tackles MCM BAD and can be looked as a shortest-path in a hypergraph. Indeed, since the adder depth is bounded, the total number of possible values that can be obtained with shift, add and subtract operations is reduced. All those values can be computed: each value is a vertex and each way to compute that value a hyperedge. First, we will give a few mathematical tools that are necessary for the comprehension of that formulation. Hypergraphs can be constructed with these notions. Second, a shortest-path solution, taking the form of an ILP based model presented after, is applied.

An A-operation is the operation that produces the fundamental of each node depending on its inputs u, v ∈ N and on an A-configuration q = (l u , l v , r, s). The integers l u , l v ∈ N correspond to the left shifts of u and v respectively, r ∈ N is the output right shift and s ∈ {0, 1} stands for addition or subtraction. This can be synthesized into a compact definition:

A q (u, v) = 2 l u u + (-1) s 2 l v v 2 -r . ( 23
)
This definition can be related to the computations of nodes of adder graphs where A q (u, v) is the fundamental computed in a node from input nodes u and v.

It is said that u, v, w ∈ N are in an A-relation if ∃q an A-configuration such that w = A q (u, v). The definition of A-configuration is adapted to only allow odd fundamentals, which simplifies its formulation without any loss for MCM problem comprehension and modeling. The statement "q is an A-configuration" will refer, somewhat imprecisely, to q = (l u , l v , r, s) is a valid configuration such that A q (u, v) is an odd integer smaller than a given c max . The statement "q is an A-configuration" is imprecise since A-configuration directly depends on the context (u, v and c max ). In our MCM problem context with d a fixed wordlength we will always assume that c max = 2 d . This definition constrains r: For l u , l v and s fixed, exactly one r gives rise to a valid A-configuration. This is direct to prove using a reductio ad absurdum argument. From Aoperation, A-sets naturally derive. This notion can be applied either on integers or on sets of integers. The A * (u, v)-set, with u, v ∈ N, contains all odd fundamentals which can be obtained from each A-operation on (u, v). It corresponds to all the integers that are in an A-relation with u and v:

A * (u, v) := A q (u, v) | q an A-configuration . (24)
Remark. Without a c max , card(A * (u, v)) = ∞, i. e. there is an infinite number of odd fundamentals that can be computed from two fundamentals.

The A * (U, V)-set, with U, V ⊂ N, contains all odd fundamentals which can be obtained from every Aoperations on each pair (u, v) ∈ U × V:

A * (U, V) := (u,v)∈U×V A * (u, v) . ( 25
)
Let X be a set of integers, concisely A * (X, X) will be written A * (X), hence:

A * (X) = A q (u, v) | u, v ∈ X, q an A-configuration . ( 26
)
Example 3. For c max = 2 4 we have: It should be noted that computing A * ({1}) with c max equal to a power of 2 is recurrent in SCM and MCM problems. This is the first step to compute

A * ({1}) = A q (1, 1) | q an A-configuration = {1, 3, 5, 7, 9, 15} .
A * • A * • • • • • A * ({1}).
Let A s be the set of odd fundamentals that can be computed in depth s of the tree hypergraph, e. g. for an adder depth that is bounded by four, which is fine in practice, A 4 = A * (A * (A * (A * ({1})))) corresponds to the nodes of the hypergraph. A s are computed using A-sets:

A 0 := {1} (27)
A s := A * A s-1 ∀s ∈ N * (28)
The size A s increases quite fast with the wordlength, until all representable odd fundamentals are produced. Kumm gave the size of A-sets for different adder depth and wordlengths [START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF]. For wordlengths less to 16 an adder depth of four is sufficient to compute all the possible odd fundamentals. These results were obtained by explicitly compute A s for different wordlengths: There is no known simple formula for card(A s ). The size of the generated hypergraph can, and should, be reduced before using the solver: S s ⊆ A s denotes the odd fundamentals that may be useful to compute the odd target constants of the given instance, T odd . The set T s denotes the set of triplets (u, v, w) such that u, v and w are in an A-relation, u ≤ v and u, v ∈ A s . The set S s and T s are both computed recursively with S the maximum possible adder depth:

S S := T odd (29) T s-1 := (u, v, w) | w = A q (u, v) , u, v ∈ A s-1 , u ≤ v, w ∈ S s ∀s ∈ N * (30) S s-1 := u, v | (u, v, w) ∈ T s-1 ∀s ∈ N * (31)
ILP Formulation 2 can be derived from these constants. Its variables are presented Table 2.

The objective is to minimize the number of adder that are used, i. e. the objective is to minimize the number of variables a s w that are equal to one: Constraints:

r S w + a S w = 1 ∀w ∈ T odd (C2.1)
r s w = 0 ∀s ∈ [[1, S -1]], w ∈ S s \ s-1 s =0 S s (C2.2) r s w -a s-1 w -r s-1 w ≤ 0 ∀s ∈ [[2, S ]], w ∈ S s (C2.3) a s w - u,v∈S s-1 (u,v,w)∈T s x s-1 (u,v) ≤ 0 ∀s ∈ [[2, S ]], w ∈ S s (C2.4) x s (u,v) -r s u -a s u ≤ 0 ∀s ∈ [[1, S -1]], u, v ∈ S s-1 , u ≤ v (C2.5) x s (u,v) -r s v -a s v ≤ 0 ∀s ∈ [[1, S -1]], u, v ∈ S s-1 , u ≤ v (C2.6)
Constraint (C2.1) ensures that, for each target constant, an adder or a register computes it. Constraint (C2.2) initializes registers. Constraint (C2.3) avoids that a register takes a value at a stage if no register or adder already contained this value at the previous stage. In order for an adder to take a specific value w, then u and v such that A * (u, v) = w have to be available at previous stage; to be available at a given stage, u and v have to be in a register or an adder: constraints (C2.4), (C2.5) and (C2.6) ensure that.

Design of filters

Design of filter is usually done in three steps [START_REF] Leich | Toolbox for the design of IIR digital filters[END_REF]: First, the filter coefficients are designed using floatingpoint arithmetic; second, the quantization process is applied, roughly it corresponds to the rounding of computed coefficients; third, the implementation minimizing the hardware cost. Volkova et al. [START_REF] Kumm | Design of Optimal Multiplierless FIR Filters[END_REF] proposed another approach that combined all the steps into a single process. Figure 5 illustrates this change of approach. For FIR filters it was recently shown that the whole block can be expressed through the form of a single ILP model. To do so, the filter design which corresponds to the filter specification constraints, is conveyed by the zero-phase frequency response [START_REF] Antoniou | Digital Filters: Analysis, Design, and Signal Processing Applications[END_REF]: Filter type Filter order Symmetry c m

H R (ω) = M-1 m=0 h m c m (ω), (33)
I Even Symmetric c m (ω) = 1 if m = 0 2 cos(ωm) if m ≥ 1 II Odd Symmetric c m (ω) = 2 cos ω m + 1 2 III Even Asymmetric c m (ω) = 2 sin(ω (m + 1)) IV Odd Asymmetric c m (ω) = 2 sin ω m + 1 2 with ω ∈ [0, π],
M the number of independent coefficients and c m trigonometric functions defined Table 3 that depend on filter order and its potential symmetry. The absolute value of the zero-phase frequency response is equal to the modulus of the transfer function on the unit circle, i. e.,

|H R (ω)| = H e iω , ω ∈ [0, π] . (34)
Coefficients h m 's are target constants to optimize, this is done using MCM through its ILP formulations.

For the MCM problem, target constants were input data. On the other hand, they are not a priori known in this one-block approach. To link together filter specifications and MCM problem, a few constraints are added. This is the glue in-between. In a way, Figure 6 represents this model. To design FIR filters using this approach an open-source tool, FIRopt, has been proposed [START_REF] Kumm | FIRopt[END_REF]. We aim to propose a similar method and tool for IIR filters which, for now, are designed using the same usual three steps process [START_REF] Leich | Toolbox for the design of IIR digital filters[END_REF].

Multiple Constant Multiplication

It's still magic even if you know how it's done.

Terry Pratchett

In this section, MCM problem is studied in several ways: ILP formulations, size of models, possibility of NLP models, etc. This allows a deeper understanding of the MCM problem and improvements of initial ILP formulations with original approaches.

Linear Models 3.1.1. Analysis of the two existing models

ILP Formulation 1 has no objective function, this has a strong impact on whether implications are efficient or not [START_REF] Klotz | Practical guidelines for solving difficult mixed integer linear programs[END_REF]. Typically, constraints with implications are dropped in the linear relaxation of the model. This affects the quality of the linear relaxation, hence, in the branch and bound process, less nodes are bounded by the linear relaxation. Note that when there is no objective function the bound process is not relevant. As we will add objective functions to the ILP Formulation 1, experimental results will show this effect of implications or big M choice.

Solving MCM problem with an adder graph or an odd adder graph is equivalent [START_REF] Dempster | Constant integer multiplication using minimum adders[END_REF]. ILP Formulation 1 uses this property by limiting possible shifts, (C1.7) and (C1.8), and ILP Formulation 2 relies on this, through A-operation. We propose to use further this property in ILP Formulation 1 by removing negative shifts from inputs and by having no more than one input with a positive shift. Hence, we will suppose that left shifts are dropped. However, to match with [19, Theorem 3] a negative shift should be allowed after the addition (or subtraction): a negative shift inside the adder. To do so, sets of new variables are needed:

• c nsh a ∈ [[0, 2 d ]], ∀a ∈ [[1, N A ]]: constant obtained in adder a before an eventual shift; • c odd a ∈ N, ∀ ∈ [[1, N A ]]
: allows forcing c a to be odd;

• Ψ a,s ∈ {0, 1}, ∀a ∈ [[1, N A ]], s ∈ [[S min , 0]]: 1 if a s shift is applied to adder a, 0 otherwise. Constraint (C1.2) is replaced by c nsh a = c
sh,sg a,l + c sh,sg a,r , then constraints c a = 2c odd a + 1 are added to force constants to be odd. Finally, a new set of constraints links c a with c nsh a and Ψ a,s :

c nsh a = 2 -s c a if Ψ a,s = 1 ∀a ∈ [[1, N A ]], s ∈ [[S min , 0]] (35) 0 s=S min Ψ a,s = 1 ∀a ∈ [[1, N A ]] (36)
or, using big M constraints,

c nsh a ≥ 2 -s c a + Ψ a,s -1 × M ∀a ∈ [[1, N A ]], s ∈ [[S min , 0]] (37)
c nsh a ≤ 2 -s c a + 1 -Ψ a,s × M ∀a ∈ [[1, N A ]], s ∈ [[S min , 0]] ( 38
) 0 s=S min Ψ a,s = 1 ∀a ∈ [[1, N A ]] (39)
This big M is bounded by 2 2d and can be adjusted to 2 d adding more variables and constraints, as detailed in Section 3.1.2.

It can be noticed that a shift inside an adder only makes sense if the sum of the inputs is even. Yet, the left input is odd since left shifts are dropped. Thus, the sum of the inputs is even if and only if the right input is odd. Hence, the sum of the inputs is even if and only if there is no right shift. From that, a constraint can be deduced to speed up the resolution:

φ a,r,0 = -1 s=S min Ψ a,s ∀a ∈ [[1, N A ]]. ( 40
)
ILP Formulation 1 modified with previous variables and constraints will be referred to as ILP Formulation 1 with odd fundamentals only, or abbreviated ILP Formulation 1odd.

Example 4. The target set {7, 19, 31} can be obtained with three adders. First, 7 and 31 are computed as a power of two (a left shift of one) minus one. Second, 19 corresponds to 7 + 31 with a negative shift inside the adder. The variable c nsh a would be equal to 38 and the inside shift Ψ a,s would be equal to one for s = -1, thus c a would be equal to 38 << -1 = 19.

Despite the combinatorial dimension, ILP Formulation 1 (and ILP Formulation 1odd) remains relatively small: the number of variables and constraints is quadratic on the number of adders, the number of outputs and the wordlength. On the other hand, ILP Formulation 2 quickly becomes huge. In a way, ILP Formulation 2 uses the whole search space of ILP Formulation 1 as its input data. Hence, the complexity of the model is passed from the search space to the size of the model. Kumm reported the size of A s for a few s and increasing wordlengths [3, Section 5.5]. For s = 4, A-sets become quickly too big and long to be computed and, for wordlengths greater than 22, A 3 already has more than a million elements. This rapidly constrains the adder depth and/or the wordlength of this formulation. However, in most practical uses a relatively small wordlength and an adder depth of three or four is enough or a material constraint.

As mentioned before, this formulation is a shortest-path in a particular finite and relatively small hypergraph while ILP Formulation 1 has a search space that is a much bigger hypergraph and which is not precomputed. Thus, for most instances, ILP Formulation 2 is supposed to perform better and we confirm this with experimental results.

Fine-tuning of big M constraints

Given a big M constraint, the better big M is the smallest one. The purpose of a big M is to inactivate a constraint, if it is too small it might not inactivate the constraint and rather interfere with values that some variables could have taken. Let a i ∈ a, a , b i ∈ b, b and c i ∈ {0, 1} be three sets of variables. Roughly, for a constraint

a i ≤ b i + c i M, we have, M = a -b. (41)
We will see that, in ILP Formulation 1, a few big M's need to take values up to 2 2d . This has a direct impact on for which wordlengths numerical instability could occur. By fine-tuning the maximum value of big M's by modifying and adding constraints we show that it allows a wider range of possible wordlengths to be used. For this analysis we will restrain shifts to positives only, s ≥ 0. This analysis naturally extends for negative shifts. First, it is straightforward that M = 2 d is big enough for constraints (C1.3a), (C1.3b), (C1.9a), (C1.9b), (C1.10a), (C1.10b), (C1.12a) and (C1.12b). However, for constraints (C1.5a) and (C1.5b) different M's are needed. For the latter M = 2 d -2 s is big enough. Per contra, for (C1.5a) M = 2 d+s is needed. Thus, for s = d, M have to take values up to 2 2d . Default integrality tolerance is 10 -5 and can be safely adjusted down to 10 -9 for CPLEX and Gurobi. This parameter will be denoted IntTol and used to deduce which maximum value the wordlength can take depending on this tolerance. Hence, for:

φ a,i,s ∈ [1 -IntTol, 1 + IntTol] , (42)
we have:

c sh a,i ∈ 2 s c a,i -IntTol × 2 2d , 2 s c a,i + IntTol × 2 2d . ( 43
)
It is clear that a problem can occur only if IntTol×2 2d ≥ 1. Thus, it is imperative that IntTol -1 > 2 2d . Hence, it is direct that d = 8 and d = 14 are the longest possible wordlengths for, respectively, default integrality tolerance, 10 -5 , and smallest integrality tolerance, 10 -9 . Those wordlengths are smalls and it would be interesting to modify big M constraints in order to increase possible wordlength values. Switching between indicators and big M constraints in a same model when M becomes too big could be a solution to avoid numerical instabilities. Nevertheless, a fully linear programming model can be obtained with lower big M. To this end, new sets of variables are needed:

• c a,i,s ∈ [[0, 2 d ]], a ∈ [[1, N A ]], i ∈ {l, r}, s ∈ [[0, S max ]]: 0 if 2 s c a > 2 d , constant c a shifted by s otherwise; • ψ a,i,s ∈ {0, 1}, a ∈ [[1, N A ]], i ∈ {l, r}, s ∈ [[1, S max ]]: 1 if 2 s c a > 2 d .
Constraints (C1.5a) and (C1.5b) can be replaced by c sh a,i ≥ c a,i,s + φ a,i,s -1 × M and c sh a,i ≤ c a,i,s + 1 -φ a,i,s × M.

In those constraints, M = 2 d is enough. To link c a,i,s and ψ a,i,s with the whole model, a few constraints are added:

2 × c a,i,s-1 ≤ 2 d + ψ a,i,s × 2 d ∀a ∈ [[1, N A ]], i ∈ {l, r} , s ∈ [[1, S max ]] (C1.5-1) c a,i,0 = c a,i ∀a ∈ [[1, N A ]], i ∈ {l, r} (C1.5-2) c a,i,s ≤ 2 × c a,i,s-1 ∀a ∈ [[1, N A ]], i ∈ {l, r} , s ∈ [[1, S max ]] (C1.5-3) c a,i,s ≥ 2 × c a,i,s-1 -ψ a,i,s × 2 d ∀a ∈ [[1, N A ]], i ∈ {l, r} , s ∈ [[1, S max ]] (C1.5-4) c a,i,s ≤ 1 -ψ a,i,s × 2 d ∀a ∈ [[1, N A ]], i ∈ {l, r} , s ∈ [[1, S max ]] (C1.5-5)
Constraint (C1.5-1) forces ψ a,i,s = 1 if 2 s c a,i > 2 d . Constraint (C1.5-5) ensures that c a,i,s = 0 if ψ a,i,s = 1. Constraints (C1.5-2), (C1.5-3) and (C1.5-4) imply that c a,i,s = 2 s c a,i if ψ a,i,s = 0. Note that the variable c a,i,s could be equal to 0 even if 2 s c a,i ≤ 2 d , yet this has no impact on the solution. This modification has a positive impact on possible wordlengths. Maximum wordlength becomes 16 for default integrality tolerance and 29 for an integrality tolerance of 10 -9 . The greater the integrality tolerance is, the faster is the solver. However, this modelization takes more variables and constraints, that might offset the gain. Hence, we use this modelization only when the wordlength is too big for the initial model.

Objective functions

As stated before, ILP Formulation 1 has no objective function. In order to minimize the number of adders, ILP Formulation 1 is solved multiple times incrementally increasing N A until a solution is found. This implies to try to solve multiple infeasible linear models, consequently to prove that for multiple values of N A the model is infeasible. This problem could be solved the other way around: use a greedy algorithm or an heuristic to find a feasible N A and solving the model by decrementing NA until a model is infeasible. Fixing N A to its upper bound and minimizing the number of used adders is a third possible approach. To do so, for each adder, a binary variable, that is activated when its adder is used, can be added. Then the objective function is to minimize the sum of those binary variables. Formally, that means adding u a ∈ {0, 1}, ∀a ∈ [[1, N A ]], to the model and the following set of constraints,

c a = 1 if u a = 0 ∀a ∈ [[1, N A ]], (45)
or, using big M,

c a ≤ u a × M + 1 ∀a ∈ [[1, N A ]], (46)
where M = 2 d can be used. This constraint ensures that c a brings no added value to the solution if u a = 0, it could be replaced by c 0 . Yet, this can be enhanced with a few constraints to break symmetries:

u a ≤ u a-1 ∀a ∈ [[2, N A ]] (47) a-1 k=1 c a,i,k ≤ 2u a ∀a ∈ [[1, N A ]], i ∈ {l, r} (48)
The first constraint forces adder a-1 to be used if adder a is. Thus, used adders are first indices. The second constraint implies that inputs of an unused adder are from the initial input. This modification of the model

c a,i,k ∈ {0, 1}, ∀a ∈ [[1, N A ]], i ∈ {l, r}, k ∈ [[0, N A -1]]: 1 if input i of adder a is adder k, 0 otherwise; s + a ∈ [[0, S max ]], ∀a ∈ [[1, N A ]]: shift of left input of adder a; s - a ∈ [[S min , 0]], ∀a ∈ [[1, N A ]]: negative shift of adder a; o a, j ∈ {0, 1}, ∀a ∈ [[1, N A ]], j ∈ [[1, N O ]]:
1 if adder a is equal to the j-th target constant, 0 otherwise. permits to give to the solver a first solution using heuristics or bounds found with a greedy algorithm. This can speed up the solver allowing it to use its full potential and not only its satisfiability part. However, this modeling requires more variables and its linear relaxation might be far from a known lower bound N A . To help the solver, u a should be fixed to 1 for all a ≤ N A . MCM problem with this objective will be referred to as MCM OBJ .

As already presented, a second objective which might be interesting to minimize is the adder depth. For each adder, a new integer variable, ad a , representing the adder depth can be added to the model. Then we just introduce an integer variable max ad, the constraints max ad ≥ ad a , ∀a ∈ [[1, N A ]], and we minimize max ad. In order to represent the adder depth, first ad 0 = 0 is fixed. Then a set of constraints is added:

ad a ≥ ad k + 1 -1 -c a,i,k × N A , ∀a ∈ [[1, N A ]], i ∈ {l, r} , k ∈ [[0, a -1]]. ( 49
)
Despite that N A can be seen as a big M here, N A remains relatively small and it would be ineffective to write this constraint using an implication. Note that for a given a, ad a might be greater than the actual adder depth of its associated adder, however it will be low enough so that the graph adder depth is actually minimized and its real adder depth could be computed from the solution. Adding this objective function permits to tackle MCM MAD . In order to handle MCM BAD it is enough to add a unique constraint, max ad ≤ AD, where AD is the wanted bound. Combining both previous objectives is possible with a priority on minimizing the number of adders. It is straightforward that minimizing N A × u a + max ad tackles this aim which allows handling the whole problem solving a single model. This corresponds to a weighted sum that ensures the first objective, which is to minimize the number of used adders, is optimized before the second objective, which is to minimize the adder depth.

Nonlinear model

No analysis of NLP models are given in the literature, however this seems to be a reasonable approach because of the nature of A-sets. We propose two models to evaluate the potential of NLP models. The first one checks for satisfiability while the second minimizes the number of used adders. Constants and variables that are used in the first model are presented Table 4.

Constraints for first NLP model:

c 0 = 1 (C3.1) 2 s - a c a = 2 s + a a-1 k=0 c k c a,l,k + (1 -2σ a ) a-1 k=0 c k c a,r,k ∀a ∈ [[1, N A ]] (C3.2) a-1 k=0 c a,i,k = 1 ∀a ∈ [[1, N A ]], i ∈ {l, r} (C3.3) s - a s + a = 0 ∀a ∈ [[1, N A ]] (C3.4) s + a -s - a ≥ 1 ∀a ∈ [[1, N A ]] (C3.5) N A a=1 c a o a, j = C j ∀ j ∈ [[1, N O ]] (C3.6) N A a=1 o a, j = 1 ∀ j ∈ [[1, N O ]] (C3.7)
Constraints (C3.2) is equivalent to the A-operation [START_REF] Kumm | FIRopt[END_REF] where a shift has been removed, where for a unique k, u = c k c a,l,k and for a unique k , v = c k c a,r,k . For a given pair of (a, i), only one k permits c k c a,i,k 0 because of constraint (C3.3). Furthermore, (-1) s has been replaced with 1 -2s because we think, a priori, that NLP solvers might handle the latter better. Constraint (C3.4) and (C3.5) ensure that either the input shift or the adder shift is equal to zero and one is not zero. Last constraints, (C3.6) and (C3.7), permit to associate adder values to output values and ensure that each output as an associated adder.

This first decision NLP model can be extended to obtain a minimization model using the same method as for the ILP approach. To do so, we start by adding variables u a ∈ {0, 1}, ∀a ∈ [[1, N A ]] and an objective to minimize: min N A a=1 u a . Then we add a few constraints that ensure that an unused adder cannot be used as an input for another adder nor associated to an output. Furthermore, the value of an unused adder is fixed to one and adders are sorted: unused adders at the tail. Formally, this is done with four constraints:

(1 -u a ) × (c a -1) = 0 ∀a ∈ [[1, N A ]] (50) u a ≤ u a-1 ∀a ∈ [[2, N A ]] (51) N O j=1 o a, j ≤ u a ∀a ∈ [[1, N A ]] (52) c a,i,k ≤ u k ∀a ∈ [[1, N A ]], k ∈ [[1, a -1]], i ∈ {l, r} (53)
This gives two NLP models to compare with. Although results are not expected to be promising, MCM problem can be combined with filter design which naturally comes with nonlinear constraints and having a NLP formulation for MCM can become convenient.

Experimental results

In this section we will compare solving times for different models. Two modeling languages were considered. First, ScaLP [START_REF] Sittel | ScaLP: A Light-Weighted (MI)LP Library[END_REF] is the modeling language used by Kumm and it was natural to use the same since many models he proposed are open source. However this language does not handle yet nonlinear terms. This limitation was a real issue since part of our objective was to compare the efficiency of nonlinear models with linear models. The other modeling language considered, which is the one we used, is JuMP [START_REF] Dunning | JuMP: A Modeling Language for Mathematical Optimization[END_REF].

gaussian 3x3 8bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 gaussian 5x5 12bit T O < 1 18 431 1285 T O 9 84 88 484 1745 highpass 5x5 8bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 highpass 9x9 10bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 highpass 15x15 12bit < 1 53 T O 1251 13 47 154 < 1 771 6 123 laplacian 3x3 8bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 lowpass 5x5 8bit < 1 13 7 3 12 T O 2 4 14 2 51 lowpass 9x9 10bit < 1 T O T O T O T O T O 40 18 OOM 473 OOM lowpass 15x15 12bit T O T O T O T O T O T O T O T O OOM T O OOM unsharp 3x3 8bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 unsharp 3x3 12bit 1301 10 5 60 21 T O 2 2 9 3 145 gaussian 3x3 8bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 gaussian 5x5 12bit 1758 244 54 756 2706 T O 4 2 368 701 44 highpass 5x5 8bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 highpass 9x9 10bit < 1 < 1 < 1 < 1 < 1 26 < 1 < 1 < 1 < 1 < 1 highpass 15x15 12bit < 1 < 1 T O T O T O T O < 1 < 1 53 2 2 laplacian 3x3 8bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 lowpass 5x5 8bit < 1 6 11 156 70 7 2 2 85 2 4 lowpass 9x9 10bit < 1 T O T O T O T O T O 19 74 T O 1169 T O lowpass 15x15 12bit T O T O T O T O T O T O T O T O T O T O T O unsharp 3x3 8bit < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 < 1 unsharp 3x3 12bit 849 56 157 345 105 T O 3 < 1 21 2 6
Table 7: Runtime comparison in seconds of ILP Formulation 2 (ILP2), ILP Formulation 1 (ILP1), 1odd (Odd) with big M constraints and the NLP models, NLP corresponds to the decision one when NLP Min corresponds to the minimization form. TO is for Time Out (3600s). Experiments were made using SCIP.

ILP2 [START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF] ILP1 [START_REF] Kumm | Optimal Constant Multiplication Using Integer Linear Programming[END_REF] Odd NLP NLP Min

gaussian 3x3 8bit < 1 < 1 < 1 3 2 gaussian 5x5 12bit T O 299 1042 T O T O highpass 5x5 8bit < 1 < 1 < 1 < 1 2 highpass 9x9 10bit < 1 < 1 < 1 < 1 155 highpass 15x15 12bit < 1 170 T O T O T O laplacian 3x3 8bit < 1 < 1 < 1 < 1 2 lowpass 5x5 8bit < 1 4 14 151 T O lowpass 9x9 10bit < 1 2486 T O T O T O lowpass 15x15 12bit T O T O T O T O T O unsharp 3x3 8bit < 1 < 1 < 1 2 3 unsharp 3x3 12bit T O 179 38 235 T O
He has both linear and nonlinear modeling possibilities and can be plugged-in with Gurobi [START_REF]Gurobi Optimizer Reference Manual[END_REF], CPLEX [START_REF] Cplex | CPLEX User's Manual[END_REF] and SCIP [START_REF] Gleixner | The SCIP Optimization Suite 6.0, software, Optimization Online[END_REF] which are solvers used in previous works. Several models are tested, both with indicator constraints and big M constraints. When big M constraints were used, integrality tolerance has been adjusted to avoid numerical instability. Models are implemented in Julia using JuMP and solved with a time limit of one hour, 3600 seconds. CPLEX 12.9 and Gurobi 9 were used for ILP formulations and SCIP 6.0.2 for the NLP models. The computer on which the experimentation was done is on Linux 64bit (Ubuntu 18.04.4 LTS), it has 16Go of DDR3 RAM and an Intel® Core™ i5-4570S, quadcore (2.90GHz). All the experiments were performed on instances from Kumm's thesis, the detail of those instances is given Table 5. No assumptions on possible shifts were made, hence the shift's range is from minus the wordlength to the wordlength.

Table 6 gives the results in seconds for ILP Formulation 1 and its variants, ILP Formulation 2 results corresponds to the ones reported by Kumm where runtimes greater than our time limit are considered timed out [3, Table 5.4]. Some instances using big M constraints raised an Out Of Memory error with CPLEX. Since models using indicator constraints and big M constraints are of the same order of magnitude this surely comes from the way the solver optimizes its branch and bound. Models with indicator constraints timed-out more often than models with big M constraints and often took more time to be solved. Globally, models with an objective function took more time than decision models. Despite its theoretical interest, the MIN model does not improve the solving and can probably be left aside in future work, unless the use of better upper bounds on the number of adders drastically increases its efficiency. MAD model, on the other hand, has showed to be efficient enough to be used when minimizing the adder depth is important. In that case, the use of big M constraints permitted faster solving time with the drawback of great memory usage. Note that in both models using an objective function, an optimal solution could have been found fast and most of the time could have been use to confirm its optimality. We left that verification aside as it is a whole experiment in itself and it is not of immediate interest for us. BAD 3 model worked just fine compared to ILP Formulation 1 and other variants. This variant is essential since it can be fairly compared to ILP Formulation 2 and reassures that the precomputation process is of great interest in solving MCM BAD problem. As expected ILP Formulation 2 seems to behave at least as good as ILP Formulation 1. Surprisingly, the gaussian 5x5 12bit and the unsharp 3x3 12bit instances have been handled better by ILP Formulation 1 BAD 3 . ILP Formulation 1odd performed as good as ILP Formulation 1 and other variants. Except for the instance lowpass 15x15 12bit, which was left unsolved by almost every model and solver, ILP Formulation 1odd solved every instance in a reasonable time.

x n x n z -1 x n-1 z -1 x n-2 y n-2 z -1 y n-1 z -1 b 0 b 1 b 2 -a 2 -a 1 y n y n
Both NLP models obtained time results are reported Table 7. ILP Formulation 2 corresponds to results reported by Kumm where runtimes greater than our time limit are considered timed out [3, Table 5.4]. On a few instances that most models handled well, NLP models shown to work too. However, NLP models timed out on instances that were solved relatively fast by ILP formulations. We conclude that ILP formulations should probably be preferred over NLP formulations. No clear result emerge from the comparison between variants and indicator/big M constraints. An experimentation on more instances and an analysis of these instances using Data Science methods seem useful to help understanding when and why a variant prevails on another. When a formulation for MCM problem will be needed in the following, we will prefer ILP Formulation 1odd with big M constraints as it seems more robust: It gave the best results on most of the instances and not too far from the best on the others, regardless of the solver.

Design of second order IIR filters

Study hard what interests you the most in the most undisciplined, irreverent and original manner possible.

Richard Feynman

We extended the design process of FIR filters to second order IIR filters. Those filters have a strongly different behavior and all the coefficients cannot be included into a single MCM problem. Indeed, in Figure 7 it can be seen that a's and b's have different inputs. Second order IIR filters are an essential block of IIR filter design since many methods rely on second order sections (SOS) to design greater order filters [START_REF] Vanuytsel | Efficient hybrid optimization of fixed-point cascaded IIR filter coefficients[END_REF][START_REF] Hourani | Hardware implementation of IIR digital filters for programmable devices[END_REF]. The transfer function of second order filters is given below:

H(z) = b(z) a(z) = b 0 + b 1 z -1 + b 2 z -2 1 + a 1 z -1 + a 2 z -2 with a 1 , a 2 , b 0 , b 1 , b 2 ∈ R, z ∈ C. (54)
In order to have the desirable design in the meaningful frequency domain, |H(z)| has to be constrained on the unit circle, i. e. z = e iω :

D(ω) ≤ H e iω ≤ D(ω) , ∀ω ∈ Ω. (55)
We obtained an explicit expression with respect to the filter parameters by expanding H e iω 2 :

H e iω 2 = b 0 2 + b 1 2 + b 2 2 + 2b 0 b 1 cos(ω) + 2b 0 b 2 cos(2ω) + 2b 1 b 2 cos(ω) 1 + a 1 2 + a 2 2 + 2a 1 cos(ω) + 2a 2 cos(2ω) + 2a 1 a 2 cos(ω) . (56)
Details are given in Appendix A. This expansion involves quadratic terms that do not fit for an ILP model. Furthermore, unlike FIR filters, IIR filters have stability issues. Those issues are addressed in the next sections and solved. It should be noted that,

b e iω 2 = b 0 2 + b 1 2 + b 2 2 + 2b 0 b 1 cos(ω) + 2b 0 b 2 cos(2ω) + 2b 1 b 2 cos(ω) , (57)

Linearize quadratic terms

In order to avoid quadratic terms as z = xy, with x, y, z ∈ N, a method was presented in 2008 by Billionnet et al. [START_REF] Billionnet | Linear Reformulations of Integer Quadratic Programs[END_REF]. Let d ∈ N be a positive integer, under the assumption of 0 ≤ x ≤ 2 d and 0 ≤ y ≤ 2 d , it suffices to introduce binary variables t y,i for i ∈ [[0, d]] and the constraint y = d i=0 2 i t y,i . Then xy can be rewrite as d multiplications between a binary variable and an integer variable which can be linearized easily. Hence z = xy will become the following set of constraints: . This means replacing x and y with, respectively, x + and y + in previous constraints and replacing (62) with the following set of inequalities where x sg and y sg are binary variables equal to 1 if x > 0 and y > 0 respectively, Figure 8: a's that ensure stability. a 1 is the abscissa and a 2 the ordinate. In blue (bottom), poles are both real valued, in green (top), poles are complex conjugates, in red a single double pole. and equal to 0 if x < 0 and y < 0 respectively:

y = d i=0 2 i t y,i (58)
u x,y,i ≤ M × t y,i ∀i ∈ [[0, d]] (59)
u x,y,i ≤ x ∀i ∈ [[0, d]] (60)
u x,y,i ≥ x -M 1 -t y,i ∀i ∈ [[0, d]] (61)
z = d i=0 2 i u x,y,i , (62)
z ≥ d i=0 2 i u x,y,i -2 × M x sg + y sg (63) z ≥ d i=0 2 i u x,y,i -2 × M 2 -x sg -y sg (64) z ≤ d i=0 2 i u x,y,i (65)
z ≤ -

d i=0 2 i u x,y,i + 2 × M 1 + x sg -y sg (66) z ≤ - d i=0 2 i u x,y,i + 2 × M 1 -x sg + y sg (67) z ≥ - d i=0 2 i u x,y,i (68)
Here, M should be replaced by the upper bound of x + times the upper bound of y + . M = 2 2d is good enough.

Stability

Stability is an issue for IIR filters. An IIR filter is stable if and only if its poles are located inside the unit circle. For a second order IIR filter there is two poles to consider, either both real, complex conjugate or one double root. To do so, three cases are considered: ∆ < 0, ∆ > 0 and ∆ = 0, with ∆ = a 1 2 -4a 2 . From those three cases, we deduced constraints on a 1 and a 2 that ensure stability.

Theorem 1 (Stability). Given F a second order IIR filter, F is stable if, and only if, a 1 ∈ ]-2, 2[ and a

2 ∈ ]|a 1 | -1, 1[.
Proof. We recall that stability ( 9) is equivalent to, for all λ, poles of the filter, |λ| < 1. For second order filters poles are the roots of the quadratic equation:

A(z) = z 2 + a 1 z + a 2 .
(69)

We have three cases to consider that depend on the sign of ∆ with:

∆ = a 1 2 -4a 2 . (70)
First, we suppose that ∆ < 0, thus that the roots of the polynomial are complex conjugates. The roots, z 1 and z 2 , of the polynomial A are equal to:

z 1 = z 2 = 1 2 -a 1 + i √ -∆ . (71)
We make the hypothesis:

∆ < 0 |z 1 | < 1 ∧ |z 2 | < 1 , (72)
which is equivalent to:

         a 1 2 4 < a 2 a 1 2 + -a 1 2 + 4a 2 2 < 2 ⇔ 0 < a 1 2 + -a 1 2 + 4a 2 < 4 ⇔ 0 < a 2 < 1 . (73)
We have obtained an equivalent and simpler set of constraints to (72): a 2 ∈ a 1 2 4 , 1 . In order to guaranty the existence of a valid value for a 2 we obtain bounds on a 1 too: a 1 ∈ ]-2, 2[. Second, we suppose that ∆ > 0, thus that the roots of the polynomial are real valued. The roots, z 1 and z 2 , of the polynomial A are equal to:

z 1 = 1 2 -a 1 + √ ∆ and z 2 = 1 2 -a 1 - √ ∆ . ( 74
)
We make the hypothesis:

∆ < 0 |z 1 | < 1 ∧ |z 2 | < 1 . (75)
By noting that for u ∈ R and v, w ∈ R + , |u ± v| < w is equivalent to |u| + v < w we deduce that:

|z 1 | < 1 ∧ |z 2 | < 1 ⇔ 1 2 |a 1 | + √ ∆ < 1. (76)
Then, System (75) is equivalent to:

       a 2 < a 1 2 4 a 1 2 -4a 2 < 2 -|a 1 | ⇔ 0 < a 1 2 -4a 2 < a 1 2 -4 |a 1 | + 4 ⇔ |a 1 | -1 < a 2 < a 1 2 4 . (77)
We have obtained an equivalent and simpler set of constraints to (75): a 2 ∈ |a 1 | -1, a 1 2 4 . Furthermore, in system (77) we had 0

< √ ∆ < 2 -|a 1 |, hence a 1 ∈ ]-2, 2[.
Third case, we suppose ∆ = 0, thus that we have a single root, z 1 , equal to:

z 1 = - 1 2 a 1 . (78)
Constraining its absolute value leads to:

∆ = 0 |z 1 | < 1 ⇔ a 1 2 -4a 2 = 0 |a 1 | < 2 ⇔      a 1 ∈ ]-2, 2[ a 2 = a 1 2 4 . (79)
Finally, we reunify the three cases into one set of constraints:

Second-order filter stability ⇔

                     a 2 ∈ |a 1 | -1, a 1 2 4 a 2 = a 1 2 4 a 2 ∈ a 1 2 4 , 1                  ⇔ a 2 ∈ ]|a 1 | -1, 1[ a 1 ∈ ]-2, 2[ . ( 80
)
Possible values for a 1 and a 2 are represented Figure 8 depending on poles. Theorem 1 gives bounds on a's and a single constraint, a 2 > |a 1 | -1, that ensure stability for a second order IIR filter. In fixed-point arithmetic the absolute value can be linearize without much effort, actually to linearize the bilinear terms, |a 1 | must be computed anyway.

Bounds on coefficients

As explained in Section 4.1, bounds on coefficients a's and b's are necessary for the linearization part. From Theorem 1, bounds on a 1 and a 2 are directly deduced: a 1 ∈ ]-2, 2[ and a 2 ∈ ]-1, 1[. Bounds on b's are left to be founded. Coefficients a's and b's satisfy the following constraint:

D(ω) 2 ≤ H e iω 2 = b e iω 2 a e iω 2 ≤ D(ω) 2 , ∀ω ∈ Ω. (81)
This links b e iω 2 and a e iω 2 together, where,

b e iω 2 = b 0 2 + b 1 2 + b 2 2 + 2b 0 b 1 cos(ω) + 2b 0 b 2 cos(2ω) + 2b 1 b 2 cos(ω) , ∀ω ∈ Ω, (82)
a e iω 2 = 1 + a 1 2 + a 2 2 + 2a 1 cos(ω) + 2a 2 cos(2ω) + 2a 1 a 2 cos(ω) , ∀ω ∈ Ω. ( 83
)
We use the bounds on a's to find bounds on a e iω 2 using (83) then we get bounds on b e iω 2 using (81). As a square, 0 is a direct lower bound for a e iω 2 . Furthermore, we found that this bound cannot be increased: lim

a 1 →0 a 2 →1 ω→ π 2 a e iω 2 = 0. ( 84
)
Studying the monotony of terms of a e iω 2 it is straightforward that a e iω 2 ≤ 16. This bound of 16 cannot be lowered since: lim

a 1 →2 a 2 →1 ω→0 a e iω 2 = 16. ( 85
)
From these bounds it can be deduced bounds on b e iω 2 :

D(ω) 2 ≤ b e iω 2 a e iω 2 ≤ D(ω) 2 , ∀ω ∈ Ω, (86)
⇒ a e iω 2 D(ω

) 2 ≤ b e iω 2 ≤ a e iω 2 D(ω) 2 , ∀ω ∈ Ω, (87)
⇒0 ≤ b e iω 2 ≤ 16D(ω) 2 , ∀ω ∈ Ω. (88)
Since D(ω) is not bounded a priori, it is not possible to derive bounds on b's independent of filter specifications. However, those bounds are mandatory to linearize quadratic terms or to situate the most significant bit (MSB) for fixed-point arithmetic. Some MILP solvers (as CPLEX or Gurobi) handle quadratic terms in case those are convex problems. Problems,

min / max b i ∀i ∈ {0, 1, 2} (89)
s.t. b e iω 2 ≤ 16D(ω) 2 ∀ω ∈ Ω d (90) b 0 , b 1 , b 2 ∈ R,
are convex, thus we can use CPLEX or Gurobi to find bounds on each b i before implementing the whole model. The following models have to be solved in order to obtain lower and upper bounds on b i 's denoted b i and b i respectively:

min / max b i ∀i ∈ {0, 1, 2} (91) s.t. b 0 2 + b 1 2 + b 2 2 + 2b 0 b 1 cos(ω) + 2b 0 b 2 cos(2ω) + 2b 1 b 2 cos(ω) ≤ 16D(ω) 2 ∀ω ∈ Ω d (92) b 0 , b 1 , b 2 ∈ R.
Knowing these bounds allows obtaining a linear model for the design of second-order IIR filters: The whole design process is summarized in the next section.

ILP Model

In fixed-point context, for this first ILP model, it is easier to suppose that variables a's and b's are in [-1, 1] and use the whole interval so that their integer counterpart is in [[-2 d , 2 d ]] where d is the wordlength. This was already mentioned in Section 2.1 with [START_REF] Prandoni | Signal Processing for Communications[END_REF]. Going forward, integer coefficients will be used more, hence real coefficients a's and b's will be denoted a 's and b 's in the following. Bounds, a i , a i , b i , b i , on a 's and b 's can be formally obtained or precomputed. Normalized integer coefficients a's and b's are deduced from a 's, b 's and their bounds:

a i = a i × 2 d 2 log 2 max{-a 1 ,-a 2 ,a 1 ,a 2 } ∀i ∈ {1, 2} (93)
b i = b i × 2 d 2 log 2 max -b 0 ,-b 1 ,-b 2 ,b 0 ,b 1 ,b 2 ∀i ∈ {0, 1, 2} (94)
D ← D × 2 log 2 max{-a 1 ,-a 2 ,a 1 ,a 2 } 2 log 2 max -b 0 ,-b 1 ,-b 2 ,b 0 ,b 1 ,b 2 (95) D ← D × 2 log 2 max{-a 1 ,-a 2 ,a 1 ,a 2 } 2 log 2 max -b 0 ,-b 1 ,-b 2 ,b 0 ,b 1 ,b 2 (96)
The model for second order IIR filter is available in Appendix B. First, there is the filter specifications constraints and the stability constraints. Second, constraints for the linearization of quadratic terms containing b's then containing a's are added. The constraints from ILP Formulation 1odd for MCM problems follow. Finally, this model has constraints that link MCM problems, one for a's and one for b's, with filter coefficients. This whole process is presented in a high-level representation Figure 9. Variables and their domain are detailed at the end of the model.

Glue constraints, detailed in Figure B.16 in Appendix B, state that the input is both available in MCM problems for a's and b's coefficients. Then a few constraints ensure that variable samemcm k,k is equal to 1 if and only if adders k and k are both available in MCM problem for a's coefficients or both available in MCM problem for b's coefficients. A constraint prevents adder k to be an input of adder a if both adders are not available for the same coefficients (a's or b's). Finally, a few constraints link outputs of MCM with filter coefficients. As already pointed out, the function b e iω 2 as two symmetries we can break. The constraint b 0 ≥ |b 2 | allows this and is easy to implement since |b 2 | is already computed in the linearization process.

As Volkova et al. [START_REF] Kumm | Design of Optimal Multiplierless FIR Filters[END_REF] did for FIR filter coefficients, we first solve the problem without MCM constraints with the objectives min / max a i and min / max b i in order to find integer bounds on coefficients. This allows reducing the coefficient range, hence the search space. Furthermore, we integrate in the final solving process an objective function which indicates to maximize the number of coefficients with a value of zero.

Experimental results

Comparison of our method with existing ones is not possible in a classical way since there are no benchmarks for second order IIR filters. Although second order IIR filters are an essential block in many methods [START_REF] Hourani | Hardware implementation of IIR digital filters for programmable devices[END_REF][START_REF] Smékal | Optimized models of IIR digital filters for fixed-point digital signal processor[END_REF][START_REF] Milić | Design of Multiplierless Elliptic IIR Filters with a Small Quantization Error[END_REF], they are part of algorithms that cannot be applied in our single block approach.

Table 9: Runtimes in seconds for the instances used to test our second order IIR filter design method. The first few columns corresponds to the name of the instance, the total time of the solving process (Total), the part of total time used by the precomputation (Pre) and the solving time without any precomputation (NoPre). Total a and Pre a correspond to the solving times with precomputation of coefficient a's then b's, Total b and Pre b correspond to the opposite. The number of adders used is denoted by NA. Then coefficient values b's and a's are given followed by adder values subscripted with a b or an a to indicate if it is used for the b or the a coefficients respectively. The symmetry breaking constraint was used (top) and inactivated (bottom). Experiments were conducted using CPLEX.

Name

Total However, lowpass and highpass filters that can be handle with second order IIR filters are proposed Table 8. Their solving times are given Table 9. In this last table we also give the number of adder used for their implementation, the value of these adders and the coefficient values we found. This experiment was done on the same computer as for MCM problems: Linux 64bit (Ubuntu 18.04.4 LTS), with 16Go of DDR3 RAM and an Intel® Core™ i5-4570S, quadcore (2.90GHz).

The instances used for those experiments are initially FIR filter instances that we manually simplified to be handled with a second order IIR filter. The obtained results are found in a reasonable time and are certified optimal for multiplierless second order IIR filters. The transfer functions we found with our method for lowpass 10bit and highpass 8bit are represented Figure 10. Our experiments show that for a filter with five coefficients we are able to obtain hardware implementations that only involve a few adders, down to a single one. Furthermore a filter coefficient was even fixed to zero for lowpass 10bit. Surprisingly the precomputation of coefficient bounds is essential for highpass 8bit while it should be avoided for lowpass 10bit. The extreme difference of utility of the precomputation does not permit to conclude anything on its usage and more research on this aspect is needed. The use of the symmetry breaking constraint allows speeding up computation, with and without precomputations. Finally, it has been shown that precomputing a coefficient ranges or b coefficient ranges first has an impact on the precomputation time. However, it is not clear that a choice prevails on the other: For lowpass 8bit and highpass 8bit a's before b's was faster with and without the symmetry breaking constraint, however for lowpass 10bit without the symmetry breaking constraint the solving time was divided by more than three. As the symmetry breaking constraint seems to give better results and should be kept in the final method we find that precomputing a coefficient ranges first is more efficient.

Conclusion and perspectives

Reserve your right to think, for even to think wrongly is better than not to think at all.

Hypatia of Alexandra

We proposed a method to design optimal second order IIR filters with respect to the number of adders used in their hardware implementation. Our method is based on ILP modeling but not only as it makes use of the quadratic capabilities of MILP solvers. Results are promising: For an instance it has been shown that only one adder suffices to implement in hardware a reasonable second order IIR filter. We tackled the stability issue by obtaining constraints that allows ensuring stability when designing second order IIR filter coefficients in fixed-point arithmetic. We limited the search space by finding bounds on coefficients and we linked the filter specifications with the hardware implementation optimization using two MCM models.

The MCM problem is an important basic block of the design of second order IIR filters. We proposed a NLP formulation as an alternative to ILP formulations and we proposed efficient variants to known ILP formulations. Our variants have shown to be competitive with the literature and future work on this would be to understand when and why a variant behave greater than the others. We left aside the complexity analysis of MCM problem, yet this should probably be tackled since the conjectured complexity might no be applicable for practical variants. Thus, it is not clear that faster solving algorithms are not possible. In either case, we think that a dedicated branch and bound for this problem as a great potential and could be of great use particularly for harder instances.

We extended a design process to second order IIR filters. Despite this is an essential step, being restricted to second order is limiting and our future work on this would be to propose an algorithm that decomposes a greater order filter into second order sections. Such algorithms already exist for the classical design of filters and we are confident that we can do the same in the context of our variant.

a i ≤

  Intervals are denoted by x ∈ [a, b], ]a, b], [a, b[, ]a, b[ and correspond to a ≤ x ≤ b, a < x ≤ b, a ≤ x < b and a < x < b respectively. The integer interval is denoted by [[a, b]] := [a, b] ∩ N, with a, b ∈ R. The notation |a| corresponds to the modulus if a ∈ C \ R and to the absolute value if a ∈ R. Sequences of u n 's with n the index which takes its values in X are denoted by (u n ) n∈X . The set of natural numbers is denoted by N. The set N * corresponds to N \ {0}.

Figure 1 :

 1 Figure 1: Example of a Low-Pass filter. In red its filter specifications moved by the gain.

Figure 2 :

 2 Figure 2: Shifts are represented with labeled arrows and adders as with an optional minus on the left or the right side.

Figure 5 :

 5 Figure 5: Filter design for hardware, classic flow (left) and proposed flow (right)

Figure 6 :

 6 Figure 6: High-level overview of FIRopt ILP model Table3: Trigonometric functions c m depend on the parity of the filter order and its symmetry. Both are deduced from the filter type.

Figure 7 :

 7 Figure 7: Direct form realization of a second-order IIR filter.

  has two symmetries: The values of the variables b 0 and b 2 could be exchanged and the values (b 0 , b 1 , b 2 ) could be sign reversed to (-b 0 , -b 1 , -b 2 ).

  with u x,y,i ∈ N for i ∈ [[0, d]] and with M equals to the upper bound of x. In Section 4.3 we will show how to get bounds a 1 , a 2 , b 0 , b 1 , b 2 , therefore M = 2 d meets the objective after normalization[START_REF] Prandoni | Signal Processing for Communications[END_REF] and this ensures that the assumptions x ≤ 2 d and y ≤ 2 d are verified in our design method. However, there is no certitude on the sign of a 1 , a 2 , b 0 , b 1 , b 2 , thus a sign variable is added for each coefficient: a sg i ∈ {0, 1} and b sg i ∈ {0, 1}. Given a variable x, the constraints x ≥ x (1x sg ) and x ≤ xx sg ensure that x sg = 0 if x < 0 and x sg = 1 if x > 0. Then variables a + i = |a i | and b + i = |b i | can be introduced using a sg i and b sg i . Consequently, quadratic products are linearized using positives a + i and b + i , then the sign is adjusted with respect to a sg i and b sg i

Figure 10 :

 10 Figure 10: Examples of obtained second order IIR filters

Figure B. 13 :

 13 Figure B.13: Constraints to linearize quadratic terms containing b's

Figure B. 14 := 1 (nsh a ≥ 2 1 Φ=0 1 Φ=0Figure B. 16 :

 14121116 Figure B.14: Constraints to linearize quadratic terms containing a's

Figure B. 17 :

 17 Figure B.17: Variables

Table 1 :

 1 Constants (top) and variables (bottom) used for ILP Formulation 1 Constants/Variables N A ∈ N: number of adders; N O ∈ N: number of outputs; C ∈ N N O : target constants; S min , S max ∈ Z: minimum and maximum shift; d ∈ N: wordlength.

Table 2 :

 2 Constants (top) and variables (bottom) used for ILP Formulation 1 1}, ∀s ∈ [[1, S ]], u, v ∈ S: 1 if both u and v are available in depth s, 0 otherwise.

	Constants/Variables
	a s w ∈ {0, 1}, ∀s ∈ [[1, S ]], w ∈ S: 1 if node w is computed by an adder at depth s, 0 otherwise. a s w corresponds to adders; r s w ∈ {0, 1}, ∀s ∈ [[1, S ]], w ∈ S: 1 if node w was computed at a depth strictly lower than s, 0 otherwise. r s w are called registers; x s (u,v) ∈ {0,

Table 4 :

 4 Constants (top) and variables (bottom) for the first NLP model

	Constants/Variables
	N A ∈ N: number of adders;
	N O ∈ N: number of outputs;
	C ∈ N N O : target constants;
	S min , S max ∈ Z: minimum and maximum shift;
	d ∈ N: wordlength.
	c a ∈ [[1, 2 d ]], ∀a ∈ [[0, N A ]]: constant of adder a;
	σ a ∈ {0, 1}, ∀a ∈ [[1, N A ]]: equivalent to s in (23), 0 if the adder a corresponds to an addition, 1 if it
	corresponds to a subtraction;

Table 5 :

 5 Instances used for the experiments

	Name	wordlength Unique odd coefficients
	gaussian 3x3 8bit	8	3, 21, 159
	gaussian 5x5 12bit	12	1, 23, 343, 1267
	highpass 5x5 8bit	8	1, 3, 5, 7, 121
	highpass 9x9 10bit	10	1, 3, 5, 7, 11, 125
	highpass 15x15 12bit	12	1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 507
	laplacian 3x3 8bit	8	5, 21, 107
	lowpass 5x5 8bit	8	11, 33, 35, 53, 103
	lowpass 9x9 10bit	10	1, 5, 7, 25, 31, 63, 65, 67, 73, 97, 117, 165, 303
	lowpass 15x15 12bit	12	1, 5, 7, 13, 17, 19, 21, 27, 41, 43, 45, 53, 61, 79,
			93, 101, 103, 113, 133, 137, 199, 331, 333, 613,
			1097, 1197
	unsharp 3x3 8bit	8	3, 11, 63
	unsharp 3x3 12bit	12	43, 171, 1109

Table 6 :

 6 Runtime comparison in seconds of ILP Formulation 2 (ILP2) and different formulations derived from ILP Formulation 1 (ILP1). Odd corresponds to ILP Formulation 1odd, MAD is the minimization of adder depth modification, BAD 3 the bounded adder depth modification for a bound of three and MIN is the minimization problem. TO is for Time Out (3600s) and OOM stands for Out Of Memory. Experiments were made using CPLEX (top) and Gurobi (bottom).

	Names	Using indicator constraints	Using big M constraints
		ILP2 [3] ILP1 [2] Odd MAD BAD 3 MIN ILP1 [2] Odd MAD BAD 3	MIN

Table 8 :

 8 Instances used for the IIR experiments This normalization feeds through to D and D, to avoid an overlayer of complexity, just as a's and b's, D and D are rewritten:

	Constraints w.r.t.		Constraints
	specifications	w.r.t. stability
		Glue in-between
	MCM A		MCM B
	Figure 9: High-level overview of second order IIR filter ILP model
	Name	wordlength passband/π stopband/π	δ
	lowpass 8bit	8	[0.0, 0.3]	[0.7, 1.0]	0.0636
	highpass 8bit	8	[0.7, 1.0]	[0.0, 0.3]	0.0636
	lowpass 10bit	10	[0.0, 0.5]	[0.9, 1.0]	0.1

  a Pre a Total b Pre b

						NoPre NA	[b 0 , b 1 , b 2 ]	[a 1 , a 2 ]	adder values
	lowpass 8bit	15	3	19	7	23	3	[56, 88, 56] /2 8	2 × [-64, 36] /2 8	[9 a , 7 b , 11 b ]
	highpass 8bit	9	3	12	6	68	3	[56, -84, 56] /2 8	2 × [64, 40] /2 8	[7 b , 5 a , 21 b ]
	lowpass 10bit	31	31	124	124	4	1	4 × [128, 144, 32] /2 10	2 × [0, 128] /2 10	[9 b ]
	lowpass 8bit	27	9	29	11	157	3	[-64, -84, -48] /2 8	2 × [-64, 40] /2 8	[3 b , 21 b , 5 a ]
	highpass 8bit	20	7	27	13	158	3	[-64, 76, -40] /2 8	2 × [80, 40] /2 8	[5 b , 19 b , 5 a ]
	lowpass 10bit	728	726	206	204	4	1	4 × [32, 144, 128] /2 10	2 × [0, 128] /2 10	[9 b ]

  2 d a ,a i , j ≤ 2 d t a i , j ∀i ∈ [[1, 2]], i ∈ [[i, 2]], j ∈ [[0, d]] (C4.32) u a i ,a i , j ≤ a + u a i ,a i , j + 2 × 2 2d 1a u a i ,a i , j + 2 × 2 2d 1 + a

		sg i				∀i ∈ [[1, 2]]	(C4.26)
	a + i ≥ a i				∀i ∈ [[1, 2]]	(C4.27)
	a + i ≤ a i + 2 × 2 d 1 -a sg i		∀i ∈ [[1, 2]]	(C4.28)
	a + i ≥ -a i				∀i ∈ [[1, 2]]	(C4.29)
	a + i ≤ -a i + 2 × 2 d a	sg i		∀i ∈ [[1, 2]]	(C4.30)
		d
	a + i =	2 j t a i , j			∀i ∈ [[1, 2]]	(C4.31)
		j=0
	u a i d
		2 j u a i ,a i , j			∀i ∈ [[1, 2]]	(C4.35)
		j=0
		d
	a i,i ≥		2 j u a sg i + a	sg i	∀i ∈ [[1, 1]], i ∈ [[i + 1, 2]]	(C4.36)
		j=0
		d
	a i,i ≥		2 j u a sg i	∀i ∈ [[1, 1]], i ∈ [[i + 1, 2]]	(C4.37)
		j=0
		d
	a i,i ≤	j=0	2 j u a sg i + a sg i	∀i ∈ [[1, 1]], i ∈ [[i + 1, 2]]	(C4.40)
		d
	a i,i ≤ -	2 sg i -a sg i	∀i ∈ [[1, 1]], i ∈ [[i + 1, 2]]	(C4.41)
		j=0

i ∀i ∈ [[1, 2]], i ∈ [[i, 2]], j ∈ [[0, d]] (C4.33) u a i ,a i , j ≥ a + i -2 d 1t a i , j ∀i ∈ [[1, 2]], i ∈ [[i, 2]], j ∈ [[0, d]] (C4.34) a i,i = i ,a i , j -2 × 2 2d a i ,a i , j -2 × 2 2d 2a sg ia i ,a i , j ∀i ∈ [[1, 1]], i ∈ [[i + 1, 2]] (C4.38) a i,i ≥ -d j=0 2 j u a i ,a i , j ∀i ∈ [[1, 1]], i ∈ [[i + 1, 2]] (C4.39) a i,i ≤ -d j=0 2 j j

In the following the term digital will be implicit, e. g. signal and filter will stand for digital signal and digital filter.

Sometimes this is referred to as Booth's recoding or Booth's encoding, e. g. in[START_REF] Lefèvre | Multiplication by an Integer Constant[END_REF][START_REF] Hewlitt | Canonical signed digit representation for FIR digital filters[END_REF]

Appendix A. Expanding H e iω 2 H e iω = b 0 (cos(ω) + i sin(ω)) 2 + b 1 (cos(ω) + i sin(ω)) + b 2 (cos(ω) + i sin(ω)) 2 + a 1 (cos(ω) + i sin(ω)) + a 2 (A.1)

⇔ H e iω = b 0 cos 2 (ω)sin 2 (ω)

2a 1 cos(2ω) cos(ω) + 2 cos(ω) sin 2 (ω) + 2a 2 cos(2ω) + 2a 1 a 2 cos(ω)

Using the following properties:

(P1) cos 2 (ω)sin 2 (ω) = cos(2ω) is used to obtain (A.3);

(P2) 2 cos(ω) sin(ω) = sin(2ω) is used to obtain (A.6);

(P3) cos 2 (ω) + sin 2 (ω) = 1 and cos(2ω) = 1 -2 sin 2 (ω) are used to obtain (A.7).

Appendix B. Second order IIR filter ILP model