最小二乗法

############################
# 最小二乗法（多項式近似）
#      coded by Y.Suganuma
############################

############################################
# 線形連立方程式を解く（逆行列を求める）
#      w : 方程式の左辺及び右辺
#      n : 方程式の数
#      m : 方程式の右辺の列の数
#      eps : 逆行列の存在を判定する規準
#      return : =0 : 正常
#               =1 : 逆行列が存在しない
#      coded by Y.Suganuma
############################################

def gauss(w, n, m, eps)

	nm  = n + m;
	ind = 0

	for i1 in 0 ... n

		y1 = 0.0
		m1 = i1 + 1
		m2 = 0
					# ピボット要素の選択
		for i2 in i1 ... n
			y2 = w[i2][i1].abs()
			if y1 < y2
				y1 = y2
				m2 = i2
			end
		end
					# 逆行列が存在しない
		if y1 < eps
			ind = 1
			break
					# 逆行列が存在する
		else   # 行の入れ替え
			for i2 in i1 ... nm
				y1        = w[i1][i2]
				w[i1][i2] = w[m2][i2]
				w[m2][i2] = y1
			end
						# 掃き出し操作
			y1 = 1.0 / w[i1][i1]

			for i2 in m1 ... nm
				w[i1][i2] *= y1
			end

			for i2 in 0 ... n
				if i2 != i1
					for i3 in m1 ... nm
						w[i2][i3] -= (w[i2][i1] * w[i1][i3])
					end
				end
			end
		end
	end

	return ind
end

#########################################
# 最小二乗法
#      m : 多項式の次数
#      n : データの数
#      x,y : データ
#      return : 多項式の係数（高次から）
#               エラーの場合はNULLを返す
#      coded by Y.Suganuma
#########################################

def least(m, n, x, y)

	m += 1
	z  = Array.new(m)
	w  = Array.new(m)
	for i1 in 0 ... m
		w[i1] = Array.new(m+1)
	end
	a = Array.new(n)
	for i1 in 0 ... n
		a[i1] = Array.new(m)
	end

	for i1 in 0 ... n
		a[i1][m-2] = x[i1]
		a[i1][m-1] = 1.0
		x1         = a[i1][m-2]
		x2         = x1
		i2         = m - 3
		while i2 > -1
			x2        *= x1
			a[i1][i2]  = x2
			i2        -= 1
		end
	end

	for i1 in 0 ... m
		for i2 in 0 ... m
			w[i1][i2] = 0.0
			for i3 in 0 ... n
				w[i1][i2] += a[i3][i1] * a[i3][i2]
			end
		end
	end

	for i1 in 0 ... m
		w[i1][m] = 0.0
		for i2 in 0 ... n
			w[i1][m] += a[i2][i1] * y[i2]
		end
	end

	sw = gauss(w, m, 1, 1.0e-10)

	if sw == 0
		for i1 in 0 ... m
			z[i1] = w[i1][m]
		end
	else
		z = Array.new(0)
	end

	return z
end

s = gets().split(" ")
m = Integer(s[0])   # 多項式の次数
n = Integer(s[1])   # データの数

x = Array.new(n)
y = Array.new(n)

for i1 in 0 ... n   # データ
	s     = gets().split()
	x[i1] = Float(s[0])
	y[i1] = Float(s[1])
end

z = least(m, n, x, y)

if z.size  > 0
	print("結果\n")
	for i1 in 0 ... m+1
		print("   " + String(m-i1) + " 次の係数 " + String(z[i1]) + "\n")
	end
else
	print("***error  逆行列を求めることができませんでした\n")
end

=begin
------------データ例１（コメント部分を除いて下さい）--------
1 10   # 多項式の次数とデータの数
0  2.450   # 以下，(x, y)
1  2.615
2  3.276
3  3.294
4  3.778
5  4.009
6  3.920
7  4.267
8  4.805
9  5.656

-------------------------データ例２-------------------------
2 6
10  28.2
15  47.0
20  44.4
26  32.8
32  20.8
40   0.8
=end