aa's homework

poly.c

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

int SolveLEG(double **CoefMatrix, int ElmtCount, double **Solve)
{
    //
    // 解多元线性方程组函数 Solve Linear Equation Group.
    //

    //
    // [in]         方程组系数矩阵  CoefMatrix[][]
    // [in]         方程组元数      ElmtCount
    // [out]        方程组的解      Solve[] *
    //
    // [return]     返回方程组是否有唯一解
    //

    // 计数器
    int i, j, _i;

    // 标志，用于判断秩
    int flag, isFullRank;

    // 初等变换中使用的比值u
    double u;

    //
    // 一、检查矩阵的行列
    //
    isFullRank = 1;
    for (i = 0; i < ElmtCount; i++) {
	flag = 0;
	for (j = 0; j < ElmtCount; j++) {
	    flag = flag || (CoefMatrix[i][j] != 0);
	    if (flag == 1)
		break;
	}
	isFullRank = isFullRank && flag;
	if (isFullRank == 0)
	    break;
    }
    for (i = 0; i < ElmtCount; i++) {
	flag = 0;
	for (j = 0; j < ElmtCount; j++) {
	    flag = flag || (CoefMatrix[j][i] != 0);
	    if (flag == 1)
		break;
	}
	isFullRank = isFullRank && flag;
	if (isFullRank == 0)
	    break;
    }
    if (!isFullRank) {
	*Solve = NULL;
	return 0;
    }
    //
    // 二、逐行进行操作
    //
    for (i = 0; i < ElmtCount; i++) {
    //
    // 若主对角线上元素是0
    //
	if (CoefMatrix[i][i] == 0) {
        //
        // 搜寻可以进行交换的行，_i 目标行标，i 需交换的行标
        //
	    for (_i = i + 1; _i < ElmtCount; _i++) {
        // 判断可交换性
		if (CoefMatrix[_i][i] == 0)
		    continue;

        // 逐列交换，j 列标
		for (j = 0; j <= ElmtCount; j++) {
		    double Temp;
		    Temp = CoefMatrix[i][j];
		    CoefMatrix[i][j] = CoefMatrix[_i][j];
		    CoefMatrix[_i][j] = Temp;
		}
	    }
	}
    //
    // 若主对角线上元素仍为0
    //
	if (CoefMatrix[i][i] == 0) {
	    *Solve = NULL;
	    return 0;
	}
	//
	// 初等变换
	// 作用：逐行进行计算，_i 为循环计数器，使第i列中除了第i行以外的元素变为0。
	//
	for (_i = 0; _i < ElmtCount; _i++) {
	    if (_i == i)
		continue;

	    // 计算比值
	    u = -(CoefMatrix[_i][i] / CoefMatrix[i][i]);

	    // 对第_i行的逐列进行初等变换计算，j 列标
	    for (j = i; j <= ElmtCount; j++)
		CoefMatrix[_i][j] += (CoefMatrix[i][j] * u);
	}
    }

    //
    // 三、变为单位矩阵
    //
    for (i = 0; i < ElmtCount; i++) {
	u = CoefMatrix[i][i];
	CoefMatrix[i][i] = 1;
	CoefMatrix[i][ElmtCount] /= u;
    }

    //
    // 四、得到方程组的解
    //

    //
    // 此时，矩阵最后一列为方程组的解

    // 为方程的解开辟内存空间
    *Solve = malloc(sizeof(double) * ElmtCount);

    // 将最后结果存入Solve数组
    for (i = 0; i < ElmtCount; i++)
	(*Solve)[i] = CoefMatrix[i][ElmtCount];

    // 返回1，方程组有唯一解
    return 1;
}





// Power 是要拟合的次数
// Known_x[]是已知的x序列
// Known_y[]是已知的y序列
// NLEGCoefMatrix[][]就是法方程组的系数矩阵了, 求到了过后就看楼上吧.
// 注意：NLEGCoefMatrix的大小是 Power+1 * Power+2 ，不要声明成数组，用二级指针

double **CalcCoef(double *Known_x, double *Known_y, int n, int Power)
{

    int i, j, k;
    double **NLEGCoefMatrix;

    NLEGCoefMatrix = malloc(sizeof(double*) * (Power+1));
    for(i = 0; i <= Power; i++) {
        NLEGCoefMatrix[i] = malloc(sizeof(double) * (Power+2));
    }
    
    for (i = 0; i <= Power; i++) {
        for (j = 0; j <= Power; j++) {
        double Sum = 0;
        for (k = 0; k < n; k++) {
            if (i + j == 0)
            Sum = Sum + 1;
            else
            Sum = Sum + pow(Known_x[k], (double) (i + j));
        }
        NLEGCoefMatrix[i][j] = Sum;
        }
    }

    for (i = 0; i <= Power; i++) {
        double Sum = 0;
        for (k = 0; k < n; k++) {
        if (i == 0)
            Sum = Sum + Known_y[k];
        else
            Sum = Sum + Known_y[k] * pow(Known_x[k], (double) i);
        }
        NLEGCoefMatrix[i][Power + 1] = Sum;
    }

    return NLEGCoefMatrix;

}

int main(int argc, char *argv[])
{
    int power, npoints, i;
    double *X, *Y;
    double *solve;

    printf("Enter the power: ");
    scanf("%d", &power);
    
    printf("Enter the number of points: ");
    scanf("%d", &npoints);
    
    X = malloc(sizeof(double) * npoints);
    Y = malloc(sizeof(double) * npoints);

    for(i=0; i<npoints; i++) {
        printf("Enter X%d: ", i+1);
        scanf("%lf", &X[i]);
        printf("Enter Y%d: ", i+1);
        scanf("%lf", &Y[i]);
    }

    SolveLEG(CalcCoef(X, Y, npoints, power), power+1, &solve);
    for(i=0; i<power+1; i++) {
        printf("%.3lf * x^%d\n", solve[i], i);
    }

    return 0;
}