最適化（ Newton 法）

　　プログラムは，Newton 法を使用して，非線形関数の最小値を求めるためのものです（プログラムの使用方法）．

/******************************/
/* Newton法による最小値の計算 */
/*      coded by Y.Suganuma   */
/******************************/
#include <stdio.h>

void dsnx1(double *, double *);
double snx1(double, double *, double *);
int hesse1(double *, double **, double);
void dsnx2(double *, double *);
double snx2(double, double *, double *);
int hesse2(double *, double **, double);
void dsnx3(double *, double *);
double snx3(double, double *, double *);
int hesse3(double *, double **, double);
int Newton(int, int, int, double, double, double *, double *, double *, double **,
          double (*)(double, double *, double *),
          void (*)(double *, double *),
          int (*)(double *, double **, double));

int main()
{
	double eps, **H, step, *x, *dx, y;
	int fun, i1, max, n, opt_1, sw = 0;
					// データの入力
	scanf("%*s %d %*s %d %*s %d %*s %d", &fun, &n, &max, &opt_1);
	scanf("%*s %lf %*s %lf", &eps, &step);
	x  = new double [n];
	dx = new double [n];
	H  = new double * [n];
	scanf("%*s");
	for (i1 = 0; i1 < n; i1++) {
		scanf("%lf", &x[i1]);
		H[i1] = new double [2*n];
	}
					// 実行
	switch (fun) {
		case 1:
			sw = Newton(opt_1, max, n, eps, step, &y, x, dx, H, snx1, dsnx1, hesse1);
			break;
		case 2:
			sw = Newton(opt_1, max, n, eps, step, &y, x, dx, H, snx2, dsnx2, hesse2);
			break;
		case 3:
			sw = Newton(opt_1, max, n, eps, step, &y, x, dx, H, snx3, dsnx3, hesse3);
			break;
	}
					// 結果の出力
	if (sw < 0) {
		printf("   収束しませんでした！");
		switch (sw) {
			case -1:
				printf("（収束回数）\n");
				break;
			case -2:
				printf("（１次元最適化の区間）\n");
				break;
			case -3:
				printf("（黄金分割法）\n");
				break;
		}
	}
	else {
		printf("   結果＝");
		for (i1 = 0; i1 < n; i1++)
			printf("%f ", x[i1]);
		printf(" 最小値＝%f  回数＝%d\n", y, sw);
	}

	return 0;
}

/********************************************************/
/* Newton法                                             */
/*      opt_1 : =0 : １次元最適化を行わない             */
/*              =1 : １次元最適化を行う                 */
/*      max : 最大繰り返し回数                          */
/*      n : 次元                                        */
/*      eps : 収束判定条件                              */
/*      step : きざみ幅                                 */
/*      y : 最小値                                      */
/*      x : x(初期値と答え)                             */
/*      dx : 関数の微分値                               */
/*      H : Hesse行列の逆行列                           */
/*      snx : 関数値を計算する関数名                    */
/*      dsnx : 関数の微分を計算する関数名(符号を変える) */
/*      hesse : Hesse行列の逆行列を計算する関数名       */
/*      return : >=0 : 正常終了(収束回数)               */
/*               =-1 : 収束せず                         */
/*               =-2 : １次元最適化の区間が求まらない   */
/*               =-3 : 黄金分割法が失敗                 */
/********************************************************/
#include <math.h>

double gold(double, double, double, double *, int *, int, double *, double *,
            double (*)(double, double *, double *));

int Newton(int opt_1, int max, int n, double eps, double step, double *y,
          double *x, double *dx, double **H, double (*snx)(double, double *, double *), 
          void (*dsnx)(double *, double *), int (*hesse)(double *, double **, double))
{
	double f1, f2, k, sp, y1, y2;
	double *wk = new double [n];
	int count = 0, i1, i2, sw = 0, sw1;

	y1 = snx(0.0, x, dx);

	while (count < max && sw == 0) {
					// 傾きの計算
		dsnx(x, wk);
					// Hesse行列の逆行列の計算
		sw1 = hesse(x, H, eps);
					// 収束していない
		if (sw1 == 0) {
						// 方向の計算
			count++;
			for (i1 = 0; i1 < n; i1++) {
				dx[i1] = 0.0;
				for (i2 = 0; i2 < n; i2++)
					dx[i1] += H[i1][n+i2] * wk[i2];
			}
						// １次元最適化を行わない
			if (opt_1 == 0) {
							// 新しい点
				for (i1 = 0; i1 < n; i1++)
					x[i1] += dx[i1];
							// 新しい関数値
				y2 = snx(0.0, x, dx);
								// 関数値の変化が大きい
				if (fabs(y2-y1) > eps)
					y1 = y2;
								// 収束（関数値の変化＜ｅｐｓ）
				else {
					sw = count;
					*y = y2;
				}
			}
						// １次元最適化を行う
			else {
							// 区間を決める
				sw1 = 0;
				f1  = y1;
				sp  = step;
				f2  = snx(sp, x, dx);
				if (f2 > f1)
					sp = -step;
				for (i1 = 0; i1 < max && sw1 == 0; i1++) {
					f2 = snx(sp, x, dx);
					if (f2 > f1)
						sw1 = 1;
					else {
						sp *= 2.0;
						f1  = f2;
					}
				}
							// 区間が求まらない
				if (sw1 == 0)
					sw = -2;
							// 区間が求まった
				else {
								// 黄金分割法
					k = gold(0.0, sp, eps, &y2, &sw1, max, x, dx, fn);
								// 黄金分割法が失敗
					if (sw1 < 0)
						sw = -3;
								// 黄金分割法が成功
					else {
									// 新しい点
						for (i1 = 0; i1 < n; i1++)
							x[i1] += k * dx[i1];
										// 関数値の変化が大きい
						if (fabs(y1-y2) > eps)
							y1 = y2;
										// 収束（関数値の変化＜ｅｐｓ）
						else {
							sw = count;
							*y = y2;
						}
					}
				}
			}
		}
					// 収束（傾き＜ｅｐｓ）
		else {
			sw = count;
			*y = y1;
		}
	}

	if (sw == 0)
		sw = -1;

	delete [] wk;

	return sw;
}

/****************************************************************/
/* 黄金分割法(与えられた方向での最小値)                         */
/*      a,b : 初期区間 a < b                                    */
/*      eps : 許容誤差                                          */
/*      val : 関数値                                            */
/*      ind : 計算状況                                          */
/*              >= 0 : 正常終了(収束回数)                       */
/*              = -1 : 収束せず                                 */
/*      max : 最大試行回数                                      */
/*      w : 位置                                                */
/*      dw : 傾きの成分                                         */
/*      snx : 関数値を計算する関数の名前                        */
/*      return : 結果(w＋y＊dwのy)                              */
/****************************************************************/
#include <math.h>

double gold(double a, double b, double eps, double *val, int *ind, int max,
            double *w, double *dw, double (*snx)(double, double *, double *))
{
	double f1, f2, fa, fb, tau = (sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
	int count = 0;
					// 初期設定
	*ind  = -1;
	x1    = b - tau * (b - a);
	x2    = a + tau * (b - a);
	f1    = snx(x1, w, dw);
	f2    = snx(x2, w, dw);
					// 計算
	while (count < max && *ind < 0) {
		count += 1;
		if (f2 > f1) {
			if (fabs(b-a) < eps) {
				*ind = 0;
				x    = x1;
				*val = f1;
			}
			else {
				b  = x2;
				x2 = x1;
				x1 = a + (1.0 - tau) * (b - a);
				f2 = f1;
				f1 = snx(x1, w, dw);
			}
		}
		else {
			if (fabs(b-a) < eps) {
				*ind = 0;
				x    = x2;
				*val = f2;
				f1   = f2;
			}
			else {
				a  = x1;
				x1 = x2;
				x2 = b - (1.0 - tau) * (b - a);
				f1 = f2;
				f2 = snx(x2, w, dw);
			}
		}
	}
					// 収束した場合の処理
	if (*ind == 0) {
		*ind = count;
		fa   = snx(a, w, dw);
		fb   = snx(b, w, dw);
		if (fb < fa) {
			a  = b;
			fa = fb;
		}
		if (fa < f1) {
			x    = a;
			*val = fa;
		}
	}

	return x;
}

/*********************************************/
/* 与えられた点xから，dx方向へk*dxだけ進んだ */
/* 点における関数値を計算する                */
/*********************************************/
					// 関数１
double snx1(double k, double *x, double *dx)
{
	double x1, y1, y, w[2];

	w[0] = x[0] + k * dx[0];
	w[1] = x[1] + k * dx[1];
	x1   = w[0] - 1.0;
	y1   = w[1] - 2.0;
	y    = x1 * x1 + y1 * y1;

	return y;
}
					// 関数２
double snx2(double k, double *x, double *dx)
{
	double x1, y1, y, w[2];

	w[0] = x[0] + k * dx[0];
	w[1] = x[1] + k * dx[1];
	x1   = w[1] - w[0] * w[0];
	y1   = 1.0 - w[0];
	y    = 100.0 * x1 * x1 + y1 * y1;

	return y;
}
					// 関数３
double snx3(double k, double *x, double *dx)
{
	double x1, y1, z1, y, w[2];

	w[0] = x[0] + k * dx[0];
	w[1] = x[1] + k * dx[1];
	x1   = 1.5 - w[0] * (1.0 - w[1]);
	y1   = 2.25 - w[0] * (1.0 - w[1] * w[1]);
	z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
	y    = x1 * x1 + y1 * y1 + z1 * z1;

	return y;
}

/********************/
/* 微係数を計算する */
/********************/
					// 関数１
void dsnx1(double *x, double *dx)
{
	dx[0] = -2.0 * (x[0] - 1.0);
	dx[1] = -2.0 * (x[1] - 2.0);
}
					// 関数２
void dsnx2(double *x, double *dx)
{
	dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0]);
	dx[1] = -200.0 * (x[1] - x[0] * x[0]);
}
					// 関数３
void dsnx3(double *x, double *dx)
{
	dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
            2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
            2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
	dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
            4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
            6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
}

/*******************************/
/* Hesse行列の逆行列を計算する */
/*******************************/
int Gauss(double **, int, int, double);
					// 関数１
int hesse1(double *x, double **H, double eps)
{
	int sw;

	H[0][0] = 2.0;
	H[0][1] = 0.0;
	H[1][0] = 0.0;
	H[1][1] = 2.0;
	H[0][2] = 1.0;
	H[0][3] = 0.0;
	H[1][2] = 0.0;
	H[1][3] = 1.0;

	sw = Gauss(H, 2, 2, eps);

	return sw;
}
					// 関数２
int hesse2(double *x, double **H, double eps)
{
	int sw;

	H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0;
	H[0][1] = -400.0 * x[0];
	H[1][0] = H[0][1];
	H[1][1] = 200.0;
	H[0][2] = 1.0;
	H[0][3] = 0.0;
	H[1][2] = 0.0;
	H[1][3] = 1.0;

	sw = Gauss(H, 2, 2, eps);

	return sw;
}
					// 関数３
