/*
 ============================================================================
 Name        : napierlog.c
 Author      : Peng Sun, sunpeng14@baidu.com
 Version     : 1.0
 Copyright   : for study only
 Description : This is a program for calculating Napier's three tables and his small table (Times Table) using the library of mpfr, written in C.
               This program makes use of the calculate process in Denis's paper "Napier’s ideal construction of the logarithms", which can be download at
               https://locomat.loria.fr/napier/napier1619construction.pdf.
 ============================================================================
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <mpfr.h>

//Table 1
mpfr_t ratio1, a[101]; //ratio1 is 0.9999999, a[0]=10^7
mpfr_t nlog_a_lower_bound[101], nlog_a_upper_bound[101], nlog_a_avg[101]; //The phrase "nlog" stands for napier's log, or Nlog(x). Yet in Denis's paper it is λn(x) at page 5
//Table 2
mpfr_t ratio2, b[51]; //ratio2 is 0.99999, b[0]=10^7
mpfr_t nlog_b_lower_bound[51], nlog_b_upper_bound[51], nlog_b_avg[51];
//Table 3
mpfr_t ratio3_row, ratio3_col, c[21][69]; //ratio3_row=0.99, ratio3_col=0.9995, c[0][0]=10^7
mpfr_t nlog_c_lower_bound[21][69], nlog_c_upper_bound[21][69], nlog_c_avg[21][69];
//Small Table (Times Table)
int idx_d[28] = {2, 4, 8, 10,
	            20, 40, 80, 100,
				200, 400, 800, 1000,
				2000, 4000, 8000, 10000,
				20000, 40000, 80000, 100000,
				200000, 400000, 800000, 1000000,
				2000000, 4000000, 8000000, 10000000};
mpfr_t nlog_d_lower_bound[28], nlog_d_upper_bound[28], nlog_d_avg[28];
/*
 * Calculate nlog value using an adjacent known nlog value, please reference to 《Mirifici Logarithmorum Canonis Constructio》RULE 39,
 * or see Denis's paper "Napier’s ideal construction of the logarithms" at page 11.
 */
