# Diophantine linear system solving

Diophantine linear system solving is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
task.

Implement the Havas-Majewski-Matthews LLL-based Hermite normal form algorithm for solving linear Diophantine systems.

clarification.

The point of this task is not comprehending the puzzles, it is implementing the LLLHnf algorithm.

The method is the result of an experimental refinement process over many iterations, terse to the point of being impenetrable, best copied verbatim from a reliable source, and hence this task mostly concerns the direct translation of existing (cryptic) code.

You may regard the test set as just random input to validate your solution, no need to delve any deeper. But to make the task a little nicer, and of course to demonstrate the power of the algorithm, the examples aren't really random.
Many are classics, with online resources abound. Others are on Rosetta Code in a different guise; some are copied from the HMM paper. Section headers like 'base cases' or 'polynomial coefficients' should be self-explanatory.
The output is deliberately left somewhat 'raw', so there's plenty of room for implementation dependent refinement. Also, to solve this task you're not obliged to click any wiki-links, but please do if you want some background information.

context.

Solving a system of linear Diophantine equations (i.e., solutions are required to be whole numbers) is a classic mathematical problem: given a coefficient matrix A (usually with less rows than columns) and vector b, the goal is to find a preferably small integer vector x that satisfies A·x = b.

How is this task different from Gaussian elimination? We now need a triangularization method that doesn't introduce fractions. The LLLHnf algorithm adapts Lenstra-Lenstra-Lovász lattice basis reduction to put the transpose of the input system into Hermite normal form, the integer analogue of the usual reduced row echelon form. In the process a unimodular transformation matrix is constructed from which (if either exists) minus the solution vector x and/or the null space basis for A can be immediately read off. These vectors are typically of small Euclidean length.

All details are found in the 1998 journal article. The relevant pseudocode (Algorithm 4) is on pages 131 and 129.

use.

LLLHnf is a versatile algorithm that will
▹find integer solutions for (nonsquare) linear systems
▹put integer matrices into Hermite normal form
▹invert unimodular integer matrices
▹compute modular inverses
▹solve Bézout's identity
▹solve Chinese remaindering problems
(not necessarily with pairwise coprime moduli)
▹solve the famous Coconut puzzle
▹give hints for packing your knapsack
▹find solutions of small norm for extended gcd problems
▹find a univariate polynomial for a given (complex) algebraic constant
In all cases: provided a solution exists.

LLLHnf puts the transpose of A into Hermite normal form, so to compute Hnf(A) one inputs AT.

The program rejects fractional matrix coefficients. Users will scale applicable rows to suitably sized integers.

limitation.

The output may be wrong if the (Gram-Schmidt) calculations produce numbers too big to be representable as integers (-larger than 253 if values are stored in 64-bit floats), though there is some tolerance.

input.

(Piped to FreeBasic stdin in this format:)
m : #rows
n : #columns
m x n coefficient matrix A (augmented with m-vector b).

```'a comment line starts with an apostrophe and contains no commas.
2
3
2 1 4| 17
-5 2 6|-13
```

To search for polynomial coefficients, follow a comment line containing the tag "const" with:
m : #digits precision
n : max. poly degree
a + b : real and imaginary parts of a complex constant.

```'constant sqrt(1 + sqrt(3)*i)
4
4
1.22474487 + .707106781
```

To quit:

```0
0
```
test vectors:
```'five base cases
'no integral solution
2
2
2 0| 1
2 1| 2
'indeterminate
2
3
1  3  5
4  6  8
'singular square
3
3
1  7  4
2  8  5
3  9  6
'overdetermined
3
2
2  1| 2
6  5| 2
7  6| 2
'square
3
3
2 -3  4| 9
5  6  7| 3
8  9 10| 3
'modular inverse
'(the smallest solution is negative
' add the modulus to make positive)
1
2
42 -2017| 1
'a classic Indian kuttaka problem
1
2
195 -221| 65
'Bachet de Méziriac "personnes en un banquet"
'(add null space vector to make positive men)
2
3
1 1 1| 41
12 9 1|120
'Malm
2
4
1   1   1   1|    80
165 235  85 389| 16324
'from the Sunzi Suanjing
3
4
1 3 0 0| 2
1 0 5 0| 3
1 0 0 7| 2
2
3
17 7  0|-1
11 0 15|-2
'from the Shushu jiuzhang
8
9
1 130   0   0   0   0   0   0   0|-60
1   0 110   0   0   0   0   0   0|-30
1   0   0 120   0   0   0   0   0|-10
1   0   0   0  60   0   0   0   0|-10
1   0   0   0   0  25   0   0   0| 10
1   0   0   0   0   0 100   0   0| 10
1   0   0   0   0   0   0  50   0| 10
1   0   0   0   0   0   0   0  20| 10
'5 sailor coconut puzzle
6
7
1 -5  0  0  0  0  0| 1
0  4 -5  0  0  0  0| 1
0  0  4 -5  0  0  0| 1
0  0  0  4 -5  0  0| 1
0  0  0  0  4 -5  0| 1
0  0  0  0  0  4 -5| 0
'unbounded knapsack with slack
3
6
3000 1800 2500 1 0 0|54500
3    2   20 0 1 0|  250
25   15    2 0 0 1|  250
'subset sum problem
1
9
575 436 1586 1030 1921 569 721 1183 1570| 6665
'HMM extended gcd (example 7.2)
1
10
763836 1066557 113192 1785102 1470060 3077752 114793 3126753 1997137 2603018| 1
'Fibonacci segment F7...F14 (example 7.3)
1
8
13 21 34 55 89 144 233 377| 1
'compute the inverse of transpose(P)
'(Fibonacci morphs into Lucas)
8
8
1  0  0  0  0  0  18 -7
1  1  0  0  0  0 -11  4
-1  1  1  0  0  0   7 -3
0 -1  1  1  0  0  -4  1
0  0 -1  1  1  0   3 -1
0  0  0 -1  1  1  -1  1
0  0  0  0 -1  1   1  0
0  0  0  0  0 -1  -1  0
'Hnf(A) with Aij = i^3 * j^2 + i + j (example 7.4)
10
10
3  11   31   69   131   223   351   521   739   1011
7  36  113  262   507   872  1381  2058  2927   4012
13  77  249  583  1133  1953  3097  4619  6573   9013
21 134  439 1032  2009  3466  5499  8204 11677  16014
31 207  683 1609  3135  5411  8587 12813 18239  25015
43 296  981 2314  4511  7788 12361 18446 26259  36016
57 401 1333 3147  6137 10597 16821 25103 35737  49017
73 522 1739 4108  8013 13838 21967 32784 46673  64018
91 659 2199 5197 10139 17511 27799 41489 59067  81019
111 812 2713 6414 12515 21616 34317 51218 72919 100020
'Gauss x*atan(1/239) + y*atan(1/57) + z*atan(1/18) = pi/4
'(fudge factor -1 to absorb round-off error
' ignore the corresponding vector entry x1)
1
4
-1 0041841 0175421 0554985| 7853982
'search for polynomial coefficients
'const sqrt(2) + i
4
4
1.41421356 + 1
'const 3^(1/3) + sqrt(2)
11
6
2.8564631326805
'some constant
12
9
-1.4172098692728
0
0
```
output.

Transformation matrix P (on the left) and the Hermite normal form of [A|b]T with messages:
inconsistent: the system is not solvable in integers.
-solution: is the negative of this particular P-vector (the 1 in the last column of P is a marker and no part of the solution).
The null space vectors in P are followed by their parenthesized lengths squared.

