二分法について、うまくいかない

//グラフェンナノリボンエネルギーバンド　spinwave susceptibility
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>

#define M  2	// 存在するサイト数   //M=24で5nmくらい
#define N  4	//サイト数×２　A+Bの個数
#define P  101	//リボン長　全電子数　P*
#define R  100

void zheev_( char*, char*, int*, double _Complex*,  int*, double*, double _Complex*, int*, double*, int* );
double f(double x, int j, int h, int k);
double g(double kai_ele[N][N]);

double u_up[P+1][N][N], u_dw[P+1][N][N];		//固有ベクトル
double u_up_1[P+1][N][N], u_dw_1[P+1][N][N];
double Energy_dw_1[P+1][N],Energy_up_1[P+1][N];  

int main(int argc, char** argv){ 
	FILE* fp;	
	FILE* fp_1;
    FILE* fp_2;
	FILE* fp_3;
	FILE* fp_4;

    double b_pi=M_PI;
	double U=1.0;		//相互作用エネルギーの平均(eV) u/t

    if ((fp = fopen("kai.txt", "w")) == NULL)
	{
		printf("Cannot open the file\n");
		exit(1);
	}
	fp_1 = fopen("u_dw.txt", "r");   
    if(fp_1 == NULL) {
        printf("ファイルを開くことが出来ませんでした．¥n");
    }
    for(int l=0; l<=P; l++){
		for(int i=0;i<N;i++){
			for(int j=0;j<N;j++){
				fscanf(fp_1, "%lf", &(u_dw[l][i][j]) );   	
			}
		} 
	}
	fp_2 = fopen("Energydw_1.txt", "r");   
    if(fp_2 == NULL) {
        printf("ファイルを開くことが出来ませんでした．¥n");
    }
    for(int l=0; l<=P; l++){
		for(int i=0;i<N;i++){
			fscanf(fp_2, "%lf", &(Energy_dw_1[l][i]) );   	
		} 
	}
    fp_3 = fopen("Energyup_1.txt", "r");   
    if(fp_3 == NULL) {
        printf("ファイルを開くことが出来ませんでした．¥n");
    }
    for(int l=0; l<=P; l++){
		for(int i=0;i<N;i++){
			fscanf(fp_3, "%lf", &(Energy_up_1[l][i]) );   	
		} 
	}
	fp_4 = fopen("u_up.txt", "r");   
    if(fp_4 == NULL) {
        printf("ファイルを開くことが出来ませんでした．¥n");
    }
    for(int l=0; l<=P; l++){
		for(int i=0;i<N;i++){
			for(int j=0;j<N;j++){
				fscanf(fp_4, "%lf", &(u_up[l][i][j]) );   	
			}
		} 
    }

	//二分法
	double x_start1 = 1.6;
    double x_start2 = 0.0;

    for(int i=0;i<(P+1)/2+3;i++){
		double err = 100.0;
    	double a = x_start1;
    	double b = x_start2;
		double c;
		// int i=8;
		while(err>0.00000001){
			c = (a+b)/2.0;
			double gc=0.0;
      		
			for(int h=0;h<N;h++){
				for(int k=0;k<N;k++){            
					double kai_ele[N][N];
					kai_ele[h][k]=f(c,i,h,k);
					gc=g(kai_ele)/(P+1);
				}				
			}
        
      		err = gc*gc;
      		if(gc > U){
        		b = c;
      		}
      		else{
        		a = c;
      		}
			// printf("result:x=%lf , q=%lf\n",c,2.0*b_pi/(P+1)*i);
			// fprintf(fp, "%2f %2f\n", c ,2.0*b_pi/(P+1)*i);
		}
        // printf("result:x=%lf , q=%lf\n",c,2.0*b_pi/(P+1)*i);
		fprintf(fp, "%lf %lf\n", c ,2.0*b_pi/(P+1)*i);
    }

    /*for(int i=0;i<R;i++){
        double e=x_start1/R*i;
        for(int j=0;j<(P+1)/2+3;j++){
			for(int h=0;h<N;h++){
				for(int k=0;k<N;k++){            
					// double fc=f(e,j);
					double kai_ele[N][N];
					kai_ele[h][k]=f(e,j,h,k);
					gc=g(kai_ele);
				}				
			}
            // fprintf(fp, "%2f %2f %2f\n",e,2*b_pi/P*j,gc/(P+1));
			if(1/U-0.005<=gc/(P+1) &&  gc/(P+1)<=1/U+0.005){
				fprintf(fp, "%2f %2f %2f\n",e,2*b_pi/P*j,gc/(P+1));
			}
        }
        // fprintf(fp, "%2f %2f\n",e,gc/(P+1));
    }*/


    fclose(fp);	
 	fclose(fp_1);
    fclose(fp_2);
 	fclose(fp_3);	
	fclose(fp_4);
    return 0;
}
double f(double x, int j, int h, int k){  	
	double kai=0;0;
	double kai_2[P+1][M]={0.0};
	double kai_1[P+1][M][M]={0.0},kai_1_dw[P+1][M][M]={0.0},kai_1_up[P+1][M][M]={0.0};
	double kai_3[P+1]={0.0};
	

	for(int v=0;v<=P;v++){	//wavenumber
		for (int t=0;t<M;t++){	//valence
			for(int s=M;s<N;s++){	//conductive

				double Omega_up[P+1][M][M], Omega_dw[P+1][M][M];
				double U_up_plus[P+1][M][M], U_up_minus[P+1][M][M];
				double U_dw_plus[P+1][M][M], U_dw_minus[P+1][M][M];	
				
				Omega_up[v][t][s]=x-Energy_dw_1[v][s]+Energy_up_1[(v+j)%(P+1)][t];
				Omega_dw[v][t][s]=x+Energy_up_1[v][s]-Energy_dw_1[(v+j)%(P+1)][t];

				U_up_plus[v][t][s]=conj(u_up[v][t][k])*u_dw[(v+j)%(P+1)][s][k];
				U_up_minus[v][t][s]=conj(u_dw[(v+j)%(P+1)][s][h])*u_up[v][t][h];
				U_dw_plus[v][t][s]=conj(u_dw[v][t][k])*u_up[(v+j)%(P+1)][s][k];
				U_dw_minus[v][t][s]=conj(u_up[(v+j)%(P+1)][s][h])*u_dw[v][t][h];

				kai_1_up[v][t][s]=U_up_plus[v][t][s]*U_up_minus[v][t][s]/Omega_up[v][t][s];
				kai_1_dw[v][t][s]=U_dw_plus[v][t][s]*U_dw_minus[v][t][s]/Omega_dw[v][t][s];
				kai_1[v][t][s]=-kai_1_up[v][t][s]+kai_1_dw[v][t][s];

				kai_2[v][t]+=kai_1[v][t][s];
				// printf("%2f\n",kai_1[v][t][s]);
			}
			kai_3[v]+=kai_2[v][t];
			// printf("%2f \n",kai_2[v][t]);
		}	
		kai+=kai_3[v];
		// printf("%2f \n",kai_3[v]);
	}										
return kai;
}