int hesse3(double *x, double **H, double eps)
{
	double x1, x2, x3;
	int sw;

	x1      = 1.0 - x[1];
	x2      = 1.0 - x[1] * x[1];
	x3      = 1.0 - x[1] * x[1] * x[1];
	H[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3;
	H[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
              4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
              6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3;
	H[1][0] = H[0][1];
	H[1][1] = 2.0 * x[0] * x[0] +
              4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
              12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
              18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1];
	H[0][2] = 1.0;
	H[0][3] = 0.0;
	H[1][2] = 0.0;
	H[1][3] = 1.0;

	sw = Gauss(H, 2, 2, eps);

	return sw;
}

/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める）              */
/*      w : 方程式の左辺及び右辺                       */
/*      n : 方程式の数                                 */
/*      m : 方程式の右辺の列の数                       */
/*      eps : 正則性を判定する規準                     */
/*      return : =0 : 正常                             */
/*               =1 : 逆行列が存在しない               */
/*******************************************************/
#include <math.h>

int Gauss(double **w, int n, int m, double eps)
{
	double y1, y2;
	int ind = 0, nm, m1, m2, i1, i2, i3;

	nm = n + m;

	for (i1 = 0; i1 < n && ind == 0; i1++) {

		y1 = .0;
		m1 = i1 + 1;
		m2 = 0;

		for (i2 = i1; i2 < n; i2++) {
			y2 = fabs(w[i2][i1]);
			if (y1 < y2) {
				y1 = y2;
				m2 = i2;
			}
		}

		if (y1 < eps)
			ind = 1;

		else {

			for (i2 = i1; i2 < nm; i2++) {
				y1        = w[i1][i2];
				w[i1][i2] = w[m2][i2];
				w[m2][i2] = y1;
			}

			y1 = 1.0 / w[i1][i1];

			for (i2 = m1; i2 < nm; i2++)
				w[i1][i2] *= y1;

			for (i2 = 0; i2 < n; i2++) {
				if (i2 != i1) {
					for (i3 = m1; i3 < nm; i3++)
						w[i2][i3] -= w[i2][i1] * w[i1][i3];
				}
			}
		}
	}

	return(ind);
}

Java

/******************************/
/* Newton法による最小値の計算 */
/*      coded by Y.Suganuma   */
/******************************/
import java.io.*;
import java.util.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double eps, H[][], step, x[], dx[], y[];
		int fun, i1, max, n, opt_1, sw = 0;
		StringTokenizer str;
		String line;
		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
					// データの入力
							// 1 行目
		line = in.readLine();
		str  = new StringTokenizer(line, " ");
		line = str.nextToken();
		line = str.nextToken();
		fun = Integer.parseInt(line);
		line = str.nextToken();
		line = str.nextToken();
		n = Integer.parseInt(line);
		line = str.nextToken();
		line = str.nextToken();
		max = Integer.parseInt(line);
		line = str.nextToken();
		line = str.nextToken();
		opt_1 = Integer.parseInt(line);
							// 2 行目
		line = in.readLine();
		str  = new StringTokenizer(line, " ");
		line = str.nextToken();
		line = str.nextToken();
		eps = Double.parseDouble(line);
		line = str.nextToken();
		line = str.nextToken();
		step = Double.parseDouble(line);
							// 3 行目
		x  = new double [n];
		dx = new double [n];
		y  = new double [1];
		H  = new double [n][2*n];
		line = in.readLine();
		str  = new StringTokenizer(line, " ");
		line = str.nextToken();
		for (i1 = 0; i1 < n; i1++) {
			line = str.nextToken();
			x[i1] = Double.parseDouble(line);
		}
					// 実行
		Kansu kn = new Kansu(fun);
		sw = kn.Newton(opt_1, max, n, eps, step, y, x, dx, H);
					// 結果の出力
		if (sw < 0) {
			System.out.print("   収束しませんでした！");
			switch (sw) {
				case -1:
					System.out.println("（収束回数）");
					break;
				case -2:
					System.out.print("（１次元最適化の区間）");
					break;
				case -3:
					System.out.print("（黄金分割法）");
					break;
			}
		}
		else {
			System.out.print("   結果＝");
			for (i1 = 0; i1 < n; i1++)
				System.out.print(x[i1] + " ");
			System.out.println(" 最小値＝" + y[0] + "  回数＝" + sw);
		}
	}
}

/************************/
/* 関数値，微係数の計算 */
/************************/
class Kansu extends Newton
{
	private int sw;
					// コンストラクタ
	Kansu (int s) {sw = s;}
					// 与えられた点xから，dx方向へk*dxだけ進んだ
					// 点における関数値を計算する
	double snx(double k, double x[], double dx[])
	{
		double x1, y1, z1, y = 0.0, w[] = new double [2];

		switch (sw) {
			case 1:
				w[0] = x[0] + k * dx[0];
				w[1] = x[1] + k * dx[1];
				x1   = w[0] - 1.0;
				y1   = w[1] - 2.0;
				y    = x1 * x1 + y1 * y1;
				break;
			case 2:
				w[0] = x[0] + k * dx[0];
				w[1] = x[1] + k * dx[1];
				x1   = w[1] - w[0] * w[0];
				y1   = 1.0 - w[0];
				y    = 100.0 * x1 * x1 + y1 * y1;
				break;
			case 3:
				w[0] = x[0] + k * dx[0];
				w[1] = x[1] + k * dx[1];
				x1   = 1.5 - w[0] * (1.0 - w[1]);
				y1   = 2.25 - w[0] * (1.0 - w[1] * w[1]);
				z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
				y    = x1 * x1 + y1 * y1 + z1 * z1;
				break;
		}

		return y;
	}
					// 関数の微分
	void dsnx(double x[], double dx[])
	{
		switch (sw) {
			case 1:
				dx[0] = -2.0 * (x[0] - 1.0);
				dx[1] = -2.0 * (x[1] - 2.0);
				break;
			case 2:
				dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0]);
				dx[1] = -200.0 * (x[1] - x[0] * x[0]);
				break;
			case 3:
				dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                        2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                        2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
				dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                        4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                        6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
				break;
		}
	}
					// Hesse行列の逆行列
	int hesse(double x[], double H[][], double eps)
	{
		double x1, x2, x3;
		int sw1 = 0;

		switch (sw) {
			case 1:
				H[0][0] = 2.0;
				H[0][1] = 0.0;
				H[1][0] = 0.0;
				H[1][1] = 2.0;
				H[0][2] = 1.0;
				H[0][3] = 0.0;
				H[1][2] = 0.0;
				H[1][3] = 1.0;
				sw1 = gauss(H, 2, 2, eps);
				break;
			case 2:
				H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0;
				H[0][1] = -400.0 * x[0];
				H[1][0] = H[0][1];
				H[1][1] = 200.0;
				H[0][2] = 1.0;
				H[0][3] = 0.0;
				H[1][2] = 0.0;
				H[1][3] = 1.0;
				sw1 = gauss(H, 2, 2, eps);
				break;
			case 3:
				x1      = 1.0 - x[1];
				x2      = 1.0 - x[1] * x[1];
				x3      = 1.0 - x[1] * x[1] * x[1];
				H[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3;
				H[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
                          4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
                          6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3;
				H[1][0] = H[0][1];
				H[1][1] = 2.0 * x[0] * x[0] +
                          4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
                          12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
                          18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1];
				H[0][2] = 1.0;
				H[0][3] = 0.0;
				H[1][2] = 0.0;
				H[1][3] = 1.0;
				sw1 = gauss(H, 2, 2, eps);
				break;
		}

		return sw1;
	}
}

abstract class Newton {

	/********************************************************/
	/* Newton法                                             */
	/*      opt_1 : =0 : １次元最適化を行わない             */
	/*              =1 : １次元最適化を行う                 */
	/*      max : 最大繰り返し回数                          */
	/*      n : 次元                                        */
	/*      eps : 収束判定条件                              */
	/*      step : きざみ幅                                 */
	/*      y : 最小値                                      */
	/*      x : x(初期値と答え)                             */
	/*      dx : 関数の微分値                               */
	/*      H : Hesse行列の逆行列                           */
	/*      return : >=0 : 正常終了(収束回数)               */
	/*               =-1 : 収束せず                         */
	/*               =-2 : １次元最適化の区間が求まらない   */
	/*               =-3 : 黄金分割法が失敗                 */
	/********************************************************/

	abstract double snx(double k, double x[], double dx[]);
	abstract void dsnx(double x[], double dx[]);
	abstract int hesse(double x[], double H[][], double eps);

	int Newton(int opt_1, int max, int n, double eps, double step, double y[], double x[], double dx[], double H[][])
	{
		double f1, f2, k, sp, y1[] = new double [1], y2[] = new double [1];
		double wk[] = new double [n];
		int count = 0, i1, i2, sw = 0, sw1[] = new int [1];

		y1[0] = snx(0.0, x, dx);

		while (count < max && sw == 0) {
					// 傾きの計算
			dsnx(x, wk);
					// Hesse行列の逆行列の計算
			sw1[0] = hesse(x, H, eps);
					// 収束していない
			if (sw1[0] == 0) {
						// 方向の計算
				count++;
				for (i1 = 0; i1 < n; i1++) {
					dx[i1] = 0.0;
					for (i2 = 0; i2 < n; i2++)
						dx[i1] += H[i1][n+i2] * wk[i2];
				}
						// １次元最適化を行わない
				if (opt_1 == 0) {
							// 新しい点
					for (i1 = 0; i1 < n; i1++)
						x[i1] += dx[i1];
							// 新しい関数値
					y2[0] = snx(0.0, x, dx);
								// 関数値の変化が大きい
					if (Math.abs(y2[0]-y1[0]) > eps)
						y1[0] = y2[0];
								// 収束（関数値の変化＜ｅｐｓ）
					else {
						sw   = count;
						y[0] = y2[0];
					}
				}
						// １次元最適化を行う
				else {
							// 区間を決める
					sw1[0] = 0;
					f1     = y1[0];
					sp     = step;
					f2     = snx(sp, x, dx);
					if (f2 > f1)
						sp = -step;
					for (i1 = 0; i1 < max && sw1[0] == 0; i1++) {
						f2 = snx(sp, x, dx);
						if (f2 > f1)
							sw1[0] = 1;
						else {
							sp *= 2.0;
							f1  = f2;
						}
					}
							// 区間が求まらない
					if (sw1[0] == 0)
						sw = -2;
							// 区間が求まった
					else {
								// 黄金分割法
						k = gold(0.0, sp, eps, y2, sw1, max, x, dx);
								// 黄金分割法が失敗
						if (sw1[0] < 0)
							sw = -3;
								// 黄金分割法が成功
						else {
									// 新しい点
							for (i1 = 0; i1 < n; i1++)
								x[i1] += k * dx[i1];
										// 関数値の変化が大きい
							if (Math.abs(y1[0]-y2[0]) > eps)
								y1[0] = y2[0];
										// 収束（関数値の変化＜ｅｐｓ）
							else {
								sw   = count;
								y[0] = y2[0];
							}
						}
					}
				}
			}
					// 収束（傾き＜ｅｐｓ）
			else {
				sw   = count;
				y[0] = y1[0];
			}
		}

		if (sw == 0)
			sw = -1;

		return sw;
	}

	/****************************************************************/
	/* 黄金分割法(与えられた方向での最小値)                         */
	/*      a,b : 初期区間 a < b                                    */
	/*      eps : 許容誤差                                          */
	/*      val : 関数値                                            */
	/*      ind : 計算状況                                          */
	/*              >= 0 : 正常終了(収束回数)                       */
	/*              = -1 : 収束せず                                 */
	/*      max : 最大試行回数                                      */
	/*      w : 位置                                                */
	/*      dw : 傾きの成分                                         */
	/*      return : 結果(w＋y＊dwのy)                              */
	/****************************************************************/
	double gold(double a, double b, double eps, double val[], int ind[], int max, double w[], double dw[])
	{
		double f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
		int count = 0;
					// 初期設定
		ind[0] = -1;
		x1     = b - tau * (b - a);
		x2     = a + tau * (b - a);
		f1     = snx(x1, w, dw);
		f2     = snx(x2, w, dw);
					// 計算
		while (count < max && ind[0] < 0) {
			count += 1;
			if (f2 > f1) {
				if (Math.abs(b-a) < eps) {
					ind[0] = 0;
					x      = x1;
					val[0] = f1;
				}
				else {
					b  = x2;
					x2 = x1;
					x1 = a + (1.0 - tau) * (b - a);
					f2 = f1;
					f1 = snx(x1, w, dw);
				}
			}
			else {
				if (Math.abs(b-a) < eps) {
					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 = snx(x2, w, dw);
				}
			}
		}
					// 収束した場合の処理
		if (ind[0] == 0) {
			ind[0] = count;
			fa     = snx(a, w, dw);
			fb     = snx(b, w, dw);
			if (fb < fa) {
				a  = b;
				fa = fb;
			}
			if (fa < f1) {
				x      = a;
				val[0] = fa;
			}
		}

		return x;
	}

	/*******************************************************/
	/* 線形連立方程式を解く（逆行列を求める）              */
	/*      w : 方程式の左辺及び右辺                       */
	/*      n : 方程式の数                                 */
	/*      m : 方程式の右辺の列の数                       */
	/*      eps : 正則性を判定する規準                     */
	/*      return : =0 : 正常                             */
	/*               =1 : 逆行列が存在しない               */
	/*******************************************************/
	int gauss(double w[][], int n, int m, double eps)
	{
		double y1, y2;
		int ind = 0, nm, m1, m2, i1, i2, i3;

		nm = n + m;

		for (i1 = 0; i1 < n && ind == 0; i1++) {

			y1 = .0;
			m1 = i1 + 1;
			m2 = 0;

			for (i2 = i1; i2 < n; i2++) {
				y2 = Math.abs(w[i2][i1]);
				if (y1 < y2) {
					y1 = y2;
					m2 = i2;
				}
			}

			if (y1 < eps)
				ind = 1;

			else {

				for (i2 = i1; i2 < nm; i2++) {
					y1        = w[i1][i2];
					w[i1][i2] = w[m2][i2];
					w[m2][i2] = y1;
				}

				y1 = 1.0 / w[i1][i1];

				for (i2 = m1; i2 < nm; i2++)
					w[i1][i2] *= y1;

				for (i2 = 0; i2 < n; i2++) {
					if (i2 != i1) {
						for (i3 = m1; i3 < nm; i3++)
							w[i2][i3] -= w[i2][i1] * w[i1][i3];
					}
				}
			}
		}

		return(ind);
	}
}

JavaScript

　　ここをクリックすると，JavaScript の仕様に適合した形で最小値を求めたい式を入力することによって，任意の関数の最小値を画面上で求めることができます．

<!DOCTYPE HTML>

<HTML>