test results:
```'five base cases
'no integral solution
P | Hnf
0 -2  1 |  1  0  1
0  1  0 |  0  1  0
-1 -2  2 | -0 -0  2   inconsistent
loop 8

'indeterminate
P | Hnf
1  0  0  0 |  1  4  0
2  1 -1 -0 | -0  6 -0
0  0  0  1 |  0  0  1
1 -2  1  0 |  0  0  0   (6)
loop 11

'singular square
P | Hnf
1  0  0  0 |  1  2  3  0
3 -1  1  0 |  0  3  6  0
0  0  0  1 |  0  0  0  1
1  1 -2 -0 | -0  0  0 -0   (6)
loop 12

'overdetermined
P | Hnf
1 -1  0 |  1  1  1  0
-1  2 -0 | -0  4  5 -0
-2  2  1 |  0  0  0  1  -solution
loop 7

'square
P | Hnf
7 -1 -4  0 |  1  1   7  0
2  0 -1 -0 | -0  3   6 -0
15 -2 -9 -0 | -0 -0  12 -0
1  1 -2  1 |  0  0   0  1  -solution
loop 15

'modular inverse
'(the smallest solution is negative
' add the modulus to make positive)
P | Hnf
-48  -1  0 |  1  0
48   1  1 |  0  1  -solution
2017  42  0 |  0  0   (4070053)
loop 7

'a classic Indian kuttaka problem
P | Hnf
8   7  0 |  13  0
-6  -5  1 |   0  1  -solution
-17 -15  0 |   0  0   (514)
loop 8

'Bachet de Méziriac "personnes en un banquet"
'(add null space vector to make positive men)
P | Hnf
-3   4   0  0 |  1  0  0
3  -4   1  0 |  0  1  0
3 -14 -30  1 |  0  0  1  -solution
-8  11  -3 -0 | -0  0 -0   (194)
loop 12

'Malm
P | Hnf
-1   2   1  -1  0 |  1  1  0
2   3  -3  -2  0 |  0  2  0
-17 -22 -25 -16  1 |  0  0  1  -solution
-4  -8   7   5  0 |  0  0  0   (154)
-15   8   7   0  0 |  0  0  0   (338)
loop 23

'from the Sunzi Suanjing
P | Hnf
-35  12   7   5  0 |  1  0  0  0
21  -7  -4  -3  0 |  0  1  0  0
15  -5  -3  -2 -0 | -0 -0  1 -0
-23   7   4   3  1 |  0  0  0  1  -solution
-105  35  21  15  0 |  0  0  0  0   (12916)
loop 23

P | Hnf
-30   73  22  0 |  1  0  0
-49  119  36  0 |  0  1  0
-23   56  17  1 |  0  0  1  -solution
105 -255 -77  0 |  0  0  0   (81979)
loop 17

'from the Shushu jiuzhang
P | Hnf
1    0    0    0     0     0    0     0     0  0 |  1   1   1   1  1    1   1   1  0
-7800   60   71   65   130   312   78   156   390  0 |  0  10   0   0  0    0   0   0  0
-35750  275  325  298   596  1430  358   715  1788  0 |  0   0  10  10  0   50   0  10  0
0    0    0    0     1     0    0     0     0  0 |  0   0   0  60  0    0   0   0  0
-34320  264  312  286   572  1373  344   687  1716  0 |  0   0   0   0  5   80  30   0  0
0    0    0    0     0     0    1     0     0  0 |  0   0   0   0  0  100   0   0  0
0    0    0    0     0     0    0     1     0  0 |  0   0   0   0  0    0  50   0  0
0    0    0    0     0     0    0     0     1  0 |  0   0   0   0  0    0   0  20  0
-3710   29   34   31    62   148   37    74   185  1 |  0   0   0   0  0    0   0   0  1  -solution
85800 -660 -780 -715 -1430 -3432 -858 -1716 -4290 -0 | -0   0   0   0  0    0   0   0 -0   (7399103669)
loop 102

'5 sailor coconut puzzle
P | Hnf
1     0     0     0     0     0     0  0 |  1  0  0  0  0  0  0
-3905  -781  -625  -500  -400  -320  -256  0 |  0  1  0  0  0  0  0
-975  -195  -156  -125  -100   -80   -64  0 |  0  0  1  0  0  0  0
-5125 -1025  -820  -656  -525  -420  -336  0 |  0  0  0  1  0  0  0
-2500  -500  -400  -320  -256  -205  -164  0 |  0  0  0  0  1  0  0
-3125  -625  -500  -400  -320  -256  -205 -0 | -0 -0 -0 -0 -0  1 -0
-3121  -624  -499  -399  -319  -255  -204  1 |  0  0  0  0  0  0  1  -solution
15625  3125  2500  2000  1600  1280  1024  0 |  0  0  0  0  0  0  0   (269403226)
loop 71

'unbounded knapsack with slack
P | Hnf
0   0   0    1    0    0  0 |  1  0  0  0
0   0   0    0    1    0  0 | -0  1 -0 -0
0   0   0    0    0    1  0 |  0  0  1  0
-6  -5 -11    0   -2   -3  1 |  0  0  0  1  -solution
3  -5   0    0    1    0  0 |  0  0  0  0   (35)
-1   3  -1  100   17  -18  0 |  0  0  0  0   (10624)
-4  15  -6    0  102 -113  0 |  0  0  0  0   (23450)
loop 55

'subset sum problem
P | Hnf
0  0  0 -1  0  1 -1  1  0  0 |  1  0
-1  0 -1 -1  0  0 -1 -1 -1  1 |  0  1  -solution
-2 -1  1  0  0  0  0  0  0  0 |  0  0   (6)
1 -1  0 -2  1  0  0  0  0  0 |  0  0   (7)
0 -1  1  1 -1  0 -2  1  0  0 |  0  0   (9)
-1  1  1 -1  0 -2  1  0  0  0 |  0  0   (9)
0  0 -1  1  1 -1  0 -2  1  0 |  0  0   (9)
0 -1  2  0  1  0 -1 -2 -1  0 |  0  0   (12)
1  1  0 -1 -2  0  1  0  2  0 |  0  0   (12)
2 -1  0  0 -1 -1  3  1 -1  0 |  0  0   (18)
loop 66

'HMM extended gcd (example 7.2)
P | Hnf
-1  0  1 -3  1  3  3 -2 -2  2  0 |  1  0
1  0 -1  3 -1 -3 -3  2  2 -2  1 |  0  1  -solution
2  0  3 -1  0  0  0  1  1 -2  0 |  0  0   (20)
0 -1  2  2 -1 -1  3 -1  1  1  0 |  0  0   (23)
-2  0  0 -1  3 -3 -1  2  1  0  0 |  0  0   (29)
0  3  2  3  2 -3  1  0  0 -1  0 |  0  0   (37)
-2  2  2  0 -1  3 -3 -2 -1  0  0 |  0  0   (36)
2  2 -2 -5 -2  1  2  1  1  0  0 |  0  0   (48)
0  2  0 -2 -4 -1 -1  4 -1  0  0 |  0  0   (43)
-3  3 -1  2 -2  1  0  1  4 -6  0 |  0  0   (81)
0 -2  1 -2  3  5  4  1 -5 -3  0 |  0  0   (94)
loop 187

'Fibonacci segment F7...F14 (example 7.3)
P | Hnf
-7   4 -3  1 -1  1  0  0  0 |  1  0
7  -4  3 -1  1 -1  0  0  1 |  0  1  -solution
-1  -1  1  0  0  0  0  0  0 |  0  0   (3)
0   0 -1 -1  1  0  0  0  0 |  0  0   (3)
0   0  0  0 -1 -1  1  0  0 |  0  0   (3)
0  -1 -1  1  0  0  0  0  0 |  0  0   (3)
0   0  0 -1 -1  1  0  0  0 |  0  0   (3)
0   0  0  0  0 -1 -1  1  0 |  0  0   (3)
18 -11  7 -4  3 -1  1 -1  0 |  0  0   (522)
loop 45

'compute the inverse of transpose(P)
'(Fibonacci morphs into Lucas)
P | Hnf
2   1    5   3  0   -5    5   13  0 |  1  0  0  0  0  0  0  0  0
3   2    8   5  0   -8    8   21  0 |  0  1  0  0  0  0  0  0  0
4   3   13   8  0  -13   13   34  0 |  0  0  1  0  0  0  0  0  0
7   4   21  13  0  -21   21   55  0 |  0  0  0  1  0  0  0  0  0
11   7   33  21  0  -34   34   89  0 |  0  0  0  0  1  0  0  0  0
18  11   54  33  0  -55   55  144  0 |  0  0  0  0  0  1  0  0  0
29  18   87  54 -1  -89   89  233  0 |  0  0  0  0  0  0  1  0  0
47  29  141  87 -1 -145  144  377  0 |  0  0  0  0  0  0  0  1  0
0   0    0   0  0    0    0    0  1 |  0  0  0  0  0  0  0  0  1
loop 84

'Hnf(A) with Aij = i^3 * j^2 + i + j (example 7.4)
P | Hnf
-10  -8 -5  1  2  3  5  3  0 -4  0 |  1  0   7  22  45   76  115  162  217  280  0
-2  -1  0  1 -1  0  1  0  1 -1  0 |  0  1   4   9  16   25   36   49   64   81  0
-15 -11 -4  0  4  5  4  3  1 -5  0 |  0  0  12  36  72  120  180  252  336  432  0
0   0  0  0  0  0  0  0  0  0  1 |  0  0   0   0   0    0    0    0    0    0  1
-1   1  1  0 -2  1  0  0  0  0  0 |  0  0   0   0   0    0    0    0    0    0  0   (8)
0  -1  1  1 -1  1 -2  1  0  0  0 |  0  0   0   0   0    0    0    0    0    0  0   (10)
-1   0  1  1  1 -2  0 -1  1  0  0 |  0  0   0   0   0    0    0    0    0    0  0   (10)
-1   0  2 -1  1 -1  1 -1 -1  1  0 |  0  0   0   0   0    0    0    0    0    0  0   (12)
1   0 -1  0 -1 -1  1  2  0 -1  0 |  0  0   0   0   0    0    0    0    0    0  0   (10)
-1   1  0  1 -1  0  0  1 -2  1  0 |  0  0   0   0   0    0    0    0    0    0  0   (10)
-1   2 -1 -1  2  0 -2  1  0  0  0 |  0  0   0   0   0    0    0    0    0    0  0   (16)
loop 99

'Gauss x*atan(1/239) + y*atan(1/57) + z*atan(1/18) = pi/4
'(fudge factor -1 to absorb round-off error
' ignore the corresponding vector entry x1)
P | Hnf
-1  -0  -0  -0 -0 |  1 -0
-1   5  -8 -12  1 |  0  1  -solution
23  42 -29   6  0 |  0  0   (3170)
-71  33 -49  13  0 |  0  0   (8700)
18 -94 -63  27  0 |  0  0   (13858)
loop 20

'search for polynomial coefficients
'const sqrt(2) + i
1.41421356 + 1*i
P | Hnf
1    0   0   0   0   0  0 |  1  0  0
0    9   0  -2   0   1  0 |  0  1  0
0    0   0   0   0   0  1 |  0  0  1
-32   17 -10   0   2   0  0 |  0  0  0   (1417)
6   33   6 -48  35  -8  0 |  0  0  0   (4754)
18  -22 -56   5  31 -20  0 |  0  0  0   (5330)
102  151 -71  91 -78  36  0 |  0  0  0   (53907)
loop 53

'const 3^(1/3) + sqrt(2)
2.8564631326805
P | Hnf
1   0   0   0    0   0   0   0  0 |  1  0  0
0   0   0   0    0   0   0   0  1 |  0  0  1
5   1 -36  12   -6  -6   0   1  0 |  0  0  0   (1539)
-27  14   0  15  -31 -23 -26  13  0 |  0  0  0   (3485)
33  20  21  22  -13 -43   1   5  0 |  0  0  0   (4458)
-18  10 -29 -49  -12 -25 -13   9  0 |  0  0  0   (4685)
53 -14  15  13  -23  35 -13   1  0 |  0  0  0   (5323)
35  78 -17 -17    8   7   0  -1  0 |  0  0  0   (8001)
10  -5 -28  33  110  11 -84  23  0 |  0  0  0   (21804)
loop 125

'some constant
-1.4172098692728
P | Hnf
1   0   0   0   0  0   0   0   0   0  0  0 |  1  0  0
0   0   0   0   0  0   0   0   0   0  0  1 |  0  0  1
0   9   0   0  -7  0  -5   0   3   0  1  0 |  0  0  0   (165)
-2   9  -2  -9   4  8   2  -2 -12  -2  4  0 |  0  0  0   (422)
-4  -1   9   2  -5 -6  -2   6   6 -11 -9  0 |  0  0  0   (441)
-5  11   6   3   4 -3   5 -12  -6   3 -1  0 |  0  0  0   (431)
-9   1   9   5   4  1 -12  -8  11   1 -5  0 |  0  0  0   (560)
-6  -2   2  17 -11 -4   3   1  -4  -5  0  0 |  0  0  0   (521)
-1   1   0   5   3  8  12   6   4   9  5  0 |  0  0  0   (402)
-7   3   0 -13  -3 -1   4   9  -4   8  9  0 |  0  0  0   (495)
-9  11   7   2  11 -7  -7   3  -8  -2  3  0 |  0  0  0   (560)
-16  11 -13  -4   8 -4  -4   3   4   1  0  0 |  0  0  0   (684)
loop 360
```

