initial exp() implementation

(refout for log(1.5) changed due to using more accurate expected result)

src/decn/decn.c

@@ -19,6 +19,8 @@
 // #define DEBUG_DIV
 #define DEBUG_LOG
 // #define DEBUG_LOG_ALL //even more verbose
+#define DEBUG_EXP
+// #define DEBUG_EXP_ALL //even more verbose
 
 #ifndef DESKTOP
 //#undef EXTRA_CHECKS
@@ -48,8 +50,8 @@ static const uint8_t num_digits_display = 16;
 
 __idata dec80 AccDecn, BDecn;
 __idata dec80 TmpDecn; //used by add_decn() and mult_decn()
-__idata dec80 Tmp2Decn; //used by div_decn() and ln_decn()
+__idata dec80 Tmp2Decn; //used by recip_decn() and ln_decn()
-__idata dec80 Tmp3Decn; //used by div_decn() and ln_decn()
+__idata dec80 Tmp3Decn; //used by recip_decn() and ln_decn()
 __idata dec80 Tmp4Decn; //used by div_decn() and ln_decn()
 
 __xdata char Buf[DECN_BUF_SIZE];
@@ -553,6 +555,8 @@ void add_decn(void){
 		copy_decn(&AccDecn, &BDecn);
 		return;
 	}
+	//save b for restoring later
+	copy_decn(&TmpDecn, &BDecn);
 	//handle cases where signs differ
 	if (AccDecn.exponent < 0 && BDecn.exponent >= 0){
 		// -acc, +x
@@ -562,23 +566,23 @@ void add_decn(void){
 			printf("|-acc| > |+x|\n");
 #endif
 			sub_mag();
-			return;
 		} else if (rel == -1){
 #ifdef DEBUG_ADD
 			printf("|-acc| < |+x|\n");
 #endif
-			copy_decn(&TmpDecn, &BDecn);
 			copy_decn(&BDecn, &AccDecn);
 			copy_decn(&AccDecn, &TmpDecn);
 			sub_mag();
-			return;
 		} else { //equal
 #ifdef DEBUG_ADD
 			printf("|-acc| == |+x|\n");
 #endif
 			set_dec80_zero(&AccDecn);
-			return;
 		}
+		//restore b
+		copy_decn(&BDecn, &TmpDecn);
+
+		return;
 	} else if (AccDecn.exponent >= 0 && BDecn.exponent < 0){
 		// +acc, -x
 		rel = compare_magn();
@@ -587,23 +591,23 @@ void add_decn(void){
 			printf("|+acc| > |-x|\n");
 #endif
 			sub_mag();
-			return;
 		} else if (rel == -1){
 #ifdef DEBUG_ADD
 			printf("|+acc| < |-x|\n");
 #endif
-			copy_decn(&TmpDecn, &BDecn);
 			copy_decn(&BDecn, &AccDecn);
 			copy_decn(&AccDecn, &TmpDecn);
 			sub_mag();
-			return;
 		} else { //equal
 #ifdef DEBUG_ADD
 			printf("|+acc| == |-x|\n");
 #endif
 			set_dec80_zero(&AccDecn);
-			return;
 		}
+		//restore b
+		copy_decn(&BDecn, &TmpDecn);
+
+		return;
 	}
 	//signs must now be the same, begin adding
 	//normalize
@@ -658,6 +662,9 @@ void add_decn(void){
 		//track sign
 		set_exponent(&AccDecn, curr_exp, rel); //rel==is_neg?
 	}
+
+	//restore b
+	copy_decn(&BDecn, &TmpDecn);
 }
 
 //AccDecn *= BDecn
@@ -752,10 +759,9 @@ void mult_decn(void){
 	remove_leading_zeros(&AccDecn);
 }
 
