A Comparative Study of Arithmetic Constraints
on Integer Intervals

arXiv:cs/0403016v1 [] 12 Mar 2004

Krzysztof R. Apt1,2 and Peter Zoeteweij1
1

CWI, P.O. Box 94079, 1090 GB Amsterdam, the Netherlands
2
University of Amsterdam, the Netherlands

Abstract. We propose here a number of approaches to implement constraint propagation for arithmetic constraints on integer intervals. To
this end we introduce integer interval arithmetic. Each approach is explained using appropriate proof rules that reduce the variable domains.
We compare these approaches using a set of benchmarks.

1
1.1

Preliminaries
Introduction

The subject of arithmetic constraints on reals has attracted a great deal of attention in the literature. For some reason arithmetic constraints on integer intervals
have not been studied even though they are supported in a number of constraint
programming systems. In fact, constraint propagation for them is present in
ECLi PSe , SICStus Prolog, GNU Prolog, ILOG Solver and undoubtedly most of
the systems that support constraint propagation for linear constraints on integer
intervals. Yet, in contrast to the case of linear constraints — see notably [5] —
we did not encounter in the literature any analysis of this form of constraint
propagation.
In this paper we study these constraints in a systematic way. It turns out
that in contrast to linear constraints on integer intervals there are a number of
natural approaches to constraint propagation for these constraints.
To define them we introduce integer interval arithmetic that is modeled after
the real interval arithmetic, see e.g., [6]. There are, however, essential differences
since we deal with integers instead of reals. For example, multiplication of two
integer intervals does not need to be an integer interval. In passing by we show
that using integer interval arithmetic we can also define succinctly the wellknown constraint propagation for linear constraints on integer intervals. In the
second part of the paper we compare the proposed approaches by means of a set
of benchmarks.
1.2

Constraint Satisfaction Problems

We review here the standard concepts of a constraint and of a constraint satisfaction problem. Consider a sequence of variables X := x1 , . . ., xn where n ≥ 0,

with respective domains D1 , . . ., Dn associated with them. So each variable xi
ranges over the domain Di . By a constraint C on X we mean a subset of
D1 × . . . × Dn . Given an element d := d1 , . . ., dn of D1 × . . . × Dn and a subsequence Y := xi1 , . . ., xiℓ of X we denote by d[Y ] the sequence di1 , . . ., diℓ . In
particular, for a variable xi from X, d[xi ] denotes di .
A constraint satisfaction problem, in short CSP, consists of a finite sequence of variables X with respective domains D, together with a finite set C of
constraints, each on a subsequence of X. We write it as hC ; x1 ∈ D1 , . . ., xn ∈
Dn i, where X := x1 , . . ., xn and D := D1 , . . ., Dn .
By a solution to hC ; x1 ∈ D1 , . . ., xn ∈ Dn i we mean an element d ∈ D1 ×
. . .×Dn such that for each constraint C ∈ C on a sequence of variables X we have
d[X] ∈ C. We call a CSP consistent if it has a solution and inconsistent if it
does not. Two CSPs with the same sequence of variables are called equivalent
if they have the same set of solutions. In what follows we consider CSPs the
constraints of which are defined in a simple language and identify the syntactic
description of a constraint with its meaning being the set of tuples that satisfy
it.
We view constraint propagation as a process of transforming CSPs that
maintains their equivalence. In what follows we define this process by means of
proof rules that act of CSPs and preserve equivalence. An interested reader can
consult [1] for a precise explanation of this approach to describing constraint
propagation.
1.3

Arithmetic Constraints

To define the arithmetic constraints use the alphabet that comprises
–
–
–
–

variables,
two constants, 0 and 1,
the unary minus function symbol ‘−’,
three binary function symbols, ‘+’,‘−’and ‘·’, all written in the infix notation.

By an arithmetic expression we mean a term formed in this alphabet and
by an arithmetic constraint a formula of the form
s op t,
where s and t are arithmetic expressions and op ∈ {<, ≤, =, 6=, ≥, >}. For example
x5 · y 2 · z 4 + 3x · y 3 · z 5 ≤ 10 + 4x4 · y 6 · z 2 − y 2 · x5 · z 4
(1)
is an arithmetic constraint. Here x5 is an abbreviation for x · x · x · x · x and
similarly with the other expressions. If ‘·’ is not used in an arithmetic constraint,
we call it a linear constraint.
By an extended arithmetic expression we mean a term formed in the
√
above alphabet extended by the unary function symbols ‘.n ’ and ‘ n .’ for each

n ≥ 1 and the binary function symbol ‘/’ written in the infix notation. For
example
p
3
(y 2 · z 4 )/(x2 · u5 )
(2)

is an extended arithmetic expression. Here, in contrast to the above x5 is a term
obtained by applying the function symbol ‘.5 ’ to the variable x. The extended
arithmetic expressions will be used only to define constraint propagation for the
arithmetic constraints.
Fix now some arbitrary linear ordering ≺ on the variables of the language.
By a monomial we mean an integer or a term of the form
a · xn1 1 · . . . · xnk k
where k > 0, x1 , . . ., xk are different variables ordered w.r.t. ≺, and a is a nonzero integer and n1 , . . ., nk are positive integers. We call then xn1 1 · . . . · xnk k the
power product of this monomial.
Next, by a polynomial we mean a term of the form
n
Σi=1
mi ,

where n > 0, at most one monomial mi is an integer, and the power products
of the monomials m1 , . . ., mn are pairwise different. Finally, by a polynomial
constraint we mean an arithmetic constraint of the form s op b, where s is a
polynomial with no monomial being an integer, op ∈ {<, ≤, =, 6=, ≥, >}, and b is
an integer. It is clear that by means of appropriate transformation rules we can
transform each arithmetic constraint to a polynomial constraint. For example,
assuming the ordering x ≺ y ≺ z on the variables, the arithmetic constraint (1)
can be transformed to the polynomial constraint
2x5 · y 2 · z 4 − 4x4 · y 6 · z 2 + 3x · y 3 · z 5 ≤ 10
So, without loss of generality, from now on we shall limit our attention to the
polynomial constraints.
Next, let us discuss the domains over which we interpret the arithmetic constraints. By an integer interval , or an interval in short, we mean an expression of the form
[a..b]
where a and b are integers; [a..b] denotes the set of all integers between a and b,
including a and b. If a > b, we call [a..b] the empty interval and denote it by
∅. Finally, by a range we mean an expression of the form
x∈I
where x is a variable and I is an interval.

2

Integer Set Arithmetic

To reason about the arithmetic constraints we employ a generalization of the
arithmetic operations to the sets of integers.

2.1

Definitions

