# An algorithm to generate the basis of solutions to homogeneous linear diophantine equations

Gérard Huet

► To cite this version:

Gérard Huet. An algorithm to generate the basis of solutions to homogeneous linear diophantine equations. [Research Report] IRIA-RR-274, IRIA. 1978, pp.12. hal-04716564

HAL Id: hal-04716564

https://inria.hal.science/hal-04716564v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dw-03678

RR274

RIA

laboria

laboratoire de recherche
en informatique
et automatique

ARCHIVES

AN ALGORITHM TO GENERATE
THE BASIS OF SOLUTIONS
TO HOMOGENEOUS LINEAR
DIOPHANTINE EQUATIONS

Gérard HUET

Rapport de Recherche N° 274

Janvier 1978

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20
# AN ALGORITHM TO GENERATE THE BASIS OF SOLUTIONS TO HOMOGENEOUS LINEAR DIOPHANTINE EQUATIONS

Gérard Huet

## Résumé :

On résout sur $\mathbb{N}$ l'équation $\sum_{i=1}^{m} a_i x_i = \sum_{j=1}^{n} b_j y_j$.

## Abstract :

We solve over $\mathbb{N}$ the equation $\sum_{i=1}^{m} a_i x_i = \sum_{j=1}^{n} b_j y_j$.
# 1. Introduction

Let $a_1, \dots, a_m, b_1, \dots, b_n$ be positive integers. We are interested in solving the equation

$$ \sum_{i=1}^{m} a_i x_i = \sum_{j=1}^{n} b_j y_j \quad (1) $$

over non-negative integers. The set of all such solutions forms a commutative monoid in $\mathbb{N}^{m+n}$. It may therefore be generated by a finite basis [1]. We say that a solution $(x,y)$ to (1) is *minimal* iff it is non-trivial (i.e. different from $(0, \dots, 0)$) and such that there is no distinct non-trivial solution $(x, y')$ to (1) such that $(x, y') \le (x,y)$ (where $\le$ is the component-wise order in $\mathbb{N}^{m+n}$). The basis consists in all the minimal solutions. We are interested in good algorithms to generate it.

This problem has applications in symbolic computation. In particular, it is a key component to unification algorithms over terms containing associative commutative function symbols [2,3].

# 2. Bounds on minimal solutions

$$ \text{Let } \max_a = \max\{a_i \mid 1 \le i \le m\}, \quad \max_b = \max\{b_j \mid 1 \le j \le n\}. $$

It is possible to bound the coordinates of minimal solutions to (1) as follows.
**Lemma 1**

Let $(x,y) \in \mathbb{N}^m \times \mathbb{N}^n$ be a minimal solution to (1). Then $\forall i \le m x_i \le \max b$ and $\forall j \le n y_j \le \max a$.

**Proof**

The problem being symmetric, we shall only prove the first half of the lemma. Let us assume that $(x,y)$ is a solution to (1) such that $\exists k \le m x_k > \max b$. We shall prove that $(x,y)$ is not minimal. We have

$$ \max b \cdot \sum_{j=1}^{n} y_j \ge \sum_{j=1}^{n} b_j y_j = \sum_{i=1}^{m} a_i x_i \ge a_k x_k > a_k \max b \quad \text{and} \\ \text{therefore} \qquad \sum_{j=1}^{n} y_j > a_k. \qquad (2) $$

Now consider $u = <u_1, \dots, u_n> \in \mathbb{N}^n$ such that

$$ \forall i \le n \quad u_i \le y_i \qquad (3) $$

$$ \text{and} \qquad 0 < \sum_{j=1}^{n} u_j \le a_k \qquad (4) $$

Because of (2), there exists an increasing chain of such $u$'s in $\mathbb{N}^n$ of length $a_k$. It follows that there exists such a vector $u$ verifying furthermore

$$ \sum_{j=1}^{n} b_j u_j \equiv 0 \pmod{a_k} \qquad (5). $$

For, otherwise, according to the pigeonhole principle, there would exist $v$ and $w$ in the chain such that

$$ \sum_{j=1}^{n} h_j v_j \equiv \sum_{j=1}^{n} b_j w_j \pmod{a_k} $$

and, assuming for instance $v<w$, the vector $u=w-v$ would fulfill conditions (3) (4) and (5), a contradiction.

Assume therefore that $u$ verifies (3) (4) and (5). By (4), we have

$$ \sum_{j=1}^{n} b_j u_j \le \max b \cdot \sum_{j=1}^{n} u_j \le \max b \cdot a_k < a_k x_k \\ \text{and thus, using (5), we get} $$
$$ \sum_{j=1}^{n} b_j u_j = a_k z_k \quad \text{with } 0 < z_k < x_k $$

