Diophantine linear system solving: Difference between revisions
Diophantine linear system solving (view source)
Revision as of 10:52, 26 November 2023
, 5 months ago→{{header|Wren}}: Changed to Wren S/H
No edit summary |
m (→{{header|Wren}}: Changed to Wren S/H) |
||
| (7 intermediate revisions by 4 users not shown) | |||
Line 532:
__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}}==
<
'Subject: Solution of an m x n linear Diophantine system
' A*x = b using LLL reduction.
Line 868 ⟶ 1,315:
print "loop"; tl
End Sub
</syntaxhighlight>
=={{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}}==
You can run this online [http://phix.x10.mx/p2js/LLLHnf.htm here].
<!--<
<span style="color: #000080;font-style:italic;">--
-- demo\rosetta\Diophantine_linear_system_solving.exw
Line 1,223 ⟶ 2,061:
<span style="color: #0000FF;">?</span><span style="color: #008000;">"done"</span>
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">wait_key</span><span style="color: #0000FF;">()</span>
<!--</
=={{header|Wren}}==
Line 1,229 ⟶ 2,067:
{{libheader|Wren-ioutil}}
{{libheader|Wren-complex}}
{{libheader|Wren-
{{libheader|Wren-seq}}
Results are the same as the FreeBASIC entry though I've just shown those for the first and last examples.
A reasonably faithful translation though I've made some slight changes to the way data is input.
<
import "./complex" for Complex
import "./
import "./seq" for Lst
Line 1,568 ⟶ 2,406:
main.call(sw)
System.print()
}</
{{out}}
Line 1,633 ⟶ 2,471:
loop 360
</pre>
| |||