For X, Y sets of integers we define the following operations:
– addition:
X + Y := {x + y | x ∈ X, y ∈ Y },
– subtraction:
X − Y := {x − y | x ∈ X, y ∈ Y },
– multiplication:
X · Y := {x · y | x ∈ X, y ∈ Y },
– division:
X/Y := {u ∈ Z | ∃x ∈ X∃y ∈ Y u · y = x},
– exponentiation:
X n := {xn | x ∈ X},
for each natural number n > 0,
– root extraction:
√
n
X := {x ∈ Z | xn ∈ X},
for each natural number n > 0.
All the operations except division are defined in the expected way. We shall
return to it at the end of Section 6. At the moment it suffices to note the division
operation is defined for all sets of integers, including Y = ∅ and Y = {0}. This
division operation corresponds to the following division operation on the sets of
reals introduced in [8]:
X
:= {u ∈ R | ∃x ∈ X∃y ∈ Y u · y = x}.
Y
For a (n integer or real) number a and op ∈ {+, −, ·, /} we identify a op X with
{a} op X and X op a with X op {a}.
To present the rules we are interested in we shall also use the addition and
division operations on the sets of real numbers. Addition is defined in the same
way as for the sets of integers, and division is defined above. In [6] it is explained
how to implement these operations.
Further, given a set A of integers or reals, we define
≤

A := {x ∈ Z | ∃a ∈ A x ≤ a},

≥

A := {x ∈ Z | ∃a ∈ A x ≥ a}.

When limiting our attention to intervals of integers the following simple observation is of importance.
Note 1. For X, Y integer intervals and a an integer the following holds:
– X ∩ Y , X + Y, X − Y are integer intervals.

–
–
–
–
–
–

X/{a} is an integer interval.
X · Y does not have to be an integer interval, even if X = {a} or Y = {a}.
X/Y does not have to be an integer interval.
For each n > 1 √
X n does not have to be an integer interval.
For odd n > 1 n√X is an integer interval.
For even n > 1 n X is an integer interval or a disjoint union of two integer
intervals.
✷
For example we have
[2..4] + [3..8] = [5..12],
[3..7] − [1..8] = [−5..6],
[3..3] · [1..2] = {3, 6},
[3..5]/[−1..2] = {−5, −4, −3, 2, 3, 4, 5},
[−3..5]/[−1..2] = Z,
[1..2]2 = {1, 4},
p
3
[−30..100] = [−3..4],
p
2
[−100..9] = [−3..3],

p
2
[1..9] = [−3.. − 1] ∪ [1..3].

To deal with the problem that non-interval domains can be produced by some
of the operations we introduce the following operation on the subsets of the set
of the integers Z:

smallest integer interval containing X if X is finite,
int(X) :=
Z
otherwise.
For example int([3..5]/[−1..2]) = [−5..5] and int([−3..5]/[−1..2]) = Z.
2.2

Implementation

To define constraint propagation for the arithmetic constraints on integer intervals we shall use the integer set arithmetic, mainly limited to the integer
intervals. This brings us to the discussion of how to implement the introduced
operations on the integer intervals. Since we are only interested in maintaining
the property that the sets remain integer intervals or the set of integers Z we
shall clarify how to implement the intersection, addition, subtraction and root
extraction operations of the integer intervals and the int(.) closure of the multiplication, division and exponentiation operations on the integer intervals. The
case when one of the intervals is empty is easy to deal with. So we assume that
we deal with non-empty intervals [a..b] and [c..d], that is a ≤ b and c ≤ d.

Intersection, addition and subtraction It is easy to see that
[a..b] ∩ [c..d] = [max(a, c)..min(b, d)],
[a..b] + [c..d] = [a + c .. b + d],
[a..b] − [c..d] = [a − d .. b − c].
So the interval intersection, addition, and subtraction are straightforward to
implement.
Root extraction The outcome of the root extraction operator applied to an integer interval will be an integer interval or a disjoint union of two integer intervals.
We shall explain in Section 4 why it is advantageous not to apply int(.) to the
outcome. This operator can be implemented by means of the following case
analysis.
Case 1. Suppose n is odd. Then
k
p
√  j√
n
n
[a..b] = [ n a ..
b ].
Case 2. Suppose n is even and b < 0. Then
p
n
[a..b] = ∅.

Case 3. Suppose n is even and b ≥ 0. Then
j √ k
l √ m
l √ m j √ k
p
n
n
n
n
n
[a..b] = [− | b| .. − | a+ | ] ∪ [ | a+ | .. | b| ]
where a+ := max (0, a).
Multiplication For the remaining operations we only need to explain how to
implement the int(.) closure of the outcome. First note that
int([a..b] · [c..d]) = [min(A)..max(A)],
where A = {a · c, a · d, b · c, b · d}.
Using an appropriate case analysis we can actually compute the bounds of
int([a..b] · [c..d]) directly in terms of the bounds of the constituent intervals.
Division In contrast, the int(.) closure of the interval division is not so straightforward to compute. The reason is that, as we shall see in a moment, we cannot
express the result in terms of some simple operations on the interval bounds.
Consider non-empty integer intervals [a..b] and [c..d]. In analyzing the outcome of int([a..b]/[c..d]) we distinguish the following cases.
Case 1. Suppose 0 ∈ [a..b] and 0 ∈ [c..d].
Then by definition int([a..b]/[c..d]) = Z. For example,
int([−1..100]/[−2..8]) = Z.

Case 2. Suppose 0 6∈ [a..b] and c = d = 0.
Then by definition int([a..b]/[c..d]) = ∅. For example,
int([10..100]/[0..0]) = ∅.
Case 3. Suppose 0 6∈ [a..b] and c < 0 and 0 < d.
It is easy to see that then
int([a..b]/[c..d]) = [−e..e],
where e = max(|a|, |b|). For example,
int([−100.. − 10]/[−2..5]) = [−100..100].
Case 4. Suppose 0 6∈ [a..b] and either c = 0 and d 6= 0 or c 6= 0 and d = 0.
Then int([a..b]/[c..d]) = int([a..b]/([c..d] − {0})). For example
int([1..100]/[−7..0]) = int([1..100]/[−7.. − 1]).
This allows us to reduce this case to Case 5 below.
Case 5. Suppose 0 6∈ [c..d].
This is the only case when we need to compute int([a..b]/[c..d]) indirectly.
First, observe that we have
int([a..b]/[c..d]) ⊆ [⌈min(A)⌉ .. ⌊max(A)⌋],
where A = {a/c, a/d, b/c, b/d}.
However, the equality does not need to hold here. Indeed, note for example
that int([155..161]/[9..11]) = [16..16], whereas for A = {155/9, 155/11, 161/9,
161/11} we have ⌈min(A)⌉ = 15 and ⌊max(A)⌋ = 17. The problem is that the
value 16 is obtained by dividing 160 by 10 and none of these two values is an
interval bound.
This complication can be solved by preprocessing the interval [c..d] so that
its bounds are actual divisors of an element of [a..b]. First, we look for the least
c′ ∈ [c..d] such that ∃x ∈ [a..b] ∃u ∈ Z u · c′ = x. Using a case analysis, the latter
property can be established without search. Suppose for example that a > 0 and
c > 0. In this case, if c′ · ⌊ cb′ ⌋ ≥ a, then c′ has the required property. Similarly,
we look for the largest d′ ∈ [c..d] for which an analogous condition holds. Now
int ([a..b]/[c..d]) = [⌈min(A)⌉..⌊max (A)⌋], where A = {a/c′ , a/d′ , b/c′ , b/d′ }.
Exponentiation The int(.) closure of the interval exponentiation is straightforward to implement by distinguishing the following cases.
Case 1. Suppose n is odd. Then
int([a..b]n ) = [an .. bn ].
Case 2. Suppose n is even and 0 ≤ a. Then
int([a..b]n ) = [an .. bn ].

