#!/bin/lua
-- simultaneous linear equation
-- n equations
-- let x be an array of n variables x[1],x[2],...,x[n]
-- let a[k] be array of n constants a[k][1],a[k][2],...,a[k][n]
-- let (a[k],x) be an inner product of a[k] and x
-- our goal is to solve:
--   (a[k],x) = b[k]  -- k=1,2,...,n
-- where b is array of n constants b[1],b[2],...,b[n]
-- givening a[k] and b

abs=math.abs
fmt=string.format
gen=require("gen")
tprint=gen.tprint
num=require("num")
frac=require("frac")
setmt=frac.setmt
dup=frac.dup
Fr=frac.Fr
fr_abs=frac.fr_abs

function iprod(a,b)
	local sum
	assert(#a==#b)
	sum = 0
	for k=1,#a do
		sum = sum + a[k]*b[k]
	end
	return sum
end

function i2frac(...)
	local a={...}
	local f={}
	for i=1,#a do f[i]=Fr.new(a[i],1) end
	return f
end

-- simultaneous linear equation
local function SLE(a,b)
	a = dup(a)
	b = dup(b)
	assert(a[1][1] and a[1][1] ~= 0)
	n = #b
	i = 1
	::Loop::
	-- making upper triangular matrix
	aii = a[i][i]
	for j=i,n do a[i][j] = a[i][j]/aii end
	b[i] = b[i]/aii
	for k=i+1,n do
		aki = a[k][i]
		for j=i,n do a[k][j] = a[k][j] - aki*a[i][j] end
		b[k] = b[k] - aki*b[i]
	end
	i = i+1
	if i<n then goto Loop end

	-- confirmation: got upper triangular matri or not?
	print("confirmation")
	for k=1,n do
		for j=1,n do printf("	%s",a[k][j]) end
		print("",b[k])
	end
	
	-- get the solution
	x = {}
	for k=n,1,-1 do
		ax = 0
		for j=k+1,n do
			ax = ax + a[k][j]*x[j]
		end
		x[k] = (b[k] - ax)/a[k][k]
	end

	return x
end

local function test1()
	--[[ Asano p.82 11
	--solve in real number field
	--	5x + 6y - 7z = 8
	--	2x - 3y + 5z = 1
	--	6x + 7y - 2z = 3
	--]]
	a = {}
	a[1] = {5,6,-7}
	a[2] = {2,-3,5}
	a[3] = {6,7,-2}
	b = {8,1,3}
	x = SLE(a,b)
	tprint("x",x)
	--[[ OK
	x
	1	1.4666666666667
	2	-1.1333333333333
	3	-1.0666666666667
	--]]
end

local function test2()
	--[[
	-- solve in rational number field
	--	5x + 6y - 7z = 8
	--	2x - 3y + 5z = 1
	--	6x + 7y - 2z = 3
	--]]
	a = {}
	a[1] = i2frac(5,6,-7)
	a[2] = i2frac(2,-3,5)
	a[3] = i2frac(6,7,-2)
	b = i2frac(8,1,3)
	x = SLE(a,b)
	tprint("x",x)
	--[[  OK
	x
	1	22/15
	2	-17/15
	3	-16/15
	--]]
end

test2()