<HEAD>

	<TITLE>最適化（Newton 法）</TITLE>
	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
	<SCRIPT TYPE="text/javascript">
		str1 = "";
		str2 = new Array();
		str3 = new Array();
		function main()
		{
					// データの設定
			let n    = parseInt(document.getElementById("n").value);
			let max  = parseInt(document.getElementById("max").value);
			let step = parseFloat(document.getElementById("step").value);
			let x    = document.getElementById("x0").value.split(/ {1,}/)
			for (let i1 = 0; i1 < n; i1++)
				x[i1] = parseFloat(x[i1]);
			let opt_1 = 0;
			if (document.getElementById("line").options[1].selected)
				opt_1 = 1;
			str1 = document.getElementById("func").value + ";";
			str2 = document.getElementById("dfunc").value.split(/\n{1,}/)
			for (let i1 = 0; i1 < n; i1++)
				str2[i1] = str2[i1] + ";";
			str3 = document.getElementById("hesse").value.split(/\n{1,}/)
			for (let i1 = 0; i1 < n*n; i1++)
				str3[i1] = str3[i1] + ";";
			let eps = 1.0e-10;
			let y   = new Array();
			let dx  = new Array();
			let H   = new Array();
			for (let i1 = 0; i1 < n; i1++) {
				dx[i1] = 0.0;
				H[i1]  = new Array();
			}
					// 実行
			sw = newton(opt_1, max, n, eps, step, y, x, dx, H, snx, dsnx);
					// 結果
			let res;
			if (sw < 0) {
				res = " 収束しませんでした！";
				switch (sw) {
					case -1:
						res += "（収束回数）";
						break;
					case -2:
						res += "（１次元最適化の区間）";
						break;
					case -3:
						res += "（黄金分割法）";
						break;
				}
				document.getElementById("res").value = res;
			}
			else {
				let c = 100000;
				res = " x ＝";
				for (let i1 = 0; i1 < n; i1++) {
					let s1 = Math.round(x[i1] * c) / c;
					res = res + " " + s1;
				}
				res += "\n";
				let s2 = Math.round(y[0] * c) / c;
				res = res + " 最小値 ＝ " + s2 + "  回数 ＝ " + sw;
				document.getElementById("res").value = res;
			}
		}

		/****************/
		/* 関数値の計算 */
		/****************/
		function snx(k, x, dx)
		{
			let y = eval(str1);
			return y;
		}

		/************************/
		/* 関数値の微分値の計算 */
		/************************/
		function dsnx(n, x, dx)
		{
			for (let i1 = 0; i1 < n; i1++)
				dx[i1] = eval(str2[i1]);
		}

		/******************************************************/
		/* Newton法                                           */
		/*      opt_1 : =0 : １次元最適化を行わない           */
		/*              =1 : １次元最適化を行う               */
		/*      max : 最大繰り返し回数                        */
		/*      n : 次元                                      */
		/*      eps : 収束判定条件                            */
		/*      step : きざみ幅                               */
		/*      y : 最小値                                    */
		/*      x : x(初期値と答え)                           */
		/*      dx : 関数の微分値                             */
		/*      H : Hesse行列の逆行列                         */
		/*      fn : 関数値を計算する関数                     */
		/*      dfn : 関数の微分値を計算する関数              */
		/*      return : >=0 : 正常終了(収束回数)             */
		/*               =-1 : 収束せず                       */
		/*               =-2 : １次元最適化の区間が求まらない */
		/*               =-3 : 黄金分割法が失敗               */
		/******************************************************/
		function newton(opt_1, max, n, eps, step, y, x, dx, H, fn, dfn)
		{
			let f1, f2, k, sp, count = 0, sw = 0;
			let y1 = new Array();
			let y2 = new Array();
			let wk = new Array();
			let sw1 = new Array();

			y1[0] = fn(0.0, x, dx);

			while (count < max && sw == 0) {
					// 傾きの計算
				dfn(n, x, wk);
					// Hesse行列の逆行列の計算
				sw1[0] = hesse(n, x, H, eps);
					// 収束していない
				if (sw1[0] == 0) {
						// 方向の計算
					count++;
					for (let i1 = 0; i1 < n; i1++) {
						dx[i1] = 0.0;
						for (let i2 = 0; i2 < n; i2++)
							dx[i1] += H[i1][n+i2] * wk[i2];
					}
						// １次元最適化を行わない
					if (opt_1 == 0) {
							// 新しい点
						for (let i1 = 0; i1 < n; i1++)
							x[i1] += dx[i1];
							// 新しい関数値
						y2[0] = fn(0.0, x, dx);
								// 関数値の変化が大きい
						if (Math.abs(y2[0]-y1[0]) > eps)
							y1[0] = y2[0];
								// 収束（関数値の変化＜ｅｐｓ）
						else {
							sw   = count;
							y[0] = y2[0];
						}
					}
						// １次元最適化を行う
					else {
							// 区間を決める
						sw1[0] = 0;
						f1     = y1[0];
						sp     = step;
						f2     = fn(sp, x, dx);
						if (f2 > f1)
							sp = -step;
						for (let i1 = 0; i1 < max && sw1[0] == 0; i1++) {
							f2 = fn(sp, x, dx);
							if (f2 > f1)
								sw1[0] = 1;
							else {
								sp *= 2.0;
								f1  = f2;
							}
						}
							// 区間が求まらない
						if (sw1[0] == 0)
							sw = -2;
							// 区間が求まった
						else {
								// 黄金分割法
							k = gold(0.0, sp, eps, y2, sw1, max, x, dx, snx);
								// 黄金分割法が失敗
							if (sw1[0] < 0)
								sw = -3;
								// 黄金分割法が成功
							else {
									// 新しい点
								for (let i1 = 0; i1 < n; i1++)
									x[i1] += k * dx[i1];
										// 関数値の変化が大きい
								if (Math.abs(y1[0]-y2[0]) > eps)
									y1[0] = y2[0];
										// 収束（関数値の変化＜ｅｐｓ）
								else {
									sw   = count;
									y[0] = y2[0];
								}
							}
						}
					}
				}
					// 収束（傾き＜ｅｐｓ）
				else {
					sw   = count;
					y[0] = y1[0];
				}
			}

			if (sw == 0)
				sw = -1;

			return sw;
		}

		/******************************************/
		/* 黄金分割法(与えられた方向での最小値)   */
		/*      a,b : 初期区間 a < b              */
		/*      eps : 許容誤差                    */
		/*      val : 関数値                      */
		/*      ind : 計算状況                    */
		/*              >= 0 : 正常終了(収束回数) */
		/*              = -1 : 収束せず           */
		/*      max : 最大試行回数                */
		/*      w : 位置                          */
		/*      dw : 傾きの成分                   */
		/*      fn : 関数値を計算する関数         */
		/*      return : 結果(w＋y＊dwのy)        */
		/******************************************/
		function gold(a, b, eps, val, ind, max, w, dw, fn)
		{
			let f1, f2, fa, fb, tau = (Math.sqrt(5.0) - 1.0) / 2.0, x = 0.0, x1, x2;
			let count = 0;
					// 初期設定
			ind[0] = -1;
			x1     = b - tau * (b - a);
			x2     = a + tau * (b - a);
			f1     = fn(x1, w, dw);
			f2     = fn(x2, w, dw);
					// 計算
			while (count < max && ind[0] < 0) {
				count += 1;
				if (f2 > f1) {
					if (Math.abs(b-a) < eps) {
						ind[0] = 0;
						x      = x1;
						val[0] = f1;
					}
					else {
						b  = x2;
						x2 = x1;
						x1 = a + (1.0 - tau) * (b - a);
						f2 = f1;
						f1 = fn(x1, w, dw);
					}
				}
				else {
					if (Math.abs(b-a) < eps) {
						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 = fn(x2, w, dw);
					}
				}
			}
					// 収束した場合の処理
			if (ind[0] == 0) {
				ind[0] = count;
				fa     = fn(a, w, dw);
				fb     = fn(b, w, dw);
				if (fb < fa) {
					a  = b;
					fa = fb;
				}
				if (fa < f1) {
					x      = a;
					val[0] = fa;
				}
			}

			return x;
		}

		/*********************/
		/* Hesse行列の逆行列 */
		/*********************/
		function hesse(n, x, H, eps)
		{
			let sw1 = 0;

			for (let i1 = 0; i1 < n; i1++) {
				for (let i2 = 0; i2 < n; i2++)
					H[i1][i2] = eval(str3[i1*n+i2]);
			}
			for (let i1 = 0; i1 < n; i1++) {
				for (let i2 = 0; i2 < n; i2++) {
					if (i1 == i2)
						H[i1][n+i2] = 1.0;
					else
						H[i1][n+i2] = 0.0;
				}
			}

			sw1 = Gauss(H, n, n, eps);

			return sw1;
		}

		/******************************************/
		/* 線形連立方程式を解く（逆行列を求める） */
		/*      w : 方程式の左辺及び右辺          */
		/*      n : 方程式の数                    */
		/*      m : 方程式の右辺の列の数          */
		/*      eps : 正則性を判定する規準        */
		/*      return : =0 : 正常                */
		/*               =1 : 逆行列が存在しない  */
		/******************************************/
		function Gauss(w, n, m, eps)
		{
			let y1, y2, ind = 0, nm, m1, m2;

			nm = n + m;

			for (let i1 = 0; i1 < n && ind == 0; i1++) {

				y1 = .0;
				m1 = i1 + 1;
				m2 = 0;

				for (let i2 = i1; i2 < n; i2++) {
					y2 = Math.abs(w[i2][i1]);
					if (y1 < y2) {
						y1 = y2;
						m2 = i2;
					}
				}

				if (y1 < eps)
					ind = 1;

				else {

					for (let i2 = i1; i2 < nm; i2++) {
						y1        = w[i1][i2];
						w[i1][i2] = w[m2][i2];
						w[m2][i2] = y1;
					}

					y1 = 1.0 / w[i1][i1];

					for (let i2 = m1; i2 < nm; i2++)
						w[i1][i2] *= y1;

					for (let i2 = 0; i2 < n; i2++) {
						if (i2 != i1) {
							for (let i3 = m1; i3 < nm; i3++)
								w[i2][i3] -= w[i2][i1] * w[i1][i3];
						}
					}
				}
			}

			return(ind);
		}
	</SCRIPT>

</HEAD>

<BODY STYLE="font-size: 130%; background-color: #eeffee;">

	<H2 STYLE="text-align:center"><B>最適化（Newton 法）</B></H2>

	<DL>
		<DT>　　テキストフィールドおよびテキストエリアには，例として，以下に示すような関数の最小値を求める場合に対する値が設定されています（他の問題を実行する場合は，それらを適切に修正してください）．
		<P STYLE="text-align:center"><IMG SRC="newton.gif"></P>
		<DT>n 変数関数 f(<B>x</B>) に対する Newton 法においては，現在の点 <B>x</B> における傾き，
		<P STYLE="text-align:center"><IMG SRC="newton1.gif"></P>
		<DT>と Hesse 行列の計算結果から方向 <B>p</B> を求め，その方向にαだけ進んだ点の関数値 f( <B>x</B> + α<B>p</B>) を計算する，という処理を繰り返します．目的関数の箇所には，この例に対する f( <B>x</B> + α<B>p</B>) の計算式，目的関数の微分の箇所には，<B>g</B> の計算式が与えられています．ただし，プログラム内では，変数 x は，この例の x と y からなる配列，変数 dx は <B>g</B> に対応し，関数 f を x 及び y で微分したものからなる配列，また，k はαに対応する変数として表現されています．また，Hesse 行列は，1 行 1 列，1 行 2 列，・・・，n 行 n 列の順で入力して下さい． なお，目的関数，目的関数の微分，及び，Hesse 行列は，JavaScript の仕様に適合した形式で記述してあることに注意してください． 
	</DL>

	<DIV STYLE="text-align:center">
		変数の数(n)：<INPUT ID="n" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="2">　
		最大繰り返し回数：<INPUT ID="max" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="200">　
		刻み幅：<INPUT ID="step" STYLE="font-size: 100%" TYPE="text" SIZE="3" VALUE="0.1"><BR>
		初期値(n個)：<INPUT ID="x0" STYLE="font-size: 100%" TYPE="text" SIZE="20" VALUE="0.0 0.0">　
		１次元最適化：<SELECT STYLE="font-size: 100%" SIZE="1" ID="line">
							<OPTION VALUE="no"> 行わない </OPTION><BR>
							<OPTION VALUE="yes"> 行う </OPTION><BR>
					  </SELECT><BR><BR>
		目的関数：　f(x+αp) = <INPUT ID="func" STYLE="font-size: 100%" TYPE="text" SIZE="70" VALUE="100.0 * (x[1] + k * dx[1] - (x[0] + k * dx[0]) * (x[0] + k * dx[0])) * (x[1] + k * dx[1] - (x[0] + k * dx[0]) * (x[0] + k * dx[0])) + (1.0 - (x[0] + k * dx[0])) * (1.0 - (x[0] + k * dx[0]))"><BR><BR>
		目的関数の微分：　f'(x) = <TEXTAREA ID="dfunc" STYLE="font-size: 110%" COLS="70" ROWS="10">
