using System;

namespace AlgorithmMain
{
 /// <summary>
 /// 最小二乘法拟合二元多次曲线
 /// </summary>
 public class MultiLineFit
 {
 ///<summary>
 /// 二元多次线性方程拟合曲线
 ///</summary>
 ///<param name="arrX">已知点的x坐标集合</param>
 ///<param name="arrY">已知点的y坐标集合</param>
 ///<param name="length">已知点的个数</param>
 ///<param name="dimension">方程的最高次数</param>
 ///<returns>拟合方程的系数矩阵</returns>
 public static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)
 {
 int n = dimension + 1; //dimension次方程需要求 dimension+1个 系数
 double[,] Guass = new double[n, n + 1]; //高斯矩阵 例如:y=a0+a1*x+a2*x*x
 for (int i = 0; i < n; i++)
 {
 int j;
 for (j = 0; j < n; j++)
 {
 Guass[i, j] = SumArr(arrX, j + i, length);
 }
 Guass[i, j] = SumArr(arrX, i, arrY, 1, length);
 }
 return ComputGauss(Guass, n);
 }

 public static double SumArr(double[] arr, int n, int length) // 求数组的元素的n次方的和
 {
 double s = 0;
 for (int i = 0; i < length; i++)
 {
 if (arr[i] != 0 || n != 0)
 s = s + Math.Pow(arr[i], n);
 else
 s = s + 1;
 }
 return s;
 }
 public static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length)
 {
 double s = 0;
 for (int i = 0; i < length; i++)
 {
 if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
 s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2);
 else
 s = s + 1;
 }
 return s;
 }

 public static double[] ComputGauss(double[,] Guass, int n)
 {
 int i, j;
 int k, m;
 double temp;
 double max;
 double s;
 double[] x = new double[n];
 for (i = 0; i < n; i++) x[i] = 0.0;//初始化
 for (j = 0; j < n; j++)
 {
 max = 0;
 k = j;
 for (i = j; i < n; i++)
 {
 if (Math.Abs(Guass[i, j]) > max)
 {
 max = Guass[i, j];
 k = i;
 }
 }
 if (k != j)
 {
 for (m = j; m < n + 1; m++)
 {
 temp = Guass[j, m];
 Guass[j, m] = Guass[k, m];
 Guass[k, m] = temp;

 }
 }
 if (0 == max)
 {
 // "此线性方程为奇异线性方程" 
 return x;
 }

 for (i = j + 1; i < n; i++)
 {
 s = Guass[i, j];
 for (m = j; m < n + 1; m++)
 {
 Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]);
 }
 }
 }//结束for (j=0;j<n;j++)
 for (i = n - 1; i >= 0; i--)
 {
 s = 0;
 for (j = i + 1; j < n; j++)
 {
 s = s + Guass[i, j] * x[j];
 }
 x[i] = (Guass[i, n] - s) / Guass[i, i];
 }
 return x;
 }//返回值是函数的系数
 }
}