To save server space, better not repeat full test results in your post. Once your solution covers all cases, a selected example will suffice.

## C

Translation of: FreeBASIC
```// Subject: Solution of an m x n linear Diophantine system
//          A*x = b using LLL reduction.
// Ref.   : G. Havas, B. Majewski, K. Matthews,
//         'Extended gcd and Hermite normal form
//          algorithms via lattice basis reduction,'
//          Experimental Mathematics 7 (1998), no.2, pp.125-136
// Code   : standard C
//          compile with (gnu compiler):
//          gcc filename.c -o diophantine -lm

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <stdbool.h>

// ---- NOTE ----
// these next few functions are useful to allocate and free the
// array d[] and the matrices la[][] and a[][]. (Numerical Recipes in C)
// useful to deal with negative array indices.

#define NR_END 1
#define FREE_ARG char*

void nrerror(char error_text[])
/* error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}

double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;

v = (double *)calloc((size_t) ((nh - nl + 1 + NR_END)), sizeof(double));
if (v) nrerror("allocation failure in dvector()");
return v - nl + NR_END;
}

double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow = nrh - nrl + 1, ncol = nch - ncl + 1;
double **m;

/* allocate pointers to rows */
m = (double **) calloc((size_t)((nrow + NR_END)), sizeof(double*));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;

/* allocate rows and set pointers to them */
m[nrl] = (double *) calloc((size_t)((nrow * ncol + NR_END)), sizeof(double));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

for(i = nrl + 1; i <= nrh; i++) m[i] = m[i - 1] + ncol;

/* return pointer to array of pointers to rows */
return m;
}

void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v + nl - NR_END));
}

void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl] + ncl - NR_END));
free((FREE_ARG) (m + nrl - NR_END));
}
// ------------------------------------------------

#define echo 1
#define cls

void swap(double *a, double *b)
{
double tmpd = *a;
*a = *b;
*b = tmpd;
}

void Main(int64_t sw);

// The complexity of the algorithm increases
// with alpha, as does the quality guarantee
// on the lattice basis vectors:
// alpha = aln / ald, 1/4 < alpha < 1
#define aln   80
#define ald   81
// rows & columns
int64_t m1 = 0, mn = 0, nx = 0, m = 0, n = 0;
// column indices
int64_t c1 = 0, c2 = 0;

// Gram-Schmidt coefficients
// mu_rs = lambda_rs / d_s
// Br = d_r / d_r-1
double **la = NULL, *d = NULL;
// work matrix
double **a = NULL;

//
// Part 1: driver, input and output
// ---------------------------------

int main()
{
char *g = alloca(1024);
int64_t i = 0, sw = 0;

while (1)
{
printf("\n");
sw = 0;
while (1)
{
printf(" rows ");
fgets(g, 1024, stdin);

if (g[0] == '\'') {
printf("%s\n", g);
} else {
break;
}
sw = strstr(g, "const") != NULL ? 1 : 0;
}
n = atof(g);
if (n < 1) break;

printf(" cols ");
fgets(g, 1024, stdin);
m = atoi(g);
if (m < 1) {
for (i = 1; i <= n; i++) {
fgets(g, 1024, stdin);
}
continue;
}

// set indices and allocate
if (sw) { sw = n - 1; n = 2; m += 2; }
m1 = m + 1; mn = m1 + n; nx = mn + 1;
la = dmatrix(0, m+1, 0, m+1);
d  = dvector(-1, m+1);
a  = dmatrix(0, m+1, 0, mn+1);

cls;
Main(sw);
printf("\n");
}

free_dmatrix(la, 0, m+1, 0, m+1);
free_dvector(d, -1, m+1);
free_dmatrix(a, 0, m+1, 0, mn+1);
}

// input complex constant, read powers into A
void Inpconst(int64_t pr)
{
int64_t r = 0, m2 = m1 + 1;
double p = 0, q = 0, t = 0, x = 0, y = 0;
char *g = alloca(1024);

printf(" a + bi: ");
fgets(g, 1024, stdin);
g = strtok(g, "+");
x = atof(g);
g = strtok(NULL , "");

printf("%lf", x);
if (g) { y = atof(g); printf(" + %lf*i", y); }
printf("\n");

// fudge factor 1
a[0][m1] = 1;
// c ^ 0
p = pow((double)(10), pr);
a[1][m1] = p;
// compute powers
for (r = 2; r <= m - 1; r++) {
t = p;
p = p * x - q * y;
q = t * y + q * x;
a[r][m1] = round(p);
a[r][m2] = round(q);
}
}

// input A and b
int64_t Inpsys()
{
int64_t r = 0, s = 0, sw = 0;
char *g = alloca(1024);

for (r = 0; r <= n - 1; r++) {
printf(" row A%ld and b%ld ", r + 1, r + 1);
fgets(g, 1024, stdin);

// reject all fractional coefficients
sw = strpbrk(g, "\\./") != NULL ? 1 : 0;

// parse row
char *token = NULL, *str = NULL;
s = 0;
for ( str = g; ; str = NULL) {
token = strtok(str, " |");
if (token == NULL)
break;
a[s][m1 + r] = atoi(token);
s ++;
}
}

if (sw) { printf("illegal input\n"); };
return sw;
}

// print row r
#define prow                                                \
for (s = 0; s <= mn; s++) {                             \
if (s == m1) { printf(" |"); }                        \
for (int64_t spc = 0; spc < p[s] - l[r][s] + 1; spc++) {  \
printf(" ");                                        \
}                                                     \
printf("% .0lf", a[r][s]);                            \
}

// print matrix A(,)
void PrntM(int64_t sw)
{
int64_t l[m+1][mn+1], p[mn+1], k = 0, r = 0, s = 0;
char *g = alloca(1024); double q = 0;

for (s = 0; s <= mn; s++) {
p[s] = 1; for (r = 0; r <= m; r++) {
// store lengths and max. length in column
// for pretty output
char *tmpg = alloca(1024);
sprintf(tmpg, "%f", fabs(a[r][s]));
l[r][s] = strlen(tmpg);
if (l[r][s] > p[s]) { p[s] = l[r][s]; }
}
}

if (sw) {
printf("P | Hnf\n");

// evaluate
for (r = 0; r <= m; r++) {
if (a[r][mn]) { k = r; break; }
}
sw = a[k][mn] == 1;
for (s = m1; s <= mn - 1; s++) {
sw &= a[k][s] == 0;
}
sw ? strcpy(g, "  -solution") : strcpy(g, "   inconsistent");
for (s = 0; s <= m - 1; s++) {
sw &= a[k][s] == 0;
}
if (sw) { g[0] = 0; } //  trivial

// Hnf and solution
for (r = m; r >= k; r--) {
prow;
printf("%s\n", r == k ? g : "");
}
// Null space with lengths squared
for (r = 0; r <= k - 1; r++) {
prow;
q = 0; for (s = 0; s <= m - 1; s++) {
q += a[r][s] * a[r][s];
}
printf("   (%.0f)\n", q);
}

} else {
printf("I | Ab~\n");

for (r = 0; r <= m; r++) {
prow;
printf("\n");
}
}
}

// ----------------------
// Part 2: HMM algorithm 4
// ------------------------

// negate rows t
void Minus(int64_t t)
{
int64_t r, s;
for (s = 0; s <= mn; s++) {
a[t][s] = -a[t][s];
}
for (r = 1; r <= m; r++) {
for (s = 0; s <= r - 1; s++) {
if (r == t || s == t) {
la[r][s] = -la[r][s];
}
}
}
}

// LLL reduce rows k
void Reduce(int64_t k, int64_t t)
{
int64_t s = 0, sx = 0;
double lk = 0, q = 0;
c1 = nx; c2 = nx;
// pivot elements Ab~ in rows t and k
for (s = m1; s <= mn; s++) {
if (a[t][s]) { c1 = s; break; }
}
for (s = m1; s <= mn; s++) {
if (a[k][s]) { c2 = s; break; }
}

q = 0;
if (c1 < nx) {
if (a[t][c1] < 0) { Minus(t); }
q = floor(a[k][c1] / a[t][c1]); // floor
} else {
lk = la[k][t];
if (2 * fabs(lk) > d[t]) {
// 2|lambda_kt| > d_t
// not LLL-reduced yet
q = round(lk / d[t]); // round;
}
}

if (q) {
sx = c1 == nx ? m : mn;
// reduce row k
for (s = 0; s <= sx; s++) {
a[k][s] -= q * a[t][s];
}
la[k][t] -= q * d[t];
for (s = 0; s <= t - 1; s++) {
la[k][s] -= q * la[t][s];
}
}
}

// exchange rows k and k-1
void Swop(int64_t k)
{
int64_t r = 0, s = 0, t = k - 1;
double db = 0, lk = 0, lr = 0;
for (s = 0; s <= mn; s++) {
swap(&a[k][s], &a[t][s]);
}
for (s = 0; s <= t - 1; s++) {
swap(&la[k][s], &la[t][s]);
}

// update Gram coefficients
// columns k, k-1 for r > k
lk = la[k][t];
db = (d[t - 1]*d[k] + lk*lk) / d[t];
for (r = k + 1; r <= m; r++) {
lr = la[r][k];
la[r][k] = (d[k] * la[r][t] - lk * lr) / d[t];
la[r][t] = (db * lr + lk * la[r][k]) / d[k];
}
d[t] = db;
}

// main limiting sequence
void Main(int64_t sw)
{
int64_t i = 0, k = 0, t = 0, tl = 0;
double db = 0, lk = 0;

if (sw) {
Inpconst(sw);
} else {
if (Inpsys()) { return; }
}
// augment Ab~ with column e_m
a[m][mn] = 1;
// prefix standard basis
for (i = 0; i <= m; i++) { a[i][i] = 1; }
// Gram sub-determinants
for (i = -1; i <= m; i++) { d[i] = 1; }

if (echo) { PrntM(0); }

k = 1;
while (k <= m)
{
t = k - 1;
// partial size reduction
Reduce(k, t);

sw = (c1 == nx && c2 == nx);
if (sw) {
// zero rows k-1, k
lk = la[k][t];
// Lovasz condition
// Bk >= (alpha - mu_kt^2) * Bt
db = d[t - 1] * d[k] + lk * lk;
// not satisfied
sw = (db * ald) < (d[t] * d[t] * aln);
}

if (sw || (c1 <= c2 && c1 < nx)) {
// test recommends a swap
Swop(k);
// decrease k
if (k > 1) { k -= 1; }
} else {
// complete size reduction
for (i = t - 1; i >= 0; i--) {
Reduce(k, i);
}
// increase k
k += 1;
}

tl += 1;
}

PrntM(-1);

printf("loop %ld\n", tl);
}
```
Output:
`The results are exactly the same as those described in the task assignment.`