400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
-200.0 * (x[1] - x[0] * x[0])
		</TEXTAREA><BR><BR>
		Hesse行列：　H = <TEXTAREA ID="hesse" STYLE="font-size: 110%" COLS="70" ROWS="10">
400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
-400.0 * x[0]
400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
200.0
		</TEXTAREA><BR><BR>
		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">実行</BUTTON><BR><BR>
		結果：<TEXTAREA ID="res" STYLE="font-size: 110%" COLS="40" ROWS="2"></TEXTAREA>
	</DIV>

</BODY>

</HTML>

PHP

<?php

/******************************/
/* Newton法による最小値の計算 */
/*      coded by Y.Suganuma   */
/******************************/

	$sw = 0;
					// データの入力
	fscanf(STDIN, "%*s %d %*s %d %*s %d %*s %d", $fun, $n, $max, $opt_1);
	fscanf(STDIN, "%*s %lf %*s %lf", $eps, $step);
	$x  = array($n);
	$dx = array(0.0, 0.0);

	$H   = array($n);
	$str = trim(fgets(STDIN));
	strtok($str, " ");
	for ($i1 = 0; $i1 < $n; $i1++) {
		$x[$i1] = floatval(strtok(" "));
		$H[$i1] = array(2*$n);
	}
					// 実行
	switch ($fun) {
		case 1:
			$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx1", "dsnx1", "hesse1");
			break;
		case 2:
			$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx2", "dsnx2", "hesse2");
			break;
		case 3:
			$sw = Newton($opt_1, $max, $n, $eps, $step, $y, $x, $dx, $H, "snx3", "dsnx3", "hesse3");
			break;
	}
					// 結果の出力
	if ($sw < 0) {
		printf("   収束しませんでした！");
		switch ($sw) {
			case -1:
				printf("（収束回数）\n");
				break;
			case -2:
				printf("（１次元最適化の区間）\n");
				break;
			case -3:
				printf("（黄金分割法）\n");
				break;
		}
	}
	else {
		printf("   結果＝");
		for ($i1 = 0; $i1 < $n; $i1++)
			printf("%f ", $x[$i1]);
		printf(" 最小値＝%f  回数＝%d\n", $y, $sw);
	}

/********************************************************/
/* Newton法                                             */
/*      opt_1 : =0 : １次元最適化を行わない             */
/*              =1 : １次元最適化を行う                 */
/*      max : 最大繰り返し回数                          */
/*      n : 次元                                        */
/*      eps : 収束判定条件                              */
/*      step : きざみ幅                                 */
/*      y : 最小値                                      */
/*      x : x(初期値と答え)                             */
/*      dx : 関数の微分値                               */
/*      H : Hesse行列の逆行列                           */
/*      snx : 関数値を計算する関数名                    */
/*      dsnx : 関数の微分を計算する関数名(符号を変える) */
/*      hesse : Hesse行列の逆行列を計算する関数名       */
/*      return : >=0 : 正常終了(収束回数)               */
/*               =-1 : 収束せず                         */
/*               =-2 : １次元最適化の区間が求まらない   */
/*               =-3 : 黄金分割法が失敗                 */
/********************************************************/

function Newton($opt_1, $max, $n, $eps, $step, &$y, &$x, $dx, $H, $snx, $dsnx, $hesse)
{
	$wk    = array($n);
	$count = 0;
	$sw    = 0;

	$y1 = $snx(0.0, $x, $dx);

	while ($count < $max && $sw == 0) {
					// 傾きの計算
		$dsnx($x, $wk);
					// Hesse行列の逆行列の計算
		$sw1 = $hesse($x, $H, $eps);
					// 収束していない
		if ($sw1 == 0) {
						// 方向の計算
			$count++;
			for ($i1 = 0; $i1 < $n; $i1++) {
				$dx[$i1] = 0.0;
				for ($i2 = 0; $i2 < $n; $i2++)
					$dx[$i1] += $H[$i1][$n+$i2] * $wk[$i2];
			}
						// １次元最適化を行わない
			if ($opt_1 == 0) {
							// 新しい点
				for ($i1 = 0; $i1 < $n; $i1++)
					$x[$i1] += $dx[$i1];
							// 新しい関数値
				$y2 = $snx(0.0, $x, $dx);
								// 関数値の変化が大きい
				if (abs($y2-$y1) > $eps)
					$y1 = $y2;
								// 収束（関数値の変化＜ｅｐｓ）
				else {
					$sw = $count;
					$y = $y2;
				}
			}
						// １次元最適化を行う
			else {
							// 区間を決める
				$sw1 = 0;
				$f1  = $y1;
				$sp  = $step;
				$f2  = $snx($sp, $x, $dx);
				if ($f2 > $f1)
					$sp = -$step;
				for ($i1 = 0; $i1 < $max && $sw1 == 0; $i1++) {
					$f2 = $snx($sp, $x, $dx);
					if ($f2 > $f1)
						$sw1 = 1;
					else {
						$sp *= 2.0;
						$f1  = $f2;
					}
				}
							// 区間が求まらない
				if ($sw1 == 0)
					$sw = -2;
							// 区間が求まった
				else {
								// 黄金分割法
					$k = gold(0.0, $sp, $eps, $y2, $sw1, $max, $x, $dx, $snx);
								// 黄金分割法が失敗
					if ($sw1 < 0)
						$sw = -3;
								// 黄金分割法が成功
					else {
									// 新しい点
						for ($i1 = 0; $i1 < $n; $i1++)
							$x[$i1] += $k * $dx[$i1];
										// 関数値の変化が大きい
						if (abs($y1-$y2) > $eps)
							$y1 = $y2;
										// 収束（関数値の変化＜ｅｐｓ）
						else {
							$sw = $count;
							$y = $y2;
						}
					}
				}
			}
		}
					// 収束（傾き＜ｅｐｓ）
		else {
			$sw = $count;
			$y = $y1;
		}
	}

	if ($sw == 0)
		$sw = -1;

	return $sw;
}

/*******************************/
/* Hesse行列の逆行列を計算する */
/*******************************/
					// 関数１
function hesse1($x, &$H, $eps)
{
	$H[0][0] = 2.0;
	$H[0][1] = 0.0;
	$H[1][0] = 0.0;
	$H[1][1] = 2.0;
	$H[0][2] = 1.0;
	$H[0][3] = 0.0;
	$H[1][2] = 0.0;
	$H[1][3] = 1.0;

	$sw = Gauss($H, 2, 2, $eps);

	return $sw;
}
					// 関数２
function hesse2($x, &$H, $eps)
{
	$H[0][0] = 400.0 * (3.0 * $x[0] * $x[0] - $x[1]) + 2.0;
	$H[0][1] = -400.0 * $x[0];
	$H[1][0] = $H[0][1];
	$H[1][1] = 200.0;
	$H[0][2] = 1.0;
	$H[0][3] = 0.0;
	$H[1][2] = 0.0;
	$H[1][3] = 1.0;

	$sw = Gauss($H, 2, 2, $eps);

	return $sw;
}
					// 関数３
function hesse3($x, &$H, $eps)
{
	$x1      = 1.0 - $x[1];
	$x2      = 1.0 - $x[1] * $x[1];
	$x3      = 1.0 - $x[1] * $x[1] * $x[1];
	$H[0][0] = 2.0 * $x1 * $x1 + 2.0 * $x2 * $x2 + 2.0 * $x3 * $x3;
	$H[0][1] = 2.0 * (1.5 - $x[0] * $x1) - 2.0 * $x[0] * $x1 +
              4.0 * $x[1] * (2.25 - $x[0] * $x2) - 4.0 * $x[0] * $x[1] * $x2 +
              6.0 * $x[1] * $x[1] * (2.625 - $x[0] * $x3) - 6.0 * $x[0] * $x[1] * $x[1] * $x3;
	$H[1][0] = $H[0][1];
	$H[1][1] = 2.0 * $x[0] * $x[0] +
              4.0 * $x[0] * (2.25 - $x[0] * $x2) + 8.0 * $x[0] * $x[0] * $x[1] * $x[1] +
              12.0 * $x[0] * $x[1] * (2.625 - $x[0] * $x3) +
              18.0 * $x[0] * $x[0] * $x[1] * $x[1] * $x[1] * $x[1];
	$H[0][2] = 1.0;
	$H[0][3] = 0.0;
	$H[1][2] = 0.0;
	$H[1][3] = 1.0;

	$sw = Gauss($H, 2, 2, $eps);

	return $sw;
}

/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める）              */
/*      w : 方程式の左辺及び右辺                       */
/*      n : 方程式の数                                 */
/*      m : 方程式の右辺の列の数                       */
/*      eps : 正則性を判定する規準                     */
/*      return : =0 : 正常                             */
/*               =1 : 逆行列が存在しない               */
/*******************************************************/
function Gauss(&$w, $n, $m, $eps)
{
	$ind = 0;
	$nm  = $n + $m;

	for ($i1 = 0; $i1 < $n && $ind == 0; $i1++) {

		$y1 = .0;
		$m1 = $i1 + 1;
		$m2 = 0;

		for ($i2 = $i1; $i2 < $n; $i2++) {
			$y2 = abs($w[$i2][$i1]);
			if ($y1 < $y2) {
				$y1 = $y2;
				$m2 = $i2;
			}
		}

		if ($y1 < $eps)
			$ind = 1;

		else {

			for ($i2 = $i1; $i2 < $nm; $i2++) {
				$y1          = $w[$i1][$i2];
				$w[$i1][$i2] = $w[$m2][$i2];
				$w[$m2][$i2] = $y1;
			}

			$y1 = 1.0 / $w[$i1][$i1];

			for ($i2 = $m1; $i2 < $nm; $i2++)
				$w[$i1][$i2] *= $y1;

			for ($i2 = 0; $i2 < $n; $i2++) {
				if ($i2 != $i1) {
					for ($i3 = $m1; $i3 < $nm; $i3++)
						$w[$i2][$i3] -= $w[$i2][$i1] * $w[$i1][$i3];
				}
			}
		}
	}

	return($ind);
}

/*********************************************/
/* 与えられた点xから，dx方向へk*dxだけ進んだ */
/* 点における関数値を計算する                */
/*********************************************/
					// 関数１
function snx1($k, $x, $dx)
{
	$w  = array(2);

	$w[0] = $x[0] + $k * $dx[0];
	$w[1] = $x[1] + $k * $dx[1];
	$x1   = $w[0] - 1.0;
	$y1   = $w[1] - 2.0;
	$y    = $x1 * $x1 + $y1 * $y1;

	return $y;
}
					// 関数２
function snx2($k, $x, $dx)
{
	$w  = array(2);

	$w[0] = $x[0] + $k * $dx[0];
	$w[1] = $x[1] + $k * $dx[1];
	$x1   = $w[1] - $w[0] * $w[0];
	$y1   = 1.0 - $w[0];
	$y    = 100.0 * $x1 * $x1 + $y1 * $y1;

	return $y;
}
					// 関数３
function snx3($k, $x, $dx)
{
	$w = array(2);

	$w[0] = $x[0] + $k * $dx[0];
	$w[1] = $x[1] + $k * $dx[1];
	$x1   = 1.5 - $w[0] * (1.0 - $w[1]);
	$y1   = 2.25 - $w[0] * (1.0 - $w[1] * $w[1]);
	$z1   = 2.625 - $w[0] * (1.0 - $w[1] * $w[1] * $w[1]);
	$y    = $x1 * $x1 + $y1 * $y1 + $z1 * $z1;

	return $y;
}

/********************/
/* 微係数を計算する */
/********************/
					// 関数１
function dsnx1($x, &$dx)
{
	$dx[0] = -2.0 * ($x[0] - 1.0);
	$dx[1] = -2.0 * ($x[1] - 2.0);
}
					// 関数２
function dsnx2($x, &$dx)
{
	$dx[0] = 400.0 * $x[0] * ($x[1] - $x[0] * $x[0]) + 2.0 * (1.0 - $x[0]);
	$dx[1] = -200.0 * ($x[1] - $x[0] * $x[0]);
}
					// 関数３
function dsnx3($x, &$dx)
{
	$dx[0] = 2.0 * (1.0 - $x[1]) * (1.5 - $x[0] * (1.0 - $x[1])) +
            2.0 * (1.0 - $x[1] * $x[1]) * (2.25 - $x[0] * (1.0 - $x[1] * $x[1])) +
            2.0 * (1.0 - $x[1] * $x[1] * $x[1]) * (2.625 - $x[0] * (1.0 - $x[1] * $x[1] * $x[1]));
	$dx[1] = -2.0 * $x[0] * (1.5 - $x[0] * (1.0 - $x[1])) -
            4.0 * $x[0] * $x[1] * (2.25 - $x[0] * (1.0 - $x[1] * $x[1])) -
            6.0 * $x[0] * $x[1] * $x[1] * (2.625 - $x[0] * (1.0 - $x[1] * $x[1] * $x[1]));
}

