// 函数功能：	求解拉格朗日三阶微分矩阵在tk处的值在第i个配置点上的分量
// 输入：point：n * 1维矩阵，是插值点的时间维度，n为插值节点数目（升序）
//       k0：为当前点的标号
//       j：以1开始的标号
// 输出：f当前点的三阶导数值
// 方法：拉格朗日微分矩阵的基本公式及其推导		
// 编写：		2021/08/16

#include "GPM.h"

double D3Lagrange_point(mat point1, int k0, int j) {

	// 初始化
	int N = (size(point1))(0);
	double b = 1.0, a_sum = 0;
	double tj, ti, a, tk, b1 ,c1;

	// 分母计算, 分母与标号L无关
	for (int L = 0; L < N; L++) {
		if (L != j) {
			tj = (point1)(L, 0);
			ti = (point1)(j, 0);
			b = b*(ti - tj);
		}
	}

	// 分子计算，分子与标号有关
	for (int i0 = 0; i0 < N; i0++) {
		if (i0 == j) {
			a = 0;
		}
		else {
			a = 0;
			for (int m = 0; m < N; m++) {
				if (m == j || m == i0) {}
				else {
					b1 = 0;
					for (int n = 0; n < N; n++) {
						if (n == j || n == i0 || n == m) {}
						else {
							c1 = 1;
							for (int L = 0; L < N; L++) {
								if (L == j || L == i0 || L == m || L == n) {}
								else {
									tj = (point1)(L, 0);
									tk = (point1)(k0, 0);
									c1 = c1*(tk - tj);
								}
							}
							b1 = b1 + c1;
						}
					}
					a = a + b1;
				}
			}
		}
		a_sum = a_sum + a;
	}

	// 相除得结果
	double f = a_sum / b;

	return f;
}