補間法（３次スプライン関数）

#*******************************/
# 補間法（３次スプライン関数） */
#      coded by Y.Suganuma     */
#*******************************/

#*********************************/
# ３次スプライン関数による補間   */
#      n : 区間の数              */
#      x, y : 点の座標           */
#      x1 : 補間値を求める値     */
#      h, b, d, g, u, r : 作業域 */
#      return : 補間値           */
#      coded by Y.Suganuma       */
#*********************************/
def spline(n, x, y, x1, h, b, d, g, u, r)
					# 区間の決定
	i = -1
	for i1 in 1 ... n
		if x1 < x[i1]
			i = i1 - 1
		end
		if i >= 0
			break
		end
	end
	if i < 0
		i = n - 1
	end
					# ステップ１
	for i1 in 0 ... n
		h[i1] = x[i1+1] - x[i1]
	end
	for i1 in 1 ... n
		b[i1] = 2.0 * (h[i1] + h[i1-1])
		d[i1] = 3.0 * ((y[i1+1] - y[i1]) / h[i1] - (y[i1] - y[i1-1]) / h[i1-1])
	end
					# ステップ２
	g[1] = h[1] / b[1]
	for i1 in 2 ... n-1
		g[i1] = h[i1] / (b[i1] - h[i1-1] * g[i1-1])
	end
	u[1] = d[1] / b[1]
	for i1 in 2 ... n
		u[i1] = (d[i1] - h[i1-1] * u[i1-1]) / (b[i1] - h[i1-1] * g[i1-1])
	end
					# ステップ３
	k      = (i > 1) ? i : 1
	r[0]   = 0.0
	r[n]   = 0.0
	r[n-1] = u[n-1]
	i1     = n - 2
	while i1 >= k
		r[i1] = u[i1] - g[i1] * r[i1+1]
		i1 -= 1;
	end
					# ステップ４
	xx = x1 - x[i]
	qi = (y[i+1] - y[i]) / h[i] - h[i] * (r[i+1] + 2.0 * r[i]) / 3.0
	si = (r[i+1] - r[i]) / (3.0 * h[i])
	y1 = y[i] + xx * (qi + xx * (r[i]  + si * xx))

	return y1
end

					# 設定
n = 5
h = Array.new(n)
b = Array.new(n)
d = Array.new(n)
g = Array.new(n)
u = Array.new(n)
r = Array.new(n+1)
x = Array[0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
y = Array[0.0, 1.0, 4.0, 9.0, 16.0, 25.0]
				# 実行と出力
x1 = 0.0
sp = 0.1
while x1 <= 5.01
	y1 = spline(n, x, y, x1, h, b, d, g, u, r)
	printf("%f %f\n", x1, y1)
	x1 += sp
end