double g(double kai_ele[N][N]){		
	double B[N];

    //感受率の対角化
                    // 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
                    char   jobz = 'V'           ;// 固有ベクトルを計算する
                    char   uplo = 'U'           ;// Aに上三角行列を格納
                    int    n    = N             ;// 対角化する正方行列のサイズ
                    double _Complex A[N*N]		;// 対角化する行列。対角化後は固有ベクトルが並ぶ

					for(int h=0;h<N;h++){
						for(int k=0;k<N;k++){
		                    A[N*h+k]=kai_ele[h][k];
						}
					}
					// fprintf(fp_1, "%2f %2f \n",A[N*k+h],kai[i][j][k][h]);
					// printf("%2f %2f\n",A[N*h+k],kai[h][k]);

                    double w[N]              	;// 実固有値が小さい順に入る
            
                    int    lda  = N	            ;// 対角化する正方行列のサイズ
                    double _Complex work[6*N]	;// 対角化する際に使用するメモリ
                    int    lwork = 6*N		    ;// workの次元
                    double rwork[3*N-2]      	;// 3*SIZE-2で固定
                    int    info                 ;// 成功すれば0、失敗すれば0以外を返す
  		            
                    zheev_( &jobz, &uplo, &n, A, &lda, w, work, &lwork, rwork, &info );
                        
                    for(int m=0;m<N;m++){
		            	B[m]=w[m];
					// printf("%2f\n",B[m]);
		            }

                double kai_acc=B[N-1];

  return kai_acc;
}