Case 3. Suppose n is even and b ≤ 0. Then
int([a..b]n ) = [bn .. an ].
Case 4. Suppose n is even and a < 0 and 0 < b. Then
int([a..b]n ) = [0..max(an , bn )].
2.3

Correctness Lemma

Given now an extended arithmetic expression s each variable of which ranges
over an integer interval, we define int(s) as the integer interval or the set Z
obtained by systematically replacing each function symbol by the application of
the int(.) operation to the correspondingp
integer set operation. For example, for
the extended arithmetic expression s := 3 (y 2 · z 4 )/(x2 · u5 ) of (2) we have
p
int(s) = int( 3 int(int(Y 2 ) · int(Z 4 ))/int(int(X 2 ) · int(U 5 ))),

where x ranges over X, etc.
The discussion in the previous subsection shows how to compute int(s) given
an extended arithmetic expression s and the integer interval domains of its variables.
The following lemma is crucial for our considerations. It is a counterpart
of the so-called ‘Fundamental Theorem of Interval Arithmetic’ established in
[7]. Because we deal here with the integer domains an additional assumption is
needed to establish the desired conclusion.
Lemma 1 (Correctness). Let s be an extended arithmetic expression with the
variables x1 , . . ., xn . Assume that each variable xi of s ranges over an integer
interval Xi . Choose ai ∈ Xi for i ∈ [1..n] and denote by s(a1 , . . ., an ) the result
of replacing in s each occurrence of a variable xi by ai .
Suppose that each subexpression of s(a1 , . . ., an ) evaluates to an integer. Then
the result of evaluating s(a1 , . . ., an ) is an element of int(s).
Proof. The proof follows by a straightforward induction on the structure of s.
✷

3

An Intermezzo: Constraint Propagation for Linear
Constraints

Even though we focus here on arithmetic constraints on integer intervals, it is
helpful to realize that the integer interval arithmetic is also useful to define
in a succinct way the well-known rules for constraint propagation for linear
n
constraints. To this end consider first a constraint Σi=1
ai · xi = b, where n ≥ 0,
a1 , . . ., an are non-zero integers, x1 , . . ., xn are different variables, and b is an
integer. To reason about it we can use the following rule parametrized by j ∈
[1..n]:

LINEAR EQUALITY
n
hΣi=1
ai
n
hΣi=1 ai

where

· xi = b ; x1 ∈ D1 , . . ., xn ∈ Dn i
· xi = b ; x1 ∈ D1′ , . . ., xn ∈ Dn′ i

– for i 6= j

Di′ := Di ,

–



Dj′ := Dj ∩ int (b − Σi∈[1..n]−{j} ai · xi )/aj .
Note that by virtue of Note 1
Dj′ = Dj ∩ (b − Σi∈[1..n]−{j} int(ai · Di ))/aj .

To see that this rule preserves equivalence suppose that for some d1 ∈
n
D1 , . . ., dn ∈ Dn we have Σi=1
ai · di = b. Then for j ∈ [1..n] we have
dj = (b − Σi∈[1..n]−{j} ai · di )/aj
which by the Correctness Lemma 1 implies that


dj ∈ int (b − Σi∈[1..n]−{j} ai · xi )/aj ,

i.e., dj ∈ Dj′ .
n
Next, consider a constraint Σi=1
ai · xi ≤ b, where a1 , . . ., an , x1 , . . ., xn and b
are as above. To reason about it we can use the following rule parametrized by
j ∈ [1..n]:
LINEAR INEQUALITY
n
hΣi=1
ai
n
hΣi=1 ai

where
– for i 6= j

· xi ≤ b ; x1 ∈ D1 , . . ., xn ∈ Dn i
· xi ≤ b ; x1 ∈ D1′ , . . ., xn ∈ Dn′ i

Di′ := Di ,

–
Dj′ := Dj ∩ (≤ int(b − Σi∈[1..n]−{j} ai · xi )/aj )
To see that this rule preserves equivalence suppose that for some d1 ∈
n
D1 , . . ., dn ∈ Dn we have Σi=1
ai · di ≤ b. Then aj · dj ≤ b − Σi∈[1..n]−{j} ai · di .
By the Correctness Lemma 1
b − Σi∈[1..n]−{j} ai · di ∈ int(b − Σi∈[1..n]−{j} ai · xi ),
so by definition
and consequently

aj · dj ∈≤ int(b − Σi∈[1..n]−{j} ai · xi )
dj ∈≤ int(b − Σi∈[1..n]−{j} ai · xi )/aj

This implies that dj ∈ Dj′ .

4

Constraint Propagation: First Approach

We now move on to a discussion of constraint propagation for the arithmetic
constraints on integer intervals. To illustrate the first approach consider the
following example. Consider the constraint
x3 y − x ≤ 40
and the ranges x ∈ [1..100] and y ∈ [1..100]. We can rewrite it as
r

40 + x
x≤ 3
y

(3)

since x assumes integer
 √values.
 The maximum value the expression on the righthand side can take is 3 140 , so we conclude x ≤ 5. By reusing (3), now with the
information that x ∈ [1..5], we conclude that the maximum
 √  value the expression
on the right-hand side of (3) can take is actually 3 45 , from which it follows
that x ≤ 3.
In the case of y we can isolate it by rewriting the original constraint as
y ≤ x403 + x12 from which it follows that y ≤ 41, since by assumption x ≥ 1. So
we could reduce the domain of x to [1..3] and the domain of y to [1..41]. This
interval reduction is optimal, since x = 1, y = 41 and x = 3, y = 1 are both
solutions to the original constraint x3 y − x ≤ 40.
m
More formally, we consider a polynomial constraint Σi=1
mi = b where m > 0,
no monomial mi is an integer, the power products of the monomials are pairwise
different and b is an integer. Suppose that x1 , . . ., xn are its variables ordered
w.r.t. ≺.
Select a non-integer monomial mℓ and assume it is of the form
a · y1n1 · . . . · yknk ,
where k > 0, y1 , . . ., yk are different variables ordered w.r.t. ≺, a is a non-zero
integer and n1 , . . ., nk are positive integers. So each yi variable equals to some
variable in {x1 , . . ., xn }. Suppose that yp equals to xj . We introduce the following
proof rule:
POLYNOMIAL EQUALITY
n
hΣi=1
mi
n
hΣi=1 mi

