黄金分割法

#********************************************/
# 黄金分割法による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
}

#*****************************************/
# 黄金分割法(関数の最小値)               */
#      a,b : 初期区間 a < b              */
#      eps : 許容誤差                    */
#      val : 間数値                      */
#      ind : 計算状況                    */
#              >= 0 : 正常終了(収束回数) */
#              = -1 : 収束せず           */
#      max : 最大試行回数                */
#      fun : 関数値を計算する関数の名前  */
#      return : 結果                     */
#*****************************************/
def gold(a, b, eps, val, ind, max, &fun)

	x      = 0.0
	tau    = (Math.sqrt(5.0) - 1.0) / 2.0
	count  = 0
	ind[0] = -1
	x1     = b - tau * (b - a)
	x2     = a + tau * (b - a)
	f1     = fun.call(x1)
	f2     = fun.call(x2)

	while count < max && ind[0] < 0
		count += 1
		if f2 > f1
			if (b-a).abs() < eps && (b-a).abs() < eps*b.abs()
				ind[0] = 0
				x      = x1
				val[0] = f1
			else
				b  = x2
				x2 = x1
				x1 = a + (1.0 - tau) * (b - a)
				f2 = f1
				f1 = fun.call(x1)
			end
		else
			if (b-a).abs() < eps && (b-a).abs() < eps*b.abs()
				ind[0] = 0
				x      = x2
				val[0] = f2
				f1     = f2
			else
				a  = x1
				x1 = x2
				x2 = b - (1.0 - tau) * (b - a)
				f1 = f2
				f2 = fun.call(x2)
			end
		end
	end

	if ind[0] == 0
		ind[0] = count
		fa     = fun.call(a)
		fb     = fun.call(b)
		if fb < fa
			a  = b
			fa = fb
		end
		if fa < f1
			x      = a
			val[0] = fa
		end
	end

	return x
end

ind = Array.new(1)
val = Array.new(1)
a   = -2.0
b   = -1.0
eps = 1.0e-10
max = 100

x = gold(a, b, eps, val, ind, max, &snx)

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