diff --git a/src/decn/decn.c b/src/decn/decn.c index b3300d1..9af3601 100644 --- a/src/decn/decn.c +++ b/src/decn/decn.c @@ -8,15 +8,17 @@ #include "../utils.h" //#define DEBUG -//#define DEBUG_ADD -#define DEBUG_MULT //#define DEBUG_COMPARE_MAGN +//#define DEBUG_ADD +//#define DEBUG_MULT +//#define DEBUG_DIV #ifndef DESKTOP #undef DEBUG -#undef DEBUG_ADD #undef DEBUG_COMPARE_MAGN +#undef DEBUG_ADD #undef DEBUG_MULT +#undef DEBUG_DIV #endif #ifdef DESKTOP @@ -672,6 +674,54 @@ void mult_decn(dec80* acc, const dec80* x){ //normalize remove_leading_zeros(acc); } + +void div_decn(dec80* acc, const dec80* x){ + dec80 tmp; //copy of x, holds current 1/x estimate + dec80 acc_copy; //holds copy of original acc + uint8_t i; + //check divide by zero + if (decn_is_zero(x)){ + set_dec80_NaN(acc); + return; + } + //store copy of acc for final multiply by 1/x + copy_decn(&acc_copy, acc); + //get initial estimate for 1/x, by negating exponent, and setting signif. to 1 + set_exponent(&tmp, -get_exponent(x), (x->exponent < 0)); + tmp.lsu[0] = 10; //1 with implicit point + for (i = 1; i < DEC80_NUM_LSU; i++){ + tmp.lsu[i] = 0; + } + copy_decn(acc, &tmp); + //do newton raphson iterations + for (i = 0; i < 20; i++){ //just fix number of iterations for now +#ifdef DEBUG_DIV + char buf[80]; + dec80_to_str(buf, &tmp); + printf("%2d: %s\n", i, buf); +#endif + mult_decn(acc, x); +#ifdef DEBUG_DIV + dec80_to_str(buf, acc); + printf(" %20s: %s\n", "recip*x", buf); +#endif + negate_decn(acc); + add_decn(acc, &DECN_1); +#ifdef DEBUG_DIV + dec80_to_str(buf, acc); + printf(" %20s: %s\n", "(1-recip*x)", buf); +#endif + mult_decn(acc, &tmp); +#ifdef DEBUG_DIV + dec80_to_str(buf, acc); + printf(" %20s: %s\n", "recip * (1-recip*x)", buf); +#endif + add_decn(acc, &tmp); + //new_est(acc) = recip + (1 - recip*x)*recip, where tmp is current recip estimate + copy_decn(&tmp, acc); + } + //acc now holds 1/x, multiply by original acc to complete division + mult_decn(acc, &acc_copy); } //buf should hold at least 18 + 4 + 5 + 1 = 28 diff --git a/src/decn/decn.h b/src/decn/decn.h index 04be408..2f1b38e 100644 --- a/src/decn/decn.h +++ b/src/decn/decn.h @@ -23,6 +23,12 @@ typedef struct { //implicit decimal point between (lsu[0]/10) and (lsu[0]%10) } dec80; +//1 constant +static const dec80 DECN_1 = { + 0, + {10, 0} +}; + //remove sign bit, and return 15 bit exponent sign-extended to 16 bits int16_t get_exponent(const dec80* x); @@ -39,6 +45,7 @@ void negate_decn(dec80* x); int8_t compare_decn(dec80* a, dec80* b); //ab: 1 void add_decn(dec80* acc, const dec80* x); void mult_decn(dec80* acc, const dec80* x); +void div_decn(dec80* acc, const dec80* x); //buf should hold at least 18 + 4 + 5 + 1 = 28 #define DECN_BUF_SIZE 28 diff --git a/src/decn/decn_test.c b/src/decn/decn_test.c index bef9385..0c1e54a 100644 --- a/src/decn/decn_test.c +++ b/src/decn/decn_test.c @@ -73,5 +73,17 @@ int main(void){ dec80_to_str(buf, &acc); printf("acc*b: %s\n\n", buf); + //new acc and b for divide test + build_dec80(&acc, "3.14", 88); + build_dec80(&b, "-1.5", -2); + dec80_to_str(buf, &acc); + printf(" acc: %s\n", buf); + dec80_to_str(buf, &b); + printf(" b: %s\n", buf); + div_decn(&acc, &b); + dec80_to_str(buf, &acc); + printf("acc/b: %s\n", buf); + + return 0; }