= b ; x1 ∈ D1 , . . ., xn ∈ Dn i
= b ; x1 ∈ D1′ , . . ., xn ∈ Dn′ i

where
– for i 6= j

Di′ := Di ,

–
Dj′
and



q

np
:= int Dj ∩
int (b − Σi∈[1..m]−{ℓ}mi )/s
n

n

p−1
p+1
s := a · y1n1 · . . . · yp−1
· yp+1
. . . · yknk .

To see that this rule preserves equivalence choose some d1 ∈ D1 , . . ., dn ∈ Dn .
To simplify the notation, given an extended arithmetic expression t denote by t′
the result of evaluating t after each occurrence of a variable xi is replaced by di .
m
Suppose that Σi=1
m′i = b. Then
n

dj p · s′ = b − Σi∈[1..m]−{ℓ} m′i ,
so by the Correctness Lemma 1 applied to b − Σi∈[1..m]−{ℓ} m′i and to s
n

dj p ∈ int(b − Σi∈[1..m]−{ℓ} mi )/int(s).
Hence
dj ∈

q
int(b − Σi∈[1..m]−{ℓ} mi )/int(s)

np

and consequently


q

np
int (b − Σi∈[1..m]−{ℓ} mi )/s
dj ∈ int Dj ∩
i.e., dj ∈ Dj′ .
Note that we do not apply int(.) to the outcome of the root extraction operation. For even np this means that the second operand of the intersection can
be a union of two intervals, instead of a single interval. To see why this is desirable, consider the constraint x2 − y = 0 in the presence of ranges x ∈ [0..10],
y ∈ [25..100]. Using the int(.) closure of the root extraction we would not be
able to update the lower bound of x to 5.
m
Next, consider a polynomial constraint Σi=1
mi ≤ b. Below we adopt the
assumptions and notation used when defining the POLYNOMIAL EQUALITY
rule. To formulate the appropriate rule we stipulate that for extended arithmetic
expressions s and t
int((≤ s)/t) := ≥ Q ∩ ≤ Q,

with Q = (≤ int(s))/int(t).
To reason about this constraint we use the following rule:
POLYNOMIAL INEQUALITY
n
hΣi=1
mi ≤ b ; x1 ∈ D1 , . . ., xn ∈ Dn i
n m ≤ b ; x ∈ D ′ , . . ., x ∈ D ′ i
hΣi=1
i
1
n
n
1

where
– for i 6= j

Di′ := Di ,

–
Dj′


q
np
int
:= int Dj ∩

≤ (b

− Σi∈[1..m]−{ℓ}mi )/s





To prove that this rule preserves equivalence choose some d1 ∈ D1 , . . ., dn ∈
Dn . As above given an extended arithmetic expression t we denote by t′ the
result of evaluating t when each occurrence of a variable xi in t is replaced by
di .
m
Suppose that Σi=1
m′i ≤ b. Then
n

dj p · s′ ≤ b − Σi∈[1..m]−{ℓ} m′i .
By the Correctness Lemma 1
b − Σi∈[1..m]−{ℓ} m′i ∈ int(b − Σi∈[1..m]−{ℓ} mi ),
so by definition
n

dj p · s′ ∈≤ int(b − Σi∈[1..m]−{ℓ} mi ).
Hence by the definition of the division operation on the sets of integers
n

dj p ∈≤ int(b − Σi∈[1..m]−{ℓ} mi )/int(s)
Consequently
dj ∈

q

np

≤ int(b −

Σi∈[1..m]−{ℓ} mi )/int(s)

This implies that dj ∈ Dj′ .
Note that the set ≤ int(b − Σi∈[1..m]−{ℓ}mi ) is not an interval. So to properly
implement this rule we need to extend the implementation of the division operation discussed in Subsection 2.2 to the case when the numerator is an extended
interval. Our implementation takes care of this.
In an optimized version of this approach we simplify the fractions of two polynomials by splitting the division over addition and subtraction and by dividing
out common powers of variables and greatest common divisors of the constant
factors. Subsequently, fractions whose denominators have identical power products are added. We used this optimization in the initial example by simplifying
40
1
40+x
x3 to x3 + x2 . The reader may check that without this simplification step we
can only deduce that y ≤ 43.
To provide details of this optimization, given two monomials s and t, we
denote by
s
[ ]
t
the result of performing this simplification operation on s and t. For example,
3
4·x3 ·y
x·y
·y
2·x3
[ 2·x
4·x2 ] equals 2 , whereas [ 2·y 2 ] equals y .
In this approach we assume that the domains of the variables y1 , . . ., yp−1 ,
yp+1 , . . ., yn of mℓ do not contain 0. (One can easily show that this restriction is
necessary here). For a monomial s involving variables ranging over the integer
intervals that do not contain 0, the set int(s) either contains only positive numbers or only negative numbers. In the first case we write sign(s) = + and in the
second case we write sign(s) = −.

The new domain of the variable xj in the POLYNOMIAL INEQUALITY
rule is defined using two sequences m′0 ...m′n and s′0 ...s′n of extended arithmetic
expressions such that
m′0
b
mi
m′
= [ ] and ′i = −[ ] for i ∈ [1..m].
′
s0
s
si
s
Let S := {s′i | i ∈ [0..m] − {ℓ}} and for an extended arithmetic expression t ∈ S
let It := {i ∈ [0..m] − {ℓ} | s′i = t}. We denote then by pt the polynomial
P
′
i∈It mi . The new domains are then defined by
Dj′

r



pt 
np
≤ int Σ
:= int Dj ∩
t∈S
t

if sign(s) = +, and by
Dj′

r



pt 
np
≥ int Σ
:= int Dj ∩
t∈S
t

if sign(s) = −. Here the int(s) notation used in the Correctness Lemma 1 is
extended to expressions involving the division operator on real intervals in the
obvious way. We define the int(.) operator applied to a bounded set of real
numbers, as produced by the division and addition operators in the above two
expressions for Dj′ , to denote the smallest interval of real numbers containing
that set.

5

Constraint Propagation: Second Approach

