prototyping of reciprocal and log functions

src/decn/proto/.gitignore

@@ -0,0 +1,3 @@
+div
+ln
+

src/decn/proto/Makefile

@@ -0,0 +1,8 @@
+all: div ln
+
+div: div_mfp.cpp
+	g++ -Wall div_mfp.cpp -o div -lgmp
+
+ln: ln_mfp.cpp
+	g++ -Wall ln_mfp.cpp -o ln -lmpfr
+

src/decn/proto/div_mfp.cpp

@@ -0,0 +1,98 @@
+//prototype for reciprocal function for division based on newton-raphson iteration
+
+#include <stdio.h>
+#include <math.h>
+#include <gmp.h>
+
+static const int BSIZE = 200;
+char buf[BSIZE];
+char pbuf[BSIZE];
+static const int PRECISION = 1000;
+
+
+int main(void){
+	static const int MULT = 1000;
+	for (int x = 1 * MULT; x < 10 * MULT; x++){
+		//build number
+		sprintf(pbuf, "%d", x);
+		//copy with decimal point
+		if (pbuf[0]){
+			buf[0] = pbuf[0];
+		} else {
+			buf[0] = '0';
+		}
+		buf[1] = '.';
+		int b_i = 2;
+		for (int pb_i = 1; pbuf[pb_i] != '\0' && b_i < BSIZE; b_i++, pb_i++){
+			buf[b_i] = pbuf[pb_i];
+		}
+		//null terminate
+		if (b_i < BSIZE){
+			buf[b_i] = '\0';
+		} else {
+			buf[BSIZE-1] = '\0';
+		}
+		//build mpf
+		mpf_t x_mpf, accum, recip, actual;
+		static const int CALC_PRECISION = 18;
+		mpf_init2(x_mpf, CALC_PRECISION); //copy of original x
+		mpf_init2(accum, CALC_PRECISION); //working register
+		mpf_init2(recip, CALC_PRECISION); //current estimate of reciprocal
+		mpf_init2(actual, PRECISION); //truth value for 1/x
+		mpf_set_str(x_mpf, buf, 10);
+		mpf_set(actual, x_mpf);
+		mpf_ui_div(actual, 1, actual); //actual = 1/x
+		//get estimate
+		if (x <= 2.0 * MULT) {
+			sprintf(buf, "%de%d", 50, -2);
+		} else if (x < 3.3 * MULT) {
+			sprintf(buf, "%de%d", 30, -2);
+		} else if (x <= 5.0 * MULT) {
+			sprintf(buf, "%de%d", 20, -2);
+		} else {
+			sprintf(buf, "%de%d", 1, -1);
+		}
+		mpf_set_str(recip, buf, 10); //recip = estimate
+		mpf_set(accum, recip);       //accum = estimate
+
+		//do division
+		for (int i = 0; i < 6; i++){ //fixed number of iterations
+			mpf_mul(accum, accum, x_mpf); //accum *= x;
+			mpf_neg(accum, accum); //accum *= -1;
+			mpf_add_ui(accum, accum, 1); //accum += 1;
+			mpf_mul(accum, accum, recip); //accum *= recip;
+			mpf_add(accum, accum, recip);//accum += recip;
+			//recip = recip + (1 - recip*x)*recip;
+			mpf_set(recip, accum); //recip = accum;
+		}
+		//calculate error
+		mpf_t calc, diff;
+		mpf_init2(calc, PRECISION);
+		mpf_init2(diff, PRECISION);
+		int cmp;
+		mpf_set(calc, accum);        //calc = accum
+		mpf_reldiff(diff, actual, calc); //diff = (actual - calc)/actual
+		cmp = mpf_cmp_d(diff, 4.5e-20);
+		mp_exp_t exponent;
+		mpf_get_str(pbuf, &exponent, 10, 20, accum);
+		mpf_get_str(buf, &exponent, 10, 20, diff);
+		//add decimal point
+		exponent--;
+		char tmp = buf[1];
+		buf[1] = '.';
+		for (int i = 2; i <= 5; i++){
+			char tmp2 = buf[i];
+			buf[i] = tmp;
+			tmp = tmp2;
+		}
+		buf[6] = '\0';
+		//print if out of tolerance
+		if (cmp > 0){
+			printf("%5d: %20s, %s(%ld)\n",
+			   	    x, pbuf, buf, exponent);
+		}
+	}
+
+	return 0;
+}
+

