Pade approxmation

pade.c

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

#define MAX_N 13


int factorial(int n){
	int i, tmp = 1;
	
	if(n != 0){
		for(i = 1; i <= n; i++) tmp *= i;
	}
	
	return (tmp);
}

int main(int argc, char *argv[])
{
	FILE *fp;
	
	char *e = NULL;
	char str[15], buf[512], csv_reco[100];
	int i, j, k, ch, row_pivot;
	int m, n, N;
	double a[MAX_N+1], p[MAX_N+1], q[MAX_N+1], b[MAX_N][MAX_N+1], x[MAX_N];
	double tmp, csv_val, b_pivot, eps = 1.0e-05;
	
	#ifdef EBUG
		int ii, jj;
	#endif
	
	fprintf(stderr, "[m, n]-type (i.e. P_m(x)/Q_n(x)) Pade approxmation. \nm = ");
	fgets(str, sizeof(str), stdin);
	fflush(stdin);
	sscanf(str, "%d", &m);
	
	fprintf(stderr, "n = ");
	fgets(str, sizeof(str), stdin);
	fflush(stdin);
	sscanf(str, "%d", &n);
	fprintf(stderr, "\n");
	
	ch = 0;
	N = m+n;
	
	/* -- Plug m+n coefficients in a[i] used by stdin or csv file -- */
	for(i = 0; i < MAX_N+1; i++) a[i] = 0.0;
	
	if(m < 1 || n < 0){
		fprintf(stderr, "Do not input a number such as m < 1 or n < 0!!\n");
		exit(1);
	}else if(N > MAX_N){
		fprintf(stderr, "The degree of your polynomial (i.e. m+n) is too large!!\n");
		exit(1);
	}
	
	fprintf(stderr, "1 : Differential Coefficients\n2 : Polynomial Coefficients of your function\n");
	do{
		fprintf(stderr, "Input 1 or 2 : ");
		fgets(str, sizeof(str), stdin);
		fflush(stdin);
		sscanf(str, "%d", &ch);
	}while(!(ch == 1 || ch == 2));
	
	fprintf(stderr, "\n");
	
	if(argc < 2){
		if(ch == 1){
			fprintf(stderr, "Input Different Coefficients of your function at x = 0.\n");
			for(i = 0; i <= N; i++){
				fprintf(stderr, "f^(%d)[x = 0] = ", i);
				fgets(str, sizeof(str), stdin);
				fflush(stdin);
				sscanf(str, "%lf", &a[i]);
				a[i] = (double)(a[i]/factorial(i));
			}
		}else if(ch == 2){
			fprintf(stderr, "Input Polynomial Coefficients of your function at x = 0.\n");
			for(i = 0; i <= N; i++){
				fprintf(stderr, "a[%d] = ", i);
				fgets(str, sizeof(str), stdin);
				fflush(stdin);
				sscanf(str, "%lf", &a[i]);
			}
		}else{
			fprintf(stderr, "Do not enter a number except for 1 or 2.\n");
			exit(1);
		}
	}else{
		if((fp = fopen(argv[1], "r" )) == NULL){
        	fprintf(stderr, "Do not open your file, %s.\n", argv[1]);
        	exit(1);
		}
		
		while(fgets(buf, sizeof(buf), fp)){
			if(strncmp(buf, "#", 1) == 0 || strcmp(buf, "\n") == 0) continue;
			
			j = 0; k = 0;
			csv_reco[0] = '0';
			csv_reco[1] = '\0';
			
			#ifdef EBUG
				fprintf(stderr, "%s", buf);
			#endif
			
			for(i = 0; i < strlen(buf); i++){
				if(buf[i] == ',' || buf[i] == '\n'){
					csv_val = strtod(csv_reco, &e);
					if(*e != '\0'){
						fprintf(stderr, "There is broken record in this csv file.\n");
						fprintf(stderr, "hoge = %lf\n", csv_val);
						fclose(fp); exit(1);
					}
					
					a[k] = csv_val;
					
					#ifdef EBUG
						fprintf(stderr, "\na[%d] = %lf\n", k, a[k]);
					#endif
					
					if(k == N) break;
					
					csv_reco[0] = '0';
					csv_reco[1] = '\0';
					k++; j = 0;
				}else{
					if(j+1 < sizeof(csv_reco)){
						csv_reco[j] = buf[i];
						csv_reco[j+1] = '\0';
						j++;
					}else if(j+1 == sizeof(csv_reco)){
						fprintf(stderr, "There is too large double record in this csv file.\n");
						fclose(fp); exit(1);
					}
				}
			}
		}
		
		fclose(fp);
		
		if(ch == 1) for(i = 0; i <= N; i++) a[i] = (double)(a[i]/factorial(i));
    }
	
	#ifdef EBUG
		for(ii = 0; ii <= N; ii++){
			fprintf(stderr, "Polynomial Coefficients of your function at x = 0.");
			fprintf(stderr, "\na[%d] = %lf", ii, a[ii]);
		}
	#endif
	
	/* -- Start to make an expand coefficient matrix B -- */
	for(i = 0; i < MAX_N; i++){
		x[i] = 0.0;
		for(j = 0; j < MAX_N+1; j++) {
			b[i][j] = 0.0;
			p[j] = 0.0;
			q[j] = 0.0;
		}
	}
	q[0] = 1.0; p[0] = a[0];
	
	for(i = 0; i < N; i++){
		b[i][N] = a[i+1];
		if(i < m) b[i][i] = 1.0;
	}
	for(j = 0; j < n; j++){
		for(i = 0; i < N; i++){
			b[i][m+j] = -a[i-j];
		}
	}
	
	#ifdef EBUG
		fprintf(stderr,"\nExpand coefficient matrix B:");
		for(ii = 0; ii < N; ii++){
			fprintf(stderr,"\n");
			for(jj = 0; jj < N+1; jj++) fprintf(stderr, "%lf ", b[ii][jj]);
		}
		fprintf(stderr,"\n");
	#endif
	
	/* -- Gaussian elimination with partial pivoting --*/
	
	/* -- Forward Elimination -- */
	for(i = m; i < N; i++){
		if(i == N-1) break;
		row_pivot = i;
		b_pivot = b[row_pivot][i];
		
		if(n > 0){
			for(k = 1; k < N-i; k++){
			if(fabs(b_pivot) < fabs(b[i+k][i])){
				b_pivot = b[i+k][i];
				row_pivot = i+k;
			}
			
			if(fabs(b_pivot) < eps){
				fprintf(stderr, "This matrix is a singular matrix.\n");
				exit(1);
			}
		}
		}
		
		#ifdef EBUG
			fprintf(stderr, "Row pivot : %d\n\n", row_pivot);
		#endif
		
		if(row_pivot != i){
			for(j = 0; j < N+1; j++){
				tmp = b[i][j];
				b[i][j] = b[row_pivot][j];
				b[row_pivot][j] = tmp;
			}
		}
		
		#ifdef EBUG
			for(ii = 0; ii < N; ii++){
				fprintf(stderr,"\n");
			for(jj = 0; jj < N+1; jj++) fprintf(stderr, "%lf ", b[ii][jj]);
			}
			fprintf(stderr,"\n");
		#endif
		
		for(k = i+1; k < N; k++){
			tmp = b[k][i]/b[i][i];
			for(j = i+1; j < N+1; j++){
				b[k][j] -= tmp*b[i][j];
				b[k][i] = 0.0;
			}
		}
	}
	
	#ifdef EBUG
		fprintf(stderr,"\nCompletly sweeped martix B:");
		for(ii = 0; ii < N; ii++){
			fprintf(stderr,"\n");
			for(jj = 0; jj < N+1; jj++) fprintf(stderr, "%lf ", b[ii][jj]);
		}
		fprintf(stderr,"\n");
	#endif
	
	if(fabs(b[N-1][N-1]) < eps){
		fprintf(stderr, "Rank is not equal to N, so this matrix is a singular matrix.\n");
		exit(1);
	}
	
	/* -- Backward Substitution -- */
	x[N-1] = b[N-1][N]/b[N-1][N-1];
	
	for(i = N-2; i >= 0; i--){
		tmp = 0.0;
		for(j = i+1; j < N; j++) tmp += b[i][j]*x[j];
		x[i] = (b[i][N]-tmp)/b[i][i];
	}
	
	for(i = 0; i < N; i++){
		if(i < m){
			p[i+1] = x[i];
		}else{
			q[i+1-m] = x[i];
		}
	}
	
	/* -- Output p[i] & q[i] -- */
	fprintf(stdout, "# a[i]: Coefficients of Polynominal f(x) (i = 0-%d)\n%lf", N, a[0]);
	for(i = 1; i <= N; i++) fprintf(stdout, ", %lf", a[i]);
	fprintf(stdout, "\n\n# p[i]: Coefficients of Polynominal P_%d(x) (i = 0-%d)\n%lf", m, m, p[0]);
	for(i = 1; i <= m; i++) fprintf(stdout, ", %lf", p[i]);
	fprintf(stdout, "\n\n# q[i]: Coefficients of Polynominal Q_%d(x) (i = 0-%d)\n%lf", n, n, q[0]);
	if(n > 0){
		for(i = 1; i <= n; i++) fprintf(stdout, ", %lf", q[i]);
	}
	
	fprintf(stdout, "\n\n# Function of [%d, %d]-type (i.e. P_%d(x)/Q_%d(x)) Pade approxmation\ny = ((%lf)", m, n, m, n, p[0]);
	for(i = 1; i <= m; i++) fprintf(stdout, "+(%lf)*x^(%d)", p[i], i);
	fprintf(stdout, ")/((%lf)", q[0]);
	if(n > 0){
		for(i = 1; i <= n; i++) fprintf(stdout, "+(%lf)*x^(%d)", q[i], i);
	}
	fprintf(stdout, ")\n");
	
	fprintf(stderr, "\n[%d, %d]-type Pade approxmation is ended without any trouble.\n", m, n);
		
	return 0;
}