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 A^T.

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 2⁵³ 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