which shows that $(x,y)$ is strictly greater than the non-trivial solution $(z,u)$, with $\forall i \neq k \ z_i = 0$. □

Now let $d_{ij} = \text{lcm}(a_i, b_j)/a_i$
and $e_{ij} = \text{lcm}(a_i, b_j)/b_j$, for $1 \le i \le m$, $1 \le j \le n$.

Let $S_{ij} = (x,y) \in \mathbb{N}^m \times \mathbb{N}^n$ such that $x_i = d_{ij}, y_j = e_{ij}$, all other coordinates 0. The $S_{ij}$'s are special solutions to (1), and furthermore it is immediate to show that they are minimal.

This gives us a second way of bounding minimal solutions to (1) :

**Lemma 2**

Let $(x,y) \in \mathbb{N}^m \times \mathbb{N}^n$ be a solution to (1) such that $\exists i \leq m \ \exists j \leq n \ x_i \geq d_{ij}$ and $y_j \geq e_{ij}$. Then $(x,y) \geq S_{ij}$.

### 3. The algorithm

We are now able to express our algorithm for the generation of the basis of (1) as a backtrack algorithm using as stopping conditions the bounds given by lemmas 1 and 2. The potential solutions are enumerated in lexicographic order.

More precisely, assume values $x_1, \dots, x_k \in \mathbb{N}$ have been chosen, for some $k \le m$. We define, for all $j \le n$ :

$$ E_j^k = \{ e_{ij} - 1 \mid \exists i \le k \ x_i \ge d_{ij} \}. $$

We compute a vector maxy of maximum values of y to consider as :

$$ \maxy_j^k(x_1, \dots, x_k) = \begin{cases} \min(E_j^k) & \text{if } E_j^k \neq \emptyset \\ \maxa & \text{if } F_j^k = \emptyset. \end{cases} \quad (j \le n) $$

(Remark that $\forall i \le m \ \forall j \le n \ e_{ij} \le a_i \le \maxa$).
**Theorem**: Let (x,y) be a minimal solution of (1) not in {$S_{ij}$}. Then $\forall k \le m$ we have :

$$a) \quad x_i \le \max_j b_i \quad \forall i \le k$$

$$b) \quad \sum_{i=1}^{k} a_i x_i \le \sum_{j=1}^{n} b_j \max_{j}^{k} (x_1, \dots, x_k)$$

and $\forall l \le n$ we have :

$$c) \quad y_j \le \max_{j}^{m} (x_1, \dots, x_m) \quad \forall j \le l$$

$$d) \quad \sum_{j=1}^{l} b_j y_j \le \sum_{i=1}^{m} a_i x_i$$

**Proof**:

a) follows from lemma 2, d) from the fact that (x,y) is a solution.

c) follows from lemmas 1 and 2. Finally

$$\sum_{i=1}^{k} a_i x_i \le \sum_{i=1}^{m} a_i x_i = \sum_{j=1}^{n} b_j y_j \quad \text{because (x,y) is solution}$$

$$\le \sum_{j=1}^{n} b_j \max_{j}^{m} (x_1, \dots, x_m) \quad \text{because of c)}$$

$$\le \sum_{j=1}^{n} b_j \max_{j}^{k} (x_1, \dots, x_m) \quad \text{because } E_j^k \subseteq E_j^m, \text{ proving b).} \quad \square$$

Our algorithm consists in generating potential solutions in increasing lexicographic order, starting with (and not including) the trivial solution $<0, \dots, 0>$. The values of the $x_i$'s are progressively bounded according to the conditions a) and b) above. When all $x_i$'s are chosen, the $y_j$'s are bounded according to the conditions c) and d).

When the algorithm generates a potential solution (x,y), it checks that 1) it is indeed a solution to equation 1)

2) it is not greater than any solution previously generated. Then it backtracks to generate further solutions. When it finally stops, all particular solutions $S_{ij}$'s are added. Remark that when a solution is generated and checked, it is indeed minimal, because any solution generated later on will either be
greater in the lexicographic order (and therefore not below it in the compo-
nent-wise order), or one of the S<sub>ij</sub>'s, which cannot be below it by construc-
tion. The S<sub>ij</sub>'s are themselves known to be minimal, which finishes the correct-
ness proof of our algorithm.

4. Actual computer implementation