## FreeBASIC

```'Subject: Solution of an m x n linear Diophantine system
'         A*x = b using LLL reduction.
'Ref.   : G. Havas, B. Majewski, K. Matthews,
'        'Extended gcd and Hermite normal form
'         algorithms via lattice basis reduction,'
'         Experimental Mathematics 7 (1998), no.2, pp.125-136
'Code   : FreeBasic 1.08.1

#define echo 1

Declare Sub Main (byval sw as integer)

'The complexity of the algorithm increases
'with alpha, as does the quality guarantee
'on the lattice basis vectors:
'alpha = aln / ald, 1/4 < alpha < 1
const aln = 80, ald = 81
'rows & columns
dim shared as integer m1, mn, nx, m, n
'column indices
dim shared as integer c1, c2

'Gram-Schmidt coefficients
'mu_rs = lambda_rs / d_s
'Br = d_r / d_r-1
dim shared as double la(any, any), d(any)
'work matrix
dim shared as double a(any, any)

'-------------------------------
'Part 1: driver, input and output
'---------------------------------

dim g as string
dim as integer i, sw

do
print
sw = 0
do
input ; " rows ", g
if instr(g, "'") then
print g
else
exit do
end if
sw or= instr(g, "const")
loop
n = val(g)
if n < 1 then exit do

input " cols ", g
m = val(g)
if m < 1 then
for i = 1 to n: input g: next
continue do
end if

'set indices and allocate
if sw then sw = n - 1: n = 2: m += 2
m1 = m + 1: mn = m1 + n: nx = mn + 1
redim as double la(m, m), d(-1 to m)
redim as double a(m, mn)

cls
Main sw
print
loop
end

'input complex constant, read powers into A
Sub Inpconst (byval pr as integer)
dim as integer r, m2 = m1 + 1
dim as double p, q = 0, t, x, y
dim g as string

input ; " a + bi:", g
g = trim(g) + " + "
r = instr(g, "+")
x = val(mid(g, 1, r - 1))
y = val(mid(g, r + 1))

print x;
if y then print " +";y;"*i";
print

'fudge factor 1
a(0, m1) = 1
'c ^ 0
p = cdbl(10) ^ pr
a(1, m1) = p
'compute powers
for r = 2 to m - 1
t = p
p = p * x - q * y
q = t * y + q * x
a(r, m1) = int(p + .5)
a(r, m2) = int(q + .5)
next r
End Sub

'input A and b
Function Inpsys () as integer
dim as integer i, j, k, r, s, t, sw = 0
dim g as string

for r = 0 to n - 1
g = str(r + 1)
print " row A"; g; " and b"; g;
input ; " ", g
'reject all fractional coefficients
sw or= instr(g, any "\./")

'parse row
i = 1: k = len(g)
for s = 0 to m
'locate next column separator (space or |)
for t = -1 to 0
j = i
for i = j to k
if (instr(" |", mid(g, i, 1)) = 0) = t then exit for
next i
next t
a(s, m1 + r) = val(mid(g, j, i - j))
next s: print
next r

if sw then print "illegal input"
Inpsys = sw
End Function

'print row r
#macro prow
for s = 0 to mn
if s = m1 then print " |";
print space(p(s) - l(r, s) + 1); a(r, s);
next s
#endmacro

'print matrix A(,)
Sub PrntM (byval sw as integer)
dim as integer l(m, mn), p(mn), k, r, s
dim g as string, q as double

for s = 0 to mn
p(s) = 1: for r = 0 to m
'store lengths and max. length in column
'for pretty output
l(r, s) = len(str(abs(a(r, s))))
if l(r, s) > p(s) then p(s) = l(r, s)
next r
next s

if sw then
print "P | Hnf"

'evaluate
for r = 0 to m
if a(r, mn) then k = r: exit for
next r
sw = a(k, mn) = 1
for s = m1 to mn - 1
sw and= a(k, s) = 0
next s
g = iif(sw,"  -solution","   inconsistent")
for s = 0 to m - 1
sw and= a(k, s) = 0
next s
if sw then g = "" ' trivial

'Hnf and solution
for r = m to k step -1
prow
print iif(r = k, g, "")
next r
'Null space with lengths squared
for r = 0 to k - 1
prow
q = 0: for s = 0 to m - 1
q += a(r, s) * a(r, s)
next s
print "   (";str(q);")"
next r

else
print "I | Ab~"

for r = 0 to m
prow
print
next r
end if
End Sub

'----------------------
'Part 2: HMM algorithm 4
'------------------------

'negate rows t
Sub Minus (byval t as integer)
dim as integer r, s
for s = 0 to mn
a(t, s) = -a(t, s)
next s
for r = 1 to m
for s = 0 to r - 1
if r = t or s = t then
la(r, s) = -la(r, s)
end if
next s
next r
End Sub

'LLL reduce rows k
Sub Reduce (byval k as integer, byval t as integer)
dim as integer s, sx
dim as double lk, q
c1 = nx: c2 = nx
'pivot elements Ab~ in rows t and k
for s = m1 to mn
if a(t, s) then c1 = s: exit for
next s
for s = m1 to mn
if a(k, s) then c2 = s: exit for
next s

q = 0
if c1 < nx then
if a(t, c1) < 0 then Minus t
q = int(a(k, c1) / a(t, c1)) ' floor
else
lk = la(k, t)
if 2 * abs(lk) > d(t) then
'2|lambda_kt| > d_t
'not LLL-reduced yet
q = int(lk / d(t) + .499) ' round
end if
end if

if q then
sx = iif(c1 = nx, m, mn)
'reduce row k
for s = 0 to sx
a(k, s) -= q * a(t, s)
next s
la(k, t) -= q * d(t)
for s = 0 to t - 1
la(k, s) -= q * la(t, s)
next s
end if
End Sub

'exchange rows k and k-1
Sub Swop (byval k as integer)
dim as integer r, s, t = k - 1
dim as double db, lk, lr
for s = 0 to mn
swap a(k, s), a(t, s)
next s
for s = 0 to t - 1
swap la(k, s), la(t, s)
next s

'update Gram coefficients
'columns k, k-1 for r > k
lk = la(k, t)
db = (d(t - 1) * d(k) + lk * lk) / d(t)
for r = k + 1 to m
lr = la(r, k)
la(r, k) = (d(k) * la(r, t) - lk * lr) / d(t)
la(r, t) = (db * lr + lk * la(r, k)) / d(k)
next r
d(t) = db
End Sub

'main limiting sequence
Sub Main (byval sw as integer)
dim as integer i, k, t, tl = 0
dim as double db, lk

if sw then
Inpconst sw
else
if Inpsys then exit sub
end if
'augment Ab~ with column e_m
a(m, mn) = 1
'prefix standard basis
for i = 0 to m: a(i, i) = 1: next
'Gram sub-determinants
for i = -1 to m: d(i) = 1: next

if echo then PrntM 0

k = 1
while k <= m
t = k - 1
'partial size reduction
Reduce k, t

sw = (c1 = nx and c2 = nx)
if sw then
'zero rows k-1, k
lk = la(k, t)
'Lovasz condition
'Bk >= (alpha - mu_kt^2) * Bt
db = d(t - 1) * d(k) + lk * lk
'not satisfied
sw = db * ald < d(t) * d(t) * aln
end if

if sw or (c1 <= c2 and c1 < nx) then
'test recommends a swap
Swop k
'decrease k
if k > 1 then k -= 1
else
'complete size reduction
for i = t - 1 to 0 step -1
Reduce k, i
next i
'increase k
k += 1
end if

tl += 1
wend

PrntM -1

print "loop"; tl
End Sub
```