In this approach we limit our attention to a special type of polynomial constraints, namely the ones of the form s op b, where s is a polynomial in which
each variable occurs at most once and where b is an integer. We call such a constraint a simple polynomial constraint. By introducing auxiliary variables that
are equated with appropriate monomials we can rewrite each polynomial constraint into a sequence of simple polynomial constraints. This allows us also to
compute the integer interval domains of the auxiliary variable from the integer
interval domains of the original variables. We apply then to the simple polynomial constraints the rules introduced in the previous section.
To see that the restriction to simple polynomial constraints can make a difference consider the constraint
100x · y − 10y · z = 212
in presence of the ranges x, y, z ∈ [1..9]. We rewrite it into the sequence
u = x · y, v = y · z, 100u − 10v = 212
where u, v are auxiliary variables, each with the domain [1..81].

It is easy to check that the POLYNOMIAL EQUALITY rule introduced
in the previous section does not yield any domain reduction when applied to
the original constraint 100x · y − 10y · z = 212. In presence of the discussed
optimization the domain of x gets reduced to [1..3].
However, if we repeatedly apply the POLYNOMIAL EQUALITY rule to
the simple polynomial constraint 100u − 10v = 212, we eventually reduce the
domain of u to the empty set (since this constraint has no integer solution in the
ranges u, v ∈ [1..81]) and consequently can conclude that the original constraint
100x · y − 10y · z = 212 has no solution in the ranges x, y, z ∈ [1..9], without
performing any search.

6

Constraint Propagation: Third Approach

In this approach we focus on a small set of ‘atomic’ arithmetic constraints. We
call an arithmetic constraint atomic if it is in one of the following two forms:
– a linear constraint,
– x · y = z.
It is easy to see that using appropriate transformation rules involving auxiliary variables we can transform each arithmetic constraint to a sequence of
atomic arithmetic constraints. In this transformation, as in the second approach,
the auxiliary variables are equated with monomials so we can easily compute
their domains.
The transformation to atomic constraints can strengthen the reduction. Consider for example the constraint
u·x·y+1= v·x·y
and ranges u ∈ [1..2], v ∈ [3..4], and x, y ∈ [1..4]. The first approach without
optimization and the second approach cannot find a solution without search.
If, as a first step in transforming this constraint into a linear constraint, we
introduce an auxiliary variable w to replace x · y, we are effectively solving the
constraint
u·w+1= v·w
with the additional range w ∈ [1..16], resulting in only one duplicate occurrence
of a variable instead of two. With variable w introduced (or using the optimized
version of the first approach) constraint propagation alone finds the solution
u = 2, v = 3, x = 1, y = 1.
We explained already in Section 3 how to reason about linear constraints.
(We omitted there the treatment of the disequalities which is routine.) Next, we
focus on the reasoning for the multiplication constraint x · y = z in presence of
the non-empty ranges x ∈ Dx , y ∈ Dy and z ∈ Dz . To this end we introduce the
following three domain reduction rules:

MULTIPLICATION 1
hx · y = z ; x ∈ Dx , y ∈ Dy , z ∈ Dz i
hx · y = z ; x ∈ Dx , y ∈ Dy , z ∈ Dz ∩ int(Dx · Dy )i
MULTIPLICATION 2
hx · y = z ; x ∈ Dx , y ∈ Dy , z ∈ Dz i
hx · y = z ; x ∈ Dx ∩ int(Dz /Dy ), y ∈ Dy , z ∈ Dz i
MULTIPLICATION 3
hx · y = z ; x ∈ Dx , y ∈ Dy , z ∈ Dz i
hx · y = z ; x ∈ Dx , y ∈ Dy ∩ int(Dz /Dx ), z ∈ Dz i
The way we defined the multiplication and the division of the integer intervals
ensures that the MULTIPLICATION rules 1,2 and 3 are equivalence preserving.
Consider for example the MULTIPLICATION 2 rule. Take some a ∈ Dx , b ∈ Dy
and c ∈ Dz such that a · b = c. Then a ∈ {x ∈ Z | ∃z ∈ Dz ∃y ∈ Dy x · y = z}, so
a ∈ Dz /Dy and a fortiori a ∈ int(Dz /Dy ). Consequently a ∈ Dx ∩ int(Dz /Dy ).
This shows that the MULTIPLICATION 2 rule is equivalence preserving.
The following example shows an interaction between all three MULTIPLICATION rules.
Example 1. Consider the CSP
hx · y = z ; x ∈ [1..20], y ∈ [9..11], z ∈ [155..161]i.

(4)

To facilitate the reading we underline the modified domains. An application
of the MULTIPLICATION 2 rule yields
hx · y = z ; x ∈ [16..16], y ∈ [9..11], z ∈ [155..161]i
since, as already noted in in Subsection 2.2, [155..161]/[9..11]) = [16..16], and
[1..20] ∩ int([16..16]) = [16..16]. Applying now the MULTIPLICATION 3 rule
we obtain
hx · y = z ; x ∈ [16..16], y ∈ [10..10], z ∈ [155..161]i
since [155..161]/[16..16] = [10..10] and [9..11] ∩ int([10..10]) = [10..10]. Next, by
the application of the MULTIPLICATION 1 rule we obtain
hx · y = z ; x ∈ [16..16], y ∈ [10..10], z ∈ [160..160]i
since [16..16] · [10..10] = [160..160] and [155..161] ∩ int([160..160]) = [160..160].
So using all three multiplication rules we could solve the CSP (4).
✷

Now let us clarify why we did not define the division of the sets of integers
Z and Y by
Z/Y := {z/y ∈ Z | y ∈ Y, z ∈ Z, y 6= 0}.

The reason is that in that case for any set of integers Z we would have Z/{0} = ∅.
Consequently, if we adopted this definition of the division of the integer intervals, the resulting MULTIPLICATION 2 and 3 rules would not be anymore
equivalence preserving. Indeed, consider the CSP
hx · y = z ; x ∈ [−2..1], y ∈ [0..0], z ∈ [−8..10]i.
Then we would have [−8..10]/[0..0] = ∅ and consequently by the MULTIPLICATION 2 rule we could conclude
hx · y = z ; x ∈ ∅, y ∈ [0..0], z ∈ [−8..10]i.
So we reached an inconsistent CSP while the original CSP is consistent.
In the remainder of the paper we will also consider variants of this third
approach that allow squaring and exponentiation as atomic constraints. For this
purpose we explain the reasoning for the constraint x = y n in presence of the
non-empty ranges x ∈ Dx and y ∈ Dy , and for n > 1. To this end we introduce
the following two rules in which to maintain the property that the domains are
intervals we use the int(.) operation of Section 2:
EXPONENTIATION
hx = y n ; x ∈ Dx , y ∈ Dy i
hx = y n ; x ∈ Dx ∩ int(Dyn ), y ∈ Dy i
ROOT EXTRACTION
hx = y n ; x ∈ Dx , y ∈ Dy i
√
hx = y n ; x ∈ Dx , y ∈ int(Dy ∩ n Dx )i