The algorithm described above has been programmed in PASCAL. The
usual program optimization techniques have been applied to push in outermost
loops as much as possible of the computation of the tests a) b) c) d) above.
The last value y_n is obtained by a division rather than by incrementations.
The complete text of the program is given in Appendix A.

Remark that an important degree of freedom is left in the formula-
tion of the problem, which is symmetric in a and b, and independent of the
order of the a_i's and b_j's. A different formulation may improve significantly
the running time of the program, since the bounds depend on the partially cons-
tructed solutions. It is our (empirical) experience that the fastest executions
are obtained in the average when :

$$
- \max b \le \max a
$$

$$
- a_1 \geq a_2 \geq \cdots \geq a_m
$$

- b_1 \le b_2 \le \dots \le b_n.

Examples of executions of the program, on an IRIS 80 computer, are
given in Appendix B.

Acknowledgements

Thanks to Gilles Kahn for his helpful comments and criticisms.
Appendix A : PASCAL program

program BASIS;
%$T-%
type TABLE[1..10]of INTEGER;
var INIMAX:TABLE;
    COEF,SOL:array[1..20]of INTEGER;
    D,E:array[1..10,1..10]of INTEGER;
    SOLS:array[0..10000]of INTEGER;
    TIM1,TIM2,LASTSOL,MAXA,MAXB,INIMAXSUM,I,J,CDIM,L,M,N,DIM:INTEGER;

function PPCM(X,Y:INTEGER):INTEGER;
var P:INTEGER;
begin
  P:=X*Y;
  while X#Y do
  begin
    if X>Y then X:=X-Y else Y:=Y-X;
    PPCM:=P div X
  end %PPCM%;

procedure INITIAL;
var I,J,P:INTEGER;
begin
  MAXA:=0;
  MAXB:=0;
  INIMAXSUM:=0;
  READ(M,N);
  DIM:=M+N;
  Writeln('M=',M:2,' N=',N:2);
  Writeln('A[I]:');
  for I:=1 to M do
    begin
      READ(P);
      WRITE(P:4);
      COEF[I]:=P;
      if MAXA<P then MAXA:=P
    end;
  Writeln;
  Writeln('B[J]:');
  for J:=1 to N do
    begin
      READ(P);
      WRITE(P:4);
      COEF[M+J]:=P;
      INIMAX[J]:=MAXA;
      INIMAXSUM:=INIMAXSUM+P;
      if MAXB<P then MAXB:=P
    end;
  INIMAXSUM:=INIMAXSUM*MAXA;
  Writeln;
  Writeln;
  LASTSOL:=-DIM;
  CDIM:=COEF[DIM];
  for I:=1 to M do
    for J:=1 to N do
      begin
        P:=PPCM(COEF[I],COEF[M+J]);
        D[I,J]:=P div COEF[I];
        E[I,J]:=P div COEF[M+J]-1
      end
end %INITIAL%;
procedure FILTER;
label
    1;
var Q,L:INTEGER;
begin
    L:=DIM %JUMP TRIVIAL SOL%;
    while L<=LASTSOL do
    begin
        %TEST WHETHER SOLS[L]<SOL%
        Q:=2 %BECAUSE ALWAYS TRUE OF 1%;
        while SOLS[L+Q]&=SOL[Q] do
        if Q=DIM then goto 1 else Q:=Q+1;
        L:=L+DIM
    end;
    LASTSOL:=LASTSOL+DIM %SOL HAS NOT BEEN FILTERED OUT%;
    for Q:=1 to DIM do SOLS[LASTSOL+Q]:=SOL[Q];
    1:
end;

procedure XLOOP(I,SUM,MAXSUM:INTEGER;MAXY:TABLE);
label
    1;
var K:INTEGER;
procedure YLOOP(I,REST:INTEGER);
label
    2;
var K:INTEGER;
begin
    if I#DIM then
        for K:=0 to MAXY[I-M] do
        begin
            SOL[I]:=K;
            YLOOP(I+1,REST);
            REST:=REST-COEF[I];
            if REST<0 then goto 2
        end
    else if REST mod CDIM=0 then
        begin
            K:=REST div CDIM;
            if K<=MAXY[N] then
                begin
                    SOL[DIM]:=K;
                    FILTER
                end
        end;
    2:
end %YLOOP%;
begin
    %BODY OF XLOOP%
    for K:=0 to MAXB do
    begin
        for L:=1 to N do
        if K=D[I,L] then
        if MAXY[L]>E[I,L] then
        begin
            MAXSUM:=MAXSUM-(MAXY[L]-E[I,L])*COEF[M+L];
            MAXY[L]:=E[I,L]
        end;
        if SUM>MAXSUM then goto 1;
        SOL[I]:=K;
        if I<M then XLOOP(I+1,SUM,MAXSUM,MAXY)else YLOOP(I+1,SUM);
        SUM:=SUM+COEF[I]
    end %FOR K%;
    1:
end %XLOOP%;
begin
INITIAL;
TIM1:=TIME('T');
XLOOP(1,0,INIMAXSUM,INIMAX);
TIM2:=TIME('T');
L:=DIM %JUMP TRIVIAL SOL%;
while L<=LASTSOL do
  begin
    for I:=1 to DIM do WRITE(SOLS[L+I]:4);
    WRITELN;
    L:=L+DIM
  end;
  for I:=1 to M do
    for J:=1 to N do
      begin
        for L:=1 to I-1 do WRITE('   0');
        WRITE(D[I,J]:4);
        for L:=I+1 to M+J-1 do WRITE('  0');
        WRITE(E[I,J]+1:4);
        for L:=J+1 to N do WRITE('   0');
        WRITELN
      end;
    WRITELN;
    WRITELN('CPU TIME: ',(TIM2-TIM1)/512:8,' SEC')
  end.
end.
## Appendix B : Computer outputs of the PASCAL program

### Example 1

```
M= 2 N= 3
A[I]:
    2    1
B[J]:
    1    1    2

    1    0    1    1    0
    1    0    2    0    0
    1    0    0    2    0
    1    0    0    0    1
    0    1    1    0    0
    0    1    0    1    0
    0    2    0    0    1

CPU TIME: 3.9E-03 SEC
```

### Example 2

```
M= 3 N= 3
A[I]:
    15    2    1
B[J]:
    3     5   10

    0    1    1    1    0    0
    0    1    3    0    1    0
    0    1    8    0    0    1
    0    2    1    0    1    0
    0    2    6    0    0    1
    0    3    4    0    0    1
    0    4    0    1    1    0
    0    4    2    0    0    1
    1    0    0    0    1    1
    1    0    1    2    0    1
    1    0    1    2    2    0
    1    0    2    4    1    0
    1    0    5    0    0    2
    1    1    0    4    1    0
    1    1    3    0    0    2
    1    2    0    3    0    1
    1    2    0    3    2    0
    1    2    1    0    0    2
    1    4    0    1    0    2
    2    1    0    4    0    2
    1    0    0    5    0    0
    1    0    0    0    3    0
    2    0    0    0    0    3
    0    3    0    2    0    0
    0    5    0    0    2    0
    0    5    0    0    0    1
    0    0    3    1    0    0
    0    0    5    0    1    0
    0    0    10   0    0    1

CPU TIME: 1.2E-01 SEC
```
### Example 3

$$M = 3 \quad N = 3$$

$$A[I]: \quad 9 \quad 5 \quad 2$$

$$B[J]: \quad 3 \quad 7 \quad 8$$

$$ \begin{matrix} 0 & 0 & 5 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 & 0 \\ 0 & 1 & 2 & 3 & 0 & 0 \end{matrix} $$

... (65 solutions generated)

$$ \begin{matrix} 6 & 0 & 1 & 0 & 0 & 7 \\ 6 & 0 & 1 & 0 & 8 & 0 \\ 7 & 0 & 0 & 0 & 1 & 7 \\ 1 & 0 & 0 & 3 & 0 & 0 \\ 7 & 0 & 0 & 0 & 9 & 0 \\ 8 & 0 & 0 & 0 & 0 & 9 \\ 0 & 3 & 0 & 5 & 0 & 0 \\ 0 & 7 & 0 & 0 & 5 & 0 \\ 0 & 8 & 0 & 0 & 0 & 5 \\ 0 & 0 & 3 & 2 & 0 & 0 \\ 0 & 0 & 7 & 0 & 2 & 0 \\ 0 & 0 & 4 & 0 & 0 & 1 \end{matrix} $$

CPU TIME: 1.4E 00 SEC

### References :

[1] Clifford A.H. and Preston G.B.
*The algebraic theory of semigroups.*
American Mathematical Society, Vol I, 1961, Vol II, 1967.

[2] Stickel M.E.
*A complete unification algorithm for associative-commutative functions.*
Proceedings of the 4th International Joint Conference on Artificial Intelligence, Tbilisi, 1975.

[3] Livesey M. and Siekmann J.
*Unification of A+C-terms (bags) and A+C+I-terms (sets).*
Interners Bericht Nr. 3/76, Institut für Informatik I, Universität Karlsruhe.