src/decn/proto/ln_mfp.cpp

@@ -0,0 +1,91 @@
+//prototype for natural log function based on HP Journal article
+// "Personal Calculator Algorithms IV: Logarithmic Functions"
+
+#include <stdio.h>
+#include <iostream>
+#include <math.h>
+#include <boost/multiprecision/mpfr.hpp>
+using namespace boost::multiprecision;
+
+using std::cout;
+using std::endl;
+
+static const int BSIZE = 200;
+char buf[BSIZE];
+char pbuf[BSIZE];
+
+static const int NUM_A_ARR = 10;
+mpfr_float a_arr[NUM_A_ARR];
+mpfr_float ln_a_arr[NUM_A_ARR];
+
+static const unsigned int CALC_PRECISION = 18;
+
+
+int main(void){
+	cout << std::scientific << std::setprecision(CALC_PRECISION);
+	//initiallize a_arr and ln_a_arr
+	for (int i = 0; i < NUM_A_ARR; i++){
+		// a[i] = (1 + 10^-i)
+		mpfr_float a(-i, CALC_PRECISION);
+		a = pow(10, a);
+		a_arr[i] = a + 1;
+		ln_a_arr[i] = log(a_arr[i]);
+		cout << i << endl;
+		cout << "     a:" << a_arr[i] << endl;
+		cout << " ln(a):" << ln_a_arr[i] << endl;
+// 		a_arr[i].assign(a_arr1000[i], CALC_PRECISION);
+// 		ln_a_arr[i].assign(ln_a_arr1000[i], CALC_PRECISION);
+	}
+	cout << endl << endl ;
+	//loop through values to test
+	for (mpfr_float x(0.0001, CALC_PRECISION); x < 1.0; x += 0.015625){
+		//build mpf values
+		mpfr_float A(x, CALC_PRECISION);
+		mpfr_float b_j(1.0 - A, CALC_PRECISION); //initial b_0 = 1 - A
+		mpfr_float accum(b_j, CALC_PRECISION);
+		mpfr_float_1000 actual;
+		actual = log(A);
+		
+		int num_times[NUM_A_ARR];
+		for (int j = 0; j < NUM_A_ARR; j++){
+			if (j != 0){
+				b_j *= 10.0;
+			}
+			accum = b_j;
+			int k_j = 0;
+			while (accum > 0){
+				b_j = accum;
+				accum *= a_arr[j];
+				accum -= 1.0;
+				k_j++;
+			}
+			num_times[j] = k_j - 1;
+// 			printf("  %d: num_times: %d, ", j, num_times[j]);
+// 			cout << b_j << endl;
+		}
+		
+		//build final value
+		accum = b_j; //remainder
+		accum *= pow(10.0, -(NUM_A_ARR - 1)); //scale remainder
+// 		cout << "remainder: " << accum << endl;
+		for (int j = NUM_A_ARR - 1; j >= 0; j--){ //sum in reverse order
+			for (int k = 0; k < num_times[j]; k++){
+				accum += ln_a_arr[j];
+			}
+		}
+		accum = -accum;
+
+		mpfr_float_1000 calc, diff;
+		calc = accum;
+		diff = (actual - calc)/actual;
+		
+		if (diff > 5e-17){
+			cout << x << ": ";
+			cout << std::setprecision(18) << accum << ", ";
+			cout << diff << endl;
+		}
+	}
+
+	return 0;
+}
+