## Nim

Translation of: freeBasic

As in Wren solution, an “i” is required for the imaginary part of a complex number. Also, numbers of the form “a - bi” are accepted. Note that, to allow arrays starting at index -1, we defined a sequence-like type with custom indexation functions.

```import std/[complex, math, parseutils, sequtils, strscans, strutils, terminal]

type
# Sequence of floats with first index at -1.
Vector = object
data: seq[float]
# Matrix of floats.
Matrix = object
data: seq[seq[float]]

func newVector(n: Natural): Vector =
## Create a vector with indices from -1 to "n".
Vector(data: newSeq[float](n + 2))

func newMatrix(m, n: Natural): Matrix =
## Create a matrix with indices from 0 to "m" and 0 to "n".
Matrix(data: newSeqWith(m + 1, newSeq[float](n + 1)))

func `[]`(v: Vector; i: int): float =
## Return the value at index "i".
v.data[i + 1]

func `[]=`(v: var Vector; i: int; value: float) =
## Set the value at index "i".
v.data[i + 1] = value

func `[]`(m: Matrix; i: int): seq[float] =
## Return the row at index "i".
m.data[i]

func `[]`(m: var Matrix; i: int): var seq[float] =
## Return the row at index "i" as a variable.
m.data[i]

proc input(prompt: string): string =
## Emit a prompt message and read a string.
stdout.write prompt
stdout.flushFile()
if not stdin.readLine(result):
quit "End of file encountered. Quitting.", QuitFailure
result = result.strip()

const Echo = true

# The complexity of the algorithm increases with alpha,
# as does the quality guarantee on the lattice basis vectors:
# alpha = aln / ald, 1/4 < alpha < 1

const
Aln = 80
Ald = 81

# Rows and columns.
var m1, mn, nx, m, n: int
# Column indices.
var c1, c2: int

# Gram-Schmidt coefficients:
# mu_rs = lambda_rs / d_s
# Br = d_r / d_r-1
var la: Matrix
var d: Vector
# Work matrix.
var a: Matrix

proc inpConst(pr: int) =
## Input complex constant, read powers into A.
let m2 = m1 + 1
var x, y: float
var op: char

var g = input(" a + bi: ").strip()

if not g.scanf("\$f\$s\$c\$s\$fi", x, op, y) or op notin "+-":
if scanf(g, "\$f\$.", x):
y = 0.0
op = '+'
else:
raise newException(ValueError, "wrong complex number")
if op == '-': y = -y
echo complex(x, y)

# Fudge factor 1.
a[0][m1] = 1
# c^0.
var p = 10.0^pr
a[1][m1] = p
var q = 0.0

# Compute powers.
for r in 2..<m:
let t = p
p = p * x - q * y
q = t * y + q * x
a[r][m1] = p.round
a[r][m2] = q.round

proc inpSys(): bool =
## Input A and b.
var sw = false

for r in 0..<n:
let g = input(" row A\$1 and b\$1 ".format(r + 1)).strip()
# Reject all fractional coefficients.
sw = sw or g.find({'.', '/'}) >= 0

# Parse row.
var s = 0
for token in g.split({' ', '|'}):
if token.len > 0:
if s > m: echo "Ignoring extra characters."
a[s][m1 + r] = token.parseFloat()
inc s

if sw: echo "illegal input"
result = sw

template prow() =
## Print row r.
for s in 0..mn:
if s == m1: stdout.write " |"
stdout.write spaces(p[s] - l[r][s] + 1)
stdout.write a[r][s].toInt

proc printM(sw: bool) =
## Print matrix A.
var l = newSeqWith(m + 1, newSeq[int](mn + 1))
var p = newSeq[int](mn + 1)
for s in 0..mn:
p[s] = 1
for r in 0..m:
# Store lengths and max. length in column for pretty output.
l[r][s] = len(\$a[r][s].toInt)
p[s] = max(p[s], l[r][s])

if sw:
echo "P | Hnf"

# Evaluate.
var k: int
for r in 0..m:
if a[r][mn] != 0:
k = r
break
var sw = a[k][mn] == 1
for s in m1..<mn:
sw = sw and a[k][s] == 0
var g = if sw: "  -solution" else: "   inconsistent"
for s in 0..<m:
sw = sw and a[k][s] == 0
if sw: g = ""   # Trivial.

# Hnf and solution.
for r in countdown(m, k):
prow()
echo if r == k: g else: ""
# Null space with lengths squared.
for r in 0..<k:
prow()
var q = 0.0
for s in 0..<m:
q += a[r][s]^2
echo "   (", q.toInt, ")"

else:
echo "I | Ab~"
for r in 0..m:
prow()
echo()

# HMM algorithm 4.

proc minus(t: int) =
## Negate rows t.
for s in 0..mn:
a[t][s] = -a[t][s]
for r in 1..m:
for s in 0..<r:
if r == t or s == t:
la[r][s] = -la[r][s]

proc reduce(k, t: int) =
## LLL reduce rows k.
c1 = nx
c2 = nx
# Pivot elements Ab~ in rows t and k.
for s in m1..mn:
if a[t][s] != 0:
c1 = s
break
for s in m1..mn:
if a[k][s] != 0:
c2 = s
break

var q = 0.0
if c1 < nx:
if a[t][c1] < 0: minus(t)
q = floor(a[k][c1] / a[t][c1])
else:
let lk = la[k][t]
if 2 * abs(lk) > d[t]:
# 2|lambda_kt| > d_t
# not LLL-reduced yet
q = round(lk / d[t])

if q != 0:
let sx = if c1 == nx: m else: mn
# Reduce row k.
for s in 0..sx:
a[k][s] -= q * a[t][s]
la[k][t] -= q * d[t]
for s in 0..<t:
la[k][s] -= q * la[t][s]

proc swop(k: int) =
## Exchange rows k and k-1.
let t = k - 1
for s in 0..mn:
swap a[k][s], a[t][s]
for s in 0..<t:
swap la[k][s], la[t][s]

# Update Gram coefficients columns k, k-1 for r > k.
let lk = la[k][t]
let db = (d[t - 1] * d[k] + lk * lk) / d[t]
for r in (k + 1)..m:
let lr = la[r][k]
la[r][k] = (d[k] * la[r][t] - lk * lr) / d[t]
la[r][t] = (db * lr + lk * la[r][k]) / d[k]
d[t] = db

proc main(sw: int) =
## Main limiting sequence.
if sw != 0: inpConst(sw)
elif inpSys(): return

# Augment Ab~ with column e_m.
a[m][mn] = 1

# Prefix standard basis.
for i in 0..m: a[i][i] = 1
# Gram sub-determinants.
for i in -1..m: d[i] = 1

if Echo: printM(false)

var k = 1
var tl = 0
while k <= m:
let t = k - 1
# Partial size reduction.
reduce(k, t)

var sw = (c1 == nx and c2 == nx)
if sw:
# Zero rows k-1, k.
let lk = la[k][t]
# Lovasz condition.
# Bk >= (alpha - mu_kt^2) * Bt.
let db = d[t - 1] * d[k] + lk * lk
# Not satisfied.
sw = db * Ald < d[t] * d[t] * Aln

if sw or c1 <= c2 and c1 < nx:
# Test recommends a swap.
swop(k)
if k > 1: dec k
else:
# Complete size reduction.
for i in countdown(t - 1, 0):
reduce(k, i)
inc k
inc tl

printM(true)
echo "loop ", tl

# Driver, input and output

var g: string

while true:
echo()
var sw = 0
while true:
g = input(" rows ")
if not g.startsWith('\''): break
echo g
sw = sw or ord("const" in g)
n = try: g.parseInt()
except ValueError: break
if n < 1: break

g = input(" cols ")
m = try: g.parseInt()
except ValueError: break
if m < 1:
for i in 1..n: g = input("")
continue

# Set indices and allocate.
if sw != 0:
sw = n - 1
n = 2
inc m, 2
m1 = m + 1
mn = m1 + n
nx = mn + 1
la = newMatrix(m, m)
d = newVector(m)
a = newMatrix(m, mn)

stdout.eraseScreen()
main(sw)
echo()
```
Output:

