Diophantine linear system solving: Difference between revisions
Created Nim solution.
m (→{{header|Wren}}: Wren-trait -> Wren-iterate) |
(Created Nim solution.) |
||
Line 1,316:
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}}==
| |||