/****************************************************************/
/* 黄金分割法(与えられた方向での最小値)                         */
/*      a,b : 初期区間 a < b                                    */
/*      eps : 許容誤差                                          */
/*      val : 関数値                                            */
/*      ind : 計算状況                                          */
/*              >= 0 : 正常終了(収束回数)                       */
/*              = -1 : 収束せず                                 */
/*      max : 最大試行回数                                      */
/*      w : 位置                                                */
/*      dw : 傾きの成分                                         */
/*      snx : 関数値を計算する関数の名前                        */
/*      return : 結果(w＋y＊dwのy)                              */
/****************************************************************/
function gold($a, $b, $eps, &$val, &$ind, $max, $w, $dw, $snx)
{
					// 初期設定
	$tau   = (sqrt(5.0) - 1.0) / 2.0;
	$x     = 0.0;
	$count = 0;
	$ind   = -1;
	$x1    = $b - $tau * ($b - $a);
	$x2    = $a + $tau * ($b - $a);
	$f1    = $snx($x1, $w, $dw);
	$f2    = $snx($x2, $w, $dw);
					// 計算
	while ($count < $max && $ind < 0) {
		$count += 1;
		if ($f2 > $f1) {
			if (abs($b-$a) < $eps) {
				$ind = 0;
				$x   = $x1;
				$val = $f1;
			}
			else {
				$b  = $x2;
				$x2 = $x1;
				$x1 = $a + (1.0 - $tau) * ($b - $a);
				$f2 = $f1;
				$f1 = $snx($x1, $w, $dw);
			}
		}
		else {
			if (abs($b-$a) < $eps) {
				$ind = 0;
				$x   = $x2;
				$val = $f2;
				$f1  = $f2;
			}
			else {
				$a = $x1;
				$x1 = $x2;
				$x2 = $b - (1.0 - $tau) * ($b - $a);
				$f1 = $f2;
				$f2 = $snx($x2, $w, $dw);
			}
		}
	}
					// 収束した場合の処理
	if ($ind == 0) {
		$ind = $count;
		$fa  = $snx($a, $w, $dw);
		$fb  = $snx($b, $w, $dw);
		if ($fb < $fa) {
			$a  = $b;
			$fa = $fb;
		}
		if ($fa < $f1) {
			$x   = $a;
			$val = $fa;
		}
	}

	return $x;
}

?>

Ruby

　　元来，関数 Gauss は，終了状態を表す変数 ind を return していました．しかし，この方法では，ブロック内で関数 Gauss を呼んだとき，ブロックに引数として渡した配列の値の変化を，ブロックを呼んだ側で知ることができません．そのため，関数 Gauss によって目的とする配列を return するように修正し，その配列の 1 行目の最後の要素に ind を記憶しています．ブロック付き関数呼び出しと C++ などにおける関数名を引数とする場合との違いが大きく出てしまいました．しかし，このままの仕様では，ブロック付き関数呼び出しにおいて，複数の種類のデータを返したいような場合は不便だと思います．

#*****************************/
# Newton法による最小値の計算 */
#      coded by Y.Suganuma   */
#*****************************/

#******************************************************/
# 線形連立方程式を解く（逆行列を求める）              */
#      w : 方程式の左辺及び右辺                       */
#      n : 方程式の数                                 */
#      m : 方程式の右辺の列の数                       */
#      eps : 正則性を判定する規準                     */
#      return : =0 : 正常                             */
#               =1 : 逆行列が存在しない               */
#******************************************************/
def Gauss(w, n, m, eps)

	ind = 0
	nm  = n + m

	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

	w[0][nm] = ind
	return w
end

#******************************************************/
# 与えられた点xから，dx方向へk*dxだけ進んだ点における */
# 関数値，微係数，及び，Hesse行列の逆行列を計算する   */
#******************************************************/
					# 関数１
snx1 = Proc.new { |sw, k, x, dx, h, eps|
	if sw == 0
		w    = Array.new(2)
		w[0] = x[0] + k * dx[0]
		w[1] = x[1] + k * dx[1]
		x1   = w[0] - 1.0
		y1   = w[1] - 2.0
		x1 * x1 + y1 * y1
	elsif sw == 1
		dx[0] = -2.0 * (x[0] - 1.0)
		dx[1] = -2.0 * (x[1] - 2.0)
	else
		h = [[2.0, 0.0, 1.0, 0.0, 0.0], [0.0, 2.0, 0.0, 1.0, 0.0]]
		Gauss(h, 2, 2, eps)
	end
}
					# 関数２
snx2 = Proc.new { |sw, k, x, dx, h, eps|
	if sw == 0
		w    = Array.new(2)
		w[0] = x[0] + k * dx[0]
		w[1] = x[1] + k * dx[1]
		x1   = w[1] - w[0] * w[0]
		y1   = 1.0 - w[0]
		100.0 * x1 * x1 + y1 * y1
	elsif sw == 1
		dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
		dx[1] = -200.0 * (x[1] - x[0] * x[0])
	else
		h[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
		h[0][1] = -400.0 * x[0]
		h[1][0] = h[0][1]
		h[1][1] = 200.0
		h[0][2] = 1.0
		h[0][3] = 0.0
		h[1][2] = 0.0
		h[1][3] = 1.0
		Gauss(h, 2, 2, eps)
	end
}
					# 関数３
snx3 = Proc.new { |sw, k, x, dx, h, eps|
	if sw == 0
		w    = Array.new(2)
		w[0] = x[0] + k * dx[0]
		w[1] = x[1] + k * dx[1]
		x1   = 1.5 - w[0] * (1.0 - w[1])
		y1   = 2.25 - w[0] * (1.0 - w[1] * w[1])
		z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1])
		x1 * x1 + y1 * y1 + z1 * z1
	elsif sw == 1
		dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
		dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
	else
		x1      = 1.0 - x[1]
		x2      = 1.0 - x[1] * x[1]
		x3      = 1.0 - x[1] * x[1] * x[1]
		h[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3
		h[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
                  4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
                  6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3
		h[1][0] = h[0][1]
		h[1][1] = 2.0 * x[0] * x[0] +
                  4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
                  12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
                  18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1]
		h[0][2] = 1.0
		h[0][3] = 0.0
		h[1][2] = 0.0
		h[1][3] = 1.0
		Gauss(h, 2, 2, eps)
	end
}

#***************************************************************/
# 黄金分割法(与えられた方向での最小値)                         */
#      a,b : 初期区間 a < b                                    */
#      eps : 許容誤差                                          */
#      val : 関数値                                            */
#      ind : 計算状況                                          */
#              >= 0 : 正常終了(収束回数)                       */
#              = -1 : 収束せず                                 */
#      max : 最大試行回数                                      */
#      w : 位置                                                */
#      dw : 傾きの成分                                         */
#      h : Hesse行列の逆行列                                   */
#      snx : 関数値を計算する関数の名前                        */
#            （微分は，その符号を変える)                       */
#      return : 結果(w＋y＊dwのy)                              */
#***************************************************************/
def gold(a, b, eps, val, ind, max, w, dw, h, &snx)
					# 初期設定
	tau    = (Math.sqrt(5.0) - 1.0) / 2.0
	x      = 0.0
	count  = 0
	ind[0] = -1
	x1     = b - tau * (b - a)
	x2     = a + tau * (b - a)
	f1     = snx.call(0, x1, w, dw, h, eps)
	f2     = snx.call(0, x2, w, dw, h, eps)
					# 計算
	while count < max && ind[0] < 0
		count += 1
		if f2 > f1
			if (b-a).abs() < eps
				ind[0] = 0
				x      = x1
				val[0] = f1
			else
				b  = x2
				x2 = x1
				x1 = a + (1.0 - tau) * (b - a)
				f2 = f1
				f1 = snx.call(0, x1, w, dw, h, eps)
			end
		else
			if (b-a).abs() < eps
				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 = snx.call(0, x2, w, dw, h, eps)
			end
		end
	end
					# 収束した場合の処理
	if ind[0] == 0
		ind[0] = count
		fa     = snx.call(0, a, w, dw, h, eps)
		fb     = snx.call(0, b, w, dw, h, eps)
		if fb < fa
			a  = b
			fa = fb
		end
		if fa < f1
			x      = a
			val[0] = fa
		end
	end

	return x
end

#*******************************************************/
# Newton法                                             */
#      opt_1 : =0 : １次元最適化を行わない             */
#              =1 : １次元最適化を行う                 */
#      max : 最大繰り返し回数                          */
#      n : 次元                                        */
#      eps : 収束判定条件                              */
#      step : きざみ幅                                 */
#      y : 最小値                                      */
#      x : x(初期値と答え)                             */
#      dx : 関数の微分値                               */
#      h : Hesse行列の逆行列                           */
#      snx : 関数値，その微分，及び，Hesse行列の逆行列 */
#      return : >=0 : 正常終了(収束回数)               */
#               =-1 : 収束せず                         */
#               =-2 : １次元最適化の区間が求まらない   */
#               =-3 : 黄金分割法が失敗                 */
#*******************************************************/
def Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx)

	sw1 = Array.new(n)
	y2  = Array.new(n)
	wk  = Array.new(n)
	for i1 in 0 ... n
		wk[i1] = 0.0
	end
	count = 0
	sw    = 0

	y1 = snx.call(0, 0.0, x, dx, h, eps)

	while count < max && sw == 0
					# 傾きの計算
		snx.call(1, 0.0, x, wk, h, eps)
					# Hesse行列の逆行列の計算
		h      = snx.call(2, 0.0, x, wk, h, eps)
		sw1[0] = h[0][2*n]   # もっと良い方法はあるか？
					# 収束していない
		if sw1[0] == 0
						# 方向の計算
			count += 1
			for i1 in 0 ... n
				dx[i1] = 0.0
				for i2 in 0 ... n
					dx[i1] += h[i1][n+i2] * wk[i2]
				end
			end
						# １次元最適化を行わない
			if opt_1 == 0
							# 新しい点
				for i1 in 0 ... n
					x[i1] += dx[i1]
				end
							# 新しい関数値
				y2[0] = snx.call(0, 0.0, x, dx, h, eps)
								# 関数値の変化が大きい
				if (y2[0]-y1).abs() > eps
					y1 = y2[0]
								# 収束（関数値の変化＜ｅｐｓ）
				else
					sw   = count
					y[0] = y2[0]
				end
						# １次元最適化を行う
			else
							# 区間を決める
				sw1[0] = 0
				f1     = y1
				sp     = step
				f2     = snx.call(0, sp, x, dx, h, eps)
				if f2 > f1
					sp = -step
				end
				for i1 in 0 ... max
					f2 = snx.call(0, sp, x, dx, h, eps)
					if f2 > f1
						sw1[0] = 1
						break
					else
						sp *= 2.0
						f1  = f2
					end
				end
							# 区間が求まらない
				if sw1[0] == 0
					sw = -2
							# 区間が求まった
				else
								# 黄金分割法
					k = gold(0.0, sp, eps, y2, sw1, max, x, dx, h, &snx)
								# 黄金分割法が失敗
					if sw1[0] < 0
						sw = -3
								# 黄金分割法が成功
					else
									# 新しい点
						for i1 in 0 ... n
							x[i1] += k * dx[i1]
						end
										# 関数値の変化が大きい
						if (y1-y2[0]).abs() > eps
							y1 = y2[0]
										# 収束（関数値の変化＜ｅｐｓ）
						else
							sw   = count
							y[0] = y2[0]
						end
					end
				end
			end
					# 収束（傾き＜ｅｐｓ）
		else
			sw = count
			y[0] = y1
		end
	end

	if sw == 0
		sw = -1
	end

	return sw
end
				# データの入力
str   = gets()
a     = str.split(" ")
fun   = Integer(a[1])
n     = Integer(a[3])
max   = Integer(a[5])
opt_1 = Integer(a[7])
str   = gets()
a     = str.split(" ")
eps   = Float(a[1])
step  = Float(a[3])
x     = Array.new(n)
dx    = Array.new(n)
y     = Array.new(1)
sw    = 0
str   = gets()
a     = str.split(" ")
h     = Array.new(n)
for i1 in 0 ... n
	x[i1]  = Float(a[i1+1])
	dx[i1] = 0.0
	h[i1]  = Array.new(2*n+1)
end
				# 実行
case fun
	when 1
		sw = Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx1)
	when 2
		sw = Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx2)
	when 3
		sw = Newton(opt_1, max, n, eps, step, y, x, dx, h, &snx3)
end
				# 結果の出力
if sw < 0
	printf("   収束しませんでした！")
	case sw
		when -1
			printf("（収束回数）\n")
		when -2
			printf("（１次元最適化の区間）\n")
		when -3
			printf("（黄金分割法）\n")
	end
else
	printf("   結果＝")
	for i1 in 0 ... n
		printf("%f ", x[i1])
	end
	printf(" 最小値＝%f  回数＝%d\n", y[0], sw)
end

Python

# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *

############################################
# 黄金分割法(与えられた方向での最小値)
#      a,b : 初期区間 a < b
#      eps : 許容誤差
#      val : 関数値
#      ind : 計算状況
#              >= 0 : 正常終了(収束回数)
#              = -1 : 収束せず
#      max : 最大試行回数
#      w : 位置
#      dw : 傾きの成分
#      snx : 関数値を計算する関数の名前
#      return : 結果(w＋y＊dwのy)
#      coded by Y.Suganuma
############################################