We provide the results for the first and last test, but all tests give with the expected result.

``` rows 'five base cases
'five base cases
rows 'no integral solution
'no integral solution
rows 2
cols 2

(on new page)
row A1 and b1  2 0| 1
row A2 and b2  2 1| 2
I | Ab~
1 0 0 | 2 2 0
0 1 0 | 0 1 0
0 0 1 | 1 2 1
P | Hnf
0 -2 1 | 1 0 1
0  1 0 | 0 1 0
-1 -2 2 | 0 0 2   inconsistent
loop 8
..........
..........
rows 'some constant
'some constant
rows 12
cols 9

(on new page)
a + bi: -1.4172098692728
(-1.4172098692728, 0.0)
I | Ab~
1 0 0 0 0 0 0 0 0 0 0 0 |              1 0 0
0 1 0 0 0 0 0 0 0 0 0 0 |   100000000000 0 0
0 0 1 0 0 0 0 0 0 0 0 0 |  -141720986927 0 0
0 0 0 1 0 0 0 0 0 0 0 0 |   200848381356 0 0
0 0 0 0 1 0 0 0 0 0 0 0 |  -284644308286 0 0
0 0 0 0 0 1 0 0 0 0 0 0 |   403400722935 0 0
0 0 0 0 0 0 1 0 0 0 0 0 |  -571703485815 0 0
0 0 0 0 0 0 0 1 0 0 0 0 |   810223822395 0 0
0 0 0 0 0 0 0 0 1 0 0 0 | -1148257197418 0 0
0 0 0 0 0 0 0 0 0 1 0 0 |  1627321432644 0 0
0 0 0 0 0 0 0 0 0 0 1 0 | -2306255994823 0 0
0 0 0 0 0 0 0 0 0 0 0 1 |              0 0 1
P | Hnf
1  0   0   0   0  0   0   0   0   0  0 0 | 1 0 0
0  0   0   0   0  0   0   0   0   0  0 1 | 0 0 1
0  9   0   0  -7  0  -5   0   3   0  1 0 | 0 0 0   (165)
-2  9  -2  -9   4  8   2  -2 -12  -2  4 0 | 0 0 0   (422)
-4 -1   9   2  -5 -6  -2   6   6 -11 -9 0 | 0 0 0   (441)
-5 11   6   3   4 -3   5 -12  -6   3 -1 0 | 0 0 0   (431)
-9  1   9   5   4  1 -12  -8  11   1 -5 0 | 0 0 0   (560)
-6 -2   2  17 -11 -4   3   1  -4  -5  0 0 | 0 0 0   (521)
-1  1   0   5   3  8  12   6   4   9  5 0 | 0 0 0   (402)
-7  3   0 -13  -3 -1   4   9  -4   8  9 0 | 0 0 0   (495)
-9 11   7   2  11 -7  -7   3  -8  -2  3 0 | 0 0 0   (560)
-16 11 -13  -4   8 -4  -4   3   4   1  0 0 | 0 0 0   (684)
loop 360
```

## Phix

You can run this online here.

