Photo by
akari_n
頑張ってラグランジュ補間多項式をC言語で①
不定期投稿で有名な関口です。
お久しぶりです。新年の祝いからしばらく経ち、私はもう通常通りの気持ちです。人によってはまだ休みが恋しいかもですね！
私の通ってる大学では後期も終盤で大学でのテストが近づいてきました。
最近、とある授業の課題でラグランジュ補間多項式から式を求める課題が出ました。
この授業ではC言語によるコーディングもあったので、ウキウキで書き始めたら…
難しい…
特に3次以上ではx1x2,x2x3,x1x3と組み合わせのように積をとるので、そこをどうするかで悩みました。プログラミングでもないのに難しすぎないか？
よく見ると手計算でした。
急にめんどいことするやん...
ただ、半分くらいコーディングしていたので、完成させたい！！
完成させましたよ...見てほしい。誰か褒めて...
２次関数専用のコード
諦めて3点から2次関数を求めるに特化した関数がまず最初に完成しました。
最初に作った良し悪しを残したかったので、次のN次版より完成度はかなり低いかな。
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#define N 3
/*プロトタイプ宣言*/
void compute_LP_Factor(int,double *x,double *f_l);
double x_product(double *x);
double x_sum(double *x);
double denominator(int num,double *x);
void output_result(double (*P)[]);
int main(void){
double x[N] = {1,2,3}; //x座標
double f[N] = {0,3,8}; //y座標
double P[N][N] = {0}; //各次数の係数
printf("\n%d次のラグランジュ補項多項式を計算します。\n",N-1);
for(int j=0;j<N;j++){
double f_l[N] = {0}; //y座標*f_l = 各項の係数
compute_LP_Factor(j,x,f_l); //ラグランジュ多項式因子(各項の係数)の計算（使いまわし、上書き）
for(int i=0;i<N;i++){
P[j][i] = f[j]*f_l[i];
}
}
for(int i=0;i<N;i++){
printf("|");
for(int j=0;j<N;j++){
if(P[i][j] >= 0) printf(" ");
printf("%.0f ",P[i][j]);
}
printf("||x^%d|\n",(N-1)-i);
}
output_result(P);
return 0;
}
///////////////////////////////////////////////////////////////////////
void compute_LP_Factor(int j,double *x,double *f_l){//係数計算
/*分子*/
//定数項はxの総積*-1
f_l[N-1] = pow(-1,N-1) * x_product(x)/x[j];
//最高次数の係数は１
f_l[0] = 1;
//1次の項はn個の組み合わせのxの総和*-1
f_l[1] = -1 * (x_sum(x) - x[j]);
/*分母に相当する値で割った値を返す*/
for(int i=0;i<N;i++){
f_l[i] /= denominator(j,x);
}
}
double x_product(double *x){
double product = 1.0;
for(int i=0;i<N;i++){
product *= x[i];
}
return product;
}
double x_sum(double *x){
double sum = 0.0;
for(int i=0;i<N;i++){
sum += x[i];
}
return sum;
}
double denominator(int num,double *x){
double product = 1.0;
for(int i = 0; i < N; i++) {
if(i != num){
product *= (x[num] - x[i]);
}
}
return product;
}
void output_result(double (*P)[N]){
double result[N] = {0}; // 最終的な係数を格納する配列
// 各基底多項式の係数の和を計算
for(int j = 0; j < N; j++){
for(int i = 0; i < N; i++){
result[j] += P[i][j];
}
}
// 結果の出力
printf("\n求まった%d次関数\nf(x)= ",N-1);
for(int i=0;i<N;i++){
if(result[i] == 0) continue;
if(i != 0){
if(result[i] > 0) printf(" + ");
else printf(" - ");
}
else if(result[i] < 0) printf("-");
if(fabs(result[i]) != 1 || i == N-1) printf("%.0f",fabs(result[i]));
if(i == 1) printf("x");
else if(i < N-1) printf("x^%d",(N-1)-i);
}
printf("\n\n");
}３次以降になるとxの組み合わせが絡むので、ゲロむずかったです。
最初に一気にN次まで対応させようとしたせいで、余計に時間かかりました。
もがきに３時間、諦めて２次限定で完成に１時間って感じです。
事前に座標を初期化する必要があるのがねー…ちょっとって感じ。
N次関数に対応！！
一回寝たあと、数時間かけ試行錯誤した結果やっとやり遂げました…
こちらが200行近いコードです。
#include <stdlib.h>
#include<stdio.h>
#include<math.h>
#define EPSILON 1e-10
#define N 3 //点の個数
/*プロトタイプ宣言*/
void compute_LP_Factor(int,double *x,double *f_l);
double x_product(double *x,int j);
double x_sum(double *x);
double x_other(int j, int num, double *x);
void generate_combination(int n,int j,int *element, int k, int start, int count,double *x,double *sum);
double combination_product(int k,double *x,int *element);
double denominator(int size,double *x);
void output_result(double (*P)[]);
int main(void){
double x[N]; //x座標
double f[N]; //y座標
double P[N][N] = {0}; //各次数の係数
printf("\n%d次のラグランジュ補項多項式を計算します。\n",N-1);
printf("\t x y\n");
/*座標入力案内所*/
for(int i=0;i<N;i++){
printf("%dつ目の座標：",i+1);
scanf("%lf %lf",&x[i],&f[i]);
}
/*ラグランジュ補間公式への処理*/
for(int j=0;j<N;j++){
double f_l[N] = {0}; //y座標*f_l = 各項の係数
compute_LP_Factor(j,x,f_l); //ラグランジュ多項式因子(各項の係数)の計算（使いまわし、上書き）
for(int i=0;i<N;i++){
P[j][i] = f[j]*f_l[i];
}
}
output_result(P);
return 0;
}
///////////////////////////////////////////////////////////////////////
void compute_LP_Factor(int j,double *x,double *f_l){//係数計算
/*分子*/
//最高次数の係数は１
f_l[0] = 1;
//n-1次の項はn個の組み合わせのxの総和*-1
f_l[1] = -1 * (x_sum(x) - x[j]);
// 3次以上の場合のみ組み合わせの計算が必要
if(N > 3){
for(int i=2;i<N-1;i++){
f_l[i] = pow(-1,i) * x_other(j,i,x); //N-i次
}
}
//定数項はxの総積*-1
f_l[N-1] = pow(-1,N-1) * x_product(x,j);
/*分母に相当する値で割った値を返す*/
for(int i=0;i<N;i++){
f_l[i] /= denominator(j,x);
}
}
/*総積の計算*/
double x_product(double *x,int j){
double product = 1.0;
for(int i=0;i<N;i++){
if(i != j){
product *= x[i];
}
}
return product;
}
/*総和の計算*/
double x_sum(double *x){
double sum = 0.0;
for(int i=0;i<N;i++){
sum += x[i];
}
return sum;
}
/*xの組み合わせの積による係数計算関数群*/
double x_other(int j,int num,double *x){ //x_j除外するため、jを除外した部分集合の要素、N-num次目の計算
double sum = 0.0;
int *element = (int *)calloc(N, sizeof(int));
if(element == NULL){
printf("メモリ確保に失敗しました\n");
return -1.0;
}
generate_combination(N,j,element,num,0,0,x,&sum);
free(element);
return sum;
}
void generate_combination(int n,int j,int *element, int k, int start, int count,double *x,double *sum){
if(count == k) {//elementは上書きされるので、その前に積をとる
*sum += combination_product(k,x,element);
return;
}
for(int i = start;i<n;i++) {
if(i != j){
element[count] = i;
generate_combination(n,j,element,k,i+1,count+1,x,sum);
}
}
}
double combination_product(int k,double *x,int *element){
double product = 1.0;
for(int i=0;i<k;i++){
product *= x[element[i]];
}
return product;
}
double denominator(int num,double *x){
double product = 1.0;
for(int i=0;i<N;i++) {
if(i != num){
product *= (x[num] - x[i]);
}
}
return product;
}
/*最終結果の出力*/
void output_result(double (*P)[N]){
double result[N] = {0}; //最終的な係数を格納する配列
/*各基底多項式の係数の和を計算*/
for(int j=0;j<N;j++){
for(int i=0;i<N;i++){
result[j] += P[i][j];// ←先ほどの行列を次数毎で和をとる
}
}
/*f_lベクトルを行列として表示*/
printf("\n\n各次数に対応する係数をベクトル表示\n");
printf("※ -0 は非整数\n");
for(int i=0;i<N;i++){
printf("\t|");
for(int j=0;j<N;j++){
if(fabs(P[i][j]) == 0) printf("%5.0f ",0.0);
else printf("%5.0f ",P[i][j]);
}
printf("||x^%d|\n",(N-1)-i);
}
/*各次数の係数*/
printf("\n各列の和 ");
for(int i=0;i<N;i++){
printf("%5.0f ",result[i]);
}
//結果の出力
printf("\n\n求まった%d次関数\nf(x)= ",N-1);
for(int i=0;i<N;i++){
if(fabs(result[i]) < EPSILON) continue;
//符号
if(i != 0){
if(result[i] > 0) printf(" + ");
else printf(" - ");
}
else if(result[i] < 0) printf("-");
// 係数（係数が1または-1の場合は省略）
if(fabs(result[i]) != 1 || i == N-1){ //定数項の場合は表示
if(0 < fabs(result[i]) && fabs(result[i]) < 1) printf("%.1f",fabs(result[i]));
else printf("%.0f",fabs(result[i]));
}
//定数項以外
if(i < N-1){
if(i == N-2) printf("x");
else printf("x^%d",(N-1)-i);
}
}
printf("\n\n");
}xの組み合わせはコンビネーションを出力する関数を最初に作り、埋め込みました。再帰を使っているのがポイントです。「単一責任の原則」が好きなので、関数が多いです。これにも欠点があり、座標の個数が事前に設定しとかないといけないので、いちいちコードを編集しなければいけないということです。もし、回避するなら全ての配列が動的確保になるので、めんどかったんです…
締め
コーディング自体は全然得意じゃなくて、周りより頭が硬い方だと自負しているのですが、書くのは好きなんですよね。自分の速度で色々し続けたいと思います。
(追記)動的確保版を作成したので、defineから設定する必要がなくなりました！！(②に掲載)