def gold(a, b, eps, val, ind, max, w, dw, snx) :
				# 初期設定
	tau    = (sqrt(5.0) - 1.0) / 2.0
	x      = 0.0
	count  = 0
	ind[0] = -1
	x1     = b - tau * (b - a)
	x2     = a + tau * (b - a)
	f1     = snx(x1, w, dw)
	f2     = snx(x2, w, dw)
				# 計算
	while count < max and ind[0] < 0 :
		count += 1
		if f2 > f1 :
			if abs(b-a) < eps :
				ind[0] = 0
				x      = x1
				val[0] = f1
			else :
				b  = x2
				x2 = x1
				x1 = a + (1.0 - tau) * (b - a)
				f2 = f1
				f1 = snx(x1, w, dw)
		else :
			if abs(b-a) < eps :
				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 = snx(x2, w, dw)
				# 収束した場合の処理
	if ind[0] == 0 :
		ind[0] = count
		fa     = snx(a, w, dw)
		fb     = snx(b, w, dw)
		if fb < fa :
			a  = b
			fa = fb
		if fa < f1 :
			x      = a
			val[0] = fa

	return x

############################################
# 線形連立方程式を解く（逆行列を求める）
#      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 range(0, n) :

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

			for i2 in range(m1, nm) :
				w[i1][i2] *= y1

			for i2 in range(0, n) :
				if i2 != i1 :
					for i3 in range(m1, nm) :
						w[i2][i3] -= (w[i2][i1] * w[i1][i3])

	return ind

############################################
# Newton法
#      opt_1 : =0 : １次元最適化を行わない
#              =1 : １次元最適化を行う
#      max : 最大繰り返し回数
#      n : 次元
#      eps : 収束判定条件
#      step : きざみ幅
#      y : 最小値
#      x : x(初期値と答え)
#      dx : 関数の微分値
#      H : Hesse行列の逆行列
#      snx : 関数値を計算する関数名
#      dsnx : 関数の微分を計算する関数名(符号を変える)
#      hesse : Hesse行列の逆行列を計算する関数名
#      return : >=0 : 正常終了(収束回数)
#               =-1 : 収束せず
#               =-2 : １次元最適化の区間が求まらない
#               =-3 : 黄金分割法が失敗
#      coded by Y.Suganuma
############################################

def Newton(opt_1, max, n, eps, step, y, x, dx, H, snx, dsnx, hesse) :

	wk    = np.empty(n, np.float)
	count = 0
	sw    = 0
	y2    = np.empty(1, np.float)
	sw1   = np.empty(1, np.int)

	y1    = snx(0.0, x, dx)

	while count < max and sw == 0 : 
			# 傾きの計算
		dsnx(x, wk)
			# Hesse行列の逆行列の計算
		sw1[0] = hesse(x, H, eps)
			# 収束していない
		if sw1[0] == 0 :
				# 方向の計算
			count += 1
			for i1 in range(0, n) :
				dx[i1] = 0.0
				for i2 in range(0, n) :
					dx[i1] += H[i1][n+i2] * wk[i2]
				# １次元最適化を行わない
			if opt_1 == 0 :
					# 新しい点
				for i1 in range(0, n) :
					x[i1] += dx[i1]
					# 新しい関数値
				y2[0] = snx(0.0, x, dx)
						# 関数値の変化が大きい
				if abs(y2[0]-y1) > eps :
					y1 = y2[0]
						# 収束（関数値の変化＜ｅｐｓ）
				else :
					sw   = count
					y[0] = y2[0]
				# １次元最適化を行う
			else :
					# 区間を決める
				sw1[0] = 0
				f1     = y1
				sp     = step
				f2     = snx(sp, x, dx)
				if f2 > f1 :
					sp = -step
				for i1 in range(0, max) :
					f2 = snx(sp, x, dx)
					if f2 > f1 :
						sw1[0] = 1
						break
					else :
						sp *= 2.0
						f1  = f2
					# 区間が求まらない
				if sw1[0] == 0 :
					sw = -2
					# 区間が求まった
				else :
						# 黄金分割法
					k = gold(0.0, sp, eps, y2, sw1, max, x, dx, snx)
						# 黄金分割法が失敗
					if sw1[0] < 0 :
						sw = -3
						# 黄金分割法が成功
					else :
							# 新しい点
						for i1 in range(0, n) :
							x[i1] += k * dx[i1]
								# 関数値の変化が大きい
						if abs(y1-y2[0]) > eps :
							y1 = y2[0]
								# 収束（関数値の変化＜ｅｐｓ）
						else :
							sw   = count
							y[0] = y2[0]
			# 収束（傾き＜ｅｐｓ）
		else :
			sw   = count
			y[0] = y1

	if sw == 0 :
		sw = -1

	return sw

############################################
# Newton法による最小値の計算
#      coded by Y.Suganuma
############################################

############################################
# 与えられた点xから，dx方向へk*dxだけ進んだ
# 点における関数値を計算する
############################################
			 # 関数１
def snx1(k, x, dx) :
	w    = np.empty(2, np.float)
	w[0] = x[0] + k * dx[0]
	w[1] = x[1] + k * dx[1]
	x1   = w[0] - 1.0
	y1   = w[1] - 2.0
	y    = x1 * x1 + y1 * y1
	return y
			# 関数２
def snx2(k, x, dx) :
	w    = np.empty(2, np.float)
	w[0] = x[0] + k * dx[0]
	w[1] = x[1] + k * dx[1]
	x1   = w[1] - w[0] * w[0]
	y1   = 1.0 - w[0]
	y    = 100.0 * x1 * x1 + y1 * y1
	return y
			# 関数３
def snx3(k, x, dx) :
	w    = np.empty(2, np.float)
	w[0] = x[0] + k * dx[0]
	w[1] = x[1] + k * dx[1]
	x1   = 1.5 - w[0] * (1.0 - w[1])
	y1   = 2.25 - w[0] * (1.0 - w[1] * w[1])
	z1   = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1])
	y    = x1 * x1 + y1 * y1 + z1 * z1
	return y

############################################
# 微係数を計算する
############################################
			# 関数１
def dsnx1(x, dx) :
	dx[0] = -2.0 * (x[0] - 1.0)
	dx[1] = -2.0 * (x[1] - 2.0)
			# 関数２
def dsnx2(x, dx) :
	dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0])
	dx[1] = -200.0 * (x[1] - x[0] * x[0])
			# 関数３
def dsnx3(x, dx) :
	dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) + 2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) + 2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))
	dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) - 4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) - 6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]))

############################################
# Hesse行列の逆行列を計算する */
############################################

			# 関数１
def hesse1(x, H, eps) :

	H[0][0] = 2.0
	H[0][1] = 0.0
	H[1][0] = 0.0
	H[1][1] = 2.0
	H[0][2] = 1.0
	H[0][3] = 0.0
	H[1][2] = 0.0
	H[1][3] = 1.0
	sw = gauss(H, 2, 2, eps)
	return sw

			# 関数２
def hesse2(x, H, eps) :

	H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0
	H[0][1] = -400.0 * x[0]
	H[1][0] = H[0][1]
	H[1][1] = 200.0
	H[0][2] = 1.0
	H[0][3] = 0.0
	H[1][2] = 0.0
	H[1][3] = 1.0
	sw = gauss(H, 2, 2, eps)
	return sw

			# 関数３
def hesse3(x, H, eps) :

	x1      = 1.0 - x[1]
	x2      = 1.0 - x[1] * x[1]
	x3      = 1.0 - x[1] * x[1] * x[1]
	H[0][0] = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3
	H[0][1] = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 + 4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 + 6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3
	H[1][0] = H[0][1]
	H[1][1] = 2.0 * x[0] * x[0] + 4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] + 12.0 * x[0] * x[1] * (2.625 - x[0] * x3) + 18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1]
	H[0][2] = 1.0
	H[0][3] = 0.0
	H[1][2] = 0.0
	H[1][3] = 1.0
	sw = gauss(H, 2, 2, eps)
	return sw

			# データの入力
sw    = 0
line  = sys.stdin.readline()
s     = line.split()
fun   = int(s[1])
n     = int(s[3])
max   = int(s[5])
opt_1 = int(s[7])

line  = sys.stdin.readline()
s     = line.split()
eps   = float(s[1])
step  = float(s[3])

x     = np.empty(n, np.float)
dx    = np.empty(n, np.float)
H     = np.empty((n, 2*n), np.float)
y     = np.empty(1, np.float)

line  = sys.stdin.readline()
s     = line.split()
for i1 in range(0, n) :
	x[i1] = float(s[i1+1])
			# 実行
if fun == 1 :
	sw = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx1, dsnx1, hesse1)
elif fun == 2 :
	sw = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx2, dsnx2, hesse2)
else :
	sw = Newton(opt_1, max, n, eps, step, y, x, dx, H, snx3, dsnx3, hesse3)
			# 結果の出力
if sw < 0 :
	print("   収束しませんでした！", end="")
	if sw == -1 :
		print("（収束回数）")
	elif sw == -2 :
		print("（１次元最適化の区間）")
	else :
		print("（黄金分割法）")
else :
	print("   結果＝", end="")
	for i1 in range(0, n) :
		print(str(x[i1]) + " ", end="")
	print(" 最小値＝" + str(y[0]) + "  回数＝" + str(sw))

/******************************/
/* Newton法による最小値の計算 */
/*      coded by Y.Suganuma   */
/******************************/
using System;

class Program
{
	static void Main()
	{
		Test1 ts = new Test1();
	}
}

class Test1
{
	public Test1()
	{
					// データの入力
							// 1 行目
		string[] str = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
		int fun      = int.Parse(str[1]);
		int n        = int.Parse(str[3]);
		int max      = int.Parse(str[5]);
		int opt_1    = int.Parse(str[7]);
							// 2 行目
		str         = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
		double eps  = double.Parse(str[1]);
		double step = double.Parse(str[3]);
							// 3 行目
		double[][] H = new double [n][];
		double[] x   = new double [n];
		double[] dx  = new double [n];
		double y     = 0.0;
		str          = Console.ReadLine().Split(new char[] {' '}, StringSplitOptions.RemoveEmptyEntries);
		for (int i1 = 0; i1 < n; i1++) {
			x[i1] = double.Parse(str[i1+1]);
			H[i1] = new double [2*n];
		}
					// 実行
		int sw;
		if (fun == 1)
			sw = Newton(opt_1, max, n, eps, step, ref y, x, dx, H, snx1, dsnx1, hesse1);
		else if (fun == 2)
			sw = Newton(opt_1, max, n, eps, step, ref y, x, dx, H, snx2, dsnx2, hesse2);
		else
			sw = Newton(opt_1, max, n, eps, step, ref y, x, dx, H, snx3, dsnx3, hesse3);
					// 結果の出力
		if (sw < 0) {
			Console.Write("   収束しませんでした！");
			switch (sw) {
				case -1:
					Console.WriteLine("（収束回数）");
					break;
				case -2:
					Console.WriteLine("（１次元最適化の区間）");
					break;
				case -3:
					Console.WriteLine("（黄金分割法）");
					break;
			}
		}
		else {
			Console.Write("   結果＝");
			for (int i1 = 0; i1 < n; i1++)
				Console.Write(x[i1] + " ");
			Console.WriteLine(" 最小値＝" + y + "  回数＝" + sw);
		}
	}
			// 関数１の値の計算
	double snx1(double k, double[] x, double[] dx)
	{
		double[] w = new double [2];
		w[0]       = x[0] + k * dx[0];
		w[1]       = x[1] + k * dx[1];
		double x1  = w[0] - 1.0;
		double y1  = w[1] - 2.0;
		double y   = x1 * x1 + y1 * y1;
		return y;
	}
			// 関数２の値の計算
	double snx2(double k, double[] x, double[] dx)
	{
		double[] w = new double [2];
		w[0]       = x[0] + k * dx[0];
		w[1]       = x[1] + k * dx[1];
		double x1  = w[1] - w[0] * w[0];
		double y1  = 1.0 - w[0];
		double y   = 100.0 * x1 * x1 + y1 * y1;
		return y;
	}
			// 関数３の値の計算
	double snx3(double k, double[] x, double[] dx)
	{
		double[] w = new double [2];
		w[0]       = x[0] + k * dx[0];
		w[1]       = x[1] + k * dx[1];
		double x1  = 1.5 - w[0] * (1.0 - w[1]);
		double y1  = 2.25 - w[0] * (1.0 - w[1] * w[1]);
		double z1  = 2.625 - w[0] * (1.0 - w[1] * w[1] * w[1]);
		double y   = x1 * x1 + y1 * y1 + z1 * z1;
		return y;
	}
					// 関数１の微分
	int dsnx1(double[] x, double[] dx)
	{
		dx[0] = -2.0 * (x[0] - 1.0);
		dx[1] = -2.0 * (x[1] - 2.0);
		return 0;
	}
					// 関数２の微分
	int dsnx2(double[] x, double[] dx)
	{
		dx[0] = 400.0 * x[0] * (x[1] - x[0] * x[0]) + 2.0 * (1.0 - x[0]);
		dx[1] = -200.0 * (x[1] - x[0] * x[0]);
		return 0;
	}
					// 関数３の微分
	int dsnx3(double[] x, double[] dx)
	{
		dx[0] = 2.0 * (1.0 - x[1]) * (1.5 - x[0] * (1.0 - x[1])) +
                2.0 * (1.0 - x[1] * x[1]) * (2.25 - x[0] * (1.0 - x[1] * x[1])) +
                2.0 * (1.0 - x[1] * x[1] * x[1]) * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
		dx[1] = -2.0 * x[0] * (1.5 - x[0] * (1.0 - x[1])) -
                4.0 * x[0] * x[1] * (2.25 - x[0] * (1.0 - x[1] * x[1])) -
                6.0 * x[0] * x[1] * x[1] * (2.625 - x[0] * (1.0 - x[1] * x[1] * x[1]));
		return 0;
	}
					// Hesse行列の逆行列１
	int hesse1(double[] x, double[][] H, double eps)
	{
		H[0][0] = 2.0;
		H[0][1] = 0.0;
		H[1][0] = 0.0;
		H[1][1] = 2.0;
		H[0][2] = 1.0;
		H[0][3] = 0.0;
		H[1][2] = 0.0;
		H[1][3] = 1.0;
		int sw1 = gauss(H, 2, 2, eps);
		return sw1;
	}
					// Hesse行列の逆行列２
	int hesse2(double[] x, double[][] H, double eps)
	{
		H[0][0] = 400.0 * (3.0 * x[0] * x[0] - x[1]) + 2.0;
		H[0][1] = -400.0 * x[0];
		H[1][0] = H[0][1];
		H[1][1] = 200.0;
		H[0][2] = 1.0;
		H[0][3] = 0.0;
		H[1][2] = 0.0;
		H[1][3] = 1.0;
		int sw1 = gauss(H, 2, 2, eps);
		return sw1;
	}
					// Hesse行列の逆行列３
	int hesse3(double[] x, double[][] H, double eps)
	{
		double x1 = 1.0 - x[1];
		double x2 = 1.0 - x[1] * x[1];
		double x3 = 1.0 - x[1] * x[1] * x[1];
		H[0][0]   = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3;
		H[0][1]   = 2.0 * (1.5 - x[0] * x1) - 2.0 * x[0] * x1 +
	                4.0 * x[1] * (2.25 - x[0] * x2) - 4.0 * x[0] * x[1] * x2 +
	                6.0 * x[1] * x[1] * (2.625 - x[0] * x3) - 6.0 * x[0] * x[1] * x[1] * x3;
		H[1][0] = H[0][1];
		H[1][1] = 2.0 * x[0] * x[0] +
	              4.0 * x[0] * (2.25 - x[0] * x2) + 8.0 * x[0] * x[0] * x[1] * x[1] +
	              12.0 * x[0] * x[1] * (2.625 - x[0] * x3) +
	              18.0 * x[0] * x[0] * x[1] * x[1] * x[1] * x[1];
		H[0][2] = 1.0;
		H[0][3] = 0.0;
		H[1][2] = 0.0;
		H[1][3] = 1.0;
		int sw1 = gauss(H, 2, 2, eps);
		return sw1;
	}