void rule39_AdjacentSolver(mpfr_t number1_known, mpfr_t number2_tobe_solved,
		mpfr_t nlog_number1_lower_bound, mpfr_t nlog_number1_upper_bound,
		mpfr_t nlog_number2_lower_bound, mpfr_t nlog_number2_upper_bound, mpfr_t nlog_number2_avg){
	mpfr_t diff_number1_2; //temporary variable storing the difference between number1_known and number2_tobe_solved
	mpfr_init2(diff_number1_2, 200);//set 200 bit
	mpfr_sub(diff_number1_2, number1_known, number2_tobe_solved, MPFR_RNDN);

	mpfr_div(nlog_number2_lower_bound, diff_number1_2, number1_known, MPFR_RNDN);
	mpfr_mul_ui(nlog_number2_lower_bound, nlog_number2_lower_bound, 10000000, MPFR_RNDN);
	mpfr_add(nlog_number2_lower_bound, nlog_number2_lower_bound, nlog_number1_lower_bound, MPFR_RNDN);

	mpfr_div(nlog_number2_upper_bound, diff_number1_2, number2_tobe_solved, MPFR_RNDN);
	mpfr_mul_ui(nlog_number2_upper_bound, nlog_number2_upper_bound, 10000000, MPFR_RNDN);
	mpfr_add(nlog_number2_upper_bound, nlog_number2_upper_bound, nlog_number1_upper_bound, MPFR_RNDN);

	mpfr_add(nlog_number2_avg, nlog_number2_lower_bound, nlog_number2_upper_bound, MPFR_RNDN);
	mpfr_div_ui(nlog_number2_avg, nlog_number2_avg, 2, MPFR_RNDN);

	mpfr_clear(diff_number1_2);
	return;
}
void generateTable1(){
	//initialize
	mpfr_init2(ratio1, 200);
	int i;
	for(i=0; i<=100; i++){
		mpfr_init2(a[i], 200);
		mpfr_init2(nlog_a_lower_bound[i], 200);
		mpfr_init2(nlog_a_upper_bound[i], 200);
		mpfr_init2(nlog_a_avg[i], 200);
	}

	/*
	 * If using mpfr_set_d (r, 0.9999999, MPFR_RNDD), it will produce: 9.9999990000000005263558477963670156896114349365234375000000000e-1,
	   yet if using mpfr_set_str (r,"0.9999999",10,GMP_RNDN), it will produce: 9.9999989999999999999999999999999999999999999999999999999999970e-1,
	   the latter will win a little bit more accuracy.
	   Thank Denis for discussing about this point.
	 */
	mpfr_set_str(ratio1, "0.9999999", 10, GMP_RNDN);
	mpfr_init_set_ui(a[0], 10000000, GMP_RNDN);

	//generate all numbers in a[], such as a[100]=10^7*(0.9999999^100)
	for(i=1; i<=100; i++){
		mpfr_pow_ui(a[i], ratio1, i, MPFR_RNDN);
		mpfr_mul(a[i], a[i], a[0], MPFR_RNDN);
	}
	//calculate the first nlog of napier's table 1, that is Nlog(a[1]),or Nlog(9999999.0000000)
	mpfr_init_set_ui(nlog_a_lower_bound[1], 1, GMP_RNDN);
	mpfr_div(nlog_a_upper_bound[1], a[0], ratio1, MPFR_RNDN);
	mpfr_sub(nlog_a_upper_bound[1], nlog_a_upper_bound[1], a[0], MPFR_RNDN);
	mpfr_add(nlog_a_avg[1], nlog_a_lower_bound[1], nlog_a_upper_bound[1], MPFR_RNDN);
	mpfr_div_ui(nlog_a_avg[1], nlog_a_avg[1], 2, MPFR_RNDN);

	//calculate the remaining nlogs of napier's talbe 1
	for(i=2; i<=100; i++){
		mpfr_mul_ui(nlog_a_lower_bound[i], nlog_a_lower_bound[1], i, MPFR_RNDN);
		mpfr_mul_ui(nlog_a_upper_bound[i], nlog_a_upper_bound[1], i, MPFR_RNDN);
		mpfr_add(nlog_a_avg[i], nlog_a_lower_bound[i], nlog_a_upper_bound[i], MPFR_RNDN);
		mpfr_div_ui(nlog_a_avg[i], nlog_a_avg[i], 2, MPFR_RNDN);
	}
	//output table 1
	printf("i,a[i],nlog_a_lower_bound,nlog_a_upper_bound,nlog_a_avg\n");
	for(i=1; i<=100; i++){
		mpfr_printf("%Pd,%.20RNf,%.20RNf,%.20RNf,%.20RNf\n", i, a[i], nlog_a_lower_bound[i], nlog_a_upper_bound[i], nlog_a_avg[i]);
	}
}
void generateTable2(){
	//initialize
	mpfr_init2(ratio2, 200);
	int i;
	for(i=0; i<=50; i++){
		mpfr_init2(b[i], 200);
		mpfr_init2(nlog_b_lower_bound[i], 200);
		mpfr_init2(nlog_b_upper_bound[i], 200);
		mpfr_init2(nlog_b_avg[i], 200);
	}
	mpfr_set_str(ratio2, "0.99999", 10, GMP_RNDN);
	mpfr_init_set_ui(b[0], 10000000, GMP_RNDN);

	//generate numbers in b[i], such as b[50]=10^7*(0.99999^50)
	for(i=1; i<=50; i++){
		mpfr_pow_ui(b[i], ratio2, i, MPFR_RNDN);
		mpfr_mul(b[i], b[i], b[0], MPFR_RNDN);
	}
	//calculate the first nlog of napier's table 2, that is Nlog(b[1]), or Nlog(9999900)
	rule39_AdjacentSolver(a[100], b[1],
			nlog_a_lower_bound[100], nlog_a_upper_bound[100],
			nlog_b_lower_bound[1], nlog_b_upper_bound[1], nlog_b_avg[1]);

	//calculate the remaining nlogs of napier's table 2
	for(i=2; i<=50; i++){
		mpfr_mul_ui(nlog_b_lower_bound[i], nlog_b_lower_bound[1], i, MPFR_RNDN);
		mpfr_mul_ui(nlog_b_upper_bound[i], nlog_b_upper_bound[1], i, MPFR_RNDN);
		mpfr_add(nlog_b_avg[i], nlog_b_lower_bound[i], nlog_b_upper_bound[i], MPFR_RNDN);
		mpfr_div_ui(nlog_b_avg[i], nlog_b_avg[i], 2, MPFR_RNDN);
	}
	//output table 2
	printf("i,b[i],nlog_b_lower_bound,nlog_b_upper_bound,nlog_b_avg\n");
	for(i=1; i<=50; i++){
		mpfr_printf("%Pd,%.20RNf,%.20RNf,%.20RNf,%.20RNf\n", i, b[i], nlog_b_lower_bound[i], nlog_b_upper_bound[i], nlog_b_avg[i]);
	}
}
void generateTable3(){
	//initialize
	mpfr_init2(ratio3_row, 200);
	mpfr_init2(ratio3_col, 200);
	int i,j;
	for(j=0; j<=68; j++){
		for(i=0; i<=20; i++){
			mpfr_init2(c[i][j], 200);
			mpfr_init2(nlog_c_lower_bound[i][j], 200);
			mpfr_init2(nlog_c_upper_bound[i][j], 200);
			mpfr_init2(nlog_c_avg[i][j], 200);
		}
	}
	mpfr_set_str(ratio3_row, "0.99", 10, GMP_RNDN);
	mpfr_set_str(ratio3_col, "0.9995", 10, GMP_RNDN);

	//1. generate numbers in c[i][j], such as c[0][j]=10^7*(0.99^j),c[i][j]=c[0][j]*(0.9995^i)
	for(j=0; j<=68; j++){
		mpfr_pow_ui(c[0][j], ratio3_row, j, MPFR_RNDN);
		mpfr_mul_ui(c[0][j], c[0][j], 10000000, MPFR_RNDN);
		for(i=1; i<=20; i++){
			mpfr_pow_ui(c[i][j], ratio3_col, i, MPFR_RNDN);
			mpfr_mul(c[i][j], c[i][j], c[0][j], MPFR_RNDN);
		}
	}
	//2. calculate the first nlog of first row of napier's table 3, that is Nlog(c[1][0]), or Nlog(9995000)
	//2.1 firstly, get x, using x/10000000=c[1][0]/b[50], see Denis's paper at page 12
	mpfr_t x, nlog_x_lower_bound, nlog_x_upper_bound, nlog_x_avg;
	mpfr_init2(x, 200);
	mpfr_init2(nlog_x_lower_bound, 200);
	mpfr_init2(nlog_x_upper_bound, 200);
	mpfr_init2(nlog_x_avg, 200);
	mpfr_div(x, c[1][0], b[50], MPFR_RNDN);
	mpfr_mul_ui(x, x, 10000000, MPFR_RNDN); //this gives x 999998.774583418771

	//2.2 then, calculate Nlog(x), using the nearest number in table 1, which is 9999999, a[1]
	rule39_AdjacentSolver(a[1], x, nlog_a_lower_bound[1], nlog_a_upper_bound[1],
			              nlog_x_lower_bound, nlog_x_upper_bound, nlog_x_avg);

	//2.3 last, calculate Nlog(c[1][0]), using Nlog(c[1][0]) = Nlog(x) + Nlog(b[50])
	mpfr_add(nlog_c_lower_bound[1][0], nlog_x_lower_bound, nlog_b_lower_bound[50], MPFR_RNDN);
	mpfr_add(nlog_c_upper_bound[1][0], nlog_x_upper_bound, nlog_b_upper_bound[50], MPFR_RNDN);
	mpfr_add(nlog_c_avg[1][0], nlog_c_lower_bound[1][0], nlog_c_upper_bound[1][0], MPFR_RNDN);
	mpfr_div_ui(nlog_c_avg[1][0], nlog_c_avg[1][0], 2, MPFR_RNDN);

	//2.4 calculate remaining nlogs of first column, c[i][0]
	mpfr_set_str(nlog_c_lower_bound[0][0], "0.0", 10, GMP_RNDN);
	mpfr_set_str(nlog_c_upper_bound[0][0], "0.0", 10, GMP_RNDN);
	mpfr_set_str(nlog_c_avg[0][0], "0.0", 10, GMP_RNDN);
	for(i=2; i<=20; i++){
		mpfr_mul_ui(nlog_c_lower_bound[i][0], nlog_c_lower_bound[1][0], i, MPFR_RNDN);
		mpfr_mul_ui(nlog_c_upper_bound[i][0], nlog_c_upper_bound[1][0], i, MPFR_RNDN);
		mpfr_add(nlog_c_avg[i][0], nlog_c_lower_bound[i][0], nlog_c_upper_bound[i][0], MPFR_RNDN);
		mpfr_div_ui(nlog_c_avg[i][0], nlog_c_avg[i][0], 2, MPFR_RNDN);
	}

	//3. calculate the second nlog of first row of napier's table 3, that is Nlog(c[0][1]), or Nlog(9900000.0000)
	//3.1 firstly, get x, using x/10000000=c[0][1]/c[20][0], see Denis's paper at page 14
	mpfr_div(x, c[0][1], c[20][0], MPFR_RNDN);
	mpfr_mul_ui(x, x, 10000000, MPFR_RNDN); //this gives x 9999521.66124220837111895696

	//3.2 then, calculate nlog(x), using the neareast number in table 2, which is 9999500.00999990000049999900, b[5]
	rule39_AdjacentSolver(b[5], x, nlog_b_lower_bound[5], nlog_b_upper_bound[5],
			nlog_c_lower_bound[0][1], nlog_c_upper_bound[0][1], nlog_c_avg[0][1]);

	//3.3 last, calculate Nlog(c[0][1]), using Nlog(c[0][1]) = Nlog(x) + Nlog(b[20][0])
	mpfr_add(nlog_c_lower_bound[0][1], nlog_c_lower_bound[0][1], nlog_c_lower_bound[20][0], MPFR_RNDN);
	mpfr_add(nlog_c_upper_bound[0][1], nlog_c_upper_bound[0][1], nlog_c_upper_bound[20][0], MPFR_RNDN);
	mpfr_add(nlog_c_avg[0][1], nlog_c_lower_bound[0][1], nlog_c_upper_bound[0][1], MPFR_RNDN);
	mpfr_div_ui(nlog_c_avg[0][1], nlog_c_avg[0][1], 2, MPFR_RNDN);

	//4. calculate the remaining nlogs of first row of napier's table 3, that is Nlog(c[0][j])
	for(j=2; j<=68; j++){
		mpfr_mul_ui(nlog_c_lower_bound[0][j], nlog_c_lower_bound[0][1], j, MPFR_RNDN);
		mpfr_mul_ui(nlog_c_upper_bound[0][j], nlog_c_upper_bound[0][1], j, MPFR_RNDN);

		mpfr_add(nlog_c_avg[0][j], nlog_c_lower_bound[0][j], nlog_c_upper_bound[0][j], MPFR_RNDN);
		mpfr_div_ui(nlog_c_avg[0][j], nlog_c_avg[0][j], 2, MPFR_RNDN);
	}
	//5. calculate all remaining nlogs of numbers in napier's table 3
	for(j=1; j<=68; j++){
		for(i=1; i<=20; i++){
			mpfr_mul_ui(nlog_c_lower_bound[i][j], nlog_c_lower_bound[1][0], i, MPFR_RNDN);
			mpfr_add(nlog_c_lower_bound[i][j], nlog_c_lower_bound[i][j], nlog_c_lower_bound[0][j], MPFR_RNDN);

			mpfr_mul_ui(nlog_c_upper_bound[i][j], nlog_c_upper_bound[1][0], i, MPFR_RNDN);
			mpfr_add(nlog_c_upper_bound[i][j], nlog_c_upper_bound[i][j], nlog_c_upper_bound[0][j], MPFR_RNDN);

			mpfr_add(nlog_c_avg[i][j], nlog_c_lower_bound[i][j], nlog_c_upper_bound[i][j], MPFR_RNDN);
			mpfr_div_ui(nlog_c_avg[i][j], nlog_c_avg[i][j], 2, MPFR_RNDN);
		}
	}
	//output table 3
	printf("i,c[i][j],nlog_c_lower_bound,nlog_c_upper_bound,nlog_c_avg\n");
	for(j=0; j<=68; j++){
		mpfr_printf("%Pd. column\n", j);
		for(i=0; i<=20; i++){
			mpfr_printf("%Pd,%Pd,%.20RNf,%.20RNf,%.20RNf,%.20RNf\n", i, j, c[i][j], nlog_c_lower_bound[i][j], nlog_c_upper_bound[i][j], nlog_c_avg[i][j]);
		}
	}
	mpfr_clear(x);
	mpfr_clear(nlog_x_lower_bound);
	mpfr_clear(nlog_x_upper_bound);
	mpfr_clear(nlog_x_avg);
}
void generateTimesTable(){
	int i;
	for(i=0; i<=27; i++){
		mpfr_init2(nlog_d_lower_bound[i], 200);
		mpfr_init2(nlog_d_upper_bound[i], 200);
		mpfr_init2(nlog_d_avg[i], 200);
	}
	//Nlog(5000000), initialization
	mpfr_t n5000000, nlog_n5000000_lower_bound, nlog_n5000000_upper_bound, nlog_n5000000_avg;
	mpfr_init2 (n5000000, 200);
	mpfr_init2 (nlog_n5000000_lower_bound, 200);
	mpfr_init2 (nlog_n5000000_upper_bound, 200);
	mpfr_init2 (nlog_n5000000_avg, 200);
	mpfr_set_str(n5000000, "5000000", 10, GMP_RNDN);
	//calculate Nlog(5000000), using the nearest number in table 3, which is c[19,68]
	rule39_AdjacentSolver(c[19][68], n5000000, nlog_c_lower_bound[19][68], nlog_c_upper_bound[19][68],
			nlog_n5000000_lower_bound, nlog_n5000000_upper_bound, nlog_n5000000_avg);

	//Nlog(8000000), initialization
	mpfr_t n8000000, nlog_n8000000_lower_bound, nlog_n8000000_upper_bound, nlog_n8000000_avg;
	mpfr_init2(n8000000, 200);
	mpfr_init2(nlog_n8000000_lower_bound, 200);
	mpfr_init2(nlog_n8000000_upper_bound, 200);
	mpfr_init2(nlog_n8000000_avg, 200);
	mpfr_set_str(n8000000, "8000000", 10, GMP_RNDN);
	//calculate Nlog(8000000), using the nearest number in table 3, which is c[4][22]
	rule39_AdjacentSolver(c[4][22], n8000000, nlog_c_lower_bound[4][22], nlog_c_upper_bound[4][22],
			nlog_n8000000_lower_bound, nlog_n8000000_upper_bound, nlog_n8000000_avg);

	//calculate ρ(2) using Nlog(5000000), the notation ρ means ρ(p/q) = Nlog(q) - Nlog(p), see Denis's paper at page 18
	mpfr_mul_si(nlog_d_lower_bound[0], nlog_n5000000_lower_bound, 1, MPFR_RNDN);
	mpfr_mul_si(nlog_d_upper_bound[0], nlog_n5000000_upper_bound, 1, MPFR_RNDN);
	mpfr_mul_si(nlog_d_avg[0], nlog_n5000000_avg, 1, MPFR_RNDN);
	//calculate ρ(4)
	mpfr_mul_si(nlog_d_lower_bound[1], nlog_d_lower_bound[0], 2, MPFR_RNDN);
	mpfr_mul_si(nlog_d_upper_bound[1], nlog_d_upper_bound[0], 2, MPFR_RNDN);
	mpfr_mul_si(nlog_d_avg[1], nlog_d_avg[0], 2, MPFR_RNDN);
	//calculate ρ(8)
	mpfr_mul_si(nlog_d_lower_bound[2], nlog_d_lower_bound[0], 3, MPFR_RNDN);
	mpfr_mul_si(nlog_d_upper_bound[2], nlog_d_upper_bound[0], 3, MPFR_RNDN);
	mpfr_mul_si(nlog_d_avg[2], nlog_d_avg[0], 3, MPFR_RNDN);

	//Nlog(1000000), initialization
	mpfr_t nlog_n1000000_lower_bound, nlog_n1000000_upper_bound, nlog_n1000000_avg;
	mpfr_init2(nlog_n1000000_lower_bound, 200);
	mpfr_init2(nlog_n1000000_upper_bound, 200);
	mpfr_init2(nlog_n1000000_avg, 200);
	//Nlog(1000000) = Nlog(8000000) + ρ(8), see Denis's paper at page 19
	mpfr_add(nlog_n1000000_lower_bound, nlog_n8000000_lower_bound, nlog_d_lower_bound[2], MPFR_RNDN);
	mpfr_add(nlog_n1000000_upper_bound, nlog_n8000000_upper_bound, nlog_d_upper_bound[2], MPFR_RNDN);
	mpfr_add(nlog_n1000000_avg, nlog_n8000000_avg, nlog_d_avg[2], MPFR_RNDN);

	//calculate ρ(10)
	mpfr_mul_si(nlog_d_lower_bound[3], nlog_n1000000_lower_bound, 1, MPFR_RNDN);
	mpfr_mul_si(nlog_d_upper_bound[3], nlog_n1000000_upper_bound, 1, MPFR_RNDN);
	mpfr_mul_si(nlog_d_avg[3], nlog_n1000000_avg, 1, MPFR_RNDN);

	/*
	 * Calculate 20,40,80,100; 200,400,800,1000; ...
	 * For convenience, the one dimension array nlog_d_lower_bound[28] is organized as 7 groups, and each group has 4 elements:
	 * "2,4,8,10","20,40,80,100", "200,400,800,1000", "2000,4000,8000,10000", "20000,40000,80000,100000",
	 * "200000,400000,800000,1000000", "200000,400000,800000,10000000", "2000000,4000000,8000000,100000000".
	 */
	int j;
	for(i = 0; i <= 5; i++){
		int group_idx = 4 + i*4;
		for(j = 0; j <= 3; j++){
			/*
			 * The formula is ρ(ab) = ρ(a) + ρ(b), see Denis's paper at page 18:
			 * ρ(20) = ρ(2) + ρ(10), ρ(200) = ρ(20) + ρ(10), ρ(2000) = ρ(200) + ρ(10), ...
			 * ρ(40) = ρ(4) + ρ(10), ρ(400) = ρ(40) + ρ(10), ρ(4000) = ρ(400) + ρ(10), ...
			 * ρ(80) = ρ(8) + ρ(10), ρ(800) = ρ(80) + ρ(10), ρ(8000) = ρ(800) + ρ(10), ...
			 * ρ(100) = ρ(10) + ρ(10), ρ(1000) = ρ(100) + ρ(10), ρ(10000) = ρ(1000) + ρ(10), ...
			 * ...
			 * that is ρ(group_idx) = ρ(group_idx-4) + ρ(10)
			 */
			mpfr_add(nlog_d_lower_bound[group_idx + j], nlog_d_lower_bound[group_idx + j - 4], nlog_d_lower_bound[3], MPFR_RNDN);
			mpfr_add(nlog_d_upper_bound[group_idx + j], nlog_d_upper_bound[group_idx + j - 4], nlog_d_upper_bound[3], MPFR_RNDN);
			mpfr_add(nlog_d_avg[group_idx + j], nlog_d_avg[group_idx + j - 4], nlog_d_avg[3], MPFR_RNDN);
		}
	}
	//output
	printf("times,nlog_avg\n");
	for(i=0; i<=27 ; i++){
		mpfr_abs(nlog_d_avg[i], nlog_d_avg[i], MPFR_RNDN);
		mpfr_printf("%Pd,%.20RNf\n", idx_d[i], nlog_d_avg[i]);
	}

	mpfr_clear(n5000000);
	mpfr_clear(nlog_n5000000_lower_bound);
	mpfr_clear(nlog_n5000000_upper_bound);
	mpfr_clear(nlog_n5000000_avg);

	mpfr_clear(n8000000);
	mpfr_clear(nlog_n8000000_lower_bound);
	mpfr_clear(nlog_n8000000_upper_bound);
	mpfr_clear(nlog_n8000000_avg);

	mpfr_clear(nlog_n1000000_lower_bound);
	mpfr_clear(nlog_n1000000_upper_bound);
	mpfr_clear(nlog_n1000000_avg);
}
int main(void){
	generateTable1();
	generateTable2();
	generateTable3();
	generateTimesTable();

	//clear mpfr variables
	mpfr_clear(ratio1);
	int i, j;
	for(i=0; i<=100; i++){
		mpfr_clear(a[i]);
		mpfr_clear(nlog_a_lower_bound[i]);
		mpfr_clear(nlog_a_upper_bound[i]);
		mpfr_clear(nlog_a_avg[i]);
	}
	mpfr_clear(ratio2);
	for(i=0; i<=50; i++){
		mpfr_clear(b[i]);
		mpfr_clear(nlog_b_lower_bound[i]);
		mpfr_clear(nlog_b_upper_bound[i]);
		mpfr_clear(nlog_b_avg[i]);
	}
	mpfr_clear(ratio3_row);
	mpfr_clear(ratio3_col);
	for(j=0; j<=68; j++){
		for(i=0; i<=20; i++){
			mpfr_clear(c[i][j]);
			mpfr_clear(nlog_c_lower_bound[i][j]);
			mpfr_clear(nlog_c_upper_bound[i][j]);
			mpfr_clear(nlog_c_avg[i][j]);
		}
	}
	for(i=0; i<=27; i++){
		mpfr_clear(nlog_d_lower_bound[i]);
		mpfr_clear(nlog_d_upper_bound[i]);
		mpfr_clear(nlog_d_avg[i]);
	}
	mpfr_free_cache ();
	return EXIT_SUCCESS;
}