# Robustness Analysis of Floating-Point Programs by Self-Composition

Authors: Liqian Chen; Jiahong Jiang; Banghu Yin; Wei Dong; Ji Wang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Robustness is a key property for critical systems that run in uncertain environments, to ensure that small input perturbations can cause only small output changes. Current critical systems often involve lots of floating-point computations which are inexact. Robustness analysis of floating-point programs needs to consider both the uncertain inputs and the inexact computation. In this paper, we propose to leverage the idea of self-composition to transform the robustness
analysis problem into a reachability problem, which enables the use of standard reachability analysis techniques such as software model checking and symbolic
execution for robustness analysis. To handle floating-point arithmetic, we employ an abstraction that encompasses the effect of rounding and that can encompass
all rounding modes. It converts floating-point expressions into linear expressions with interval coefficients in exact real arithmetic. On this basis, we employ interval
linear programming to compute the maximum output change or maximum allowed input perturbation for the abstracted programs. Preliminary experimental
results of our prototype implementation are encouraging.


---

## Body

## 1. Introduction

Uncertainty and inexactness in computing have attracted much attention in computer science. In Cyber Physical Systems (CPS), the discrete world of computation is integrated with the continuous world of physical processes. Moreover, CPS run in the open environmental context and thus have to deal with uncertain data which may come from noisy sensor data or approximate computation. Hence, inputs for programs in CPS are of intrinsic uncertainty. On the other hand, due to finite precision on computers, physical values are truncated into digital ones. In modern computers, real numbers are approximated by a finite set of floating-point numbers. Due to the pervasive rounding errors, numerical computation using floating-point arithmetic is not exact. Since many safety-critical CPS systems (such as aircrafts, automobiles, and medical devices) often involve lots of numerical computations, there is a great need to ensure that these programs arerobust with respect to the uncertain input as well as the inexact computation.Althoughrobustness is long known as a standard correctness property for control systems [1], considering the robustness of programs is quite recent [2–5]. Intuitively, robustness of a program means that small input perturbations of the program can cause only small output changes. Much existing work on analyzing robustness of programs assumes that the analyzed program is in exact real arithmetic, although floating-point computation is pervasive in practical applications. This paper targets the analysis of robustness properties of floating-point programs.A program using floating-point arithmetic often exhibits more robustness issues than that using exact real arithmetic, due to the misunderstandings and nonintuitive behaviors of floating-point semantics. Although floating-point arithmetic is quite different from the exact real arithmetic, most developers of floating-point programs will write programs as if computations were done in exact arithmetic. For the same input, the control flow of the program using floating-point arithmetic can be different from the one that would be taken assuming exact real arithmetic. Similarly, for two inputs whose values are close to each other, the resulting two control flows of the same program can be different (even when following exact real arithmetic). Two different control flows may lead to very large difference in outputs.We illustrate the robustness problem due to floating-point computation using a motivating example shown in Figure1, which is a “classroom” example of a robustness problem frequently used in the field of geometric computations [6]. The program  Orientation implements the 2D orientation test that decides whether a point r lies to the left of, to the right of, or on the line p
q
→ defined by the 2 points p
,
q, by evaluating the sign of a determinant  det which is expressed in terms of the coordinates of the input points. Due to rounding errors, the floating-point computation of the determinant  det may lead to a wrong result when the true determinant (via exact real arithmetic) is close to zero. From the robustness point of view, for this program, even a very small input perturbation may lead to an output change of 1 or 2. If the rounding modes for the floating-point operations are not determinate in the program, the output change can be 2 even when there is no perturbation in the inputs (by running the program in different rounding modes). This misinformation may then lead to a failure of a computational geometry application (e.g., crash or not terminate) or produce wrong results [6].Figure 1
Floating-point implementation for orientation test of 2D points.Analyzing robustness of floating-point programs is more challenging than analyzing programs assuming exact real arithmetic, since besides the input perturbations, we need to consider also the inexactness of floating-point computation. The floating-point program itself acts as if inputs were perturbated due to the pervasive rounding errors or nondeterminate rounding modes. There exist a few known pitfalls of analyzing and verifying floating-point programs [7].In this paper, we present a robustness analysis method for floating-point programs. The key idea is to leverage the self-composition technique from the field of secure information flow to transform the robustness analysis problem into a reachability (safety) problem. Then we use standard rechability analysis techniques such as software model checking and symbolic execution to analyze the self-composed programs. To cope with floating-point arithmetic, we utilize a rounding mode insensitive abstraction method to abstract floating-point expressions into linear expressions with interval coefficients in the field of reals. On this basis, we use interval linear programming to compute the maximum output change (when given the input perturbation) or the maximum allowed input perturbation (when given the output change) for the abstracted programs. The preliminary experimental results are promising on benchmark programs.The rest of the paper is organized as follows. Section2 reviews the IEEE 754 floating-point arithmetic and the basic theory of interval linear systems as well as interval linear programming. Section 3 presents the robustness analysis approach via self-composition for programs (that assume exact real arithmetic). Section 4 presents the techniques to handle floating-point arithmetic. Section 5 presents our prototype implementation together with preliminary experimental results. Section 6 discusses some related work before Section 7 concludes.


## 2. Preliminaries

In this section, we briefly provide the background on the IEEE 754 floating-point arithmetic and the basic theory on interval arithmetics as well as interval linear programming.


### 2.1. The IEEE 754 Floating-Point Arithmetic

A digital computer cannot represent all possible real numbers in mathematics exactly. In computing, floating-point numbers provide an approach to represent a finite subset of the real numbers. In this paper, we focus on analyzing programs with respect to the binary formats of the IEEE 754 floating-point standard [8] which is the most commonly used floating-point representation and is followed by almost all modern computers.In the IEEE 754 standard, the binary representation of a floating-point numberx can be described as x
=
(
-
1
)
S
×
M
×
2
E, where(i)
Sis the 1-bitsign of x, which represents that x is positive (when S
=
0) or negative (when S
=
1);
(ii)
E
=
e
-
b
i
a
s is called theexponent, where e is a biased e-bit unsigned integer and b
i
a
s
=
2
e
-
1
-
1;
(iii)
M
=
m
0
.
m
1
m
2
…
m
p is called thesignificand, where f
=
.
m
1
m
2
…
m
p represents a p-bit fraction and m
0 is the hidden bit without need of storage.The values ofe
,
b
i
a
s
,
p depend on the floating-point formats. The IEEE 754 standard supports several formats, among which the basic formats include(i)
32-bit single-precision format, wheree
=
8 (and thus b
i
a
s
=
127), p
=
23;
(ii)
64-bit double-precision format, wheree
=
11 (and thus b
i
a
s
=
1023), p
=
52.According to the value ofe, the floating-point numbers can be divided into the following categories:(i)
normalized number (
-
1
)
S
×
1
.
f
×
2
e
-
b
i
a
s, when 1
≤
e
≤
2
e
-
2;
(ii)
denormalized number (
-
1
)
S
×
0
.
f
×
2
1
-
b
i
a
s, when e
=
0 and f
≠
0;
(iii)
+
0 or -
0, when e
=
0 and f
=
0;
(iv)
+
∞ or -
∞, when e
=
2
e
-
1 and f
=
0;
(v)
NaN (Not a Number), when e
=
2
e
-
1 and f
≠
0.Let F be the set of floating-point formats. For each f
∈
F, we define(i)
m
f
f
=
def
2
1
-
b
i
a
s
-
p, the smallest nonzero positive floating-point number;
(ii)
M
f
f
=
def
(
2
-
2
-
p
)
2
2
e
-
b
i
a
s
-
2, the largest noninfinity floating-point number.In general, the result of a floating-point operation may not be exactly representable in the floating-point representation, and thus the result needs to be rounded into a floating-point number. The IEEE 754 standard supports four rounding modes: toward nearest, toward+
∞, toward -
∞, and toward zero. In this paper, in order to distinguish floating-point arithmetic operations from exact real arithmetic ones, we introduce additional notations. As usual, {
+
,
-
,
×
,
/
} are used as exact rational arithmetic operations. The corresponding floating-point operations are denoted by {
⊕
f
,
r
,
⊖
f
,
r
,
⊗
f
,
r
,
⊘
f
,
r
}, tagged with a floating-point format f
∈
F and a rounding mode r
∈
{
+
∞
,
-
∞
,
0
,
n
} (n representing rounding to nearest). We also use ? to denote arbitrary rounding mode.Due to rounding errors, many well-known algebraic properties (such as associativity and distributivity) over the reals do not hold for floating-point arithmetic.Example 1.
Consider the following expressions in the 32-bit single-precision floating-point arithmetic:(1)
(
2
24
⊕
32
,
?
-
2
24
)
⊕
32
,
?
1
=
1
(
2
24
⊕
32
,
-
∞
1
)
⊕
32
,
-
∞
-
2
24
=
0
(
2
24
⊕
32
,
+
∞
1
)
⊕
32
,
+
∞
-
2
24
=
2
.
Note that in the 32-bit single-precision format, the significand is M
=
m
0
.
m
1
m
2
…
m
23. However, to represent the exact result of 2
24
+
1 over the reals, we need one more bit for the significand M (say m
24). Hence, rounding happens. 2
24
⊕
32
,
-
∞
1 will result in 2
24, while 2
24
⊕
32
,
+
∞
1 will result in (
1
+
2
-
23
)
×
2
24.


### 2.2. Interval Linear Systems and Interval Linear Programming

