Diophantine linear system solving
- task.
Implement the Havas-Majewski-Matthews LLL-based Hermite normal form algorithm for solving linear Diophantine systems.
- 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.
- clarification.
The point of this task is not comprehending the above puzzles, it is implementing the LLLHnf algorithm.
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.
FreeBASIC
<lang 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 </lang>
Phix
A mechanical and (no disrespect) somewhat laborious translation of FreeBasic, with a bit of help where needed from the Wren entry.
(Admittedly made harder by the need to xlate 0 and -1 based idx to 1 based.)
-- -- demo\rosetta\Diophantine_linear_system_solving.exw -- ================================================== -- without js -- (due to file reading) constant echo = false, intext = get_text("Diophantine_linear_system_problems.txt",GT_LF_STRIPPED), outtxt = get_text("Diophantine_linear_system_solution.txt",GT_LF_STRIPPED) integer nxi = 1, nxo = 1 function input(string prompt) string in = intext[nxi] nxi += 1 return in 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() 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 InputABC(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 r = find('+',g) string line = trim(g[1..r-1]) atom x = to_number(line), y = iff(r?to_integer(trim(g[r+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, line 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 -- Hnf and solution for r=m+1 to k by -1 do output(print_row(r,fmts)&iif(r=k?g:"")) end for -- Null space with lengths squared for r=1 to k-1 do line = print_row(r,fmts) q = 0 for s=1 to m+1 do q += a[r,s]*a[r,s] end for output(line&sprintf(" (%d)",q)) end for else printf(1,"I | Ab~\n") for r=1 to m+1 do printf(1,print_row(r,fmts)&"\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) integer s, sx atom lk, q 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 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 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 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 r, s, t = k-1 atom db, lk, lr {a[k], a[t]} = {a[t], a[k]} {lambda[k][1..t-1], lambda[t][1..t-1]} = {lambda[t][1..t-1], lambda[k][1..t-1]} -- update Gram coefficients -- columns k, k-1 for r > k lk = lambda[k,t] db = (d[t]*d[k+1]+lk*lk)/d[t+1] for r=k+1 to m+1 do lr = lambda[r,k] lambda[r,k] = (d[k+1]*lambda[r,t]-lk*lr)/d[t+1] lambda[r,t] = (db*lr+lk*lambda[r,k])/d[k+1] end for d[t+1] = db end procedure integer problem_no = 0 -- main limiting sequence procedure Main(integer sw) problem_no += 1 printf(1,"problem #%d\n",problem_no) integer i, k, t, tl = 0 atom db, lk InputABC(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 k = 1 while k<=m do 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 lk = lambda[k+1,t+1] -- Lovasz condition -- Bk >= (alpha - mu_kt^2) * Bt 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
Note that for problem 16 (HMM extended gcd (example 7.2)), the signs of the (20) and (37) rows in the null space are flipped, somehow, otherwise output matches that of FreeBasic (and is auto-checked). (I suspect the issue is FreeBasic distinguishing 0 and -0 in a way that Phix does not) Eagle eyed viewers may have spotted my blase use of null space there as if I knew what I was talking about (tee hee).
Wren
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. <lang ecmascript>import "./ioutil" for Input import "./complex" for Complex import "./trait" 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()
}</lang>
- 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