To prove that these rules are equivalence preserving suppose that for some
a ∈ Dx and b ∈ Dy we have a = bn . Then a ∈ Dyn , so a ∈ int(Dyn ) and conse√
√
quently a ∈ Dx ∩ int(Dyn ). Also b ∈ n Dx , so b ∈ Dy ∩ n Dx , and consequently
√
b ∈ int(Dy ∩ n Dx ).

7

Implementation Details

In this section we describe the benchmark experiments that were performed
to compare the proposed approaches. These experiments were performed using
a single solver of the DICE (DIstributed Constraint Environment) framework.
DICE [10] is a framework for solver cooperation, implemented using techniques
from coordination programming. It is developed around an experimental constraint solver, called OpenSolver, which is particularly suited for coordination.
The coordination and cooperation aspects are irrelevant from the point of view
of this paper. Relevant aspects of the OpenSolver are:

– It implements a branch-and-infer tree search algorithm for constraint solving. The inference stage corresponds to constraint propagation and is performed by repeated application of domain reduction functions (DRFs) that
correspond to the domain reduction rules associated with the considered
constraints.
– This algorithm is abstract in the sense that the actual functionality is determined by software plug-ins in a number of predefined categories. These categories correspond to various aspects of the abstract branch-and-infer tree
search algorithm. Relevant categories are: variable domain types, domain
reduction functions, schedulers that control the application of the DRFs,
branching strategies that split the search tree after constraint propagation
has terminated, and several categories corresponding to different aspects of
a search strategy that determine how to traverse a search tree.
All experiments were performed using the IntegerInterval variable domain type plug-in. Domains of this type consist of an indication of the type of
the interval (bounded, unbounded, left/right-bounded, or empty), and a pair
of arbitrary precision integer bounds. This plug-in, and the interval arithmetic
operations on it are built using the GNU MP library [4].
The branching strategy that we used selects variables using the chronological
ordering in which the auxiliary variables come last. The domain of the selected
variable is split into two subdomains using bisection, so the resulting search trees
are binary trees. In all experiments we searched for all solutions, traversing the
entire search tree by means of depth-first leftmost-first chronological backtracking.
For the experiments in this paper a DRF plug-in has been developed that
implements the domain reduction rules discussed in the previous sections. The
scheduler plug-in used in the benchmarks keeps cycling through the sequence of
DRFs, applying DRFs that have been scheduled for execution. When a DRF is
applied, and some variable domain is modified, all DRFs that depend on these
changes are scheduled for execution, including possibly the one that has just been
applied. The cycling stops when no more DRFs are scheduled for execution, or
when the domain of a variable becomes empty.
As an alternative to cycling, the scheduler can be supplied with a schedule:
a sequence of indices into the sequence of DRFs. The scheduler will then cycle
through this schedule instead, and consider DRFs for application in the specified
order. This is used in combination with the second and third approach, where we
distinguish user constraints from the constraints that are introduced to define
the values of auxiliary variables. Before considering for execution a DRF f that
is part of the implementation of a user constraint, we make sure that all auxiliary
variables that f relies on are updated. For this purpose, the indices of the DRFs
that update these variables precede the index of f in the schedule. If f can
change the value of an auxiliary variable, its index is followed by the indices of
the DRFs that propagate back these changes to the variables that define the
value of this auxiliary variable.

For the third approach, there can be hierarchical dependencies between auxiliary variables. Much like the HC4 algorithm of [2], the schedule specifies a
bottom-up traversal of this hierarchy in a forward evaluation phase and a topdown traversal in a backward propagation phase before and after applying a
DRF of a user constraint, respectively. In the forward evaluation phase, the
DRFs that are executed correspond to rules MULTIPLICATION 1 and EXPONENTIATION. The DRFs of the backward propagation phase correspond to
MULTIPLICATION 2 and 3 , and ROOT EXTRACTION. It is easy to construct examples showing that the use of hierarchical schedules can be beneficial
compared to cycling through the rules.
The proposed approaches were implemented by first rewriting arithmetic
constraints to polynomial constraints, and then to a sequence of DRFs that
correspond with the rules of the approach used. We considered the following
methods:
1a the first approach, discussed in Section 4,
1b the optimization of the first approach discussed at the end of Section 4 that
involves dividing out common powers of variables,
2a the second approach, discussed in Section 5. The conversion to simple polynomial constraints is implemented by introducing an auxiliary variable for
every non-linear monomial. This procedure may introduce more auxiliary
variables than necessary.
2b an optimized version of approach 2a, where we stop introducing auxiliary
variables as soon as the constraints contain no more duplicate occurrences
of variables.
3a the third approach, discussed in Section 6, allowing only linear constraints
and multiplication as atomic constraints.
3b idem, but also allowing x = y 2 as an atomic constraint.
3c idem, allowing x = y n for all n > 1 as an atomic constraint.
Approaches 2 and 3 involve an extra rewrite step, where the auxiliary variables are introduced. The resulting CSP is then rewritten according to approach
1a. During the first rewrite step the hierarchical relations between the auxiliary
variables are recorded and the schedules are generated as a part of the second
rewrite step. For approaches 2b and 3 the question of which auxiliary variables
to introduce is an optimization problem in itself. Some choices result in more
auxiliary variables than others. We have not treated this issue as an optimization problem but relied on heuristics. We are confident that these yield a realistic
implementation. In our experiments we used the following benchmarks.
Cubes The problem is to find all natural numbers n ≤ 1000 that are a sum of
four different cubes, for example
13 + 23 + 33 + 43 = 100.
This problem is modeled as follows:
h1 ≤ x1 , x1 ≤ x2 − 1, x2 ≤ x3 − 1, x3 ≤ x4 − 1, x4 ≤ n,
x31 + x32 + x33 + x34 = n; n ∈ [1..1000], x1 , x2 , x3 , x4 ∈ Zi

Opt We are interested in finding a solution to the constraint x3 + y 2 = z 3 in the
integer interval [1..100000] for which the value of 2x · y − z is maximal.
Fractions This problem is taken from [9]: find distinct nonzero digits such that
the following equation holds:
A
D
G
+
+
=1
BC
EF
HI
There is a variable for each letter. The initial domains are [1..9]. To avoid symmetric solutions an ordering is imposed:
D
G
A
≥
≥
BC
EF
HI
Also two redundant constraints are added:
3

A
≥1
BC

and

3

G
≤1
HI

