Implementation

A system of $m$ distance functions has a specific description. First consider distance function that involve only points of constant. Such function may be written as:

$$\sum (x_i-x_j)^2+L_i =0$$

where the maximal number of square term is $n$ and is described by a row of 2 matrices APOW, ACONS and a vector LI. APOW is an integer matrix with $2n$ columns and $m$ rows that contain the unknown number of each term in the function and in which a 0 means that the unknown is a constant. The value of theses constants are given in the matrix of real ACONS of size $m \times 2n$. Finally the value of the constant $L_i$ is indicated in the vector LI of size $m$. For example consider a system involving the 4 unknowns $x_1, x_2, y_1, y_2$ numbered from 1 to 4:

$$\begin{eqnarray*} &&x_1^2+(y_1-4)^2-6 =0\\ &&(x_1-x_2)^2+(y_1-y_2)^2-12=0 \end{eqnarray*}$$

In that case APOW, ACONS, LI will be:

$${\tt APOW}= \left( \begin{tabular}{cccc} 1& 0& 3& 0\\ 1& 2 ... ...I}= \left( \begin{tabular}{c} -6\\ -12\\ \end{tabular} \right)$$

The two first elements of the first row of APOW (1, 0) describes the first square term of the first equation and state that it is $(x_1-a)^2$, the value of $a$ being given in the ACONS(1,2). Note that each square term must be written as (unknown-unknown or constant$)^2$ and not as (constant -unknown$)^2$.

Consider now function involving virtual points. Each square term may thus be written as:

$$(\sum \lambda_j x_j -X_k)^2$$

where $X_k$ may be either a constant, an unknown or the coordinate of a virtual point. Let $k$ be the number of term of the form $\sum \lambda_j x_j$ existing in the system. Such equation is described by a matrix AVARV with $k$ rows and a number of columns equal to the number of unknowns. Each term $\sum \lambda_j x_j$ is numbered from 1 to k and the row $i$ of AVARV will contain in position $j$ the value of $\lambda_j$. The existence of the coordinate of a virtual point in a distance function will be indicated in APOW by a negative number whose opposite is the number of the virtual coordinates. Hence if we add to the previous system the third equation:

$$(0.1 x_1-0.2 x_2-3)^2+(0.1 y_1-0.2y_2-1)^2-10 =0$$

the third row of APOW will be (-1 0 -2 0), the third row of ACONS will be (0 3 0 1) and the third row of LI will be -10. AVARV will have 2 row and 4 columns and is given as:

$${\tt AVARV}= \left( \begin{tabular}{cccc} 0.1 & -0.2 & 0 & 0\\ 0 & 0 & 0.1 & -0.2\\ \end{tabular} \right)$$

Among the $m$ equations there will be $p$ equations involving virtual points. The system must be written in such way that first are defined the $m-p$ equations not involving virtual points and then the $p$ equations. Note also that in the current implementation inequalities are handled although with less efficiency than equations. The algorithm is implemented as:

 
int Solve_Distance(int DimVar,int DimEq,
                   INTEGER_VECTOR &Type_Eq,
                   INTEGER_MATRIX &APOW,MATRIX &ACONS,VECTOR &LI,
                   int p,int k,MATRIX &AVARV,
                   INTERVAL_VECTOR & TheDomain, 
                   int M,
                   double epsilonf,
                   int Stop,
                   INTERVAL_MATRIX & Solution,int Nb,
                   int (* Simp_Proc)(INTERVAL_VECTOR &))

The arguments are:

• DimVar: number of unknowns
• DimEq: number of equations and inequalities in the system
• Type_Eq: an integer array of dimension DimEq where the j-th element indicates if the j-th function is an equation (value =0), an inequality $\le 0$ (value = -1) or an inequality $\ge 0$ (value = 1)
• p: number of functions involving virtual points
• k: number of virtual points
• TheDomain: box in which we are looking for solution of the system
• M: the maximum number of boxes which may be stored. In the algorithm we use the reverse storage mode except if the global integer ALIAS_Parallel_Slave is set to 1 (its default value is 0)
• epsilonf: either the maximal width of the function intervals for a solution if it is not determined by using Newton scheme or the accuracy used in the Newton scheme
• Stop: the possible values are 0,1,2
  • 0: the algorithm will look for every solution in TheDomain
  • 1: the algorithm will stop as soon as 1 solution has been found
  • 2: the algorithm will stop as soon as Nb solutions have been found
• Solution: an interval matrix of size (Nb,m) which will contained the solution intervals. This list is sorted using the order specified by Order
• Nb: the maximal number of solution which will be returned by the algorithm
• Simp_Proc: a user-supplied procedure that take as input the current box and proceed to some further reduction of the width of the box or even determine that there is no solution for this box, in which case it should return -1 (see note 2.3.3). Remember also that you may use the 3B method to improve the efficiency of this algorithm (see section 2.3.2).

The argument Simp_Proc in this procedure may be omitted.

The following variables play also a role in the computation:

• ALIAS_Diam_Max_Kraw: if the maximal width of the ranges in a box is lower than this value, then the Krawczyk operator will be used to determine if there is a unique solution in the box (default value: 1.e10)
• ALIAS_Diam_Max_Newton: if the maximal width of the ranges in a box is lower than this value, then the interval Newton method will be used either to try to reduce the width of the box or to to ensure that there is no solution of the system in the box (default value: 1.e10)
• ALIAS_Permute_List: if the value of this flag is $n$ the algorithm will permute the current list with the largest box in the list of boxes to process (as the algorithm uses systematically the Newton scheme with as initial guess the center of the current box permutation may allow to find quickly new solutions)