Diophantine linear system solving: Difference between revisions
mNo edit summary |
m →{{header|Wren}}: Changed to Wren S/H |
||
(10 intermediate revisions by 5 users not shown) | |||
Line 1: | Line 1: | ||
{{draft task|matrices}} |
{{draft task|matrices}} |
||
{{Difficulty}} |
|||
;task. |
;task. |
||
Implement the Havas-Majewski-Matthews |
Implement the Havas-Majewski-Matthews |
||
Line 6: | Line 5: | ||
[[wp:Hermite_normal_form| Hermite normal form]] algorithm for solving |
[[wp:Hermite_normal_form| Hermite normal form]] algorithm for solving |
||
[[wp:Diophantine_equation#Linear_Diophantine_equations| linear Diophantine systems]]. |
[[wp:Diophantine_equation#Linear_Diophantine_equations| 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.<br/> |
|||
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.<br/> |
|||
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. |
;context. |
||
Line 513: | Line 529: | ||
Once your solution covers all cases, a selected example will suffice. |
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. |
|||
__TOC__ |
|||
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.<br/> |
|||
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.<br/> |
|||
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. |
|||
__TOC__ |
|||
=={{header|C}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="c"> |
|||
// 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); |
|||
} |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>The results are exactly the same as those described in the task assignment.</pre> |
|||
=={{header|FreeBASIC}}== |
=={{header|FreeBASIC}}== |
||
< |
<syntaxhighlight lang="freebasic"> |
||
'Subject: Solution of an m x n linear Diophantine system |
'Subject: Solution of an m x n linear Diophantine system |
||
' A*x = b using LLL reduction. |
' A*x = b using LLL reduction. |
||
Line 865: | Line 1,315: | ||
print "loop"; tl |
print "loop"; tl |
||
End Sub |
End Sub |
||
</syntaxhighlight> |
|||
</lang> |
|||
=={{header|Nim}}== |
|||
{{trans|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. |
|||
<syntaxhighlight lang="Nim">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() |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
We provide the results for the first and last test, but all tests give with the expected result. |
|||
<pre> 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 |
|||
</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
You can run this online [http://phix.x10.mx/p2js/LLLHnf.htm here]. |
|||
A mechanical and (no disrespect) somewhat laborious translation of FreeBasic, with a bit of help where needed from the Wren entry.<br> |
|||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
(Admittedly made harder by the need to xlate 0 and -1 based idx to 1 based.) |
|||
<!--<lang Phix>(notonline)--> |
|||
<span style="color: #000080;font-style:italic;">-- |
<span style="color: #000080;font-style:italic;">-- |
||
-- demo\rosetta\Diophantine_linear_system_solving.exw |
-- 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. |
|||
--</span> |
--</span> |
||
<span style="color: #008080;"> |
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> <span style="color: #000080;font-style:italic;">-- (using an include instead of file i/o)</span> |
||
<span style="color: #008080;">include</span> <span style="color: #000000;">Diophantine_linear_system_constants</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> <span style="color: #000080;font-style:italic;">-- DLS_PROBS/SOLNS</span> |
|||
<span style="color: #008080;">constant</span> <span style="color: #000000;">echo</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">false</span><span style="color: #0000FF;">,</span> |
<span style="color: #008080;">constant</span> <span style="color: #000000;">echo</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">false</span><span style="color: #0000FF;">,</span> |
||
<span style="color: #000000;">intext</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;"> |
<span style="color: #000000;">intext</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">split</span><span style="color: #0000FF;">(</span><span style="color: #000000;">DLS_PROBS</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">),</span> |
||
<span style="color: #000000;">outtxt</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;"> |
<span style="color: #000000;">outtxt</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">split</span><span style="color: #0000FF;">(</span><span style="color: #000000;">DLS_SOLNS</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">,</span><span style="color: #004600;">false</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">integer</span> <span style="color: #000000;">nxi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nxo</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> |
<span style="color: #004080;">integer</span> <span style="color: #000000;">nxi</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nxo</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> |
||
<span style="color: #008080;">function</span> <span style="color: #000000;">input</span><span style="color: #0000FF;">(</span><span style="color: #004080;">string</span> <span style="color: # |
<span style="color: #008080;">function</span> <span style="color: #000000;">input</span><span style="color: #0000FF;">(</span><span style="color: #004080;">string</span> <span style="color: #000080;font-style:italic;">/*prompt*/</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">string</span> <span style="color: # |
<span style="color: #004080;">string</span> <span style="color: #000000;">in_line</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">intext</span><span style="color: #0000FF;">[</span><span style="color: #000000;">nxi</span><span style="color: #0000FF;">]</span> |
||
<span style="color: #000000;">nxi</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
<span style="color: #000000;">nxi</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
||
<span style="color: #008080;">return</span> <span style="color: # |
<span style="color: #008080;">return</span> <span style="color: #000000;">in_line</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
||
Line 892: | Line 1,739: | ||
<span style="color: #008080;">if</span> <span style="color: #000000;">out_line</span> <span style="color: #0000FF;">!=</span> <span style="color: #000000;">expected</span> <span style="color: #008080;">then</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">out_line</span> <span style="color: #0000FF;">!=</span> <span style="color: #000000;">expected</span> <span style="color: #008080;">then</span> |
||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s <<=== expected ***ERROR***\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">expected</span><span style="color: #0000FF;">})</span> |
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s <<=== expected ***ERROR***\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">expected</span><span style="color: #0000FF;">})</span> |
||
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">wait_key</span><span style="color: #0000FF;">()</span> |
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">wait_key</span><span style="color: #0000FF;">()</span> <span style="color: #000080;font-style:italic;">-- (nb does nowt in a browser)</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
<span style="color: #000000;">nxo</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
<span style="color: #000000;">nxo</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
||
Line 914: | Line 1,761: | ||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span> |
<span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span> |
||
<span style="color: #008080;">procedure</span> <span style="color: #000000;"> |
<span style="color: #008080;">procedure</span> <span style="color: #000000;">InputAb_or_c</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">pr</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #000080;font-style:italic;">-- input A and b, or a complex constant and compute powers into A</span> |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">pr</span> <span style="color: #008080;">then</span> <span style="color: #000080;font-style:italic;">-- (complex constant)</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">pr</span> <span style="color: #008080;">then</span> <span style="color: #000080;font-style:italic;">-- (complex constant)</span> |
||
<span style="color: #004080;">integer</span> <span style="color: #000000;">m2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> |
<span style="color: #004080;">integer</span> <span style="color: #000000;">m2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">m1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> |
||
Line 921: | Line 1,768: | ||
<span style="color: #004080;">string</span> <span style="color: #000000;">g</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">input</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" a + bi:"</span><span style="color: #0000FF;">)</span> |
<span style="color: #004080;">string</span> <span style="color: #000000;">g</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">input</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" a + bi:"</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">integer</span> <span style="color: #000000;"> |
<span style="color: #004080;">integer</span> <span style="color: #000000;">plus</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">find</span><span style="color: #0000FF;">(</span><span style="color: #008000;">'+'</span><span style="color: #0000FF;">,</span><span style="color: #000000;">g</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">string</span> <span style="color: #000000;">line</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">trim</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;"> |
<span style="color: #004080;">string</span> <span style="color: #000000;">line</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">trim</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">plus</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span> |
||
<span style="color: #004080;">atom</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">to_number</span><span style="color: #0000FF;">(</span><span style="color: #000000;">line</span><span style="color: #0000FF;">),</span> |
<span style="color: #004080;">atom</span> <span style="color: #000000;">x</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">to_number</span><span style="color: #0000FF;">(</span><span style="color: #000000;">line</span><span style="color: #0000FF;">),</span> |
||
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;"> |
<span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">plus</span><span style="color: #0000FF;">?</span><span style="color: #7060A8;">to_integer</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">trim</span><span style="color: #0000FF;">(</span><span style="color: #000000;">g</span><span style="color: #0000FF;">[</span><span style="color: #000000;">plus</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..$])):</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #008080;">if</span> <span style="color: #000000;">y</span> <span style="color: #008080;">then</span> <span style="color: #000000;">line</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" + %g*i"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">y</span> <span style="color: #008080;">then</span> <span style="color: #000000;">line</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" + %g*i"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
<span style="color: #000000;">output</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" "</span><span style="color: #0000FF;">&</span><span style="color: #000000;">line</span><span style="color: #0000FF;">)</span> |
<span style="color: #000000;">output</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" "</span><span style="color: #0000FF;">&</span><span style="color: #000000;">line</span><span style="color: #0000FF;">)</span> |
||
Line 970: | Line 1,817: | ||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">PrntM</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">sw</span><span style="color: #0000FF;">)</span> |
<span style="color: #008080;">procedure</span> <span style="color: #000000;">PrntM</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">sw</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">integer</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span> |
<span style="color: #004080;">integer</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">s</span> |
||
<span style="color: #004080;">string</span> <span style="color: #000000;">g |
<span style="color: #004080;">string</span> <span style="color: #000000;">g</span> |
||
<span style="color: #004080;">atom</span> <span style="color: #000000;">q</span> |
<span style="color: #004080;">atom</span> <span style="color: #000000;">q</span> |
||
Line 994: | Line 1,841: | ||
<span style="color: #000000;">sw</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sw</span> <span style="color: #008080;">and</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">0</span> |
<span style="color: #000000;">sw</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sw</span> <span style="color: #008080;">and</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">0</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<span style="color: #008080;">if</span> <span style="color: #000000;">sw</span> <span style="color: #008080;">then</span> <span style="color: #000000;">g</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">""</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #000080;font-style:italic;">-- trivial |
<span style="color: #008080;">if</span> <span style="color: #000000;">sw</span> <span style="color: #008080;">then</span> <span style="color: #000000;">g</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">""</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #000080;font-style:italic;">-- trivial</span> |
||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">lensq</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #008000;">""</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
-- Hnf and solution</span> |
|||
<span style="color: #000000;">lensq</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">g</span> |
|||
<span style="color: #000080;font-style:italic;">-- Calculate lengths squared</span> |
|||
<span style="color: #000000;">output</span><span style="color: #0000FF;">(</span><span style="color: #000000;">print_row</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmts</span><span style="color: #0000FF;">)&</span><span style="color: #008080;">iif</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">k</span><span style="color: #0000FF;">?</span><span style="color: #000000;">g</span><span style="color: #0000FF;">:</span><span style="color: #008000;">""</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #000080;font-style:italic;">-- Null space with lengths squared</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
||
<span style="color: #000000;">line</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">print_row</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmts</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
||
<span style="color: #008080;">for</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
<span style="color: #008080;">for</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
||
<span style="color: #000000;">q</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">]</span> |
<span style="color: #000000;">q</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">s</span><span style="color: #0000FF;">]</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<span style="color: #000000;"> |
<span style="color: #000000;">lensq</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" (%d)"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">q</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
<span style="color: #000080;font-style:italic;">-- Hnf and solution and null space</span> |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">order</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)&</span><span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">rr</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">order</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #000000;">output</span><span style="color: #0000FF;">(</span><span style="color: #000000;">print_row</span><span style="color: #0000FF;">(</span><span style="color: #000000;">rr</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmts</span><span style="color: #0000FF;">)&</span><span style="color: #000000;">lensq</span><span style="color: #0000FF;">[</span><span style="color: #000000;">rr</span><span style="color: #0000FF;">])</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
Line 1,014: | Line 1,865: | ||
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">print_row</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmts</span><span style="color: #0000FF;">) |
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">print_row</span><span style="color: #0000FF;">(</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmts</span><span style="color: #0000FF;">))</span> |
||
<span style="color: #000080;font-style:italic;">-- (not particularly helpful:) |
|||
-- printf(1,", L:") |
|||
-- for s=1 to m+1 do |
|||
-- printf(1," %.20g",lambda[r,s]) |
|||
-- end for</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
Line 1,037: | Line 1,894: | ||
<span style="color: #000080;font-style:italic;">-- LLL reduce rows k</span> |
<span style="color: #000080;font-style:italic;">-- LLL reduce rows k</span> |
||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">Reduce</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
<span style="color: #008080;">procedure</span> <span style="color: #000000;">Reduce</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">integer</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sx</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">lk</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">q</span> |
|||
<span style="color: #000000;">col1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nx</span> |
<span style="color: #000000;">col1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nx</span> |
||
<span style="color: #000000;">col2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nx</span> |
<span style="color: #000000;">col2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nx</span> |
||
Line 1,054: | Line 1,909: | ||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
<span style="color: #004080;">atom</span> <span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
||
<span style="color: #008080;">if</span> <span style="color: #000000;">col1</span><span style="color: #0000FF;"><</span><span style="color: #000000;">nx</span> <span style="color: #008080;">then</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">col1</span><span style="color: #0000FF;"><</span><span style="color: #000000;">nx</span> <span style="color: #008080;">then</span> |
||
<span style="color: #008080;">if</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">col1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]<</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #000000;">Minus</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">col1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]<</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #000000;">Minus</span><span style="color: #0000FF;">(</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">col1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]/</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">col1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span> |
<span style="color: #000000;">q</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">col1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]/</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">,</span><span style="color: #000000;">col1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span> |
||
<span style="color: #008080;">else</span> |
<span style="color: #008080;">else</span> |
||
<span style="color: #000000;">lk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]</span> |
<span style="color: #004080;">atom</span> <span style="color: #000000;">lk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]</span> |
||
<span style="color: #008080;">if</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">)></span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">)></span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span> |
||
<span style="color: #000080;font-style:italic;">-- 2|lambda_kt| > d_t |
<span style="color: #000080;font-style:italic;">-- 2|lambda_kt| > d_t |
||
Line 1,068: | Line 1,923: | ||
<span style="color: #008080;">if</span> <span style="color: #000000;">q</span> <span style="color: #008080;">then</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">q</span> <span style="color: #008080;">then</span> |
||
<span style="color: #000000;">sx</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iif</span><span style="color: #0000FF;">(</span><span style="color: #000000;">col1</span><span style="color: #0000FF;">=</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">?</span><span style="color: #000000;">m</span><span style="color: #0000FF;">:</span><span style="color: #000000;">mn</span><span style="color: #0000FF;">)</span> |
<span style="color: #004080;">integer</span> <span style="color: #000000;">sx</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iif</span><span style="color: #0000FF;">(</span><span style="color: #000000;">col1</span><span style="color: #0000FF;">=</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">?</span><span style="color: #000000;">m</span><span style="color: #0000FF;">:</span><span style="color: #000000;">mn</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #000080;font-style:italic;">-- reduce row k</span> |
<span style="color: #000080;font-style:italic;">-- reduce row k</span> |
||
<span style="color: #008080;">for</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">sx</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
<span style="color: #008080;">for</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">sx</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
||
Line 1,082: | Line 1,937: | ||
<span style="color: #000080;font-style:italic;">-- exchange rows k and k-1</span> |
<span style="color: #000080;font-style:italic;">-- exchange rows k and k-1</span> |
||
<span style="color: #008080;">procedure</span> <span style="color: #000000;">Swop</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
<span style="color: #008080;">procedure</span> <span style="color: #000000;">Swop</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">integer</span> <span style="color: #000000;"> |
<span style="color: #004080;">integer</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">tm1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> |
||
<span style="color: #004080;">atom</span> <span style="color: #000000;">db</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lk</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lr</span> |
|||
<span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]}</span> |
<span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]}</span> |
||
<span style="color: # |
<span style="color: #004080;">object</span> <span style="color: #000000;">tmp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">tm1</span><span style="color: #0000FF;">]</span> |
||
<span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">tm1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">tm1</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">1</span><span style="color: #0000FF;">..</span><span style="color: #000000;">tm1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">tmp</span><span style="color: #0000FF;">;</span> <span style="color: #000000;">tmp</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
|||
<span style="color: #000080;font-style:italic;">-- update Gram coefficients |
<span style="color: #000080;font-style:italic;">-- update Gram coefficients |
||
-- columns k, k-1 for r > k</span> |
-- columns k, k-1 for r > k</span> |
||
<span style="color: #000000;">lk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]</span> |
<span style="color: #004080;">atom</span> <span style="color: #000000;">lk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">],</span> |
||
<span style="color: #000000;">db</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]+</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;"> |
<span style="color: #000000;">db</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]+</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> |
||
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span> |
||
<span style="color: #000000;">lr</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> |
<span style="color: #004080;">atom</span> <span style="color: #000000;">lr</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span> |
||
<span style="color: #000000;">dk1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
|||
<span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;"> |
<span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">dk1</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]-</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lr</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> |
||
<span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">db</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lr</span><span style="color: #0000FF;">+</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">])/</span><span style="color: #000000;">dk1</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;"> |
<span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">db</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span> |
||
Line 1,105: | Line 1,962: | ||
<span style="color: #000000;">problem_no</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
<span style="color: #000000;">problem_no</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span> |
||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"problem #%d\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">problem_no</span><span style="color: #0000FF;">)</span> |
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"problem #%d\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">problem_no</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #004080;">integer</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">tl</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">db</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lk</span> |
|||
<span style="color: #000000;"> |
<span style="color: #000000;">InputAb_or_c</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sw</span><span style="color: #0000FF;">)</span> |
||
<span style="color: #000080;font-style:italic;">-- augment Ab~ with column e_m</span> |
<span style="color: #000080;font-style:italic;">-- augment Ab~ with column e_m</span> |
||
Line 1,119: | Line 1,974: | ||
<span style="color: #008080;">if</span> <span style="color: #000000;">echo</span> <span style="color: #008080;">then</span> <span style="color: #000000;">PrntM</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">echo</span> <span style="color: #008080;">then</span> <span style="color: #000000;">PrntM</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
<span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> |
<span style="color: #004080;">integer</span> <span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> |
||
<span style="color: #000000;">tl</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">k</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">m</span> <span style="color: #008080;">do</span> |
<span style="color: #008080;">while</span> <span style="color: #000000;">k</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">m</span> <span style="color: #008080;">do</span> |
||
<span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> |
<span style="color: #004080;">integer</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> |
||
<span style="color: #000080;font-style:italic;">-- partial size reduction</span> |
<span style="color: #000080;font-style:italic;">-- partial size reduction</span> |
||
<span style="color: #000000;">Reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
<span style="color: #000000;">Reduce</span><span style="color: #0000FF;">(</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
||
Line 1,128: | Line 1,984: | ||
<span style="color: #008080;">if</span> <span style="color: #000000;">sw</span> <span style="color: #008080;">then</span> |
<span style="color: #008080;">if</span> <span style="color: #000000;">sw</span> <span style="color: #008080;">then</span> |
||
<span style="color: #000080;font-style:italic;">-- zero rows k-1, k</span> |
<span style="color: #000080;font-style:italic;">-- zero rows k-1, k</span> |
||
<span style="color: #000000;">lk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
<span style="color: #004080;">atom</span> <span style="color: #000000;">lk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lambda</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> |
||
<span style="color: #000080;font-style:italic;">-- Lovasz condition |
<span style="color: #000080;font-style:italic;">-- Lovasz condition |
||
-- Bk >= (alpha - mu_kt^2) * Bt</span> |
-- Bk >= (alpha - mu_kt^2) * Bt</span> |
||
<span style="color: #000000;">db</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]+</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lk</span> |
<span style="color: #004080;">atom</span> <span style="color: #000000;">db</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]+</span><span style="color: #000000;">lk</span><span style="color: #0000FF;">*</span><span style="color: #000000;">lk</span> |
||
<span style="color: #000080;font-style:italic;">-- not satisfied</span> |
<span style="color: #000080;font-style:italic;">-- not satisfied</span> |
||
<span style="color: #000000;">sw</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">db</span><span style="color: #0000FF;">*</span><span style="color: #000000;">ald</span><span style="color: #0000FF;"><</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">aln</span> |
<span style="color: #000000;">sw</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">db</span><span style="color: #0000FF;">*</span><span style="color: #000000;">ald</span><span style="color: #0000FF;"><</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">d</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">aln</span> |
||
Line 1,202: | Line 2,058: | ||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
||
<!--</lang>--> |
|||
<span style="color: #0000FF;">?</span><span style="color: #008000;">"done"</span> |
|||
Note that for problem 16 (HMM extended gcd (example 7.2)), the signs of |
|||
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">wait_key</span><span style="color: #0000FF;">()</span> |
|||
the (20) and (37) rows in the null space are flipped, somehow, otherwise |
|||
<!--</syntaxhighlight>--> |
|||
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). |
|||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
Line 1,214: | Line 2,067: | ||
{{libheader|Wren-ioutil}} |
{{libheader|Wren-ioutil}} |
||
{{libheader|Wren-complex}} |
{{libheader|Wren-complex}} |
||
{{libheader|Wren- |
{{libheader|Wren-iterate}} |
||
{{libheader|Wren-seq}} |
{{libheader|Wren-seq}} |
||
Results are the same as the FreeBASIC entry though I've just shown those for the first and last examples. |
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. |
A reasonably faithful translation though I've made some slight changes to the way data is input. |
||
< |
<syntaxhighlight lang="wren">import "./ioutil" for Input |
||
import "./complex" for Complex |
import "./complex" for Complex |
||
import "./ |
import "./iterate" for Stepped |
||
import "./seq" for Lst |
import "./seq" for Lst |
||
Line 1,553: | Line 2,406: | ||
main.call(sw) |
main.call(sw) |
||
System.print() |
System.print() |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
Latest revision as of 10:52, 26 November 2023
- task.
Implement the Havas-Majewski-Matthews LLL-based Hermite normal form algorithm for solving linear Diophantine systems.
- clarification.
The point of this task is not comprehending the puzzles, it is implementing the LLLHnf algorithm.
The method is the result of an experimental refinement process over many iterations, terse to the point of being impenetrable, best copied verbatim from a reliable source, and hence this task mostly concerns the direct translation of existing (cryptic) code.
You may regard the test set as just random input to validate your solution, no need to delve any deeper.
But to make the task a little nicer, and of course to demonstrate the power of the algorithm,
the examples aren't really random.
Many are classics, with online resources abound. Others are on Rosetta Code in a different guise;
some are copied from the HMM paper. Section headers like 'base cases' or 'polynomial coefficients'
should be self-explanatory.
The output is deliberately left somewhat 'raw', so there's plenty of room for implementation
dependent refinement. Also, to solve this task you're not obliged to click any wiki-links,
but please do if you want some background information.
- context.
Solving a system of linear Diophantine equations (i.e., solutions are required to be whole numbers) is a classic mathematical problem: given a coefficient matrix A (usually with less rows than columns) and vector b, the goal is to find a preferably small integer vector x that satisfies A·x = b.
How is this task different from Gaussian elimination? We now need a triangularization method that doesn't introduce fractions. The LLLHnf algorithm adapts Lenstra-Lenstra-Lovász lattice basis reduction to put the transpose of the input system into Hermite normal form, the integer analogue of the usual reduced row echelon form. In the process a unimodular transformation matrix is constructed from which (if either exists) minus the solution vector x and/or the null space basis for A can be immediately read off. These vectors are typically of small Euclidean length.
All details are found in the 1998 journal article. The relevant pseudocode (Algorithm 4) is on pages 131 and 129.
- use.
LLLHnf is a versatile algorithm that will
▹find integer solutions for (nonsquare) linear systems
▹put integer matrices into Hermite normal form
▹invert unimodular integer matrices
▹compute modular inverses
▹solve Bézout's identity
▹solve Chinese remaindering problems
(not necessarily with pairwise coprime moduli)
▹solve the famous Coconut puzzle
▹give hints for packing your knapsack
▹find solutions of small norm for extended gcd problems
▹find a univariate polynomial for a given (complex) algebraic constant
In all cases: provided a solution exists.
LLLHnf puts the transpose of A into Hermite normal form, so to compute Hnf(A) one inputs AT.
The program rejects fractional matrix coefficients. Users will scale applicable rows to suitably sized integers.
- limitation.
The output may be wrong if the (Gram-Schmidt) calculations produce numbers too big to be representable as integers (-larger than 253 if values are stored in 64-bit floats), though there is some tolerance.
- input.
(Piped to FreeBasic stdin in this format:)
m : #rows
n : #columns
m x n coefficient matrix A (augmented with m-vector b).
'a comment line starts with an apostrophe and contains no commas. 2 3 2 1 4| 17 -5 2 6|-13
To search for polynomial coefficients,
follow a comment line containing the tag "const" with:
m : #digits precision
n : max. poly degree
a + b : real and imaginary parts of a complex constant.
'constant sqrt(1 + sqrt(3)*i) 4 4 1.22474487 + .707106781
To quit:
0 0
- test vectors:
'five base cases 'no integral solution 2 2 2 0| 1 2 1| 2 'indeterminate 2 3 1 3 5 4 6 8 'singular square 3 3 1 7 4 2 8 5 3 9 6 'overdetermined 3 2 2 1| 2 6 5| 2 7 6| 2 'square 3 3 2 -3 4| 9 5 6 7| 3 8 9 10| 3 'modular inverse '(the smallest solution is negative ' add the modulus to make positive) 1 2 42 -2017| 1 'a classic Indian kuttaka problem 1 2 195 -221| 65 'Bachet de Méziriac "personnes en un banquet" '(add null space vector to make positive men) 2 3 1 1 1| 41 12 9 1|120 'Malm 2 4 1 1 1 1| 80 165 235 85 389| 16324 'from the Sunzi Suanjing 3 4 1 3 0 0| 2 1 0 5 0| 3 1 0 0 7| 2 2 3 17 7 0|-1 11 0 15|-2 'from the Shushu jiuzhang 8 9 1 130 0 0 0 0 0 0 0|-60 1 0 110 0 0 0 0 0 0|-30 1 0 0 120 0 0 0 0 0|-10 1 0 0 0 60 0 0 0 0|-10 1 0 0 0 0 25 0 0 0| 10 1 0 0 0 0 0 100 0 0| 10 1 0 0 0 0 0 0 50 0| 10 1 0 0 0 0 0 0 0 20| 10 '5 sailor coconut puzzle 6 7 1 -5 0 0 0 0 0| 1 0 4 -5 0 0 0 0| 1 0 0 4 -5 0 0 0| 1 0 0 0 4 -5 0 0| 1 0 0 0 0 4 -5 0| 1 0 0 0 0 0 4 -5| 0 'unbounded knapsack with slack 3 6 3000 1800 2500 1 0 0|54500 3 2 20 0 1 0| 250 25 15 2 0 0 1| 250 'subset sum problem 1 9 575 436 1586 1030 1921 569 721 1183 1570| 6665 'HMM extended gcd (example 7.2) 1 10 763836 1066557 113192 1785102 1470060 3077752 114793 3126753 1997137 2603018| 1 'Fibonacci segment F7...F14 (example 7.3) 1 8 13 21 34 55 89 144 233 377| 1 'compute the inverse of transpose(P) '(Fibonacci morphs into Lucas) 8 8 1 0 0 0 0 0 18 -7 1 1 0 0 0 0 -11 4 -1 1 1 0 0 0 7 -3 0 -1 1 1 0 0 -4 1 0 0 -1 1 1 0 3 -1 0 0 0 -1 1 1 -1 1 0 0 0 0 -1 1 1 0 0 0 0 0 0 -1 -1 0 'Hnf(A) with Aij = i^3 * j^2 + i + j (example 7.4) 10 10 3 11 31 69 131 223 351 521 739 1011 7 36 113 262 507 872 1381 2058 2927 4012 13 77 249 583 1133 1953 3097 4619 6573 9013 21 134 439 1032 2009 3466 5499 8204 11677 16014 31 207 683 1609 3135 5411 8587 12813 18239 25015 43 296 981 2314 4511 7788 12361 18446 26259 36016 57 401 1333 3147 6137 10597 16821 25103 35737 49017 73 522 1739 4108 8013 13838 21967 32784 46673 64018 91 659 2199 5197 10139 17511 27799 41489 59067 81019 111 812 2713 6414 12515 21616 34317 51218 72919 100020 'Gauss x*atan(1/239) + y*atan(1/57) + z*atan(1/18) = pi/4 '(fudge factor -1 to absorb round-off error ' ignore the corresponding vector entry x1) 1 4 -1 0041841 0175421 0554985| 7853982 'search for polynomial coefficients 'const sqrt(2) + i 4 4 1.41421356 + 1 'const 3^(1/3) + sqrt(2) 11 6 2.8564631326805 'some constant 12 9 -1.4172098692728 0 0
- output.
Transformation matrix P (on the left) and the Hermite normal form of
[A|b]T with messages:
inconsistent: the system is not solvable in integers.
-solution: is the negative of this particular P-vector
(the 1 in the last column of P is a marker and no part of the solution).
The null space vectors in P are followed by their parenthesized lengths squared.
- test results:
'five base cases 'no integral solution P | Hnf 0 -2 1 | 1 0 1 0 1 0 | 0 1 0 -1 -2 2 | -0 -0 2 inconsistent loop 8 'indeterminate P | Hnf 1 0 0 0 | 1 4 0 2 1 -1 -0 | -0 6 -0 0 0 0 1 | 0 0 1 1 -2 1 0 | 0 0 0 (6) loop 11 'singular square P | Hnf 1 0 0 0 | 1 2 3 0 3 -1 1 0 | 0 3 6 0 0 0 0 1 | 0 0 0 1 1 1 -2 -0 | -0 0 0 -0 (6) loop 12 'overdetermined P | Hnf 1 -1 0 | 1 1 1 0 -1 2 -0 | -0 4 5 -0 -2 2 1 | 0 0 0 1 -solution loop 7 'square P | Hnf 7 -1 -4 0 | 1 1 7 0 2 0 -1 -0 | -0 3 6 -0 15 -2 -9 -0 | -0 -0 12 -0 1 1 -2 1 | 0 0 0 1 -solution loop 15 'modular inverse '(the smallest solution is negative ' add the modulus to make positive) P | Hnf -48 -1 0 | 1 0 48 1 1 | 0 1 -solution 2017 42 0 | 0 0 (4070053) loop 7 'a classic Indian kuttaka problem P | Hnf 8 7 0 | 13 0 -6 -5 1 | 0 1 -solution -17 -15 0 | 0 0 (514) loop 8 'Bachet de Méziriac "personnes en un banquet" '(add null space vector to make positive men) P | Hnf -3 4 0 0 | 1 0 0 3 -4 1 0 | 0 1 0 3 -14 -30 1 | 0 0 1 -solution -8 11 -3 -0 | -0 0 -0 (194) loop 12 'Malm P | Hnf -1 2 1 -1 0 | 1 1 0 2 3 -3 -2 0 | 0 2 0 -17 -22 -25 -16 1 | 0 0 1 -solution -4 -8 7 5 0 | 0 0 0 (154) -15 8 7 0 0 | 0 0 0 (338) loop 23 'from the Sunzi Suanjing P | Hnf -35 12 7 5 0 | 1 0 0 0 21 -7 -4 -3 0 | 0 1 0 0 15 -5 -3 -2 -0 | -0 -0 1 -0 -23 7 4 3 1 | 0 0 0 1 -solution -105 35 21 15 0 | 0 0 0 0 (12916) loop 23 P | Hnf -30 73 22 0 | 1 0 0 -49 119 36 0 | 0 1 0 -23 56 17 1 | 0 0 1 -solution 105 -255 -77 0 | 0 0 0 (81979) loop 17 'from the Shushu jiuzhang P | Hnf 1 0 0 0 0 0 0 0 0 0 | 1 1 1 1 1 1 1 1 0 -7800 60 71 65 130 312 78 156 390 0 | 0 10 0 0 0 0 0 0 0 -35750 275 325 298 596 1430 358 715 1788 0 | 0 0 10 10 0 50 0 10 0 0 0 0 0 1 0 0 0 0 0 | 0 0 0 60 0 0 0 0 0 -34320 264 312 286 572 1373 344 687 1716 0 | 0 0 0 0 5 80 30 0 0 0 0 0 0 0 0 1 0 0 0 | 0 0 0 0 0 100 0 0 0 0 0 0 0 0 0 0 1 0 0 | 0 0 0 0 0 0 50 0 0 0 0 0 0 0 0 0 0 1 0 | 0 0 0 0 0 0 0 20 0 -3710 29 34 31 62 148 37 74 185 1 | 0 0 0 0 0 0 0 0 1 -solution 85800 -660 -780 -715 -1430 -3432 -858 -1716 -4290 -0 | -0 0 0 0 0 0 0 0 -0 (7399103669) loop 102 '5 sailor coconut puzzle P | Hnf 1 0 0 0 0 0 0 0 | 1 0 0 0 0 0 0 -3905 -781 -625 -500 -400 -320 -256 0 | 0 1 0 0 0 0 0 -975 -195 -156 -125 -100 -80 -64 0 | 0 0 1 0 0 0 0 -5125 -1025 -820 -656 -525 -420 -336 0 | 0 0 0 1 0 0 0 -2500 -500 -400 -320 -256 -205 -164 0 | 0 0 0 0 1 0 0 -3125 -625 -500 -400 -320 -256 -205 -0 | -0 -0 -0 -0 -0 1 -0 -3121 -624 -499 -399 -319 -255 -204 1 | 0 0 0 0 0 0 1 -solution 15625 3125 2500 2000 1600 1280 1024 0 | 0 0 0 0 0 0 0 (269403226) loop 71 'unbounded knapsack with slack P | Hnf 0 0 0 1 0 0 0 | 1 0 0 0 0 0 0 0 1 0 0 | -0 1 -0 -0 0 0 0 0 0 1 0 | 0 0 1 0 -6 -5 -11 0 -2 -3 1 | 0 0 0 1 -solution 3 -5 0 0 1 0 0 | 0 0 0 0 (35) -1 3 -1 100 17 -18 0 | 0 0 0 0 (10624) -4 15 -6 0 102 -113 0 | 0 0 0 0 (23450) loop 55 'subset sum problem P | Hnf 0 0 0 -1 0 1 -1 1 0 0 | 1 0 -1 0 -1 -1 0 0 -1 -1 -1 1 | 0 1 -solution -2 -1 1 0 0 0 0 0 0 0 | 0 0 (6) 1 -1 0 -2 1 0 0 0 0 0 | 0 0 (7) 0 -1 1 1 -1 0 -2 1 0 0 | 0 0 (9) -1 1 1 -1 0 -2 1 0 0 0 | 0 0 (9) 0 0 -1 1 1 -1 0 -2 1 0 | 0 0 (9) 0 -1 2 0 1 0 -1 -2 -1 0 | 0 0 (12) 1 1 0 -1 -2 0 1 0 2 0 | 0 0 (12) 2 -1 0 0 -1 -1 3 1 -1 0 | 0 0 (18) loop 66 'HMM extended gcd (example 7.2) P | Hnf -1 0 1 -3 1 3 3 -2 -2 2 0 | 1 0 1 0 -1 3 -1 -3 -3 2 2 -2 1 | 0 1 -solution 2 0 3 -1 0 0 0 1 1 -2 0 | 0 0 (20) 0 -1 2 2 -1 -1 3 -1 1 1 0 | 0 0 (23) -2 0 0 -1 3 -3 -1 2 1 0 0 | 0 0 (29) 0 3 2 3 2 -3 1 0 0 -1 0 | 0 0 (37) -2 2 2 0 -1 3 -3 -2 -1 0 0 | 0 0 (36) 2 2 -2 -5 -2 1 2 1 1 0 0 | 0 0 (48) 0 2 0 -2 -4 -1 -1 4 -1 0 0 | 0 0 (43) -3 3 -1 2 -2 1 0 1 4 -6 0 | 0 0 (81) 0 -2 1 -2 3 5 4 1 -5 -3 0 | 0 0 (94) loop 187 'Fibonacci segment F7...F14 (example 7.3) P | Hnf -7 4 -3 1 -1 1 0 0 0 | 1 0 7 -4 3 -1 1 -1 0 0 1 | 0 1 -solution -1 -1 1 0 0 0 0 0 0 | 0 0 (3) 0 0 -1 -1 1 0 0 0 0 | 0 0 (3) 0 0 0 0 -1 -1 1 0 0 | 0 0 (3) 0 -1 -1 1 0 0 0 0 0 | 0 0 (3) 0 0 0 -1 -1 1 0 0 0 | 0 0 (3) 0 0 0 0 0 -1 -1 1 0 | 0 0 (3) 18 -11 7 -4 3 -1 1 -1 0 | 0 0 (522) loop 45 'compute the inverse of transpose(P) '(Fibonacci morphs into Lucas) P | Hnf 2 1 5 3 0 -5 5 13 0 | 1 0 0 0 0 0 0 0 0 3 2 8 5 0 -8 8 21 0 | 0 1 0 0 0 0 0 0 0 4 3 13 8 0 -13 13 34 0 | 0 0 1 0 0 0 0 0 0 7 4 21 13 0 -21 21 55 0 | 0 0 0 1 0 0 0 0 0 11 7 33 21 0 -34 34 89 0 | 0 0 0 0 1 0 0 0 0 18 11 54 33 0 -55 55 144 0 | 0 0 0 0 0 1 0 0 0 29 18 87 54 -1 -89 89 233 0 | 0 0 0 0 0 0 1 0 0 47 29 141 87 -1 -145 144 377 0 | 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 | 0 0 0 0 0 0 0 0 1 loop 84 'Hnf(A) with Aij = i^3 * j^2 + i + j (example 7.4) P | Hnf -10 -8 -5 1 2 3 5 3 0 -4 0 | 1 0 7 22 45 76 115 162 217 280 0 -2 -1 0 1 -1 0 1 0 1 -1 0 | 0 1 4 9 16 25 36 49 64 81 0 -15 -11 -4 0 4 5 4 3 1 -5 0 | 0 0 12 36 72 120 180 252 336 432 0 0 0 0 0 0 0 0 0 0 0 1 | 0 0 0 0 0 0 0 0 0 0 1 -1 1 1 0 -2 1 0 0 0 0 0 | 0 0 0 0 0 0 0 0 0 0 0 (8) 0 -1 1 1 -1 1 -2 1 0 0 0 | 0 0 0 0 0 0 0 0 0 0 0 (10) -1 0 1 1 1 -2 0 -1 1 0 0 | 0 0 0 0 0 0 0 0 0 0 0 (10) -1 0 2 -1 1 -1 1 -1 -1 1 0 | 0 0 0 0 0 0 0 0 0 0 0 (12) 1 0 -1 0 -1 -1 1 2 0 -1 0 | 0 0 0 0 0 0 0 0 0 0 0 (10) -1 1 0 1 -1 0 0 1 -2 1 0 | 0 0 0 0 0 0 0 0 0 0 0 (10) -1 2 -1 -1 2 0 -2 1 0 0 0 | 0 0 0 0 0 0 0 0 0 0 0 (16) loop 99 'Gauss x*atan(1/239) + y*atan(1/57) + z*atan(1/18) = pi/4 '(fudge factor -1 to absorb round-off error ' ignore the corresponding vector entry x1) P | Hnf -1 -0 -0 -0 -0 | 1 -0 -1 5 -8 -12 1 | 0 1 -solution 23 42 -29 6 0 | 0 0 (3170) -71 33 -49 13 0 | 0 0 (8700) 18 -94 -63 27 0 | 0 0 (13858) loop 20 'search for polynomial coefficients 'const sqrt(2) + i 1.41421356 + 1*i P | Hnf 1 0 0 0 0 0 0 | 1 0 0 0 9 0 -2 0 1 0 | 0 1 0 0 0 0 0 0 0 1 | 0 0 1 -32 17 -10 0 2 0 0 | 0 0 0 (1417) 6 33 6 -48 35 -8 0 | 0 0 0 (4754) 18 -22 -56 5 31 -20 0 | 0 0 0 (5330) 102 151 -71 91 -78 36 0 | 0 0 0 (53907) loop 53 'const 3^(1/3) + sqrt(2) 2.8564631326805 P | Hnf 1 0 0 0 0 0 0 0 0 | 1 0 0 0 0 0 0 0 0 0 0 1 | 0 0 1 5 1 -36 12 -6 -6 0 1 0 | 0 0 0 (1539) -27 14 0 15 -31 -23 -26 13 0 | 0 0 0 (3485) 33 20 21 22 -13 -43 1 5 0 | 0 0 0 (4458) -18 10 -29 -49 -12 -25 -13 9 0 | 0 0 0 (4685) 53 -14 15 13 -23 35 -13 1 0 | 0 0 0 (5323) 35 78 -17 -17 8 7 0 -1 0 | 0 0 0 (8001) 10 -5 -28 33 110 11 -84 23 0 | 0 0 0 (21804) loop 125 'some constant -1.4172098692728 P | Hnf 1 0 0 0 0 0 0 0 0 0 0 0 | 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 | 0 0 1 0 9 0 0 -7 0 -5 0 3 0 1 0 | 0 0 0 (165) -2 9 -2 -9 4 8 2 -2 -12 -2 4 0 | 0 0 0 (422) -4 -1 9 2 -5 -6 -2 6 6 -11 -9 0 | 0 0 0 (441) -5 11 6 3 4 -3 5 -12 -6 3 -1 0 | 0 0 0 (431) -9 1 9 5 4 1 -12 -8 11 1 -5 0 | 0 0 0 (560) -6 -2 2 17 -11 -4 3 1 -4 -5 0 0 | 0 0 0 (521) -1 1 0 5 3 8 12 6 4 9 5 0 | 0 0 0 (402) -7 3 0 -13 -3 -1 4 9 -4 8 9 0 | 0 0 0 (495) -9 11 7 2 11 -7 -7 3 -8 -2 3 0 | 0 0 0 (560) -16 11 -13 -4 8 -4 -4 3 4 1 0 0 | 0 0 0 (684) loop 360
To save server space, better not repeat full test results in your post. Once your solution covers all cases, a selected example will suffice.
C
// 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
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
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