/********************************/ /* 最小二乗法(多項式近似) */ /* coded by Y.Suganuma */ /********************************/ import java.io.*; public class Test { public static void main(String args[]) throws IOException { double x[], y[], z[]; int i1, m, n, sw; String str; Data dt; BufferedReader in = new BufferedReader(new InputStreamReader(System.in)); // 多項式の次数とデータの数 str = in.readLine(); dt = new Data(str); dt.next(); m = Integer.parseInt(dt.data); dt.next(); n = Integer.parseInt(dt.data); x = new double [n]; y = new double [n]; z = new double [m+1]; // データ for (i1 = 0; i1 < n; i1++) { str = in.readLine(); dt = new Data(str); dt.next(); x[i1] = Double.parseDouble(dt.data); dt.next(); y[i1] = Double.parseDouble(dt.data); } sw = App.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 逆行列を求めることができませんでした"); } } /**************************** /* 科学技術系算用の手法 */ /****************************/ class App { /*****************************************/ /* 最初2乗法 */ /* 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 = App.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); } } -------------------------class Data------------------------- /**********************************************/ /* 1行のデータから個々のデータを取り出す */ /* coded by Y.Suganuma */ /**********************************************/ public class Data { private int len, start, end; private String line; String data; // コンストラクタ Data (String l) { line = l; len = line.length(); start = 0; end = 0; data = new String (); // 次のデータ } // 次のデータの取り出し // =0 : 成功,=-1 : 失敗 int next() { int k = 0; if (start < len) { int sw = 0; while (sw == 0) { end = line.indexOf(" ", start); if (end >= 0) { if ((end - start) > 0) { data = line.substring(start, end); sw = 1; } } else { end = len; sw = 1; if ((end - start) > 0) data = line.substring(start, end); else k = -1; } start = end + 1; } } else k = -1; return k; } } -------------------------データ例1------------------------- 1 10 0 2.450 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------------------------- 2 6 10 28.2 15 47.0 20 44.4 26 32.8 32 20.8 40 0.8