import std/random import std/os import std/parseutils from mmath/rational import Rational, newRational, `$`, zero, abs, `<`, `*`, `-`, `/=`, `+`, `+=` from mmath/hmatrix import newHMatrix, clone, lup, lu_solve, `*`, `-`, `$` from mmath/hvector import newHvector, `-`, norm from mmath/bigint import BigInt, newBigInt, `of`, one, zero, abs, `*`, `mod`, `div`, `-`, `+`, `$`, `==` let argv = commandLineParams() let size = block: var n : int = 0 if argv.len() > 0: discard parseInt(argv[0], n) else: n = 3 n proc nextInt(r : var Rand, min : int, max : int) : int = min + r.rand(max - min) var rand = initRand(101325) let mtx = block: let valueGenerator = proc(i : int, j : int): Rational[BigInt] = return newRational[BigInt](BigInt.of(rand.nextInt(-1000 * size, 1000 * size).int64), BigInt.of((1000 * size).int64)) newHMatrix[Rational[BigInt]](size, size, valueGenerator) var lu = mtx.clone() let pivot = lu.lup() let b = block: let generator = proc(i : int) : Rational[BigInt] = newRational(BigInt.of(rand.nextInt(0, size)), BigInt.of(size)) newHVector[Rational[BigInt]](size, generator) let x = lu.lu_solve(b, pivot) let error = mtx * x - b echo $norm(error)