Because division is not present in the arithmetic expressions, the above constraints are multiplied by the denominators of the fractions to obtain arithmetic
constraints.
Two representations for this problem were studied:
– fractions1 in which five constraints are used: one equality and four inequalities for the ordering and the redundant constraints,
– fractions2, used in [9], in which three auxiliary variables, BC, EF and HI,
are introduced to simplify the arithmetic constraints: BC = 10B + C, EF =
10E + F , and HI = 10H + I.
Additionally, in both representations, 36 disequalities A 6= B, A 6= C, ..., H 6= I
are used.
Kyoto The problem1 is to find the number n such that the alphanumeric equation
KYOTO
KYOTO
+KYOTO
TOKYO
has a solution in the base-n number system. Our model uses a variable for each
letter and one variable for the base number. The variables K and T may not
be zero. There is one large constraint for the addition, 6 disequalities K 6= Y
... T 6= O and four constraints stating that the individual digits K, Y, O, T ,
are smaller than the base number. To spend some CPU time, we searched base
numbers 2..100.
1

V. Dubrovsky and A. Shvetsov. Quantum
http://www.nsta.org/quantum/kyotoarc.asp

cyberteaser:

May/June

1995,

8

Results

Tables 1 and 2 compare the proposed approaches on the problems defined in the
previous section. The first two columns of table 1 list the number of variables and
the DRFs that were used. Column nodes lists the size of the search tree, including
failures and solutions. The next two columns list the number of times that a
DRF was executed, and the percentage of these activations that the domain of a
variable was actually modified. For the opt problem, the DRF that implements
the optimization is not counted, and its activation is not taken into account. The
elapsed times in the last column are the minimum times (in seconds) recorded
for 5 runs on a 1200 MHz Athlon CPU.

cubes

1a
2a
3a
3b
opt
1a
2a
3a
3b
fractions1 1a
1b
2a
2b
3
fractions2 1a
1b
2a
2b
3
kyoto
1a
1b
2a
2b
3a
3b
3c

nvar nDRF
nodes
activated %effective elapsed
5
14
167
1880
13.03
0.013
9
22
167
2370
22.15
0.014
13
34
359
4442
26.23
0.024
13
34
227
3759
29.24
0.021
4
7
115,469
5,186,968
42.16 22.037
8
15
115,469
9,799,967
60.00 23.544
10
21
?
?
?
?
10
21 5,065,137 156,903,869
46.49 518.898
9
154
11,289
1,193,579
3.65 16.586
9
154
7,879
734,980
3.45 17.811
37
210
11,289
1,410,436
23.27
5.575
32
200
11,289
1,385,933
21.65
5.957
43
208
11,131
1,426,186
27.76
5.635
12
105
2,449
270,833
9.72
0.660
12
105
989
94,894
9.12
0.538
20
121
2,449
350,380
22.19
0.597
15
111
2,449
301,855
17.51
0.547
22
123
1,525
293,038
27.33
0.509
5
37
87,085
3,299,736
6.09 23.680
5
37
87,085
3,288,461
5.94 45.406
13
53
87,085
3,781,414
23.03 11.784
12
51
87,085
3,622,361
21.45 12.138
16
60
87,087
4,275,930
26.70 22.538
16
60
87,085
4,275,821
26.70 22.530
16
59
87,085
3,746,532
23.26 10.466

E
+
+
=
+
+
=
=
=
=
=
=
=
=
=
=
=
+
=
=
=
=
=

I
=
=
+
+
=
=
=
=
=
=
=
=
=
=
=
+
=
=
=
=
=

Table 1. Statistics and comparison with other solvers

Table 2 lists measured numbers of basic interval operations. Note that for
approach 1b, there are two versions of the division and addition operations:
one for integer intervals, and one for intervals of reals of which the bounds are
rational numbers (marked Q). Columns multI and multF list the numbers of

multiplications of two integer intervals, and of an integer interval and an integer
factor, respectively. These are different operations in our implementation.

cubes

1a
2a
3a
3b
opt
1a
2a
3a
3b
fractions1 1a
1b
2a
2b
3
fractions2 1a
1b

kyoto

2a
2b
3
1a
1b
2a
2b
3a
3b
3c

root
exp
div
multI multF
sum
total
1
4
0
0
5
4
14
< 0.5 < 0.5
0
0
5
4
9
0
0
1
1
6
5
13
< 0.5 < 0.5
1
< 0.5
5
5
11
2,299 4,599 1,443
1,444 11,064 5,187
26,037
1,636 1,538 2,150
738
8,138 4,445
18,645
?
?
?
?
?
?
?
21,066 18,105 54,171 18,284 106,651 57,469 275,747
0
0
868 28,916 14,238 13,444
57,466
0
0
51 11,892
8,010 6,727
29,584
1,550 Q
1,355 Q
0
0
734
933
4,736 4,669
11,071
0
0
776
1,509
5,292 5,147
12,725
0
0
693
339
4,835 4,769
10,636
0
0
142
690
304
212
1,347
0
0
19
127
59
26
344
65 Q
49 Q
0
0
124
149
138
94
505
0
0
124
206
210
118
658
0
0
114
46
142
101
403
735 11,040 1,963 13,852 10,852 13,946
52,388
735 8,146
218
8,955 12,516 10,592
48,749
4,310 Q
3,277 Q
383
759 1,590
484
5,324 7,504
16,044
383
759 1,597
1,360
5,756 8,008
17,863
0
0 1,991
578
5,324 7,505
15,397
< 0.5 < 0.5 1,990
578
5,324 7,504
15,397
1
1 1,554
484
5,324 7,504
14,868

Table 2. Measured numbers (thousands) of interval operations

For the cubes and opt problems, the constraints are already in simple form, so
approaches 1a, 1b and 2b are identical. Also all non-linear terms involve either a
multiplication or an exponentiation, so also approaches 2a and 3c are the same.
The results of these experiments clearly show the disadvantage of implementing
exponentiation by means of multiplication: the search space grows because we
increase the number of variable occurrences and lose the information that it is
the same number that is being multiplied. For opt and approach 3a, the run did
not complete within reasonable time and was aborted.