```--
-- demo\rosetta\Diophantine_linear_system_solving.exw
-- ==================================================
--
--  Translation of FreeBasic, with some help from Wren, admittedly made harder
--  by the need to xlate 0 and -1 based idx to 1 based.
--
--  Note that for problem 16 (HMM extended gcd (example 7.2)), the signs of
--  the (20) and (37) rows are flipped, which I'm told is OK.
--
with javascript_semantics -- (using an include instead of file i/o)
include Diophantine_linear_system_constants.e -- DLS_PROBS/SOLNS
constant echo = false,
intext = split(DLS_PROBS,"\n"),
outtxt = split(DLS_SOLNS,"\n",false)

integer nxi = 1, nxo = 1
function input(string /*prompt*/)
string in_line = intext[nxi]
nxi += 1
return in_line
end function

procedure output(string out_line)
printf(1,"%s\n",{out_line})
string expected = outtxt[nxo]
if out_line != expected then
printf(1,"%s <<=== expected ***ERROR***\n",{expected})
{} = wait_key() -- (nb does nowt in a browser)
end if
nxo += 1
end procedure

-- The complexity of the algorithm increases
-- with alpha, as does the quality guarantee
-- on the lattice basis vectors:
-- alpha = aln / ald, 1/4 < alpha < 1
constant aln = 80, ald = 81
-- rows & columns
integer m1, mn, nx, m, n
-- column indices
integer col1, col2

-- Gram-Schmidt coefficients
-- mu_rs = lambda_rs / d_s
-- Br = d_r / d_r-1
sequence lambda, d
-- work matrix
sequence a

procedure InputAb_or_c(integer pr)
-- input A and b, or a complex constant and compute powers into A
if pr then -- (complex constant)
integer m2 = m1+1
atom p, q = 0, t
string g = input(" a + bi:")

integer plus = find('+',g)
string line = trim(g[1..plus-1])
atom x = to_number(line),
y = iff(plus?to_integer(trim(g[plus+1..\$])):0)
if y then line &= sprintf(" + %g*i",y) end if
output(" "&line)

-- fudge factor 1
a[1,m1+1] = 1
-- c ^ 0
p = power(10,pr)
a[2,m1+1] = p
-- compute powers
for r=3 to m do
t = p
p = p*x-q*y
q = t*y+q*x
a[r,m1+1] = round(p)
a[r,m2+1] = round(q)
end for
else -- (input A and b)
for r=1 to n do
--          printf(1," row A%d and b%d\n",r+1)
string g = input(" ")
-- reject all fractional coefficients
assert(not find_any(`\./`,g))
sequence gi = apply(split(substitute(g,'|',' ')),to_integer)
for s=1 to length(gi) do
a[s,m1+r] = gi[s]
end for
end for
end if
end procedure

function cf(sequence c) return apply(apply(c,sprint),length) end function
function col_formats(sequence c) return apply(true,sprintf,{{" %%%dd"},apply(apply(columnize(c),cf),maxsq)}) end function

function print_row(integer r, sequence fmts)
string line = ""
for s=1 to mn+1 do
if s=m1+1 then line &= " |" end if
line &= sprintf(fmts[s],a[r,s])
end for
return line
end function

-- print matrix a[,]
procedure PrntM(integer sw)
integer k, r, s
string g
atom q

sequence fmts = col_formats(a)

if sw then
output("P | Hnf")

-- evaluate
k = 1
for r=1 to m+1 do
if a[r,mn+1] then
k = r
exit
end if
end for
sw = a[k,mn+1]=1
for s=m1+1 to mn do
sw = sw and a[k,s]=0
end for
g = iif(sw,"  -solution","   inconsistent")
for s=1 to m do
sw = sw and a[k,s]=0
end for
if sw then g = "" end if -- trivial

sequence lensq = repeat("",m+1)
lensq[k] = g
-- Calculate lengths squared
for r=1 to k-1 do
q = 0
for s=1 to m+1 do
q += a[r,s]*a[r,s]
end for
lensq[r] &= sprintf("   (%d)",q)
end for

-- Hnf and solution and null space
sequence order = tagset(k,m+1,-1)&tagset(k-1)
for r=1 to m+1 do
integer rr = order[r]
output(print_row(rr,fmts)&lensq[rr])
end for

else
printf(1,"I | Ab~\n")

for r=1 to m+1 do
printf(1,print_row(r,fmts))
-- (not particularly helpful:)
--          printf(1,", L:")
--          for s=1 to m+1 do
--              printf(1," %.20g",lambda[r,s])
--          end for
printf(1,"\n")
end for
end if
end procedure

-- ----------------------
-- Part 2: HMM algorithm 4
-- ------------------------

-- negate rows t
procedure Minus(integer t)
a[t] = sq_uminus(a[t])
for r=1+1 to m+1 do
for s=1 to r+1 do
if r=t or s=t then
lambda[r,s] = -lambda[r,s]
end if
end for
end for
end procedure

-- LLL reduce rows k
procedure Reduce(integer k, t)
col1 = nx
col2 = nx
-- pivot elements Ab~ in rows t and k
for s=m1+1 to mn+1 do
if a[t,s] then
col1 = s-1
exit
end if
end for
for s=m1+1 to mn+1 do
if a[k,s] then
col2 = s-1
exit
end if
end for
atom q = 0
if col1<nx then
if a[t,col1+1]<0 then Minus(t) end if
q = floor(a[k,col1+1]/a[t,col1+1])
else
atom lk = lambda[k,t]
if 2*abs(lk)>d[t+1] then
-- 2|lambda_kt| > d_t
-- not LLL-reduced yet
q = round(lk/d[t+1])
end if
end if

if q then
integer sx = iif(col1=nx?m:mn)
-- reduce row k
for s=1 to sx+1 do
a[k,s] -= q*a[t,s]
end for
lambda[k,t] -= q*d[t+1]
for s=1 to t+1 do
lambda[k,s] -= q*lambda[t,s]
end for
end if
end procedure

-- exchange rows k and k-1
procedure Swop(integer k)
integer t = k-1, tm1 = t-1
{a[k], a[t]} = {a[t], a[k]}
object tmp = lambda[t][1..tm1]
lambda[t][1..tm1] = lambda[k][1..tm1]
lambda[k][1..tm1] = tmp; tmp = 0

-- update Gram coefficients
-- columns k, k-1 for r > k
atom lk = lambda[k,t],
db = (d[t]*d[k+1]+lk*lk)/d[k]
for r=k+1 to m+1 do
atom lr = lambda[r,k],
dk1 = d[k+1]
lambda[r,k] = (dk1*lambda[r,t]-lk*lr)/d[k]
lambda[r,t] = (db*lr+lk*lambda[r,k])/dk1
end for
d[k] = db
end procedure

integer problem_no = 0

-- main limiting sequence
procedure Main(integer sw)
problem_no += 1
printf(1,"problem #%d\n",problem_no)

InputAb_or_c(sw)

-- augment Ab~ with column e_m
a[m+1,mn+1] = 1
-- prefix standard basis
for i=1 to m+1 do a[i,i] = 1 end for
-- Gram sub-determinants
d = repeat(1,m+2)

if echo then PrntM(0) end if

integer k = 1,
tl = 0
while k<=m do
integer t = k-1
-- partial size reduction
Reduce(k+1, t+1)

sw = (col1=nx and col2=nx)
if sw then
-- zero rows k-1, k
atom lk = lambda[k+1,t+1]
-- Lovasz condition
-- Bk >= (alpha - mu_kt^2) * Bt
atom db = d[t+1]*d[k+2]+lk*lk
-- not satisfied
sw = db*ald<d[t+2]*d[t+2]*aln
end if

if sw or (col1<=col2 and col1<nx) then
-- test recommends a swap
Swop(k+1)
-- decrease k
if k>1 then k -= 1 end if
else
-- complete size reduction
for i=t-1 to 0 by -1 do
Reduce(k+1, i+1)
end for
-- increase k
k += 1
end if

tl += 1
end while

PrntM(-1)

output(sprintf("loop %d",tl))
end procedure

-- -------------------------------
-- Part 1: driver, input and output
-- ---------------------------------

while true do
printf(1,"\n")
integer sw = 0
string g
while true do
g = input(" rows ")
if match("'",g) then
output(g)
else
exit
end if
sw = sw or match("const",g)
end while
n = to_integer(trim(g))
if n<1 then exit end if

g = input(" cols ")
m = to_integer(trim(g))
if m<1 then
for i=1 to n do
g = input("")
end for
else
-- set indices and allocate
if sw then
sw = n-1
n = 2
m += 2
end if
m1 = m+1
mn = m1+n
nx = mn+1
lambda = repeat(repeat(0,m+1),m+1)
a = repeat(repeat(0,mn+1),m+1)

Main(sw)
output("")
output("")
end if
end while

?"done"
{} = wait_key()
```

## Wren

Translation of: FreeBASIC
Library: Wren-ioutil
Library: Wren-complex
Library: Wren-iterate
Library: Wren-seq

Results are the same as the FreeBASIC entry though I've just shown those for the first and last examples.

A reasonably faithful translation though I've made some slight changes to the way data is input.

