I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# 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.

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 & columnsint64_t m1 = 0, mn = 0, nx = 0, m = 0, n = 0;// column indicesint64_t c1 = 0, c2 = 0; // Gram-Schmidt coefficients// mu_rs = lambda_rs / d_s// Br = d_r / d_r-1double **la = NULL, *d = NULL;// work matrixdouble **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 Avoid 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 bint64_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 tvoid 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 kvoid 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-1void 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 sequencevoid 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 < 1const aln = 80, ald = 81'rows & columnsdim shared as integer m1, mn, nx, m, n'column indicesdim shared as integer c1, c2 'Gram-Schmidt coefficients'mu_rs = lambda_rs / d_s'Br = d_r / d_r-1dim shared as double la(any, any), d(any)'work matrixdim shared as double a(any, any) '-------------------------------'Part 1: driver, input and output'--------------------------------- dim g as stringdim 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   printloopend  'input complex constant, read powers into ASub Inpconst (byval pr as integer)dim as integer r, m2 = m1 + 1dim as double p, q = 0, t, x, ydim 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 rEnd Sub 'input A and bFunction Inpsys () as integerdim as integer i, j, k, r, s, t, sw = 0dim 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 = swEnd 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, sdim 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 rnext 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 rend ifEnd Sub '----------------------'Part 2: HMM algorithm 4'------------------------ 'negate rows tSub 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 rEnd Sub 'LLL reduce rows kSub Reduce (byval k as integer, byval t as integer)dim as integer s, sxdim 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 ifEnd Sub 'exchange rows k and k-1Sub Swop (byval k as integer)dim as integer r, s, t = k - 1dim 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) = dbEnd Sub 'main limiting sequenceSub Main (byval sw as integer)dim as integer i, k, t, tl = 0dim as double db, lk if sw then   Inpconst swelse   if Inpsys then exit subend if'augment Ab~ with column e_ma(m, mn) = 1'prefix standard basisfor i = 0 to m: a(i, i) = 1: next'Gram sub-determinantsfor i = -1 to m: d(i) = 1: next if echo then PrntM 0 k = 1while 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 += 1wend PrntM -1 print "loop"; tlEnd Sub `

## 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))
--          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-trait
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 Inputimport "./complex" for Compleximport "./trait" for Steppedimport "./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 < 1var aln = 80var ald = 81 // rows and columnsvar m1 = 0var mn = 0var nx = 0var m  = 0var n  = 0 // column indicesvar c1 = 0var c2 = 0 // Gram-Schmidt coefficients// mu_rs = lambda_rs / d_s// Br = d_r / d_r-1var la = []var d  = {} // use map instead of list to deal with an index of -1 // work matrixvar 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 bvar 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 rvar 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 Avar 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 tvar 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 kvar 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 - 1var 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 sequencevar 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 = 0while (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
```