	/********************************************************/
	/* Newton法                                             */
	/*      opt_1 : =0 : １次元最適化を行わない             */
	/*              =1 : １次元最適化を行う                 */
	/*      max : 最大繰り返し回数                          */
	/*      n : 次元                                        */
	/*      eps : 収束判定条件                              */
	/*      step : きざみ幅                                 */
	/*      y : 最小値                                      */
	/*      x : x(初期値と答え)                             */
	/*      dx : 関数の微分値                               */
	/*      H : Hesse行列の逆行列                           */
	/*      fn : 関数値を計算する関数                       */
	/*      dfn : 関数の微分値を計算する関数                */
	/*      hesse : Hesse行列の逆行列を計算する関数         */
	/*      return : >=0 : 正常終了(収束回数)               */
	/*               =-1 : 収束せず                         */
	/*               =-2 : １次元最適化の区間が求まらない   */
	/*               =-3 : 黄金分割法が失敗                 */
	/********************************************************/
	int Newton(int opt_1, int max, int n, double eps, double step, ref double y,
	           double[] x, double[] dx, double[][] H,
	           Func<double, double[], double[], double> fn,
	           Func<double[], double[], int> dfn,
	           Func<double[], double[][], double, int> hesse)
	{
		double[] wk = new double [n];
		double y1   = fn(0.0, x, dx);
		double y2   = 0.0;
		int count   = 0;
		int sw      = 0;

		while (count < max && sw == 0) {
					// 傾きの計算
			dfn(x, wk);
					// Hesse行列の逆行列の計算
			int sw1 = hesse(x, H, eps);
					// 収束していない
			if (sw1 == 0) {
						// 方向の計算
				count++;
				for (int i1 = 0; i1 < n; i1++) {
					dx[i1] = 0.0;
					for (int i2 = 0; i2 < n; i2++)
						dx[i1] += H[i1][n+i2] * wk[i2];
				}
						// １次元最適化を行わない
				if (opt_1 == 0) {
							// 新しい点
					for (int i1 = 0; i1 < n; i1++)
						x[i1] += dx[i1];
							// 新しい関数値
					y2 = fn(0.0, x, dx);
								// 関数値の変化が大きい
					if (Math.Abs(y2-y1) > eps)
						y1 = y2;
								// 収束（関数値の変化＜ｅｐｓ）
					else {
						sw = count;
						y  = y2;
					}
				}
						// １次元最適化を行う
				else {
							// 区間を決める
					sw1       = 0;
					double f1 = y1;
					double sp = step;
					double f2 = fn(sp, x, dx);
					if (f2 > f1)
						sp = -step;
					for (int i1 = 0; i1 < max && sw1 == 0; i1++) {
						f2 = fn(sp, x, dx);
						if (f2 > f1)
							sw1 = 1;
						else {
							sp *= 2.0;
							f1  = f2;
						}
					}
							// 区間が求まらない
					if (sw1 == 0)
						sw = -2;
							// 区間が求まった
					else {
								// 黄金分割法
						double k = gold(0.0, sp, eps, ref y2, ref sw1, max, x, dx, fn);
								// 黄金分割法が失敗
						if (sw1 < 0)
							sw = -3;
								// 黄金分割法が成功
						else {
									// 新しい点
							for (int i1 = 0; i1 < n; i1++)
								x[i1] += k * dx[i1];
										// 関数値の変化が大きい
							if (Math.Abs(y1-y2) > eps)
								y1 = y2;
										// 収束（関数値の変化＜ｅｐｓ）
							else {
								sw = count;
								y  = y2;
							}
						}
					}
				}
			}
					// 収束（傾き＜ｅｐｓ）
			else {
				sw = count;
				y  = y1;
			}
		}

		if (sw == 0)
			sw = -1;

		return sw;
	}

	/******************************************/
	/* 黄金分割法(与えられた方向での最小値)   */
	/*      a,b : 初期区間 a < b              */
	/*      eps : 許容誤差                    */
	/*      val : 関数値                      */
	/*      ind : 計算状況                    */
	/*              >= 0 : 正常終了(収束回数) */
	/*              = -1 : 収束せず           */
	/*      max : 最大試行回数                */
	/*      w : 位置                          */
	/*      dw : 傾きの成分                   */
	/*      fn : 関数値を計算する関数         */
	/*      return : 結果(w＋y＊dwのy)        */
	/******************************************/
	double gold(double a, double b, double eps, ref double val, ref int ind,
	            int max, double[] w, double[] dw,
	            Func<double, double[], double[], double> fn)
	{
					// 初期設定
		double tau = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0;
		double x1  = b - tau * (b - a);
		double x2  = a + tau * (b - a);
		double f1  = fn(x1, w, dw);
		double f2  = fn(x2, w, dw);
		int count  = 0;
		ind        = -1;
					// 計算
		while (count < max && ind < 0) {
			count += 1;
			if (f2 > f1) {
				if (Math.Abs(b-a) < eps) {
					ind = 0;
					x   = x1;
					val = f1;
				}
				else {
					b  = x2;
					x2 = x1;
					x1 = a + (1.0 - tau) * (b - a);
					f2 = f1;
					f1 = fn(x1, w, dw);
				}
			}
			else {
				if (Math.Abs(b-a) < eps) {
					ind = 0;
					x   = x2;
					val = f2;
					f1  = f2;
				}
				else {
					a  = x1;
					x1 = x2;
					x2 = b - (1.0 - tau) * (b - a);
					f1 = f2;
					f2 = fn(x2, w, dw);
				}
			}
		}
					// 収束した場合の処理
		if (ind == 0) {
			ind       = count;
			double fa = fn(a, w, dw);
			double fb = fn(b, w, dw);
			if (fb < fa) {
				a  = b;
				fa = fb;
			}
			if (fa < f1) {
				x   = a;
				val = fa;
			}
		}

		return x;
	}

	/*******************************************************/
	/* 線形連立方程式を解く（逆行列を求める）              */
	/*      w : 方程式の左辺及び右辺                       */
	/*      n : 方程式の数                                 */
	/*      m : 方程式の右辺の列の数                       */
	/*      eps : 逆行列の存在を判定する規準               */
	/*      return : =0 : 正常                             */
	/*               =1 : 逆行列が存在しない               */
	/*******************************************************/
	int gauss(double[][] w, int n, int m, double eps)
	{
		int ind = 0;
		int nm  = n + m;

		for (int i1 = 0; i1 < n && ind == 0; i1++) {

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

				for (int i2 = m1; i2 < nm; i2++)
					w[i1][i2] *= y1;

				for (int i2 = 0; i2 < n; i2++) {
					if (i2 != i1) {
						for (int i3 = m1; i3 < nm; i3++)
							w[i2][i3] -= w[i2][i1] * w[i1][i3];
					}
				}
			}
		}

		return ind;
	}
}

''''''''''''''''''''''''''''''
' Newton法による最小値の計算 '
'      coded by Y.Suganuma   '
''''''''''''''''''''''''''''''
Imports System.Text.RegularExpressions