LetA
_
,
A
¯
∈
R
m
×
n be two matrices with A
_
≤
A
¯, where comparison operators are defined element-wise; then the set of matrices A
∈
I
R
m
×
n defined by
(2)
A
=
[
A
_
,
A
¯
]
=
{
A
∈
R
m
×
n
:
A
_
≤
A
≤
A
¯
}
is called aninterval matrix, and the matrices A
_
,
A
¯ are called its bounds. Let us define thecenter matrix of A as A
c
=
(
1
/
2
)
(
A
_
+
A
¯
) and theradius matrix as Δ
A
=
(
1
/
2
)
(
A
¯
-
A
_
). Then, A
=
[
A
_
,
A
¯
]
=
[
A
c
-
Δ
A
,
A
c
+
Δ
A
]. Aninterval vector is a one-column interval matrix d
=
[
d
_
,
d
¯
]
=
{
d
∈
R
m
:
d
_
≤
d
≤
d
¯
}, where d
_
,
d
¯
∈
R
m and d
_
≤
d
¯.LetA be an m
×
n interval matrix and b be a vector of size m. The following system of interval linear inequalities
(3)
A
x
≤
b
denotes aninterval linear system, that is, thefamily of all systems of linear inequalities A
x
≤
b such that A
∈
A.Definition 2 (weak solution).
A vectorx
∈
R
n is called aweak solution of the interval linear system A
x
≤
b, if it satisfies A
x
≤
b for some A
∈
A. Furthermore, the set
(4)
Σ
∃
(
A
,
b
)
=
{
x
∈
R
n
:
∃
A
∈
A
,
A
x
≤
b
}
is said to be theweak solution set of the system A
x
≤
b.The weak solution set of an interval linear system is characterized by the following theorem [9].Theorem 3.
A vectorx
∈
R
n is a weak solution of A
x
≤
b if and only if it satisfies A
c
x
-
Δ
A
|
x
|
≤
b.LetA
∈
I
R
m
×
n be an m
×
n interval matrix, b
∈
R
m be an m-dimensional vector, and c
∈
I
R
n be an n-dimensional interval vector. Thefamily of linear programming (LP) problems
(5)
f
(
A
,
b
,
c
)
=
max
⁡
{
c
T
x
:
A
x
≤
b
}
with data satisfying
(6)
A
∈
A
,
c
∈
c
is called an interval linear programming (ILP) problem.In this paper, we are only interested in computing the upper boundf
¯
(
A
,
b
,
c
)
=
sup
⁡
{
f
(
A
,
b
,
c
)
:
A
∈
A
,
c
∈
c
}. In general, according to Theorem 3, to compute the exact f
¯
(
A
,
b
,
c
), in the worst case up to 2
n LP problems have to be solved, one for each orthant. Recall that a (closed) orthant is one of the 2
n subsets of an n-dimensional Euclidean space defined by constraining each Cartesian coordinate axis to be either nonnegative or nonpositive. In each orthant, we consider the following LP problem:
(7)
max
⁡
∑
j
=
1
n
‍
c
j
′
x
j
s
.
t
.
⋀
0
≤
i
≤
m
∑
j
=
1
n
‍
A
i
j
′
x
j
≤
b
i
,
where
(8)
c
j
′
=
{
c
¯
j
if
x
j
≥
0
c
_
j
if
x
j
<
0
A
i
j
′
=
{
A
_
i
j
if
x
j
≥
0
A
¯
i
j
if
x
j
<
0
.
And f
¯
(
A
,
b
,
c
) will be the the maximum over all the optimal values of the 2
n LP problems with one per each orthant.


## 2.1. The IEEE 754 Floating-Point Arithmetic

A digital computer cannot represent all possible real numbers in mathematics exactly. In computing, floating-point numbers provide an approach to represent a finite subset of the real numbers. In this paper, we focus on analyzing programs with respect to the binary formats of the IEEE 754 floating-point standard [8] which is the most commonly used floating-point representation and is followed by almost all modern computers.In the IEEE 754 standard, the binary representation of a floating-point numberx can be described as x
=
(
-
1
)
S
×
M
×
2
E, where(i)
Sis the 1-bitsign of x, which represents that x is positive (when S
=
0) or negative (when S
=
1);
(ii)
E
=
e
-
b
i
a
s is called theexponent, where e is a biased e-bit unsigned integer and b
i
a
s
=
2
e
-
1
-
1;
(iii)
M
=
m
0
.
m
1
m
2
…
m
p is called thesignificand, where f
=
.
m
1
m
2
…
m
p represents a p-bit fraction and m
0 is the hidden bit without need of storage.The values ofe
,
b
i
a
s
,
p depend on the floating-point formats. The IEEE 754 standard supports several formats, among which the basic formats include(i)
32-bit single-precision format, wheree
=
8 (and thus b
i
a
s
=
127), p
=
23;
(ii)
64-bit double-precision format, wheree
=
11 (and thus b
i
a
s
=
1023), p
=
52.According to the value ofe, the floating-point numbers can be divided into the following categories:(i)
normalized number (
-
1
)
S
×
1
.
f
×
2
e
-
b
i
a
s, when 1
≤
e
≤
2
e
-
2;
(ii)
denormalized number (
-
1
)
S
×
0
.
f
×
2
1
-
b
i
a
s, when e
=
0 and f
≠
0;
(iii)
+
0 or -
0, when e
=
0 and f
=
0;
(iv)
+
∞ or -
∞, when e
=
2
e
-
1 and f
=
0;
(v)
NaN (Not a Number), when e
=
2
e
-
1 and f
≠
0.Let F be the set of floating-point formats. For each f
∈
F, we define(i)
m
f
f
=
def
2
1
-
b
i
a
s
-
p, the smallest nonzero positive floating-point number;
(ii)
M
f
f
=
def
(
2
-
2
-
p
)
2
2
e
-
b
i
a
s
-
2, the largest noninfinity floating-point number.In general, the result of a floating-point operation may not be exactly representable in the floating-point representation, and thus the result needs to be rounded into a floating-point number. The IEEE 754 standard supports four rounding modes: toward nearest, toward+
∞, toward -
∞, and toward zero. In this paper, in order to distinguish floating-point arithmetic operations from exact real arithmetic ones, we introduce additional notations. As usual, {
+
,
-
,
×
,
/
} are used as exact rational arithmetic operations. The corresponding floating-point operations are denoted by {
⊕
f
,
r
,
⊖
f
,
r
,
⊗
f
,
r
,
⊘
f
,
r
}, tagged with a floating-point format f
∈
F and a rounding mode r
∈
{
+
∞
,
-
∞
,
0
,
n
} (n representing rounding to nearest). We also use ? to denote arbitrary rounding mode.Due to rounding errors, many well-known algebraic properties (such as associativity and distributivity) over the reals do not hold for floating-point arithmetic.Example 1.
Consider the following expressions in the 32-bit single-precision floating-point arithmetic:(1)
(
2
24
⊕
32
,
?
-
2
24
)
⊕
32
,
?
1
=
1
(
2
24
⊕
32
,
-
∞
1
)
⊕
32
,
-
∞
-
2
24
=
0
(
2
24
⊕
32
,
+
∞
1
)
⊕
32
,
+
∞
-
2
24
=
2
.
Note that in the 32-bit single-precision format, the significand is M
=
m
0
.
m
1
m
2
…
m
23. However, to represent the exact result of 2
24
+
1 over the reals, we need one more bit for the significand M (say m
24). Hence, rounding happens. 2
24
⊕
32
,
-
∞
1 will result in 2
24, while 2
24
⊕
32
,
+
∞
1 will result in (
1
+
2
-
23
)
×
2
24.


## 2.2. Interval Linear Systems and Interval Linear Programming

LetA
_
,
A
¯
∈
R
m
×
n be two matrices with A
_
≤
A
¯, where comparison operators are defined element-wise; then the set of matrices A
∈
I
R
m
×
n defined by
(2)
A
=
[
A
_
,
A
¯
]
=
{
A
∈
R
m
×
n
:
A
_
≤
A
≤
A
¯
}
is called aninterval matrix, and the matrices A
_
,
A
¯ are called its bounds. Let us define thecenter matrix of A as A
c
=
(
1
/
2
)
(
A
_
+
A
¯
) and theradius matrix as Δ
A
=
(
1
/
2
)
(
A
¯
-
A
_
). Then, A
=
[
A
_
,
A
¯
]
=
[
A
c
-
Δ
A
,
A
c
+
Δ
A
]. Aninterval vector is a one-column interval matrix d
=
[
d
_
,
d
¯
]
=
{
d
∈
R
m
:
d
_
≤
d
≤
d
¯
}, where d
_
,
d
¯
∈
R
m and d
_
≤
d
¯.LetA be an m
×
n interval matrix and b be a vector of size m. The following system of interval linear inequalities
(3)
A
x
≤
b
denotes aninterval linear system, that is, thefamily of all systems of linear inequalities A
x
≤
b such that A
∈
A.Definition 2 (weak solution).
A vectorx
∈
R
n is called aweak solution of the interval linear system A
x
≤
b, if it satisfies A
x
≤
b for some A
∈
A. Furthermore, the set
(4)
Σ
∃
(
A
,
b
)
=
{
x
∈
R
n
:
∃
A
∈
A
,
A
x
≤
b
}
is said to be theweak solution set of the system A
x
≤
b.The weak solution set of an interval linear system is characterized by the following theorem [9].Theorem 3.
A vectorx
∈
R
n is a weak solution of A
x
≤
b if and only if it satisfies A
c
x
-
Δ
A
|
x
|
≤
b.LetA
∈
I
R
m
×
n be an m
×
n interval matrix, b
∈
R
m be an m-dimensional vector, and c
∈
I
R
n be an n-dimensional interval vector. Thefamily of linear programming (LP) problems
(5)
f
(
A
,
b
,
c
)
=
max
⁡
{
c
T
x
:
A
x
≤
b
}
with data satisfying
(6)
A
∈
A
,
c
∈
c
is called an interval linear programming (ILP) problem.In this paper, we are only interested in computing the upper boundf
¯
(
A
,
b
,
c
)
=
sup
⁡
{
f
(
A
,
b
,
c
)
:
A
∈
A
,
c
∈
c
}. In general, according to Theorem 3, to compute the exact f
¯
(
A
,
b
,
c
), in the worst case up to 2
n LP problems have to be solved, one for each orthant. Recall that a (closed) orthant is one of the 2
n subsets of an n-dimensional Euclidean space defined by constraining each Cartesian coordinate axis to be either nonnegative or nonpositive. In each orthant, we consider the following LP problem:
(7)
max
⁡
∑
j
=
1
n
‍
c
j
′
x
j
s
.
t
.
⋀
0
≤
i
≤
m
∑
j
=
1
n
‍
A
i
j
′
x
j
≤
b
i
,
where
(8)
c
j
′
=
{
c
¯
j
if
x
j
≥
0
c
_
j
if
x
j
<
0
A
i
j
′
=
{
A
_
i
j
if
x
j
≥
0
A
¯
i
j
if
x
j
<
0
.
And f
¯
(
A
,
b
,
c
) will be the the maximum over all the optimal values of the 2
n LP problems with one per each orthant.


## 3. Robustness Analysis via Self-Composition

### 3.1. Robustness of Programs

In this paper, we follow the definition for robustness of programs used by Majumdar and Saha [5]. Let f be a function with inputs x
1
,
…
,
x
n and output y; that is, y
=
f
(
x
1
,
…
,
x
n
). The function f is said to be (
δ
,
ϵ
)-r
o
b
u
s
t in the ith input x
i if a perturbation of at most δ in the input x
i can only cause a change of at most ϵ in the output; that is,
(9)
∀
x
i
,
x
i
′
·
|
x
i
-
x
i
′
|
≤
δ
⟹
|
f
(
x
1
,
…
,
x
i
,
…
,
x
n
)
-
f
(
x
1
,
…
,
x
i
′
,
…
,
x
n
)
|
≤
ϵ
,
where δ
,
ϵ
∈
R are nonnegative constant parameters specified by users. Recall that we consider the perturbation over only one input at a time while assume that there is no perturbation over all other inputs at the same time.Moreover, in practice, users may be interested in the maximum output change ofy with respect to x
i and δ; that is,
(10)
ϵ
¯
δ
=
def
max
⁡
x
,
x
′
{
|
y
-
y
′
|
∣
y
=
f
(
x
1
,
…
,
x
i
,
…
,
x
n
)
y
′
=
f
(
x
1
,
…
,
x
i
′
,
…
,
x
n
)
|
x
i
-
x
i
′
|
≤
δ
}
.
Similarly, users may be interested in the maximum input perturbation allowed over x
i with respect to y and ϵ; that is,
(11)
δ
¯
ϵ
=
def
max
⁡
y
,
y
′
{
|
x
i
-
x
i
′
|
∣
y
=
f
(
x
1
,
…
,
x
i
,
…
,
x
n
)
y
′
=
f
(
x
1
,
…
,
x
i
′
,
…
,
x
n
)
|
y
-
y
′
|
≤
ϵ
}
.Example 4.
Consider the program shown in Figure2, which implements a piece-wise linear function. When x
=
1.001, the two branches give the same result y
=
1002.001 in exact real arithmetic (assuming floats are reals). It is easy to see that in exact real arithmetic, this program is (
0.1
,
ϵ
0
)-r
o
b
u
s
t for all ϵ
0
≥
100.1 but is not (
0.1
,
ϵ
1
)-r
o
b
u
s
t for all ϵ
1
<
100.1. This can be deduced by observing that in exact real arithmetic, given the input perturbation δ
=
0.1, the maximum output change of y is ϵ
¯
δ
=
100.1; given the output change ϵ
=
100.1, the maximum input perturbation allowed over x is δ
¯
ϵ
=
0.1.Figure 2
A floating-point programpiecewise_linear.


### 3.2. Self-Composition

The idea of self-composition is firstly used in the field of secure information flow [10, 11] to characterize noninterference. Let P be a program and P
′ be a copy of P with each variable x in P replaced by a fresh variable x
′. Using Hoare triples, noninterference can be characterized as
(12)
{
L
=
L
′
}
P
;
P
′
{
L
=
L
′
}
,
where L denotes low-security variables. In other words, it requires that running two instances of the same program with equal low-security values and arbitrary high-security values results in equal low-security values. Hence, via self-composition, a secure information flow property of P reduces to a reachability property over single program executions of the program P
;
P
′.In this paper, we would like to leverage the idea of self-composition to reduce the robustness problem of a programP into an equivalent reachability problem over P
;
P
′. Assume that program P has n input variables x
1
,
…
,
x
n and an output variable y. Similarly, using Hoare triples, the (
δ
,
ϵ
)-r
o
b
u
s
t
n
e
s
s of program P over the ith input x
i can be characterized as
(13)
{
|
x
i
-
x
i
′
|
≤
δ
∧
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
}
P
;
P
′
{
|
y
-
y
′
|
≤
ϵ
}
.Example 5.
Consider again the programpiecewise_linear in Figure 2. The self-composition of the function body is shown in Figure 3. To express the robustness property, we add the assumption |
x
-
x
′
|
≤
δ as a precondition at the beginning of the self-composed program and add an assertion |
y
-
y
′
|
≤
ϵ as a postcondition at the end.Figure 3
Self-composition of thepiecewise_linear program for robustness analysis.Essentially, the copied programP
′ has the same program code as P but uses variables with different initial values. Hence, there exists inherent symmetry and redundancy in the self-composed programs. In order to make the following analysis and verification process for self-composed programs easier, program transformations can be used to optimize the self-composed programs. In the field of secure information flow analysis, Terauchi and Aiken [12] proposed type-directed transformation to improve self-composition. The main idea of type-directed transformation is not to self-compose branch (or loop) statements when the branch (or loop) condition is only dependent on the values of low-security variables. In addition, for an assignment statement {
x
:
=
e
;
}, when the right-hand expression e is only dependent on the values of low-security variables, its self-composition is simplified as {
x
:
=
e
;
x
′
:
=
x
;
}.With respect to robustness, a similar transformation can be applied. Intuitively, we could consider the perturbed input variablex
i as a high-security variable and all other input variables x
j's as low-security variables where j
≠
i. Hence, similarly to type directed transformation, we do not self-compose branch (or loop) statements when the branch (or loop) condition is not dependent on the values of perturbed input variables. For an assignment statement {
x
:
=
e
;
}, when the right-hand expression e is not dependent on the values of perturbed input variables, its self-composition is simplified as {
x
:
=
e
;
x
′
:
=
x
;
}.Example 6.
Consider the functionm
i
n
_
p
l
u
s
1 shown in Figure 4, which implements min
⁡
(
x
+
1
,
y
) by adding 0.1 to x ten times. The optimized self-composition result of the function body after applying transformation is given in Figure 5, when we consider the perturbation over the input variable x (while assuming no perturbation over y). More specifically, since the loop condition i
<
10 in the original program is not dependent on the value of the perturbed input variable x, we do not self-compose the loop statement and thus there is only one loop in the transformed resulting self-composed program.Figure 4
A floating-point programmin_plus1.Figure 5
Transformed self-composition of themin_plus1 program for robustness analysis.


### 3.3. Robustness Analysis of Self-Composed Programs

Via self-composition, the robustness analysis problem can be reduced to solving a standard reachability (safety) problem. The recent success of automatic analysis and verification tools (such as SLAM [13], CBMC [14], and ASTRÉE [15]) aiming at checking reachability properties in programs makes this approach promising. In the following, we will present two popular reachability analysis approaches that fit for analyzing robustness, that is, software model checking and symbolic execution.


#### 3.3.1. Checking Robustness by Software Model Checking

Software model checking [16] provides an automatic approach to check whether a program satisfies a property by exploring the state space of the program. For the robustness analysis problem, the property to be checked is an assertion at the end of the self-composed programs stating that the output change is bounded by ϵ, that is,  assert
(
-
ϵ
≤
y
-
y
′
≤
ϵ
). A main advantage of using software model checking is that it will generate a counterexample when the robustness property does not hold. The counterexample shows an execution trace which violates the robustness property. A counterexample is very helpful for the users to identify the source of nonrobustness.


#### 3.3.2. Finding Maximum Output Change (or Input Perturbation) by Symbolic Execution

Symbolic execution [17, 18] is a technique to analyze a program by executing the program with symbolic rather than concrete values as program inputs. The process of symbolic execution essentially generates and explores a symbolic execution tree which represents all execution paths followed during the process. Each tree node represents a symbolic execution state, while each edge represents a program transition between the states. At any tree node, the symbolic execution state includes a program counter, apath condition (PC) that encodes the constraints on the symbolic inputs to reach that node, apath function (PF) that represents the current values of the program variables as function of symbolic inputs when the path condition holds true. The path condition is a boolean expression over the symbolic inputs. The path function describes the expected result of the program, under the given path condition. Due to conditional branches and loops in a program, the symbolic execution of a program will result in a set of paths, each of which is described by a pair 〈
PC
,
PF
〉 of the path condition PC and the associated path function PF.We now show how to use symbolic execution to conduct a robustness analysis of programP with inputs x
1
,
…
,
x
n and output y. First, the analysis algorithm performs symbolic execution on the self-composed program P
;
P
′. Assume that the algorithm collects, at the end of the self-composed program, a set S of pairs 〈
PC
,
PF
〉 of the path condition PC and the associated path function PF. Then for each s
≜
〈
PC
,
PF
〉
∈
S, we compute the maximum output change ϵ
¯
δ
s:
(14)
max
⁡
|
P
F
y
-
P
F
y
′
|
h
s
.
t
.
(
|
x
i
-
x
i
′
|
≤
δ
∧
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
)
∧
PC
.
Here, P
F
y and P
F
y
′ denote the symbolic expressions that the path function PF maps the variables y and y
′ to, respectively. Let ϵ
¯
δ be the maximum element of {
ϵ
¯
δ
s
∣
s
∈
S
}; that is, ϵ
¯
δ
=
max
⁡
(
{
ϵ
¯
δ
s
∣
s
∈
S
}
). If ϵ
¯
δ
≤
ϵ, then the original program P is (
δ
,
ϵ
)-r
o
b
u
s
t.Similarly, given the bound of output changeϵ, computing the maximum allowed input perturbation is reduced to solving a series of the following optimization problems for each s
≜
〈
PC
,
PF
〉
∈
S:
(15)
max
⁡
|
x
i
-
x
i
′
|
h
s
.
t
.
(
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
)
∧
PC
∧
|
P
F
y
-
P
F
y
′
|
≤
ϵ
.
And δ
¯
ϵ will be the maximum element of {
δ
¯
ϵ
s
∣
s
∈
S
}; that is, δ
¯
ϵ
=
max
⁡
(
{
δ
¯
ϵ
s
∣
s
∈
S
}
).


## 3.1. Robustness of Programs

In this paper, we follow the definition for robustness of programs used by Majumdar and Saha [5]. Let f be a function with inputs x
1
,
…
,
x
n and output y; that is, y
=
f
(
x
1
,
…
,
x
n
). The function f is said to be (
δ
,
ϵ
)-r
o
b
u
s
t in the ith input x
i if a perturbation of at most δ in the input x
i can only cause a change of at most ϵ in the output; that is,
(9)
∀
x
i
,
x
i
′
·
|
x
i
-
x
i
′
|
≤
δ
⟹
|
f
(
x
1
,
…
,
x
i
,
…
,
x
n
)
-
f
(
x
1
,
…
,
x
i
′
,
…
,
x
n
)
|
≤
ϵ
,
where δ
,
ϵ
∈
R are nonnegative constant parameters specified by users. Recall that we consider the perturbation over only one input at a time while assume that there is no perturbation over all other inputs at the same time.Moreover, in practice, users may be interested in the maximum output change ofy with respect to x
i and δ; that is,
(10)
ϵ
¯
δ
=
def
max
⁡
x
,
x
′
{
|
y
-
y
′
|
∣
y
=
f
(
x
1
,
…
,
x
i
,
…
,
x
n
)
y
′
=
f
(
x
1
,
…
,
x
i
′
,
…
,
x
n
)
|
x
i
-
x
i
′
|
≤
δ
}
.
Similarly, users may be interested in the maximum input perturbation allowed over x
i with respect to y and ϵ; that is,
(11)
δ
¯
ϵ
=
def
max
⁡
y
,
y
′
{
|
x
i
-
x
i
′
|
∣
y
=
f
(
x
1
,
…
,
x
i
,
…
,
x
n
)
y
′
=
f
(
x
1
,
…
,
x
i
′
,
…
,
x
n
)
|
y
-
y
′
|
≤
ϵ
}
.Example 4.
Consider the program shown in Figure2, which implements a piece-wise linear function. When x
=
1.001, the two branches give the same result y
=
1002.001 in exact real arithmetic (assuming floats are reals). It is easy to see that in exact real arithmetic, this program is (
0.1
,
ϵ
0
)-r
o
b
u
s
t for all ϵ
0
≥
100.1 but is not (
0.1
,
ϵ
1
)-r
o
b
u
s
t for all ϵ
1
<
100.1. This can be deduced by observing that in exact real arithmetic, given the input perturbation δ
=
0.1, the maximum output change of y is ϵ
¯
δ
=
100.1; given the output change ϵ
=
100.1, the maximum input perturbation allowed over x is δ
¯
ϵ
=
0.1.Figure 2
A floating-point programpiecewise_linear.


## 3.2. Self-Composition

The idea of self-composition is firstly used in the field of secure information flow [10, 11] to characterize noninterference. Let P be a program and P
′ be a copy of P with each variable x in P replaced by a fresh variable x
′. Using Hoare triples, noninterference can be characterized as
(12)
{
L
=
L
′
}
P
;
P
′
{
L
=
L
′
}
,
where L denotes low-security variables. In other words, it requires that running two instances of the same program with equal low-security values and arbitrary high-security values results in equal low-security values. Hence, via self-composition, a secure information flow property of P reduces to a reachability property over single program executions of the program P
;
P
′.In this paper, we would like to leverage the idea of self-composition to reduce the robustness problem of a programP into an equivalent reachability problem over P
;
P
′. Assume that program P has n input variables x
1
,
…
,
x
n and an output variable y. Similarly, using Hoare triples, the (
δ
,
ϵ
)-r
o
b
u
s
t
n
e
s
s of program P over the ith input x
i can be characterized as
(13)
{
|
x
i
-
x
i
′
|
≤
δ
∧
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
}
P
;
P
′
{
|
y
-
y
′
|
≤
ϵ
}
.Example 5.
Consider again the programpiecewise_linear in Figure 2. The self-composition of the function body is shown in Figure 3. To express the robustness property, we add the assumption |
x
-
x
′
|
≤
δ as a precondition at the beginning of the self-composed program and add an assertion |
y
-
y
′
|
≤
ϵ as a postcondition at the end.Figure 3
Self-composition of thepiecewise_linear program for robustness analysis.Essentially, the copied programP
′ has the same program code as P but uses variables with different initial values. Hence, there exists inherent symmetry and redundancy in the self-composed programs. In order to make the following analysis and verification process for self-composed programs easier, program transformations can be used to optimize the self-composed programs. In the field of secure information flow analysis, Terauchi and Aiken [12] proposed type-directed transformation to improve self-composition. The main idea of type-directed transformation is not to self-compose branch (or loop) statements when the branch (or loop) condition is only dependent on the values of low-security variables. In addition, for an assignment statement {
x
:
=
e
;
}, when the right-hand expression e is only dependent on the values of low-security variables, its self-composition is simplified as {
x
:
=
e
;
x
′
:
=
x
;
}.With respect to robustness, a similar transformation can be applied. Intuitively, we could consider the perturbed input variablex
i as a high-security variable and all other input variables x
j's as low-security variables where j
≠
i. Hence, similarly to type directed transformation, we do not self-compose branch (or loop) statements when the branch (or loop) condition is not dependent on the values of perturbed input variables. For an assignment statement {
x
:
=
e
;
}, when the right-hand expression e is not dependent on the values of perturbed input variables, its self-composition is simplified as {
x
:
=
e
;
x
′
:
=
x
;
}.Example 6.
Consider the functionm
i
n
_
p
l
u
s
1 shown in Figure 4, which implements min
⁡
(
x
+
1
,
y
) by adding 0.1 to x ten times. The optimized self-composition result of the function body after applying transformation is given in Figure 5, when we consider the perturbation over the input variable x (while assuming no perturbation over y). More specifically, since the loop condition i
<
10 in the original program is not dependent on the value of the perturbed input variable x, we do not self-compose the loop statement and thus there is only one loop in the transformed resulting self-composed program.Figure 4
A floating-point programmin_plus1.Figure 5
Transformed self-composition of themin_plus1 program for robustness analysis.


## 3.3. Robustness Analysis of Self-Composed Programs

Via self-composition, the robustness analysis problem can be reduced to solving a standard reachability (safety) problem. The recent success of automatic analysis and verification tools (such as SLAM [13], CBMC [14], and ASTRÉE [15]) aiming at checking reachability properties in programs makes this approach promising. In the following, we will present two popular reachability analysis approaches that fit for analyzing robustness, that is, software model checking and symbolic execution.


### 3.3.1. Checking Robustness by Software Model Checking

Software model checking [16] provides an automatic approach to check whether a program satisfies a property by exploring the state space of the program. For the robustness analysis problem, the property to be checked is an assertion at the end of the self-composed programs stating that the output change is bounded by ϵ, that is,  assert
(
-
ϵ
≤
y
-
y
′
≤
ϵ
). A main advantage of using software model checking is that it will generate a counterexample when the robustness property does not hold. The counterexample shows an execution trace which violates the robustness property. A counterexample is very helpful for the users to identify the source of nonrobustness.


### 3.3.2. Finding Maximum Output Change (or Input Perturbation) by Symbolic Execution

Symbolic execution [17, 18] is a technique to analyze a program by executing the program with symbolic rather than concrete values as program inputs. The process of symbolic execution essentially generates and explores a symbolic execution tree which represents all execution paths followed during the process. Each tree node represents a symbolic execution state, while each edge represents a program transition between the states. At any tree node, the symbolic execution state includes a program counter, apath condition (PC) that encodes the constraints on the symbolic inputs to reach that node, apath function (PF) that represents the current values of the program variables as function of symbolic inputs when the path condition holds true. The path condition is a boolean expression over the symbolic inputs. The path function describes the expected result of the program, under the given path condition. Due to conditional branches and loops in a program, the symbolic execution of a program will result in a set of paths, each of which is described by a pair 〈
PC
,
PF
〉 of the path condition PC and the associated path function PF.We now show how to use symbolic execution to conduct a robustness analysis of programP with inputs x
1
,
…
,
x
n and output y. First, the analysis algorithm performs symbolic execution on the self-composed program P
;
P
′. Assume that the algorithm collects, at the end of the self-composed program, a set S of pairs 〈
PC
,
PF
〉 of the path condition PC and the associated path function PF. Then for each s
≜
〈
PC
,
PF
〉
∈
S, we compute the maximum output change ϵ
¯
δ
s:
(14)
max
⁡
|
P
F
y
-
P
F
y
′
|
h
s
.
t
.
(
|
x
i
-
x
i
′
|
≤
δ
∧
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
)
∧
PC
.
Here, P
F
y and P
F
y
′ denote the symbolic expressions that the path function PF maps the variables y and y
′ to, respectively. Let ϵ
¯
δ be the maximum element of {
ϵ
¯
δ
s
∣
s
∈
S
}; that is, ϵ
¯
δ
=
max
⁡
(
{
ϵ
¯
δ
s
∣
s
∈
S
}
). If ϵ
¯
δ
≤
ϵ, then the original program P is (
δ
,
ϵ
)-r
o
b
u
s
t.Similarly, given the bound of output changeϵ, computing the maximum allowed input perturbation is reduced to solving a series of the following optimization problems for each s
≜
〈
PC
,
PF
〉
∈
S:
(15)
max
⁡
|
x
i
-
x
i
′
|
h
s
.
t
.
(
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
)
∧
PC
∧
|
P
F
y
-
P
F
y
′
|
≤
ϵ
.
And δ
¯
ϵ will be the maximum element of {
δ
¯
ϵ
s
∣
s
∈
S
}; that is, δ
¯
ϵ
=
max
⁡
(
{
δ
¯
ϵ
s
∣
s
∈
S
}
).


## 3.3.1. Checking Robustness by Software Model Checking

Software model checking [16] provides an automatic approach to check whether a program satisfies a property by exploring the state space of the program. For the robustness analysis problem, the property to be checked is an assertion at the end of the self-composed programs stating that the output change is bounded by ϵ, that is,  assert
(
-
ϵ
≤
y
-
y
′
≤
ϵ
). A main advantage of using software model checking is that it will generate a counterexample when the robustness property does not hold. The counterexample shows an execution trace which violates the robustness property. A counterexample is very helpful for the users to identify the source of nonrobustness.


## 3.3.2. Finding Maximum Output Change (or Input Perturbation) by Symbolic Execution

Symbolic execution [17, 18] is a technique to analyze a program by executing the program with symbolic rather than concrete values as program inputs. The process of symbolic execution essentially generates and explores a symbolic execution tree which represents all execution paths followed during the process. Each tree node represents a symbolic execution state, while each edge represents a program transition between the states. At any tree node, the symbolic execution state includes a program counter, apath condition (PC) that encodes the constraints on the symbolic inputs to reach that node, apath function (PF) that represents the current values of the program variables as function of symbolic inputs when the path condition holds true. The path condition is a boolean expression over the symbolic inputs. The path function describes the expected result of the program, under the given path condition. Due to conditional branches and loops in a program, the symbolic execution of a program will result in a set of paths, each of which is described by a pair 〈
PC
,
PF
〉 of the path condition PC and the associated path function PF.We now show how to use symbolic execution to conduct a robustness analysis of programP with inputs x
1
,
…
,
x
n and output y. First, the analysis algorithm performs symbolic execution on the self-composed program P
;
P
′. Assume that the algorithm collects, at the end of the self-composed program, a set S of pairs 〈
PC
,
PF
〉 of the path condition PC and the associated path function PF. Then for each s
≜
〈
PC
,
PF
〉
∈
S, we compute the maximum output change ϵ
¯
δ
s:
(14)
max
⁡
|
P
F
y
-
P
F
y
′
|
h
s
.
t
.
(
|
x
i
-
x
i
′
|
≤
δ
∧
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
)
∧
PC
.
Here, P
F
y and P
F
y
′ denote the symbolic expressions that the path function PF maps the variables y and y
′ to, respectively. Let ϵ
¯
δ be the maximum element of {
ϵ
¯
δ
s
∣
s
∈
S
}; that is, ϵ
¯
δ
=
max
⁡
(
{
ϵ
¯
δ
s
∣
s
∈
S
}
). If ϵ
¯
δ
≤
ϵ, then the original program P is (
δ
,
ϵ
)-r
o
b
u
s
t.Similarly, given the bound of output changeϵ, computing the maximum allowed input perturbation is reduced to solving a series of the following optimization problems for each s
≜
〈
PC
,
PF
〉
∈
S:
(15)
max
⁡
|
x
i
-
x
i
′
|
h
s
.
t
.
(
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
)
∧
PC
∧
|
P
F
y
-
P
F
y
′
|
≤
ϵ
.
And δ
¯
ϵ will be the maximum element of {
δ
¯
ϵ
s
∣
s
∈
S
}; that is, δ
¯
ϵ
=
max
⁡
(
{
δ
¯
ϵ
s
∣
s
∈
S
}
).


## 4. Robustness Analysis of Floating-Point Programs

In this section, we consider the robustness analysis problem of floating-point programs. In Section3.3, we propose to utilize software model checking and symbolic execution to perform robustness analysis of self-composed programs (in exact real arithmetic). However, most existing software model checkers and symbolic execution tools can not be directly applied to floating-point programs, since they rely on constraint solvers that often assume good algebraic properties such as associativity and distributivity over the reals which do not hold for floating-point arithmetic. To handle floating-point arithmetic, we have to resort to bit-precise modeling of floating-point arithmetic or abstracting floating-point arithmetic to real number arithmetic.CBMC (C Bounded Model Checker) [14] is one of the few software model checkers that have considered floating-point arithmetic. CBMC employs a sound and complete decision procedure for floating-point arithmetic [19, 20]. It precisely encodes floating-point operations as functions on bit-vectors. Each floating-point operation is further modeled as a formula in propositional logic. The formula is then handled by a SAT-solver in the backend to check for satisfiability.When we consider symbolic execution of floating-point programs, both the path condition and the path function will involve floating-point expressions. Hence, to compute the maximum output change (or maximum allowed input perturbation), we need optimization methods supporting floating-point constraints. However, as far as we know, even for linear programming, there is no available sound solver supporting floating-point constraints. To this end, in this paper, we abstract the optimization problem with floating-point constraints into an interval linear programming problem (i.e., linear programming problem with interval coefficients) over the reals. The main idea is to use the so-calledfloating-point linearization technique [21, 22] to abstract floating-point expressions into linear real number expressions with interval coefficients (in the form of Σ
i
[
a
i
,
b
i
]
x
i).


### 4.1. Floating-Point Abstraction

In this subsection, we will explain how to abstract floating-point expressions into interval linear expressions over the reals.First, let us consider the upper bound on rounding errors due to one floating-point operation. LetR
f
,
r
(
x
) denote the floating-point rounding function that maps a real number x to a floating-point number (or a runtime error due to, for example, overflows) with respect to the floating-point format f and the rounding mode r. The amount of the rounding error due to R
f
,
r
(
x
) depends on the category of x. (i)
Ifx is in the range of normalized numbers, then |
R
f
,
r
(
x
)
-
x
|
≤
ε
rel
⁡
·
|
x
| where ε
rel
⁡
=
2
-
p (wherein p is the number of bits of fraction in the significand of the floating-point format f). In this case we consider the relative rounding error ε
rel
⁡.
(ii)
Ifx is in the range of denormalized number, then |
R
f
,
r
(
x
)
-
x
|
≤
ε
abs, where ε
abs
=
m
f
f (wherein m
f
f is the smallest nonzero positive denormalized floating-point number in the floating-point format f, which is also the gap between two neighboring denormalized numbers). In this case, we consider the absolute rounding error ε
abs.The rounding errors of these two cases can be unified as(16)
|
R
f
,
r
(
x
)
-
x
|
≤
max
⁡
(
ε
rel
⁡
·
|
x
|
,
ε
abs
)
.
Since max
⁡ is not a linear operation, we derive an overapproximation
(17)
|
R
f
,
r
(
x
)
-
x
|
≤
ε
rel
⁡
·
|
x
|
+
ε
abs
.
Furthermore, when b
≥
0, |
y
|
≤
b is equivalent to y
=
[
-
1,1
]
×
b. Hence,
(18)
R
f
,
r
(
x
)
-
x
=
[
-
1,1
]
(
ε
rel
⁡
·
|
x
|
+
ε
abs
)
;
that is,
(19)
R
f
,
r
(
x
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
-
ε
abs
,
ε
abs
]
.
In general, we could abstract floating-point operations into interval linear expressions in real number semantics. For example,
(20)
x
⊕
f
,
r
y
,
that is,
(21)
R
f
,
r
(
x
+
y
)
can be abstracted into
(22)
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
(
x
+
y
)
+
[
-
ε
abs
,
ε
abs
]
;
that is,
(23)
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
y
+
[
-
ε
abs
,
ε
abs
]
.The advantage of this kind of rounding mode insensitive floating-point abstractions is that the result is sound with respect to arbitrary rounding modes, sinceR
f
,
r
(
x
) always satisfies R
f
,
-
∞
(
x
)
≤
R
f
,
r
(
x
)
≤
R
f
,
+
∞
(
x
) while |
R
f
,
r
(
x
)
-
x
|
≤
ε
rel
⁡
·
|
x
|
+
ε
abs has already taken into account the extreme cases of r
=
-
∞ and r
=
+
∞. This is of practical importance, since we may not know the exact rounding mode for each floating-point operation. For example, C99 provides thefesetround() function to set the current rounding mode. Of course, when we know the exact rounding mode for the floating-point operation, we could make the floating-point abstraction more precise. For example, if the current rounding mode is toward nearest, then
(24)
R
f
,
n
(
x
)
=
[
1
-
ε
rel
⁡
2
,
1
+
ε
rel
⁡
2
]
×
x
+
[
-
ε
abs
2
,
ε
abs
2
]
.
In addition, if we know the range of x, we may also define more precise floating-point abstractions. For example, if we know that x is in the range of denormalized numbers, then
(25)
R
f
,
r
(
x
)
=
x
+
[
-
ε
abs
,
ε
abs
]
.For the sake of generality, in this paper, we use the following rounding mode insensitive floating-point abstraction:(26)
R
f
,
?
#
(
x
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
-
ε
abs
,
ε
abs
]
,
where we assume |
x
|
<
M
f
f.More clearly, we use the following abstraction for floating-point arithmetic:(27)
R
f
,
?
#
(
x
⊕
f
,
?
y
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
y
+
[
-
ε
abs
,
ε
abs
]
R
f
,
?
#
(
x
⊖
f
,
?
y
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
-
1
-
ε
rel
⁡
,
-
1
+
ε
rel
⁡
]
×
y
+
[
-
ε
abs
,
ε
abs
]
R
f
,
?
#
(
x
⊗
f
,
?
y
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
×
y
+
[
-
ε
abs
,
ε
abs
]
R
f
,
?
#
(
x
⊘
f
,
?
y
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
y
+
[
-
ε
abs
,
ε
abs
]
.
Specially, for a constant number c that appears in the source code, we use the following abstraction:
(28)
R
f
,
?
#
(
c
)
=
[
R
f
,
-
∞
#
(
c
)
,
R
f
,
+
∞
#
(
c
)
]
.


### 4.2. Symbolic Execution of Abstracted Floating-Point Programs

From Section4.1, we see that floating-point expressions can be soundly abstracted into real number expressions with interval coefficients. Since the multiplication x
×
y and division x
/
y are not linear expressions when both x and y are not constant numbers, in order to obtain linear expressions with interval coefficients, we replace y with its interval range denoted as [
y
_
,
y
¯
]. In symbolic execution, y is always an expression over the symbolic input values. We assume users provide the interval ranges for those symbolic input values. Then, all floating-point expressions can be abstracted as interval linear expressions. Therefore, the resulting path conditions of symbolic execution consist of interval linear constraints while the resulting path functions consist of interval linear expressions.Finally, the problems of computing the maximum output change and the maximum allowed input perturbation are reduced to solving a series of interval linear programming problems. For example, computing the maximum output change requires the solutions of the following interval linear programming problems:(29)
max
⁡
Σ
i
[
a
_
i
,
a
¯
i
]
×
x
i
+
Σ
i
[
a
_
i
′
,
a
¯
i
′
]
×
x
i
′
+
b
h
s
.
t
.
(
|
x
i
-
x
i
′
|
≤
δ
∧
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
)
hs
.
t
.
∧
⋀
k
Σ
i
[
A
_
k
i
,
A
¯
k
i
]
×
x
i
+
Σ
i
[
A
_
k
i
′
,
A
¯
k
i
′
]
×
x
i
′
≤
c
k
,
where Σ
i
[
a
_
i
,
a
¯
i
]
×
x
i
+
Σ
i
[
a
_
i
′
,
a
¯
i
′
]
×
x
i
′
+
b denotes the abstracted output change P
F
y
-
P
F
y
′ (or P
F
y
′
-
P
F
y) while ⋀
k
Σ
i
[
A
_
k
i
,
A
¯
k
i
]
×
x
i
+
Σ
i
[
A
_
k
i
′
,
A
¯
k
i
′
]
×
x
i
′
≤
c
k denotes the abstracted PC.Example 7.
Consider the self-composed programpiecewise_linear in Example 5. Suppose we would like to compute the maximum output change, given the input perturbation δ
=
0.1 over x. The self-composed program includes four paths overall. Let us consider for example the path that takes theelse branch in both the unprimed program P and the primed program P
′. Since x
,
y are offloat type, ε
rel
⁡
=
2
-
23 and ε
abs
=
2
-
149 for the 32-bit single precision floating-point format. We will have
(30)
l
PC
:
x
≥
R
f
,
-
∞
#
(
1.001
)
∧
x
′
≥
R
f
,
-
∞
#
(
1.001
)
P
F
y
:
[
1
-
2
-
23
,
1
+
2
-
23
]
P
F
y
:
×
[
R
f
,
-
∞
#
(
1001.0
)
,
R
f
,
+
∞
#
(
1001.0
)
]
×
x
P
F
y
:
+
[
-
2
-
149
,
2
-
149
]
P
F
y
′
:
[
1
-
2
-
23
,
1
+
2
-
23
]
P
F
y
:
×
[
R
f
,
-
∞
#
(
1001.0
)
,
R
f
,
+
∞
#
(
1001.0
)
]
×
x
′
P
F
y
:
+
[
-
2
-
149
,
2
-
149
]
.
Thus, we get the following interval linear programming problem:
(31)
max
⁡
[
1
-
2
-
23
,
1
+
2
-
23
]
×
[
R
f
,
-
∞
#
(
1001.0
)
,
R
f
,
+
∞
#
(
1001.0
)
]
×
x
+
[
-
1
-
2
-
23
,
-
1
+
2
-
23
]
×
[
R
f
,
-
∞
#
(
1001.0
)
,
R
f
,
+
∞
#
(
1001.0
)
]
×
x
′
+
[
-
2
-
148
,
2
-
148
]
(32)
h
s
.
t
.
x
-
x
′
≤
0.1
∧
-
x
+
x
′
≤
0.1
∧
-
x
≤
-
R
f
,
-
∞
#
(
1.001
)
∧
-
x
′
≤
-
R
f
,
-
∞
#
(
1.001
)
.
Solving the above interval linear programming problem by the method described in Section 2.2 will give us 100.10023889571252. After we deal with all other paths in the same way, we will find that 100.10023889571252 is the maximum output change with respect to the given input perturbation 0.1. Hence, the programpiecewise_linear in floating-point arithmetic is at least (
0.1,100.10023889571252
)-r
o
b
u
s
t.


## 4.1. Floating-Point Abstraction

In this subsection, we will explain how to abstract floating-point expressions into interval linear expressions over the reals.First, let us consider the upper bound on rounding errors due to one floating-point operation. LetR
f
,
r
(
x
) denote the floating-point rounding function that maps a real number x to a floating-point number (or a runtime error due to, for example, overflows) with respect to the floating-point format f and the rounding mode r. The amount of the rounding error due to R
f
,
r
(
x
) depends on the category of x. (i)
Ifx is in the range of normalized numbers, then |
R
f
,
r
(
x
)
-
x
|
≤
ε
rel
⁡
·
|
x
| where ε
rel
⁡
=
2
-
p (wherein p is the number of bits of fraction in the significand of the floating-point format f). In this case we consider the relative rounding error ε
rel
⁡.
(ii)
Ifx is in the range of denormalized number, then |
R
f
,
r
(
x
)
-
x
|
≤
ε
abs, where ε
abs
=
m
f
f (wherein m
f
f is the smallest nonzero positive denormalized floating-point number in the floating-point format f, which is also the gap between two neighboring denormalized numbers). In this case, we consider the absolute rounding error ε
abs.The rounding errors of these two cases can be unified as(16)
|
R
f
,
r
(
x
)
-
x
|
≤
max
⁡
(
ε
rel
⁡
·
|
x
|
,
ε
abs
)
.
Since max
⁡ is not a linear operation, we derive an overapproximation
(17)
|
R
f
,
r
(
x
)
-
x
|
≤
ε
rel
⁡
·
|
x
|
+
ε
abs
.
Furthermore, when b
≥
0, |
y
|
≤
b is equivalent to y
=
[
-
1,1
]
×
b. Hence,
(18)
R
f
,
r
(
x
)
-
x
=
[
-
1,1
]
(
ε
rel
⁡
·
|
x
|
+
ε
abs
)
;
that is,
(19)
R
f
,
r
(
x
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
-
ε
abs
,
ε
abs
]
.
In general, we could abstract floating-point operations into interval linear expressions in real number semantics. For example,
(20)
x
⊕
f
,
r
y
,
that is,
(21)
R
f
,
r
(
x
+
y
)
can be abstracted into
(22)
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
(
x
+
y
)
+
[
-
ε
abs
,
ε
abs
]
;
that is,
(23)
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
y
+
[
-
ε
abs
,
ε
abs
]
.The advantage of this kind of rounding mode insensitive floating-point abstractions is that the result is sound with respect to arbitrary rounding modes, sinceR
f
,
r
(
x
) always satisfies R
f
,
-
∞
(
x
)
≤
R
f
,
r
(
x
)
≤
R
f
,
+
∞
(
x
) while |
R
f
,
r
(
x
)
-
x
|
≤
ε
rel
⁡
·
|
x
|
+
ε
abs has already taken into account the extreme cases of r
=
-
∞ and r
=
+
∞. This is of practical importance, since we may not know the exact rounding mode for each floating-point operation. For example, C99 provides thefesetround() function to set the current rounding mode. Of course, when we know the exact rounding mode for the floating-point operation, we could make the floating-point abstraction more precise. For example, if the current rounding mode is toward nearest, then
(24)
R
f
,
n
(
x
)
=
[
1
-
ε
rel
⁡
2
,
1
+
ε
rel
⁡
2
]
×
x
+
[
-
ε
abs
2
,
ε
abs
2
]
.
In addition, if we know the range of x, we may also define more precise floating-point abstractions. For example, if we know that x is in the range of denormalized numbers, then
(25)
R
f
,
r
(
x
)
=
x
+
[
-
ε
abs
,
ε
abs
]
.For the sake of generality, in this paper, we use the following rounding mode insensitive floating-point abstraction:(26)
R
f
,
?
#
(
x
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
-
ε
abs
,
ε
abs
]
,
where we assume |
x
|
<
M
f
f.More clearly, we use the following abstraction for floating-point arithmetic:(27)
R
f
,
?
#
(
x
⊕
f
,
?
y
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
y
+
[
-
ε
abs
,
ε
abs
]
R
f
,
?
#
(
x
⊖
f
,
?
y
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
+
[
-
1
-
ε
rel
⁡
,
-
1
+
ε
rel
⁡
]
×
y
+
[
-
ε
abs
,
ε
abs
]
R
f
,
?
#
(
x
⊗
f
,
?
y
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
×
y
+
[
-
ε
abs
,
ε
abs
]
R
f
,
?
#
(
x
⊘
f
,
?
y
)
=
[
1
-
ε
rel
⁡
,
1
+
ε
rel
⁡
]
×
x
y
+
[
-
ε
abs
,
ε
abs
]
.
Specially, for a constant number c that appears in the source code, we use the following abstraction:
(28)
R
f
,
?
#
(
c
)
=
[
R
f
,
-
∞
#
(
c
)
,
R
f
,
+
∞
#
(
c
)
]
.


## 4.2. Symbolic Execution of Abstracted Floating-Point Programs

From Section4.1, we see that floating-point expressions can be soundly abstracted into real number expressions with interval coefficients. Since the multiplication x
×
y and division x
/
y are not linear expressions when both x and y are not constant numbers, in order to obtain linear expressions with interval coefficients, we replace y with its interval range denoted as [
y
_
,
y
¯
]. In symbolic execution, y is always an expression over the symbolic input values. We assume users provide the interval ranges for those symbolic input values. Then, all floating-point expressions can be abstracted as interval linear expressions. Therefore, the resulting path conditions of symbolic execution consist of interval linear constraints while the resulting path functions consist of interval linear expressions.Finally, the problems of computing the maximum output change and the maximum allowed input perturbation are reduced to solving a series of interval linear programming problems. For example, computing the maximum output change requires the solutions of the following interval linear programming problems:(29)
max
⁡
Σ
i
[
a
_
i
,
a
¯
i
]
×
x
i
+
Σ
i
[
a
_
i
′
,
a
¯
i
′
]
×
x
i
′
+
b
h
s
.
t
.
(
|
x
i
-
x
i
′
|
≤
δ
∧
⋀
1
≤
j
≤
n
,
j
≠
i
x
j
=
x
j
′
)
hs
.
t
.
∧
⋀
k
Σ
i
[
A
_
k
i
,
A
¯
k
i
]
×
x
i
+
Σ
i
[
A
_
k
i
′
,
A
¯
k
i
′
]
×
x
i
′
≤
c
k
,
where Σ
i
[
a
_
i
,
a
¯
i
]
×
x
i
+
Σ
i
[
a
_
i
′
,
a
¯
i
′
]
×
x
i
′
+
b denotes the abstracted output change P
F
y
-
P
F
y
′ (or P
F
y
′
-
P
F
y) while ⋀
k
Σ
i
[
A
_
k
i
,
A
¯
k
i
]
×
x
i
+
Σ
i
[
A
_
k
i
′
,
A
¯
k
i
′
]
×
x
i
′
≤
c
k denotes the abstracted PC.Example 7.
Consider the self-composed programpiecewise_linear in Example 5. Suppose we would like to compute the maximum output change, given the input perturbation δ
=
0.1 over x. The self-composed program includes four paths overall. Let us consider for example the path that takes theelse branch in both the unprimed program P and the primed program P
′. Since x
,
y are offloat type, ε
rel
⁡
=
2
-
23 and ε
abs
=
2
-
149 for the 32-bit single precision floating-point format. We will have
(30)
l
PC
:
x
≥
R
f
,
-
∞
#
(
1.001
)
∧
x
′
≥
R
f
,
-
∞
#
(
1.001
)
P
F
y
:
[
1
-
2
-
23
,
1
+
2
-
23
]
P
F
y
:
×
[
R
f
,
-
∞
#
(
1001.0
)
,
R
f
,
+
∞
#
(
1001.0
)
]
×
x
P
F
y
:
+
[
-
2
-
149
,
2
-
149
]
P
F
y
′
:
[
1
-
2
-
23
,
1
+
2
-
23
]
P
F
y
:
×
[
R
f
,
-
∞
#
(
1001.0
)
,
R
f
,
+
∞
#
(
1001.0
)
]
×
x
′
P
F
y
:
+
[
-
2
-
149
,
2
-
149
]
.
Thus, we get the following interval linear programming problem:
(31)
max
⁡
[
1
-
2
-
23
,
1
+
2
-
23
]
×
[
R
f
,
-
∞
#
(
1001.0
)
,
R
f
,
+
∞
#
(
1001.0
)
]
×
x
+
[
-
1
-
2
-
23
,
-
1
+
2
-
23
]
×
[
R
f
,
-
∞
#
(
1001.0
)
,
R
f
,
+
∞
#
(
1001.0
)
]
×
x
′
+
[
-
2
-
148
,
2
-
148
]
(32)
h
s
.
t
.
x
-
x
′
≤
0.1
∧
-
x
+
x
′
≤
0.1
∧
-
x
≤
-
R
f
,
-
∞
#
(
1.001
)
∧
-
x
′
≤
-
R
f
,
-
∞
#
(
1.001
)
.
Solving the above interval linear programming problem by the method described in Section 2.2 will give us 100.10023889571252. After we deal with all other paths in the same way, we will find that 100.10023889571252 is the maximum output change with respect to the given input perturbation 0.1. Hence, the programpiecewise_linear in floating-point arithmetic is at least (
0.1,100.10023889571252
)-r
o
b
u
s
t.


## 5. Implementation and Experimental Results

We have implemented a robustness analysis tool RAFP, based on the symbolic execution and floating-point abstraction techniques presented in Section4. Given an input perturbation δ over one input variable of the program, RAFP can compute the maximum output change. Furthermore, if the user also provides a candidate output change ϵ and would like to check whether the program is (
δ
,
ϵ
)-r
o
b
u
s
t, RAFP will check this property during the process of computing maximum output change and will stop once one path violating the property is found. Also, given an output change ϵ, RAFP can compute the maximum allowed input perturbation for floating-point programs. RAFP is built on top of Symbolic PathFinder (SPF) [23] which is a symbolic execution engine for Java programs. We use SPF to extract the path conditions together with the associated path functions. For linear programming, RAFP makes use of the Java Binding for GLPK (GNU Linear programming kit) called GLPK-Java [24].To conduct experiments on checking robustness properties of floating-point programs via software model checking, we choose CBMC (C Bounded Model Checker) [14] which implements bounded model checking for ANSI-C programs using SAT/SMT solvers. CBMC utilizes a bit-precise modeling for floating-point operations and employs a sound and complete decision procedure for floating-point arithmetic. CBMC provides an option- -floatbv to use IEEE floating point arithmetic and options for choosing rounding modes. However, CBMC does not support to use different rounding modes for the floating-point operations in the same program. In other words, all floating-point operations in a program are of the same rounding mode during the analysis. We use the default rounding mode- -round-to-nearest during our experiments. Moreover, CBMC provides_  _CPROVER_assume() and_  _CPROVER_assert() statements, which are needed for robustness analysis of self-composed programs. Both statements take Boolean conditions. The_  _CPROVER_assume() statement restricts that the program traces should satisfy the assumed condition. For the_  _CPROVER_assert() statement, CBMC will check whether the asserted condition holds true for all runs of the program.We have conducted experiments on a selection of benchmark examples using both RAFP and CBMC. Table1 shows the comparison of performance and the resulting output changes. The column “δ
in” shows the considered input perturbation over one input variable of the program. The column “ϵ
max
⁡” shows the resulting maximum output change computed by RAFP with respect to the given input perturbation. The column “ϵ
unr” gives the largest possible output change that we have tried with CBMC such that the program is not (
δ
in
,
ϵ
unr
)-r
o
b
u
s
t with respect to the given input perturbation. The column “ϵ
r” gives the smallest output change that we have tried with CBMC such that the program is (
δ
in
,
ϵ
r
)-r
o
b
u
s
t with respect to the given input perturbation (Note that CBMC can be used only to check whether a program is (
δ
,
ϵ
)-r
o
b
u
s
t and can not be used to compute the amount of output change with respect to the given input perturbation. During our experiments, we try CBMC with different candidate values of ϵ to find ϵ
unr and ϵ
r). Since CBMC uses the same rounding mode for all floating-point operations in the same program during the analysis, the output change is always 0 when the given input perturbation is 0. Hence, for those rows that specify input perturbation as 0, we do not need to run CBMC and thus we mark the table entry with ⋆ in this case. Our tool RAFP utilizes rounding mode insensitive floating-point abstraction and thus in principle it holds that ϵ
max
⁡
≥
ϵ
r
≥
ϵ
unr, which is confirmed by the experimental results.Table 1
Experimental results for benchmark examples.
Program
δ
in
RAFP
CBMC
ϵ
max
⁡
t (ms)
ϵ
unr
t (ms)
ϵ
r
t (ms)
p
iecewise_linear
0
2.3889571220452006
e
-
4
29
⋆
⋆
⋆
⋆
0.01
10.010238895712442
33
?
>1 h
?
>1 h
0.1
100.10023889571252
34
?
>1 h
?
>1 h
Max1
0
1.4012987984203268
e
-
45
44
⋆
⋆
⋆
⋆
0.01
0.010000000000000007
55
0.0099
215
0.01
16066
0.1
0.10000000000000006
54
0.099
227
0.1
16262
MorePaths
0
1.0000000027939686
59
⋆
⋆
⋆
⋆
0.01
1.0000000027939686
79
0.99
279
1.0
379
0.1
1.0000000027939686
86
0.99
314
1.0
471
Orientation
0
2.0
55
⋆
⋆
⋆
⋆
0.1
e
-
8
*
E
2.0
68
0
1840
1.0
6845
0.1
e
-
3
*
E
2.0
73
0
1338
1.0
13327
0.1
e
-
2
*
E
2.0
80
1.0
14165
2.0
413
Filtered_Orientation
0
1.0
54
⋆
⋆
⋆
⋆
0.1
e
-
8
*
E
1.0
70
0
1044
1.0
29641
0.1
e
-
2
*
E
1.0
73
0
898
1.0
30261
0.1
*
E
2.0
75
0
719
1.0
26463
E
2.0
68
1.0
13206
2.0
654The programpiecewise_linear corresponds to the program shown in Example 4.   Max1,MorePaths come from JPF Continuity [25].Max1 is a floating-point program that implements max
⁡
(
x
,
y
), and thus it is (
δ
in
,
δ
in
)-r
o
b
u
s
t.MorePaths is a floating-point program that involves both a step function and a max
⁡ function, and thus it is (
δ
in
,
1.0
)-r
o
b
u
s
t for all δ
in
≤
1.0. Orientation (which corresponds to the program shown in Figure 1) together withFiltered_Orientation are extracted from the computational geometry algorithms library CGAL [26] and address robust geometric computation.Filtered_Orientation is an improved version of  Orientation via static filter technique. The approximate result of computing the sign of a determinant is compared with a given positive filter bound E (rather than compared with zero). When the approximate result is in the interval [
-
E
,
E
],Filtered_Orientation gives 0. During our experiments, we set E
=
1.5
e
-
5 (and for the sake of comparison, we express the input perturbation in terms of E also forOrientation although here E does not appear). The outputs ofOrientation andFiltered_Orientation are always −1 (negative), 0 (zero), or 1 (positive). Hence, in Table 1, the resulting output changes for these two programs are always 0, 1.0, or 2.0. From Table 1, we could find thatFiltered_Orientation is more robust thanOrientation. For example, given the input perturbation δ
=
0.1
e
-
2
*
E, CBMC finds that for ϵ
=
1.0,Filtered_Orientation is robust whileOrientation is not. Similarly, given the input perturbation δ
=
0.1
e
-
2
*
E, RAFP gives ϵ
max
⁡
=
1.0 fors we mark the table entry withFiltered_Orientation but gives ϵ
max
⁡
=
2.0 forOrientation.The column “t
(
ms
)” presents the analysis times in milliseconds when the analyzers run on a 2.5 GHz PC with 4 GB of RAM running Windows 7. (RAFP runs further on a Java Virtual Machine (JVM) while CBMC runs further on a virtual machine VMWare running Fedora 12.) From Table 1, we could see that RAFP outperforms CBMC in time efficiency. Especially forpiecewise_linear, CBMC could not even finish the analysis process in 1 hour. The low efficiency of CBMC is because that CBMC uses a sound and complete decision procedure for floating-point arithmetic. Especially, the multiplication and division floating-point operations may generate formulae that are expensive to decide and quite hard for SAT solvers to solve [27]. Hence, checking robustness properties of floating-point programs via CBMC may have limitations in scalability due to the current expensive decision procedures for floating-point logic. During our experiments, the approach via symbolic execution of abstracted floating-point programs is much more efficient. In principle, symbolic execution may suffer from the path explosion problem. However, the recent success of symbolic execution tools such as KLEE [28] on analyzing large-scale programs [17] makes this approach promising.


## 6. Related Work

### 6.1. Robustness Analysis of Programs

Robustness is a standard correctness property for control systems [1]. Robustness analysis of programs has received increasing attention in the recent years. Majumdar and Saha [5] took a first step toward analyzing the robustness of programs in control systems. They also utilized symbolic execution and optimization techniques to compute the maximum difference in program outputs with respect to the given input perturbation. However, they assumed exact real arithmetic in the program. Continuity as one aspect of robustness for software was firstly considered in [29]. Recently, Chaudhuri et al. presented logic-based mostly automated methods to determine whether a program is continuous [2] or Lipschitz continuous [3, 4], and more recently to determine whether a decision-making program is consistent under uncertainty [30]. Quite recently, Shahrokni and Feldt [31] conducted a systematic review of software robustness. However, much existing work on robustness analysis does not handle floating-point arithmetic in the program. Bushnell [25] presented a symbolic execution based approach to identify continuties and discontinuties associated with path condition boundaries for floating-point software, but it did not consider the true floating-point semantics. Besides, Gazeau et al. [32] presented a nonlocal method for proving the robustness of floating-point programs but which needs much manual work. Recently, Goubault and Putot [33] proposed an abstract interpretation based robustness analysis method for finite precision implementations.


### 6.2. Safety Analysis of Floating Point Programs

Monniaux [7] described common pitfalls in analyzing and verifying floating-point programs. Abstract interpretation [34] based techniques have shown quite successful on analysis of floating-point programs. In [35], Goubault analyzed the origin of the loss of precision in floating-point programs based on abstract interpretation. Following this direction, a static analyzer FLUCTUAT [36] was developed. The abstract interpretation based static analyzer ASTRÉE [15] checks for floating-point run-time errors based on the computed set of reachable values for floating-point variables. As in ASTRÉE, we rely on the floating-point abstraction technique of [21] to soundly abstract floating-point expressions into ones over the field of reals. Chen et al. [37, 38] utilized interval linear constraints to design numerical abstract domains and to construct sound floating-point implementations [39]. Ivančić et al. [40] used bounded model checking based on SMT solvers to detect numerical instabilities in floating-point programs, based on a mixed integer-real model for floating-point variables and operations. Brain et al. [41] recently improved the bit-precise decision procedure for the theory of floating-point arithmetic based on a strict lifting of the conflict-driven clause learning algorithm in modern SAT solvers to abstract domains. Barr et al. [42] presented a method to automatically detect the floating-point exception through symbolic execution.


### 6.3. Self-Composition

The idea of self-composition is firstly used in the field of secure information flow [10, 11], to characterize noninterference. Terauchi and Aiken [12] proposed the type-directed transformation approach to make self-composition work in practice with off-the-shelf automatic safety analysis tools. Recently, Barthe et al. [43] proposed a general notion of product program that is beneficial to relational verification, which could be considered as the generalization of self-composition. Kovacs et al. [44] presented a general method to analyze 2-hypersafety properties by applying abstract interpretation on the self-compositions of the control flow graphs of programs.


## 6.1. Robustness Analysis of Programs

Robustness is a standard correctness property for control systems [1]. Robustness analysis of programs has received increasing attention in the recent years. Majumdar and Saha [5] took a first step toward analyzing the robustness of programs in control systems. They also utilized symbolic execution and optimization techniques to compute the maximum difference in program outputs with respect to the given input perturbation. However, they assumed exact real arithmetic in the program. Continuity as one aspect of robustness for software was firstly considered in [29]. Recently, Chaudhuri et al. presented logic-based mostly automated methods to determine whether a program is continuous [2] or Lipschitz continuous [3, 4], and more recently to determine whether a decision-making program is consistent under uncertainty [30]. Quite recently, Shahrokni and Feldt [31] conducted a systematic review of software robustness. However, much existing work on robustness analysis does not handle floating-point arithmetic in the program. Bushnell [25] presented a symbolic execution based approach to identify continuties and discontinuties associated with path condition boundaries for floating-point software, but it did not consider the true floating-point semantics. Besides, Gazeau et al. [32] presented a nonlocal method for proving the robustness of floating-point programs but which needs much manual work. Recently, Goubault and Putot [33] proposed an abstract interpretation based robustness analysis method for finite precision implementations.


## 6.2. Safety Analysis of Floating Point Programs

Monniaux [7] described common pitfalls in analyzing and verifying floating-point programs. Abstract interpretation [34] based techniques have shown quite successful on analysis of floating-point programs. In [35], Goubault analyzed the origin of the loss of precision in floating-point programs based on abstract interpretation. Following this direction, a static analyzer FLUCTUAT [36] was developed. The abstract interpretation based static analyzer ASTRÉE [15] checks for floating-point run-time errors based on the computed set of reachable values for floating-point variables. As in ASTRÉE, we rely on the floating-point abstraction technique of [21] to soundly abstract floating-point expressions into ones over the field of reals. Chen et al. [37, 38] utilized interval linear constraints to design numerical abstract domains and to construct sound floating-point implementations [39]. Ivančić et al. [40] used bounded model checking based on SMT solvers to detect numerical instabilities in floating-point programs, based on a mixed integer-real model for floating-point variables and operations. Brain et al. [41] recently improved the bit-precise decision procedure for the theory of floating-point arithmetic based on a strict lifting of the conflict-driven clause learning algorithm in modern SAT solvers to abstract domains. Barr et al. [42] presented a method to automatically detect the floating-point exception through symbolic execution.


## 6.3. Self-Composition

The idea of self-composition is firstly used in the field of secure information flow [10, 11], to characterize noninterference. Terauchi and Aiken [12] proposed the type-directed transformation approach to make self-composition work in practice with off-the-shelf automatic safety analysis tools. Recently, Barthe et al. [43] proposed a general notion of product program that is beneficial to relational verification, which could be considered as the generalization of self-composition. Kovacs et al. [44] presented a general method to analyze 2-hypersafety properties by applying abstract interpretation on the self-compositions of the control flow graphs of programs.


## 7. Conclusion

We have proposed a self-composition based approach for robustness analysis of programs, which enables making use of off-the-shelf automatic reachability analysis tools to analyze robustness properties of programs. Then, we have shown how to use software model checking and symbolic execution techniques on self-composed programs to analyze program robustness properties. In particular, we have considered the robustness analysis problem of floating-point programs. To deal with floating-point arithmetic during symbolic execution, we have utilized a rounding mode insensitive floating-point abstraction to abstract floating-point expressions into interval linear expressions in exact real arithmetic. On this basis, the maximum output change (when given the input perturbation) or maximum allowed input perturbation (when given the input perturbation) are computed based on symbolic execution and interval linear programming for abstracted floating-point programs. Experimental results of our prototype implementation are encouraging.It remains for future work to exploit the intrinsic symmetry of self-composed programs to reduce the number of considered paths during robustness analysis. We also plan to improve the prototype implementation and to conduct more experiments on larger realistic floating-point programs.



---

*Source: 789213-2014-05-20.xml*