最小二乗法（多項式近似）

/****************************/
/* 最小二乗法（多項式近似） */
/*      coded by Y.Suganuma */
/****************************/
import java.io.*;
import java.util.StringTokenizer;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double x[], y[], z[];
		int i1, m, n, sw;
		StringTokenizer str;
		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
					// 多項式の次数とデータの数
		str = new StringTokenizer(in.readLine(), " ");
		m   = Integer.parseInt(str.nextToken());
		n   = Integer.parseInt(str.nextToken());

		x = new double [n];
		y = new double [n];
		z = new double [m+1];
					// データ
		for (i1 = 0; i1 < n; i1++) {
			str = new StringTokenizer(in.readLine(), " ");
			x[i1] = Double.parseDouble(str.nextToken());
			y[i1] = Double.parseDouble(str.nextToken());
		}

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

		if (sw == 0) {
			System.out.println("結果");
			for (i1 = 0; i1 < m+1; i1++)
				System.out.println("   " + (m-i1) + " 次の係数 " + z[i1]);
		}
		else
			System.out.println("***error  逆行列を求めることができませんでした");
	}

	/*************************************/
	/* 最初２乗法                        */
	/*      m : 多項式の次数             */
	/*      n : データの数               */
	/*      x,y : データ                 */
	/*      z : 多項式の係数（高次から） */
	/*      return : =0 : 正常           */
	/*               =1 : エラー         */
	/*************************************/
	static int least(int m, int n, double x[], double y[], double z[])
	{
		double A[][], w[][], x1, x2;
		int i1, i2, i3, sw = 0;

		m++;
		w = new double [m][m+1];
		A = new double [n][m];

		for (i1 = 0; i1 < n; i1++) {
			A[i1][m-2] = x[i1];
			A[i1][m-1] = 1.0;
			x1         = A[i1][m-2];
			x2         = x1;
			for (i2 = m-3; i2 >= 0; i2--) {
				x2        *= x1;
				A[i1][i2]  = x2;
			}
		}

		for (i1 = 0; i1 < m; i1++) {
			for (i2 = 0; i2 < m; i2++) {
				w[i1][i2] = 0.0;
				for (i3 = 0; i3 < n; i3++)
					w[i1][i2] += A[i3][i1] * A[i3][i2];
			}
		}

		for (i1 = 0; i1 < m; i1++) {
			w[i1][m] = 0.0;
			for (i2 = 0; i2 < n; i2++)
				w[i1][m] += A[i2][i1] * y[i2];
		}

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

		if (sw == 0) {
			for (i1 = 0; i1 < m; i1++)
				z[i1] = w[i1][m];
		}
		else
			sw = 1;

		return sw;
	}

	/*******************************************************/
	/* 線形連立方程式を解く（逆行列を求める）              */
	/*      w : 方程式の左辺及び右辺                       */
	/*      n : 方程式の数                                 */
	/*      m : 方程式の右辺の列の数                       */
	/*      eps : 逆行列の存在を判定する規準               */
	/*      return : =0 : 正常                             */
	/*               =1 : 逆行列が存在しない               */
	/*******************************************************/
	static 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);
	}
}

-----------データ例１（コメント部分を除いて下さい）---------
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