-void div_decn(void){
+void recip_decn(void){
 #define CURR_RECIP Tmp2Decn //copy of x, holds current 1/x estimate
 #define X_COPY     Tmp3Decn //holds copy of original x
-#define ACC_COPY   Tmp4Decn //holds copy of original acc
 	uint8_t i;
 	exp_t initial_exp;
 	//check divide by zero
@@ -771,8 +777,6 @@ void div_decn(void){
 	//normalize
 	remove_leading_zeros(&AccDecn);
 	remove_leading_zeros(&BDecn);
-	//store copy of acc for final multiply by 1/x
-	copy_decn(&ACC_COPY, &AccDecn);
 	//store copy of x
 	copy_decn(&X_COPY, &BDecn);
 	//get initial exponent of estimate for 1/x
@@ -834,30 +838,41 @@ void div_decn(void){
 		//new_est(Accum) = recip + (1 - recip*x)*recip, where recip is current recip estimate
 		copy_decn(&CURR_RECIP, &AccDecn);
 	}
+
+//try not to pollute namespace
+#undef CURR_RECIP
+#undef X_COPY
+}
+
+inline void div_decn(void){
+#define ACC_COPY   Tmp4Decn //holds copy of original acc
+	//store copy of acc for final multiply by 1/x
+	copy_decn(&ACC_COPY, &AccDecn);
+	recip_decn();
 	//Accum now holds 1/x, multiply by original acc to complete division
 	copy_decn(&BDecn, &ACC_COPY);
 	mult_decn();
 
 //try not to pollute namespace
-#undef CURR_RECIP
-#undef X_COPY
 #undef ACC_COPY
 }
 
-void ln_decn(void){
+
+//constants used for ln(x) and exp(x)
 #define NUM_A_ARR 9
-	static const dec80 LN_A_ARR[NUM_A_ARR] = {
+static const dec80 LN_A_ARR[NUM_A_ARR] = {
-		{0, {69, 31, 47, 18,  5, 59, 94, 53,  9}},
+	{-1 & 0x7fff, {69, 31, 47, 18,  5, 59, 94, 53,  9}},
-		{0, {95, 31,  1, 79, 80, 43, 24, 86,  0}},
+	{-2 & 0x7fff, {95, 31,  1, 79, 80, 43, 24, 86,  0}},
-		{0, {99, 50, 33,  8, 53, 16, 80, 82, 84}},
+	{-3 & 0x7fff, {99, 50, 33,  8, 53, 16, 80, 82, 84}},
-		{0, {99, 95,  0, 33, 30, 83, 53, 31, 67}},
+	{-4 & 0x7fff, {99, 95,  0, 33, 30, 83, 53, 31, 67}},
-		{0, {99, 99, 50,  0, 33, 33,  8, 33, 33}},
+	{-5 & 0x7fff, {99, 99, 50,  0, 33, 33,  8, 33, 33}},
-		{0, {99, 99, 95,  0,  0, 33, 33, 29, 95}},
+	{-6 & 0x7fff, {99, 99, 95,  0,  0, 33, 33, 29, 95}},
-		{0, {99, 99, 99, 50,  0,  0, 33,  5, 35}},
+	{-7 & 0x7fff, {99, 99, 99, 50,  0,  0, 33,  5, 35}},
-		{0, {99, 99, 99, 95,  0,  0,  2, 76, 40}},
+	{-8 & 0x7fff, {99, 99, 99, 95,  0,  0,  2, 76, 40}},
-		{0, {99, 99, 99, 99, 50,  0, 15, 98, 65}},
+	{-9 & 0x7fff, {99, 99, 99, 99, 50,  0, 15, 98, 65}},
-// 		{0, {99, 99, 99, 99, 95,  2, 19, 27, 11}},
+};
-	};
+
+void ln_decn(void){
 	uint8_t j, k;
 	exp_t initial_exp;
 #define B_j Tmp2Decn
@@ -946,6 +961,7 @@ void ln_decn(void){
 		for (k = 0; k < NUM_TIMES.lsu[j]; k++){
 			//accum += ln_a_arr[j];
 			copy_decn(&BDecn, &LN_A_ARR[j]);
+			BDecn.exponent = 0; //tracking exponent through shifts/saved initial exponent
 			add_decn();
 		}
 		shift_right(&AccDecn);
@@ -1007,7 +1023,6 @@ void ln_decn(void){
 //try not to pollute namespace
 #undef B_j
 #undef NUM_TIMES
-#undef NUM_A_ARR
 }
 
 //inline void log10_decn(void){
@@ -1016,6 +1031,158 @@ void ln_decn(void){
 //	div_decn();
 //}
 
+
+
+void exp_decn(void){
+	uint8_t j, k;
+	uint8_t need_recip;
+	#define SAVED Tmp2Decn
+	#define NUM_TIMES Tmp3Decn
+
+	//check if negative
+	if (AccDecn.exponent < 0){
+		negate_decn(&AccDecn);
+		need_recip = 1;
+	}
+
+	//check if in range
+	copy_decn(&SAVED, &AccDecn); //save = accum
+	set_dec80_zero(&BDecn);
+	BDecn.lsu[0] = 23;
+	BDecn.lsu[1] = 02;
+	BDecn.lsu[2] = 58;
+	BDecn.lsu[2] = 51;
+	BDecn.exponent = 2; //b = 230.25851
+	negate_decn(&BDecn);
+	add_decn(); //accum = x - 230.25851 (should be negative if in range)
+	if (!(AccDecn.exponent < 0)){ //if not negative
+		set_dec80_NaN(&AccDecn);
+		return;
+	}
+	copy_decn(&AccDecn, &SAVED); //restore
+
+// 	//normalize
+// 	remove_leading_zeros(&AccDecn);
+
+	//initial b = -10*ln(10)
+	copy_decn(&BDecn, &DECN_LN_10); //b=ln(10)
+	copy_decn(&SAVED, &AccDecn); //save = accum
+	copy_decn(&AccDecn, &DECN_1);
+	AccDecn.exponent = 1; //accum = 10
+	mult_decn();             //accum =  10*ln(10)
+	copy_decn(&BDecn, &AccDecn); //b =  10*ln(10)
+	negate_decn(&BDecn);         //b = -10*ln(10)
+	copy_decn(&AccDecn, &SAVED); //restore
+	//track number of times 10*ln(10) can be subtracted
+	k = 0;
+	while (!(AccDecn.exponent < 0)){ //while not negative
+		copy_decn(&SAVED, &AccDecn); //save = accum
+		//accum -= 10*ln10
+		add_decn();
+		k++;
+	}
+	//subtracted 1 time too many
+	NUM_TIMES.exponent = (k - 1) * 10; //use exp to store number of times ln(10) subtracted
+	copy_decn(&AccDecn, &SAVED); //restore
+#ifdef DEBUG_EXP
+	decn_to_str_complete(&AccDecn);
+	printf("exp() num_times for 10*ln(10): %s (%d)\n", Buf, NUM_TIMES.exponent);
+#endif
+
+	//load b with -ln(10)
+	copy_decn(&BDecn, &DECN_LN_10); //b =  ln(10)
+	negate_decn(&BDecn);            //b = -ln(10)
+	//track number of times ln(10) and then (1 + 10^-j) can be subtracted
+	j = UINT8_MAX; //becomes 0 after incrementing to start (1 + 10^-j) series
+	do {
+		k = 0;
+		while (!(AccDecn.exponent < 0)){ //while not negative
+			copy_decn(&SAVED, &AccDecn); //save = accum
+			//accum -= b (ln10, then ln(1 + 10^-j) series)
+			add_decn();
+#ifdef DEBUG_EXP_ALL
+			decn_to_str_complete(&AccDecn);
+			printf("    %u: %s\n", k, Buf);
+// 			decn_to_str_complete(&BDecn);
+// 			printf("(%s)\n", Buf);
+#endif
+			k++;
+#ifdef DEBUG_EXP
+			assert(k != 255);
+#endif
+		}
+		//subtracted 1 time too many:
+		if (j == UINT8_MAX){
+			NUM_TIMES.exponent += k - 1;
+		} else {
+			NUM_TIMES.lsu[j] = k - 1;
+		}
+		copy_decn(&AccDecn, &SAVED); //restore
+#ifdef DEBUG_EXP
+		decn_to_str_complete(&AccDecn);
+		printf("exp() num_times for %d: %s (%d)\n", j, Buf, k-1);
+#endif
+		//next j
+		j++;
+		if (j < NUM_A_ARR){
+			//get next ln(1 + 10^-j) for subtraction
+			copy_decn(&BDecn, &LN_A_ARR[j]);
+			negate_decn(&BDecn);
+		} else {
+			break;
+		}
+	} while (1);
+
+	//build final value
+	// (currently accum = save = remainder)
+	// calculate 1+remainder
+	copy_decn(&BDecn, &DECN_1);
+	add_decn();
+	//get initial multiplier (10) for ln(10)
+	copy_decn(&BDecn, &DECN_1);
+	BDecn.exponent = 1; //BDecn = 10
+	//do multiplies
+	j = UINT8_MAX; //becomes 0 after incrementing to start (1 + 10^-j) series
+	do {
+		for (k = 0; k < (j==UINT8_MAX ? NUM_TIMES.exponent : NUM_TIMES.lsu[j]); k++){
+			mult_decn();
+		}
+#ifdef DEBUG_EXP
+		decn_to_str_complete(&AccDecn);
+		printf("exp() current val for %d: %s\n", j, Buf);
+#endif
+		//next j
+		j++;
+		if (j < NUM_A_ARR){
+			//get next multiplier (1 + 10^-j) for ln(1 + 10^-j)
+			if (j == 0){
+				//set to 2
+				BDecn.lsu[0] = 20;
+				BDecn.exponent = 0;
+			} else if (j == 1) {
+				//set to 1.1
+				BDecn.lsu[0] = 11;
+				//exponent is already 0
+			} else {
+				//get next (1 + 10^-j)
+				shift_right(&BDecn);
+				BDecn.lsu[0] = 10;
+			}
+		} else {
+			break;
+		}
+	} while (1);
+
+	//take reciprocal if exp was negative
+	if (need_recip){
+		recip_decn();
+	}
+
+//try not to pollute namespace
+#undef SAVED
+#undef NUM_TIMES
+}
+
 static void set_str_error(void){
 	Buf[0] = 'E';
 	Buf[1] = 'r';

src/decn/decn.h

@@ -64,18 +64,26 @@ void set_dec80_NaN(dec80* dest);
 uint8_t decn_is_nan(const dec80* x);
 
 void negate_decn(dec80* x);
-void add_decn(void);
+void add_decn(void);   //calculate AccDecn -= BDecn (BDecn is preserved)
-void mult_decn(void);
+void mult_decn(void);  //calculate AccDecn *= BDecn (BDecn is preserved)
+void recip_decn(void);
 void div_decn(void);
 
 void ln_decn(void);
 //inline void log10_decn(void);
-#define log10_decn() do {               \
+#define log10_decn() do {           \
-		ln_decn();                      \
+	ln_decn();                      \
-		copy_decn(&BDecn, &DECN_LN_10); \
+	copy_decn(&BDecn, &DECN_LN_10); \
-		div_decn();                     \
+	div_decn();                     \
-	} while(0);
+} while(0);
 
+void exp_decn(void);
+//exp10_decn() = exp_decn(AccDecn * ln(10))
+#define exp10_decn() do {           \
+	copy_decn(&BDecn, &DECN_LN_10); \
+	mult_decn();                    \
+	exp_decn();                     \
+} while(0);
 
 //Buf should hold at least 18 + 4 + 5 + 1 = 28
 #define DECN_BUF_SIZE 28

src/decn/decn_test.c

@@ -68,16 +68,32 @@ static void log10_test(
 {
 	build_dec80(x_str, x_exp);
 	decn_to_str_complete(&AccDecn);
-	printf("  a  : %s\n", Buf);
+	printf("   a  : %s\n", Buf);
 	log10_decn();
 	decn_to_str_complete(&AccDecn);
-	printf("ln(a): %s\n", Buf);
+	printf("log(a): %s\n", Buf);
 	build_decn_at(&diff, res_str, res_exp);
 	take_diff();
 	decn_to_str_complete(&diff);
 	printf("     : %s\n\n", Buf);
 }
 
+static void exp_test(
+	const char* x_str, int x_exp,
+	const char* res_str, int res_exp)
+{
+	build_dec80(x_str, x_exp);
+	decn_to_str_complete(&AccDecn);
+	printf("  a  : %s\n", Buf);
+	exp_decn();
+	decn_to_str_complete(&AccDecn);
+	printf("exp(a): %s\n", Buf);
+	build_decn_at(&diff, res_str, res_exp);
+	take_diff();
+	decn_to_str_complete(&diff);
+	printf("      : %s\n\n", Buf);
+}
+
 int main(void){
 //	dec80 acc, b;
 
@@ -245,7 +261,29 @@ int main(void){
 	//new acc for log test
 	log10_test(
 		"1.5", 0,
-		"0.176091259", 0
+		"0.176091259055681242", 0
 	);
+
+	//new acc for exp test
+	exp_test(
+		"4.4", 0,
+		"81.4508686649681174", 0
+	);
+
+	//get new acc for exp test (calculate 3^201)
+// 	build_dec80("3", 0);
+// 	ln_decn();
+	log_test(
+		"3", 0,
+		"1.09861228866810969", 0
+	);
+	build_decn_at(&BDecn, "201", 0);
+	mult_decn(); //acc = 201*ln(3)
+	int8_t exp_test_exp = decn_to_str(&AccDecn);
+	exp_test(
+		Buf, exp_test_exp,
+		"7.96841966627624308", 95
+	);
+
 	return 0;
 }

src/decn/decn_test_refout.txt

@@ -172,7 +172,7 @@ ln() exponent from initial: -9.
 ln(a): -22.3227534185273433
      : 1.E-16
 
-  a  : 1.5
+   a  : 1.5
 ln() accum scaled between 1,10: 1.5
 ln() initial accum: 8.5 (1)
   0: num_times: 2, 4.
@@ -196,6 +196,83 @@ ln() remainder: 0.953729349636862
 0: ln() remainder: 1.8971199848858812
 ln() accum after summing: 1.8971199848858812
 ln() exponent from initial: 1.
-ln(a): 0.176091259055681285
+log(a): 0.176091259055681285
-     : 5.5681285E-11
+     : 4.3E-17
+
+  a  : 4.4
+exp() num_times for 10*ln(10): 4.4 (0)
+exp() num_times for 255: 2.09741490700595432 (1)
+exp() num_times for 0: 0.017973365326118411 (3)
+exp() num_times for 1: 0.017973365326118411 (0)
+exp() num_times for 2: 0.0080230344729503282 (1)
+exp() num_times for 3: 2.703180828206292E-5 (8)
+exp() num_times for 4: 2.703180828206292E-5 (0)
+exp() num_times for 5: 7.0319082813962602E-6 (2)
+exp() num_times for 6: 3.191178139394649E-8 (7)
+exp() num_times for 7: 3.191178139394649E-8 (0)
+exp() num_times for 8: 1.9117815439416942E-9 (3)
+exp() current val for 255: 10.0000000191178154
+exp() current val for 0: 80.0000001529425232
+exp() current val for 1: 80.0000001529425232
+exp() current val for 2: 80.8000001544719485
+exp() current val for 3: 81.4486670861725846
+exp() current val for 4: 81.4486670861725846
+exp() current val for 5: 81.4502960676591747
+exp() current val for 6: 81.450866221442107
+exp() current val for 7: 81.450866221442107
+exp() current val for 8: 81.450868664968118
+exp(a): 81.450868664968118
+      : 6.E-16
+
+  a  : 3.
+ln() accum scaled between 1,10: 3.
+ln() initial accum: 7. (1)
+  0: num_times: 1, 4.
+  1: num_times: 5, 3.3694
+  2: num_times: 3, 4.413961894
+  3: num_times: 4, 4.2564024200762553
+  4: num_times: 4, 2.5750519644445702
+  5: num_times: 2, 5.7509346574136428
+  6: num_times: 5, 7.5095341213443924
+  7: num_times: 7, 5.0953727801950431
+  8: num_times: 5, 0.953729349636862
+ln() remainder: 0.953729349636862
+8: ln() remainder: 5.0953729099644855
+7: ln() remainder: 7.509536940996642
+6: ln() remainder: 5.7509511941013168
+5: ln() remainder: 2.5750851194767976
+4: ln() remainder: 4.257308525280013
+3: ln() remainder: 4.4237321848621339
+2: ln() remainder: 3.4274724744366382
+1: ln() remainder: 5.1082562376599068
+0: ln() remainder: 1.2039728043259359
+ln() accum after summing: 1.2039728043259359
+ln() exponent from initial: 1.
+ln(a): 1.09861228866810978
+     : 9.E-17
+
+  a  : 220.821070022290065
+exp() num_times for 10*ln(10): 13.5884116528259586 (90)
+exp() num_times for 255: 2.07548618785573036 (5)
+exp() num_times for 0: 0.68919182673583976 (2)
+exp() num_times for 1: 0.02202056810556574 (7)
+exp() num_times for 2: 0.0021199063992295744 (2)
+exp() num_times for 3: 1.2090573306250808E-4 (2)
+exp() num_times for 4: 2.0910732729199747E-5 (1)
+exp() num_times for 5: 9.108327285330872E-7 (2)
+exp() num_times for 6: 9.108327285330872E-7 (0)
+exp() num_times for 7: 1.0832773533062324E-8 (9)
+exp() num_times for 8: 8.327735830607254E-10 (1)
+exp() current val for 255: 1.00000000083277358E95
+exp() current val for 0: 4.00000000333109432E95
+exp() current val for 1: 7.79486840649136042E95
+exp() current val for 2: 7.95154526146183677E95
+exp() current val for 3: 7.96745630353002189E95
+exp() current val for 4: 7.96825304916037489E95
+exp() current val for 5: 7.96841241501818339E95
+exp() current val for 6: 7.96841241501818339E95
+exp() current val for 7: 7.9684195865922255E95
+exp() current val for 8: 7.96841966627642136E95
+exp(a): 7.96841966627642136E95
+      : 1.7828E82
 

src/decn/proto/.gitignore

@@ -1,3 +1,4 @@
 div
 ln
+exp
 

src/decn/proto/Makefile

@@ -1,4 +1,4 @@
-all: div ln
+all: div ln exp
 
 div: div_mfp.cpp
 	g++ -Wall div_mfp.cpp -o div -lgmp
@@ -6,3 +6,5 @@ div: div_mfp.cpp
 ln: ln_mfp.cpp
 	g++ -Wall ln_mfp.cpp -o ln -lmpfr
 
+exp: exp.cpp
+	g++ -Wall exp.cpp -o exp -lmpfr

src/decn/proto/exp.cpp

@@ -0,0 +1,130 @@
+//prototype for exponential function
+// complement of ln() function based on HP Journal article
+// "Personal Calculator Algorithms IV: Logarithmic Functions"
+// (uses same set of constants)
+
+#include <stdio.h>
+#include <iostream>
+#include <math.h>
+#include <boost/multiprecision/mpfr.hpp>
+using namespace boost::multiprecision;
+
+// #define DEBUG
+// #define DEBUG_ALL
+
+using std::cout;
+using std::endl;
+
+static const int NUM_A_ARR = 9;
+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);
+	mpfr_float ln10 (10, CALC_PRECISION);
+	ln10 = log( ln10 );
+	cout << "ln(10):" << ln10 << endl;
+	//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;
+	}
+	cout << endl << endl ;
+	//loop through values to test
+	for (mpfr_float x(0.0001, CALC_PRECISION); x < 230.0; x += 0.015625){
+// 	for (mpfr_float x(230, CALC_PRECISION); x < 230.01; x += 0.015625){
+		//build mpf values
+		mpfr_float accum(x, CALC_PRECISION);
+		mpfr_float save (0, CALC_PRECISION);
+		mpfr_float_1000 actual;
+		actual = exp(x);
+
+		int num_times_ln10;
+		int num_times[NUM_A_ARR];
+		int k_j = 0;
+		//track number of times 10*ln(10) subtracted
+		while (accum > 0){
+			save = accum;
+			accum -= 10*ln10;
+			k_j++;
+		}
+		num_times_ln10 = (k_j - 1) * 10;
+		accum = save; //restore
+		//track number of times ln(10) subtracted
+		k_j = 0;
+		while (accum > 0){
+			save = accum;
+			accum -= ln10;
+			k_j++;
+		}
+		num_times_ln10 += k_j - 1;
+		accum = save; //restore
+#ifdef DEBUG
+		printf("  10: num_times: %d, ", num_times_ln10);
+		cout << accum << endl;
+#endif
+		//track number of times ln_a_arr[j] subtracted
+		for (int j = 0; j < NUM_A_ARR; j++){
+			k_j = 0;
+			while (accum > 0){
+				save = accum;
+				accum -= ln_a_arr[j];
+#ifdef DEBUG_ALL
+				printf("     %d: ", k_j);
+				cout << accum << endl;
+#endif
+				k_j++;
+			}
+			num_times[j] = k_j - 1;
+			accum = save; //restore
+#ifdef DEBUG
+			printf("  %d: num_times: %d, ", j, num_times[j]);
+			cout << accum << endl;
+#endif
+		}
+
+		//build final value
+		accum = 1+save; //initialize product to 1+remainder
+		//multiplies for ln(10)
+		for (int k = 0; k < num_times_ln10; k++){
+			accum *= 10;
+		}
+		//multiplies for ln(1 + 10^-j)
+		for (int j = 0; j < NUM_A_ARR; j++){
+			for (int k = 0; k < num_times[j]; k++){
+				accum *= a_arr[j];
+			}
+		}
+
+
+		//calculate relative error
+		mpfr_float_1000 calc, diff;
+		calc = accum;
+		diff = abs((actual - calc)/actual);
+
+#ifdef DEBUG_ALL
+		if (1){
+#else
+		if (diff > 5e-17){
+#endif
+			cout << x << ": ";
+			cout << std::setprecision(18) << accum << ", ";
+#ifdef DEBUG_ALL
+			cout << actual << ", ";
+#endif
+			cout << diff << endl;
+		}
+	}
+
+	return 0;
+}
+

src/decn/proto/ln_mfp.cpp

@@ -80,7 +80,7 @@ int main(void){
 		//calculate relative error
 		mpfr_float_1000 calc, diff;
 		calc = accum;
-		diff = (actual - calc)/actual;
+		diff = abs((actual - calc)/actual);
 
 		if (diff > 5e-17){
 			cout << x << ": ";