多項式近似(関数の最小値)

#********************************************/
# 多項式近似によるy=x^4+3x^3+2x^2+1の最小値 */
#      coded by Y.Suganuma                  */
#********************************************/

#***************/
# 関数値の計算 */
#***************/
snx = Proc.new { |x|
	x * x * x * x + 3.0 * x * x * x + 2.0 * x * x + 1.0
}

#*****************************************/
# 多項式近似(関数の最小値)               */
#      x0  : 初期値                      */
#      d0  : 初期ステップ                */
#      eps : 許容誤差                    */
#      val : 関数値                      */
#      ind : 計算状況                    */
#              >= 0 : 正常終了(収束回数) */
#              = -1 : 収束せず           */
#      max : 最大試行回数                */
#      fun : 関数値を計算する関数の名前  */
#      return : 結果                     */
#*****************************************/
def approx(x0, d0, eps, val, ind, max, &fun)

	f      = Array.new(4)
	x      = Array.new(4)
	xx     = 0.0
	d      = d0
	x[1]   = x0
	f[1]   = fun.call(x0)
	k      = 0
	count  = 0
	ind[0] = -1

	while count < max && ind[0] < 0
		x[3] = x[1] + d
		f[3] = fun.call(x[3])
		while k < max && f[3] <= f[1]
			k   += 1
			d   *= 2.0
			x[0] = x[1]
			f[0] = f[1]
			x[1] = x[3]
			f[1] = f[3]
			x[3] = x[1] + d
			f[3] = fun.call(x[3])
		end
					# 初期値が不適当
		if k >= max
			count = max
		else
					# ３点の選択
			sw = 0
			if k > 0
				x[2] = x[3] - 0.5 * d
				f[2] = fun.call(x[2])
				min  = -1
				for i1 in 0 ... 4
					if min < 0 || f[i1] < f[min]
						min = i1
					end
				end
				if min >= 2
					for i1 in 0 ... 3
						x[i1] = x[i1+1]
						f[i1] = f[i1+1]
					end
				end
				sw = 1
			else
				x[0] = x[1] - d0
				f[0] = fun.call(x[0])
				if f[0] > f[1]
					x[2] = x[3]
					f[2] = f[3]
					sw   = 1
				else
					x[1] = x[0]
					f[1] = f[0]
					d0   = -d0
					d    = 2.0 * d0
					k    = 1
				end
			end
					# 収束？
			if sw > 0
				count += 1
				dl     = 0.5 * d * (f[2] - f[0]) / (f[0] - 2.0 * f[1] + f[2])
				xx     = x[1] - dl
				val[0] = fun.call(xx)
				if dl.abs() < eps
					ind[0] = count
				else
					k  = 0
					d0 = 0.5 * d
					d  = d0
					if val[0] < f[1]
						x[1] = xx
						f[1] = val[0]
					end
				end
			end
		end
	end

	return xx
end

ind = Array.new(1)
val = Array.new(1)
x   = -3.0
d   = 0.1
eps = 1.0e-7
max = 100

x = approx(x, d, eps, val, ind, max, &snx)

printf("x %f val %f ind %d\n", x, val[0], ind[0])