連立方程式のプログラムがうまく動かない

前提・実現したいこと

LU分解を用いた連立方程式のプログラムの作成。

発生している問題・エラーメッセージ

xの値がおかしい

正しいプログラムではx[0]=1,x[1]=2,x[2]=3,x[3]=4となります。

該当のソースコード

c
1#include <stdio.h>
2
3int main(){
4 double a[N][N]={{4,-6,-10,-2},//係数行列
5                {-6,8,21,-1},
6                {-10,21,-20,23},
7                {-2,-1,23,-18}};
8 double b[N] ={-46,69,64,-7};//右辺ベクトル
9 double x[N],y[N];//中間ベクトル
10 double fct; //正規化比率
11 double sum; //積和
12 int i,j,k;//係数
13
14 // 連立方程式の行列
15 printf("a[%d][%d]b[%d]:\n",N,N,N);
16 for(i=0;i<N;i++){
17 for(j=0;j<N;j++){
18 printf("%9f,",a[i][j]);
19 }printf(" : %9f\n",b[i]);
20}
21//解法
22printf("\nStart LU-decomposition:\n");
23  for(k=0;k<N;++k){
24    if(a[k][k]==0){
25      printf("a[k][k]==0：confirm equations！\n");
26      return 1;
27    }
28    for(i=k+1;i<N;++i){
29      fct = a[i][k]/a[k][k];
30      for(j=k+1;j<N;++j){
31        a[i][j] -= fct*a[k][j];
32      }
33    }
34  }
35  printf("\ny[] vector:\n");
36  for(i=0;i<N;++i){
37    sum = 0.0;
38    for(j=0;j<i;++j){
39      sum += a[i][j]*y[j];
40    }
41    y[i] = b[i]-sum;
42    printf("  y[%d] = %9f\n",i,y[i]);
43  }
44  printf("\nx[] vector:\n");
45  for(i=N-1;i>=0;--i){
46    sum = 0.0;
47    for(j=i+1;j<N;++j){
48      sum += a[i][j]*x[j];
49    }
50    x[i] = (y[i]-sum)/a[i][i];
51     printf("  x[%d] = %9f\n",i,x[i]);
52  }
53  printf("\nSolution x[%d]:\n",N);
54  for(i=0;i<N;i++){
55    printf("x[%d] = %20.17f\n",i,x[i]);
56  }
57  return 0;
58}
59