Module Test

	Sub Main()
					' データの入力
							' 1 行目
		Dim MS As Regex      = New Regex("\s+") 
		Dim str() As String  = MS.Split(Console.ReadLine().Trim())
		Dim fun As Integer   = Integer.Parse(str(1))
		Dim n As Integer     = Integer.Parse(str(3))
		Dim max As Integer   = Integer.Parse(str(5))
		Dim opt_1 As Integer = Integer.Parse(str(7))
							' 2 行目
		str = MS.Split(Console.ReadLine().Trim())
		Dim eps As Double   = Double.Parse(str(1))
		Dim step1 As Double = Double.Parse(str(3))
							' 3 行目
		Dim H(n,2*n) As Double
		Dim x(n) As Double
		Dim dx(n) As Double
		Dim y As Double = 0.0
		str             = MS.Split(Console.ReadLine().Trim())
		For i1 As Integer = 0 To n-1
			x(i1) = Double.Parse(str(i1+1))
		Next
					' 関数１の値の計算（ラムダ式）
		Dim snx1 = Function(k0, x0, dx0) As Double
			Dim w(2) As Double
			w(0) = x0(0) + k0 * dx0(0)
			w(1) = x0(1) + k0 * dx0(1)
			Dim x1 As Double = w(0) - 1.0
			Dim y1 As Double = w(1) - 2.0
			Dim y0 As Double = x1 * x1 + y1 * y1
			Return y0
		End Function
					' 関数２の値の計算（ラムダ式）
		Dim snx2 = Function(k0, x0, dx0) As Double
			Dim w(2) As Double
			w(0) = x0(0) + k0 * dx0(0)
			w(1) = x0(1) + k0 * dx0(1)
			Dim x1 As Double = w(1) - w(0) * w(0)
			Dim y1 As Double = 1.0 - w(0)
			Dim y0 As Double = 100.0 * x1 * x1 + y1 * y1
			Return y0
		End Function
					' 関数３の値の計算（ラムダ式）
		Dim snx3 = Function(k0, x0, dx0) As Double
			Dim w(2) As Double
			w(0) = x0(0) + k0 * dx0(0)
			w(1) = x0(1) + k0 * dx0(1)
			Dim x1 As Double = 1.5 - w(0) * (1.0 - w(1))
			Dim y1 As Double = 2.25 - w(0) * (1.0 - w(1) * w(1))
			Dim z1 As Double = 2.625 - w(0) * (1.0 - w(1) * w(1) * w(1))
			Dim y0 As Double = x1 * x1 + y1 * y1 + z1 * z1
			Return y0
		End Function
					' 関数１の微分（ラムダ式）
		Dim dsnx1 = Function(x0, dx0) As Integer
			dx0(0) = -2.0 * (x0(0) - 1.0)
			dx0(1) = -2.0 * (x0(1) - 2.0)
			Return 0
		End Function
					' 関数２の微分（ラムダ式）
		Dim dsnx2 = Function(x0, dx0) As Integer
			dx0(0) = 400.0 * x0(0) * (x0(1) - x0(0) * x0(0)) + 2.0 * (1.0 - x0(0))
			dx0(1) = -200.0 * (x0(1) - x0(0) * x0(0))
			Return 0
		End Function
					' 関数３の微分（ラムダ式）
		Dim dsnx3 = Function(x0, dx0) As Integer
			dx0(0) = 2.0 * (1.0 - x0(1)) * (1.5 - x0(0) * (1.0 - x0(1))) +
	                 2.0 * (1.0 - x0(1) * x0(1)) * (2.25 - x0(0) * (1.0 - x0(1) * x0(1))) +
	                 2.0 * (1.0 - x0(1) * x0(1) * x0(1)) * (2.625 - x0(0) * (1.0 - x0(1) * x0(1) * x0(1)))
			dx0(1) = -2.0 * x0(0) * (1.5 - x0(0) * (1.0 - x0(1))) -
	                 4.0 * x0(0) * x0(1) * (2.25 - x0(0) * (1.0 - x0(1) * x0(1))) -
	                 6.0 * x0(0) * x0(1) * x0(1) * (2.625 - x0(0) * (1.0 - x0(1) * x0(1) * x0(1)))
			Return 0
		End Function
					' Hesse行列の逆行列１（ラムダ式）
		Dim hesse1 = Function(x0, H0, eps0)
			H0(0,0) = 2.0
			H0(0,1) = 0.0
			H0(1,0) = 0.0
			H0(1,1) = 2.0
			H0(0,2) = 1.0
			H0(0,3) = 0.0
			H0(1,2) = 0.0
			H0(1,3) = 1.0
			Dim sw1 As Integer = gauss(H0, 2, 2, eps0)
			Return sw1
		End Function
					' Hesse行列の逆行列２（ラムダ式）
		Dim hesse2 = Function(x0, H0, eps0)
			H0(0,0) = 400.0 * (3.0 * x0(0) * x0(0) - x0(1)) + 2.0
			H0(0,1) = -400.0 * x0(0)
			H0(1,0) = H0(0,1)
			H0(1,1) = 200.0
			H0(0,2) = 1.0
			H0(0,3) = 0.0
			H0(1,2) = 0.0
			H0(1,3) = 1.0
			Dim sw1 As Integer = gauss(H0, 2, 2, eps0)
			Return sw1
		End Function
					' Hesse行列の逆行列３（ラムダ式）
		Dim hesse3 = Function(x0, H0, eps0)
			Dim x1 As Double = 1.0 - x0(1)
			Dim x2 As Double = 1.0 - x0(1) * x0(1)
			Dim x3 As Double = 1.0 - x0(1) * x0(1) * x0(1)
			H0(0,0) = 2.0 * x1 * x1 + 2.0 * x2 * x2 + 2.0 * x3 * x3
			H0(0,1) = 2.0 * (1.5 - x0(0) * x1) - 2.0 * x0(0) * x1 +
			            4.0 * x0(1) * (2.25 - x0(0) * x2) - 4.0 * x0(0) * x0(1) * x2 +
			            6.0 * x0(1) * x0(1) * (2.625 - x0(0) * x3) - 6.0 * x0(0) * x0(1) * x0(1) * x3
			H0(1,0) = H0(0,1)
			H0(1,1) = 2.0 * x0(0) * x0(0) +
			          4.0 * x0(0) * (2.25 - x0(0) * x2) + 8.0 * x0(0) * x0(0) * x0(1) * x0(1) +
			          12.0 * x0(0) * x0(1) * (2.625 - x0(0) * x3) +
			          18.0 * x0(0) * x0(0) * x0(1) * x0(1) * x0(1) * x0(1)
			H0(0,2) = 1.0
			H0(0,3) = 0.0
			H0(1,2) = 0.0
			H0(1,3) = 1.0
			Dim sw1 As Integer = gauss(H0, 2, 2, eps0)
			Return sw1
		End Function
					' 実行
		Dim sw As Integer
		If fun = 1
			sw = Newton(opt_1, max, n, eps, step1, y, x, dx, H, snx1, dsnx1, hesse1)
		ElseIf fun = 2
			sw = Newton(opt_1, max, n, eps, step1, y, x, dx, H, snx2, dsnx2, hesse2)
		Else
			sw = Newton(opt_1, max, n, eps, step1, y, x, dx, H, snx3, dsnx3, hesse3)
		End If
					' 結果の出力
		If sw < 0
			Console.Write("   収束しませんでした！")
			If sw = -1
				Console.WriteLine("（収束回数）")
			ElseIf sw = -2
				Console.WriteLine("（１次元最適化の区間）")
			ElseIf sw = -3
				Console.WriteLine("（黄金分割法）")
			End If
		Else
			Console.Write("   結果＝")
			For i1 As Integer = 0 To n-1
				Console.Write(x(i1) & " ")
			Next
			Console.WriteLine(" 最小値＝" & y & "  回数＝" & sw)
		End If

	End Sub

	''''''''''''''''''''''''''''''''''''''''''''''''''''''''
	' Newton法                                             '
	'      opt_1 : =0 : １次元最適化を行わない             '
	'              =1 : １次元最適化を行う                 '
	'      max : 最大繰り返し回数                          '
	'      n : 次元                                        '
	'      eps : 収束判定条件                              '
	'      step1 : きざみ幅                                '
	'      y : 最小値                                      '
	'      x : x(初期値と答え)                             '
	'      dx : 関数の微分値                               '
	'      H : Hesse行列の逆行列                           '
	'      fn : 関数値を計算する関数                       '
	'      dfn : 関数の微分値を計算する関数                '
	'      hesse : Hesse行列の逆行列を計算する関数         '
	'      return : >=0 : 正常終了(収束回数)               '
	'               =-1 : 収束せず                         '
	'               =-2 : １次元最適化の区間が求まらない   '
	'               =-3 : 黄金分割法が失敗                 '
	''''''''''''''''''''''''''''''''''''''''''''''''''''''''
	Function Newton(opt_1 As Integer, max As Integer, n As Integer, eps As Double,
	                step1 As Double, ByRef y As Double, x() As Double, dx() As Double,
	                H(,) As Double, fn As Func(Of Double, Double(), Double(), Double),
	                dfn As Func(Of Double(), Double(), Integer),
	                hesse As Func(Of Double(), Double(,), Double, Integer))

		Dim wk(n) As Double
		Dim y1 As Double     = fn(0.0, x, dx)
		Dim y2 As Double     = 0.0
		Dim count As Integer = 0
		Dim sw As Integer    = 0

		Do while count < max and sw = 0
					' 傾きの計算
			dfn(x, wk)
					' Hesse行列の逆行列の計算
			Dim sw1 As Integer = hesse(x, H, eps)
					' 収束していない
			If sw1 = 0
						' 方向の計算
				count += 1
				For i1 As Integer = 0 To n-1
					dx(i1) = 0.0
					For i2 As Integer = 0 To n-1
						dx(i1) += H(i1,n+i2) * wk(i2)
					Next
				Next
						' １次元最適化を行わない
				If opt_1 = 0
							' 新しい点
					For i1 As Integer = 0 To n-1
						x(i1) += dx(i1)
					Next
							' 新しい関数値
					y2 = fn(0.0, x, dx)
								' 関数値の変化が大きい
					If Math.Abs(y2-y1) > eps
						y1 = y2
								' 収束（関数値の変化＜ｅｐｓ）
					Else
						sw = count
						y  = y2
					End If
						' １次元最適化を行う
				Else
							' 区間を決める
					sw1 = 0
					Dim f1 As Double = y1
					Dim sp As Double = step1
					Dim f2 As Double = fn(sp, x, dx)
					If f2 > f1
						sp = -step1
					End If
					Dim ii As Integer = 0
					Do While ii < max and sw1 = 0
						f2 = fn(sp, x, dx)
						If f2 > f1
							sw1 = 1
						Else
							sp *= 2.0
							f1  = f2
						End If
						ii += 1
					Loop
							' 区間が求まらない
					If sw1 = 0
						sw = -2
							' 区間が求まった
					Else
								' 黄金分割法
						Dim k As Double = gold(0.0, sp, eps, y2, sw1, max, x, dx, fn)
								' 黄金分割法が失敗
						If sw1 < 0
							sw = -3
								' 黄金分割法が成功
						Else
									' 新しい点
							For i1 As Integer = 0 To n-1
								x(i1) += k * dx(i1)
							Next
										' 関数値の変化が大きい
							If Math.Abs(y1-y2) > eps
								y1 = y2
										' 収束（関数値の変化＜ｅｐｓ）
							Else
								sw = count
								y  = y2
							End If
						End If
					End If
				End If
					' 収束（傾き＜ｅｐｓ）
			Else
				sw = count
				y  = y1
			End If
		Loop

		If sw = 0
			sw = -1
		End If

		Return sw

	End Function

	''''''''''''''''''''''''''''''''''''''''''
	' 黄金分割法(与えられた方向での最小値)   '
	'      a,b : 初期区間 a < b              '
	'      eps : 許容誤差                    '
	'      val : 関数値                      '
	'      ind : 計算状況                    '
	'              >= 0 : 正常終了(収束回数) '
	'              = -1 : 収束せず           '
	'      max : 最大試行回数                '
	'      w : 位置                          '
	'      dw : 傾きの成分                   '
	'      fn : 関数値を計算する関数         '
	'      return : 結果(w＋y＊dwのy)        '
	''''''''''''''''''''''''''''''''''''''''''
	Function gold(a As Double, b As Double, eps As Double, ByRef val As Double,
	              ByRef ind As Integer, max As Integer, w() As Double, dw() As Double,
	              fn As Func(Of Double, Double(), Double(), Double))

					' 初期設定
		Dim tau As Double    = (Math.Sqrt(5.0) - 1.0) / 2.0, x = 0.0
		Dim x1 As Double     = b - tau * (b - a)
		Dim x2 As Double     = a + tau * (b - a)
		Dim f1 As Double     = fn(x1, w, dw)
		Dim f2 As Double     = fn(x2, w, dw)
		Dim count As Integer = 0
		ind = -1
					' 計算
		Do While count < max and ind < 0
			count += 1
			If f2 > f1
				If Math.Abs(b-a) < eps
					ind = 0
					x   = x1
					val = f1
				Else
					b  = x2
					x2 = x1
					x1 = a + (1.0 - tau) * (b - a)
					f2 = f1
					f1 = fn(x1, w, dw)
				End If
			Else
				If Math.Abs(b-a) < eps
					ind = 0
					x   = x2
					val = f2
					f1  = f2
				Else
					a  = x1
					x1 = x2
					x2 = b - (1.0 - tau) * (b - a)
					f1 = f2
					f2 = fn(x2, w, dw)
				End If
			End If
		Loop
					' 収束した場合の処理
		If ind = 0
			ind = count
			Dim fa As Double = fn(a, w, dw)
			Dim fb As Double = fn(b, w, dw)
			If fb < fa
				a  = b
				fa = fb
			End If
			If fa < f1
				x   = a
				val = fa
			End If
		End If

		Return x

	End Function

	''''''''''''''''''''''''''''''''''''''''''
	' 線形連立方程式を解く（逆行列を求める） '
	'      w : 方程式の左辺及び右辺          '
	'      n : 方程式の数                    '
	'      m : 方程式の右辺の列の数          '
	'      eps : 逆行列の存在を判定する規準  '
	'      return : =0 : 正常                '
	'               =1 : 逆行列が存在しない  '
	''''''''''''''''''''''''''''''''''''''''''
	Function gauss(w(,) As Double, n As Integer, m As Integer, eps As Double) As Integer

		Dim ind As Integer = 0
		Dim nm As Integer  = n + m

		Dim i1 As Integer = 0
		Do While i1 < n and ind = 0

			Dim y1 As Double  = 0.0
			Dim m1 As Integer = i1 + 1
			Dim m2 As Integer = 0
						' ピボット要素の選択
			For i2 As Integer = i1 To n-1
				Dim y2 As Double = Math.Abs(w(i2,i1))
				If y1 < y2
					y1 = y2
					m2 = i2
				End If
			Next
						' 逆行列が存在しない
			If y1 < eps
				ind = 1
						' 逆行列が存在する
			Else
							' 行の入れ替え
				For i2 As Integer = i1 To nm-1
					y1       = w(i1,i2)
					w(i1,i2) = w(m2,i2)
					w(m2,i2) = y1
				Next
							' 掃き出し操作
				y1 = 1.0 / w(i1,i1)

				For i2 As Integer = m1 To nm-1
					w(i1,i2) *= y1
				Next

				For i2 As Integer = 0 To n-1
					If i2 <> i1
						For i3 As Integer = m1 To nm-1
							w(i2,i3) -= w(i2,i1) * w(i1,i3)
						Next
					End If
				Next
			End If
			i1 += 1
		Loop

		Return ind

	End Function

End Module