```import "./ioutil" for Input
import "./complex" for Complex
import "./iterate" for Stepped
import "./seq" for Lst

var echo = true

// The complexity of the algorithm increases
// with alpha, as does the quality guarantee
// on the lattice basis vectors:
// alpha = aln / ald, 1/4 < alpha < 1
var aln = 80
var ald = 81

// rows and columns
var m1 = 0
var mn = 0
var nx = 0
var m  = 0
var n  = 0

// column indices
var c1 = 0
var c2 = 0

// Gram-Schmidt coefficients
// mu_rs = lambda_rs / d_s
// Br = d_r / d_r-1
var la = []
var d  = {} // use map instead of list to deal with an index of -1

// work matrix
var a = []

// input complex constant, read powers into 'a'
var inpConst = Fn.new { |pr|
var m2 = m1 + 1
var q = 0
var g = Input.text(" a + bi: ").trim()  // unlike FB, requires the 'i' for any imaginary part
var cmplx = Complex.fromString(g)
Complex.showAsReal = true
System.print(cmplx)
var x = cmplx.real
var y = cmplx.imag

// fudge factor 1
a[0][m1] = 1
// c ^ 0
var p = 10.pow(pr)
a[1][m1] = p

// compute powers
for (r in Stepped.ascend(2...m)) {
var t = p
p = p * x - q * y
q = t * y + q * x
a[r][m1] = p.round
a[r][m2] = q.round
}
}

// input A and b
var inpSys = Fn.new {
var sw = 0
var g = ""
for (r in 0...n) {
g = Input.text(" row A%(r+1) and b%(r+1) ")
// reject all fractional coefficients
sw = sw | (Lst.indexOfAny(g.toList, ["/", "."]) + 1)

// parse row
var i = 0
var k = g.count
for (s in 0..m) {
// locate next column separator (space or |)
var j
for (t in -1..0) {
j = i
while (i < k) {
if (((" |".indexOf(g[i]) == -1) ? -1 : 0) == t) break
i = i + 1
}
}
var e = Num.fromString(g[j...i])
a[s][m1+r] = e ? e : 0
}
}
if (sw != 0) System.print("illegal input")
return sw
}

// print row r
var prow = Fn.new { |r, l, p|
for (s in 0..mn) {
if (s == m1) System.write(" |")
System.write(" " * (p[s] - l[r][s] + 1))
System.write(a[r][s])
}
}

// print matrix A
var printM = Fn.new { |sw|
var l = List.filled(m+1, null)
for (i in 0..m) l[i] = List.filled(mn+1, 0)
var p = List.filled(mn+1, 0)
for (s in 0..mn) {
p[s] = 1
for (r in 0..m) {
// store lengths and max. length in column
// for pretty output
l[r][s] = a[r][s].toString.count
if (l[r][s] > p[s]) p[s] = l[r][s]
}
}

if (sw != 0) {
System.print("P | Hnf")

// evaluate
var k = 0
for (r in 0..m) {
if (a[r][mn] != 0) {
k = r
break
}
}
sw = (a[k][mn] == 1) ? -1 : 0
for (s in m1...mn) sw = sw & ((a[k][s] == 0) ? -1: 0)
var g = (sw != 0) ? "  -solution" : "   inconsistent"
for (s in 0...m) sw = sw & ((a[k][s] == 0) ? -1: 0)
if (sw != 0) g = "" // trivial

// Hnf and solution
for (r in Stepped.descend(m..k)) {
prow.call(r, l, p)
System.print((r == k) ? g : "")
}

// Null space with lengths squared
for (r in 0...k) {
prow.call(r, l, p)
var q = 0
for (s in 0...m) q = q + a[r][s] * a[r][s]
System.print("   (%(q))")
}
} else {
System.print("I | Ab~")
for (r in 0..m) {
prow.call(r, l, p)
System.print()
}
}
}

/* HMM algorithm 4 */

// negate rows t
var minus = Fn.new { |t|
for (s in 0..mn) a[t][s] = -a[t][s]
for (r in 1..m) {
for (s in 0...r) {
if (r == t || s == t) la[r][s] = -la[r][s]
}
}
}

// LLL reduce rows k
var reduce = Fn.new { |k, t|
c1 = nx
c2 = nx

// pivot elements Ab~ in rows t and k
for (s in m1..mn) {
if (a[t][s] != 0) {
c1 = s
break
}
}
for (s in m1..mn) {
if (a[k][s] != 0) {
c2 = s
break
}
}

var q = 0
if (c1 < nx) {
if (a[t][c1] < 0) minus.call(t)
q = (a[k][c1] / a[t][c1]).floor
} else {
var lk = la[k][t]
if (2 * lk.abs > d[t]) {
// 2|lambda_kt| > d_t
// not LLL-reduced yet
q = (lk/d[t]).round
}
}

if (q != 0) {
var sx = (c1 == nx) ? m : mn

// reduce row k
for (s in 0..sx) a[k][s] = a[k][s] - q * a[t][s]
la[k][t] = la[k][t] - q * d[t]
for (s in 0...t) la[k][s] = la[k][s] - q * la[t][s]
}
}

// exchange rows k and k - 1
var swop = Fn.new { |k|
var t = k - 1
for (s in 0..mn) {
var tmp = a[k][s]
a[k][s] = a[t][s]
a[t][s] = tmp
}
for (s in 0...t) {
var tmp = la[k][s]
la[k][s] = la[t][s]
la[t][s] = tmp
}

// update Gram coefficients
// columns k, k-1 for r > k
var lk = la[k][t]
var db = (d[t-1] * d[k] + lk * lk) / d[t]
for (r in Stepped.ascend(k+1..m)) {
var lr = la[r][k]
la[r][k] = (d[k] * la[r][t] - lk * lr) / d[t]
la[r][t] = (db * lr + lk * la[r][k]) / d[k]
}
d[t] = db
}

// main limiting sequence
var main = Fn.new { |sw|
if (sw != 0) {
inpConst.call(sw)
} else if (inpSys.call() != 0) {
return
}
// augment Ab~ with column e_m
a[m][mn] = 1

// prefix standard basis
for (i in 0..m) a[i][i] = 1

// Gram sum-determinants
for (i in -1..m) d[i] = 1

if (echo) printM.call(0)
var tl = 0
var k = 1
while (k <= m) {
var t = k - 1

// partial size reduction
reduce.call(k, t)

sw = (c1 == nx && c2 == nx) ? -1 : 0
if (sw != 0) {
// zero rows k-1, k
var lk = la[k][t]

// Lovasz condition
// Bk >= (alpha - mu_kt^2) * Bt
var db = d[t-1] * d[k] + lk * lk

// not satisfied
sw = (db * ald < d[t] * d[t] * aln) ? -1 : 0
}
if (sw != 0 || (c1 <= c2 && c1 < nx)) {
// test recommends a swap
swop.call(k)

// decrease k
if (k > 1) k = k - 1
} else {
// complete size reduction
for (i in Stepped.descend(t-1..0)) reduce.call(k, i)

// increase k
k = k + 1
}
tl = tl + 1
}
printM.call(-1)
System.print("loop %(tl)")
}

/* driver, input and output */

var g  = ""
var sw = 0
while (true) {
System.print()
sw = 0
while (true) {
g = Input.text(" rows ")
if (g.indexOf("'") >= 0) {
System.print(g)
} else {
break
}
sw = sw | (g.indexOf("const") + 1)
}
n = Num.fromString(g)
if (!n || n < 1) break

g = Input.text(" cols ")
m = Num.fromString(g)
if (!m || m < 1) {
for (i in 1..n) g = Input.text("")
continue
}

// set indices and allocate
if (sw != 0) {
sw = n - 1
n = 2
m = m + 2
}
m1 = m + 1
mn = m1 + n
nx = mn + 1
la = List.filled(m+1, null)
for (i in 0..m) la[i] = List.filled(m+1, 0)
for (i in -1..m) d[i] = 0
a = List.filled(m+1, null)
for (i in 0..m) a[i] = List.filled(mn+1, 0)
System.write("\e[2J")   // clear the terminal
System.write("\e[0;0H") // home the cursor
main.call(sw)
System.print()
}
```
Output:
```(first example)
rows 'five base cases
'five base cases
rows 'no integral solution
'no integral solution
rows 2
cols 2

(on new page)
row A1 and b1 2 0| 1
row A2 and b2 2 1| 2
I | Ab~
1 0 0 | 2 2 0
0 1 0 | 0 1 0
0 0 1 | 1 2 1
P | Hnf
0 -2 1 |  1  0 1
0  1 0 |  0  1 0
-1 -2 2 | -0 -0 2   inconsistent
loop 8

...
...

(last example)
rows 'some constant
'some constant
rows 12
cols 9

(on new page)
a + bi: -1.4172098692728
-1.4172098692728
I | Ab~
1 0 0 0 0 0 0 0 0 0 0 0 |              1  0 0
0 1 0 0 0 0 0 0 0 0 0 0 |   100000000000  0 0
0 0 1 0 0 0 0 0 0 0 0 0 |  -141720986927  0 0
0 0 0 1 0 0 0 0 0 0 0 0 |   200848381356 -0 0
0 0 0 0 1 0 0 0 0 0 0 0 |  -284644308286  0 0
0 0 0 0 0 1 0 0 0 0 0 0 |   403400722935 -0 0
0 0 0 0 0 0 1 0 0 0 0 0 |  -571703485815  0 0
0 0 0 0 0 0 0 1 0 0 0 0 |   810223822395 -0 0
0 0 0 0 0 0 0 0 1 0 0 0 | -1148257197418  0 0
0 0 0 0 0 0 0 0 0 1 0 0 |  1627321432644 -0 0
0 0 0 0 0 0 0 0 0 0 1 0 | -2306255994823  0 0
0 0 0 0 0 0 0 0 0 0 0 1 |              0  0 1
P | Hnf
1  0   0   0   0  0   0   0   0   0  0 0 | 1  0 0
0  0   0   0   0  0   0   0   0   0  0 1 | 0  0 1
0  9   0   0  -7  0  -5   0   3   0  1 0 | 0  0 0   (165)
-2  9  -2  -9   4  8   2  -2 -12  -2  4 0 | 0  0 0   (422)
-4 -1   9   2  -5 -6  -2   6   6 -11 -9 0 | 0 -0 0   (441)
-5 11   6   3   4 -3   5 -12  -6   3 -1 0 | 0  0 0   (431)
-9  1   9   5   4  1 -12  -8  11   1 -5 0 | 0 -0 0   (560)
-6 -2   2  17 -11 -4   3   1  -4  -5  0 0 | 0  0 0   (521)
-1  1   0   5   3  8  12   6   4   9  5 0 | 0  0 0   (402)
-7  3   0 -13  -3 -1   4   9  -4   8  9 0 | 0 -0 0   (495)
-9 11   7   2  11 -7  -7   3  -8  -2  3 0 | 0  0 0   (560)
-16 11 -13  -4   8 -4  -4   3   4   1  0 0 | 0 -0 0   (684)
loop 360
```