Columns E and I of table 1 compare the propagation achieved by our approaches with two other systems, respectively ECLi PSe Version 5.62 using the
ic library, and ILOG Solver 5.13 using type ILOINT. For this purpose we ran
the test problems without search, and compared the results of constraint propagation. A mark ‘=’ means that the computed domains are the same, ‘+’ that
our approach achieved stronger propagation than the solver that we compare
with, and ‘-’ that propagation is weaker. For cubes, ECLi PSe computes the same
domains as those computed according to approach 3b, so here the reduction is
stronger than for 3a, but weaker than for the other approaches. For opt ECLi PSe
and ILOG Solver compute the same domains. These domains are narrower than
those computed according to approaches 3a and 3b, but the other approaches
achieve stronger reduction. In all other cases except for kyoto and approach 1b
the results of all three solvers are the same.
For both representations for the fractions puzzle, the symbolic manipulation
of approach 1b is able to achieve a significant reduction of the search tree, but
this is not reflected in the timings. For fractions1 the elapsed time even increases.
The reason is that computing the domain updates involves adding intervals of
real numbers. The arithmetic operations on such intervals are more expensive
than their counterparts on integer intervals, because the bounds have to be maintained as rational numbers. Arithmetic operations on rational numbers are more
expensive because they involve the computation of greatest common divisors.
For kyoto the symbolic manipulation did not reduce the size of the search tree,
so the effect is even more severe.
In general, the introduction of auxiliary variables leads to a reduction of
the number of interval operations compared to approach 1a. The reason is that
auxiliary variables prevent the evaluation of subexpressions that did not change.
This effect is strongest for fractions1, where the main constraint contains a
large number of different power products. Without auxiliary variables all power
products are evaluated for every POLYNOMIAL EQUALITY rule defined by
this constraint, even those power products the variable domains of which did
not change. With auxiliary variables the intervals for such unmodified terms are
available immediately, which leads to a significant reduction of the number of
interval multiplications.
The effect that stronger reduction is achieved as a result of introducing auxiliary variables, mentioned in Section 6, is seen for both representations of the
fractions benchmark. The effect described in Section 5 is not demonstrated by
these experiments.
If we don’t consider the symbolic manipulation of approach 1b, approach
3c leads to the smallest total number of interval operations in all cases, but
the scheduling mechanism discussed in Section 7 is essential for a consistent
good performance. If for example the schedule is omitted for opt, the number
2

3

ECLi PSe
Constraint
Logic
icparc.doc.ic.ac.uk/eclipse
See http://www.ilog.com

Programming

System.

See

http://www-

of interval operations almost triples, and performance of approach 2a and 3c is
then much worse than that of approach 1a.
The total numbers of interval operations in table 2 do not fully explain all
differences in elapsed times. One of the reasons is that different interval operations have different costs. Especially the preprocessing of the numerator interval
for integer interval division, discussed in Subsection 2.2, is potentially expensive, which may explain why for opt , approach 1a runs faster than approach 2a,
even though the total number of interval operations is higher. Among the many
other factors that may be of influence, some overhead is involved in applying a
DRF, so if the number of applications differs significantly for two experiments,
this probably influences the elapsed times as well (cubes, 1a, 2a, opt , 1a, 2a,
fractions2 , 2a, 2b). The elapsed times are not the only measure that is subject
to implementation details. For example, we implemented division by a constant
interval [−1.. − 1] as multiplication by a constant, which is more efficient in our
implementation. Such decisions are reflected in the numbers reported in table 2.

9

Discussion

In this paper we discussed a number of approaches to constraint propagation
for arithmetic constraints on integer intervals. To assess them we implemented
them using the DICE (DIstributed Constraint Environment) framework of [10],
and compared their performance on a number of benchmark problems. We can
conclude that:
– Implementation of exponentiation by multiplication gives weak reduction.
In our third approach x = y n should be an atomic constraint.
– The optimization of the first approach, where common powers of variables
are divided out, can significantly reduce the size of the search tree, but
the resulting reduction steps rely heavily on the division and addition of
rational numbers. These operations can be expected to be more expensive
than their integer counterparts, because they involve the computation of
greatest common divisors.
– Introducing auxiliary variables can be beneficial in two ways: it may strengthen the propagation, as discussed in Sections 5 and 6, and it may prevent the
evaluation of subexpressions the variable domains of which did not change.
– As a result, given a proper scheduling of the rules, the second and third
approach perform better than the first approach without the optimization,
in terms of numbers of interval operations. Actual performance depends on
many implementation aspects. However for our test problems the results of
variants 2a, 2b and 3c do not differ significantly.
In general, our implementation is slow compared to, for example, ILOG
Solver. A likely cause is that we use arbitrary precision integers. We chose this
representation to avoid having to deal with overflow, but an additional benefit
is that large numbers can be represented exactly.

A different approach would be to use floating-point arithmetic and then round
intervals inwards to the largest enclosed integer interval. This was suggested in [3]
and implemented in for example RealPaver4. A benefit of this inward rounding
approach is that all algorithms that were developed for constraints on the reals
are immediately available. A disadvantage is that for large numbers no precise
representation exists, i.e., the interval defined by two consecutive floating-point
numbers contains more than one integer. But it is debatable whether an exact
representation is required for such large numbers.
We realize that the current set of test problems is rather limited. In addition
to puzzles, some more complex non-linear integer optimization problems should
be studied. We plan to further evaluate the proposed approaches on non-linear
integer models for the SAT problem. Also we would like to study the relationship
with the local consistency notions that have been defined for constraints on the
reals and give a proper characterization of the local consistencies computed by
our reduction rules.
Note This work was performed during the first author’s stay at the School of
Computing of the National University of Singapore. The work of the second
author was supported by NWO, The Netherlands Organization for Scientific
Research, under project number 612.069.003.

References
1. K. R. Apt. A proof theoretic view of constraint programming. Fundamenta Informaticae, 33(3):263–293, 1998. Available via http://arXiv.org/archive/cs/.
2. F. Benhamou, F. Goualard, L. Granvilliers, and J.-F. Puget. Revising hull and
box consistency. In Proceedings of the 16th International Conference on Logic
Programming (ICLP’99), pages 230–244. The MIT Press, 1999.
3. F. Benhamou and W. Older. Applying interval arithmetic to real, integer and
boolean constraints. Journal of Logic Programming, 32(1):1–24.
4. T. Granlund. GNU MP, The GNU Multiple Precision Arithmetic Library, Edition
4.1.2. Swox AB, December 2002.
5. W. Harvey and P. J. Stuckey. Improving linear constraint propagation by changing
constraint representation. Constraints, 8(2):173–207, 2003.
6. T. J. Hickey, Q. Ju, and M. H. van Emden. Interval arithmetic: from principles to
implementation. Journal of the ACM, 48(5):1038–1068, 2001.
7. R. E. Moore. Interval Analysis. Prentice-Hall, Englewood Cliffs, NJ, 1966.
8. D. Ratz. Inclusion isotone extended interval arithmetic. Technical report, University of Karlsruhe, 1996. Report No. D-76128 Karlsruhe.
9. C. Schulte and G. Smolka. Finite domain constraint programming in Oz. A tutorial, August 2002. Version 1.2.4 (20020829). Available via http://www.mozartoz.org/documentation/fdt/index.html.
10. P. Zoeteweij. Coordination-based distributed constraint solving in DICE. In Brajendra Panda, editor, Proceedings of the 2003 ACM Symposium on Applied Computing, pages 360–366, 2003.
4

http://www.sciences.univ-nantes.fr/info/perso/permanents/granvil/realpaver/main.html

