From 2903579efd17a59d950ea04af5b3639ddf646c2d Mon Sep 17 00:00:00 2001 From: Mirko Scholz Date: Thu, 10 Sep 2020 18:49:47 +0200 Subject: [PATCH] added arcsin and arccos added a stack for temporaries, otherwise it became to tedious to keep track which function uses which temporary for later --- qt_gui/main.qml | 2 +- src/calc.c | 27 ++++------ src/decn/decn.c | 98 ++++++++++++++++++++++++++++++------ src/decn/decn.h | 9 +++- src/decn/decn_tests_trig.cpp | 34 +++++++++++++ 5 files changed, 135 insertions(+), 35 deletions(-) diff --git a/qt_gui/main.qml b/qt_gui/main.qml index e3bb33c..98374b3 100644 --- a/qt_gui/main.qml +++ b/qt_gui/main.qml @@ -180,7 +180,7 @@ ApplicationWindow ["", "", "", ""], ["", "", "", ""], ["R▲", "", "", ""], - ["", "", "tan−1", "►deg"], + ["sin−1", "cos−1", "tan−1", "►deg"], ["", "", "", ""] ] diff --git a/src/calc.c b/src/calc.c index 69387be..bb45688 100644 --- a/src/calc.c +++ b/src/calc.c @@ -23,6 +23,10 @@ #include "calc.h" #include "stack_debug.h" +#ifdef DESKTOP +#include +#endif + __xdata dec80 StoredDecn; __xdata dec80 LastX; @@ -163,21 +167,7 @@ void process_cmd(char cmd){ case '<':{ //use as +/- and sqrt if (IsShiftedUp){ //take sqrt IsShiftedUp = 0; - if (decn_is_zero(&stack(STACK_X))){ - //sqrt(0) = 0 - } else if (!decn_is_nan(&stack(STACK_X))){ - copy_decn(&LastX, &stack(STACK_X)); //save LastX - copy_decn(&AccDecn, &stack(STACK_X)); - if (AccDecn.exponent < 0){ //negative - set_dec80_NaN(&stack(STACK_X)); - break; - } - //b = 0.5 - set_dec80_zero(&BDecn); - BDecn.lsu[0] = 5; - pow_decn(); - copy_decn(&stack(STACK_X), &AccDecn); - } + do_unary_op(sqrt_decn); } else { // +/- if (!decn_is_nan(&stack(STACK_X))){ negate_decn(&stack(STACK_X)); @@ -215,7 +205,7 @@ void process_cmd(char cmd){ if (IsShiftedUp){ do_unary_op(sin_decn); } else if (IsShiftedDown){ - // do_unary_op(arcsin_decn); + do_unary_op(arcsin_decn); } } break; ////////// @@ -223,7 +213,7 @@ void process_cmd(char cmd){ if (IsShiftedUp){ do_unary_op(cos_decn); } else if (IsShiftedDown){ - // do_unary_op(arccos_decn); + do_unary_op(arccos_decn); } } break; ////////// @@ -283,6 +273,9 @@ void process_cmd(char cmd){ } //switch(cmd) IsShiftedUp = 0; IsShiftedDown = 0; +#ifdef DESKTOP + assert(TmpStackPtr == 0); // there should be no items on the temporaries stack after one global operation +#endif } diff --git a/src/decn/decn.c b/src/decn/decn.c index 6419b84..68360b3 100644 --- a/src/decn/decn.c +++ b/src/decn/decn.c @@ -67,10 +67,14 @@ static const uint8_t num_digits_display = 16; dec80 AccDecn; __idata dec80 BDecn; -__idata dec80 TmpDecn; //used by add_decn() and mult_decn() -__idata dec80 Tmp2Decn; //used by recip_decn() and ln_decn() and sincos_decn() -__xdata dec80 Tmp3Decn; //used by recip_decn() and ln_decn() and sincos_decn() -__xdata dec80 Tmp4Decn; //used by div_decn() and pow_decn() and sincos_decn() +__idata dec80 TmpDecn; //used by add_decn() and mult_decn() and sqrt_decn() +__idata dec80 Tmp2Decn; //used by recip_decn() and ln_decn() and sincos_decn() and sqrt_decn() +__xdata dec80 Tmp3Decn; //used by recip_decn() and ln_decn() and sincos_decn() and sqrt_decn() +__xdata dec80 Tmp4Decn; //used by sincos_decn() + +__xdata dec80 TmpStackDecn[4]; +#define TMP_STACK_SIZE (sizeof TmpStackDecn / sizeof TmpStackDecn[0]) +__idata uint8_t TmpStackPtr; __xdata char Buf[DECN_BUF_SIZE]; @@ -89,11 +93,36 @@ const dec80 DECN_2PI = { 0, {62, 83, 18, 53, 7, 17, 95, 86, 48} }; +// pi/2 +const dec80 DECN_PI2 = { + 0, {15, 70, 79, 63, 26, 79, 48, 96, 62} +}; + // 180/pi = 1rad in degree const dec80 DECN_1RAD = { 1, {57, 29, 57, 79, 51, 30, 82, 32, 9} }; +void st_push_decn(const dec80 * const src) +{ + copy_decn(&TmpStackDecn[TmpStackPtr], src); + TmpStackPtr++; + assert(TmpStackPtr < TMP_STACK_SIZE); +} + +void st_pop_decn(dec80 * const dst) +{ + assert(TmpStackPtr >= 1); + TmpStackPtr--; + if (dst) copy_decn(dst, &TmpStackDecn[TmpStackPtr]); +} + +void st_load_decn(dec80 * const dst) +{ + assert(TmpStackPtr >= 1); + copy_decn(dst, &TmpStackDecn[TmpStackPtr - 1]); +} + void copy_decn(dec80* const dest, const dec80* const src){ uint8_t i; @@ -599,6 +628,7 @@ void add_decn(void){ return; } //save b for restoring later + //n.b. don't use TmpStackDecn here, it is called quite often. So you'd need to increase TMP_STACK_SIZE copy_decn(&TmpDecn, &BDecn); //handle cases where signs differ if (AccDecn.exponent < 0 && BDecn.exponent >= 0){ @@ -805,7 +835,6 @@ void mult_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 uint8_t i; exp_t initial_exp; //check divide by zero @@ -821,7 +850,7 @@ void recip_decn(void){ //normalize remove_leading_zeros(&AccDecn); //store copy of x - copy_decn(&X_COPY, &AccDecn); + st_push_decn(&AccDecn); //get initial exponent of estimate for 1/x initial_exp = get_exponent(&AccDecn); #ifdef DEBUG_DIV @@ -854,7 +883,7 @@ void recip_decn(void){ printf("%2d: %s\n", i, Buf); #endif //Accum *= x_copy - copy_decn(&BDecn, &X_COPY); + st_load_decn(&BDecn); mult_decn(); #ifdef DEBUG_DIV decn_to_str_complete(&AccDecn); @@ -881,24 +910,20 @@ void recip_decn(void){ //new_est(Accum) = recip + (1 - recip*x)*recip, where recip is current recip estimate copy_decn(&CURR_RECIP, &AccDecn); } + st_pop_decn(0); //try not to pollute namespace #undef CURR_RECIP -#undef X_COPY } 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); + st_push_decn(&AccDecn); copy_decn(&AccDecn, &BDecn); recip_decn(); //Accum now holds 1/x, multiply by original acc to complete division - copy_decn(&BDecn, &ACC_COPY); + st_pop_decn(&BDecn); mult_decn(); - -//try not to pollute namespace -#undef ACC_COPY } @@ -1256,9 +1281,9 @@ void pow_decn(void) { return; } //calculate AccDecn = AccDecn ^ BDecn - copy_decn(&Tmp4Decn, &BDecn); //save b + st_push_decn(&BDecn); ln_decn(); - copy_decn(&BDecn, &Tmp4Decn); //restore b + st_pop_decn(&BDecn); mult_decn(); //accum = b*ln(accum) exp_decn(); } @@ -1379,6 +1404,29 @@ void tan_decn(void) { void arctan_decn(void) { sincos_decn(1); } + +// see W.E. Egbert, "Personal Calculator Algorithms III: Inverse Trigonometric Functions" +void arcsin_decn(void) { + st_push_decn(&AccDecn); + copy_decn(&BDecn, &AccDecn); + mult_decn(); + negate_decn(&AccDecn); + copy_decn(&BDecn, &DECN_1); + add_decn(); + sqrt_decn(); + recip_decn(); + st_pop_decn(&BDecn); + mult_decn(); + + sincos_decn(1); +} + +void arccos_decn(void) { + arcsin_decn(); + negate_decn(&AccDecn); + copy_decn(&BDecn, &DECN_PI2); + add_decn(); +} #undef SIN #undef COS #undef THETA @@ -1400,6 +1448,24 @@ void pi_decn(void) { mult_decn(); } +void sqrt_decn(void) { + if (decn_is_zero(&AccDecn)) { + return; + } + if (decn_is_nan(&AccDecn)) { + return; + } + if (AccDecn.exponent < 0){ //negative + set_dec80_NaN(&AccDecn); + return; + } + st_push_decn(&BDecn); // sqrt should behave like an unary operation + //b = 0.5 + set_dec80_zero(&BDecn); + BDecn.lsu[0] = 5; + pow_decn(); + st_pop_decn(&BDecn); +} static void set_str_error(void){ Buf[0] = 'E'; diff --git a/src/decn/decn.h b/src/decn/decn.h index cdffaa0..ae04654 100644 --- a/src/decn/decn.h +++ b/src/decn/decn.h @@ -64,7 +64,7 @@ void copy_decn(dec80* const dest, const dec80* const src); extern dec80 AccDecn; extern __idata dec80 BDecn; -extern __xdata dec80 Tmp4Decn; +extern __idata uint8_t TmpStackPtr; void build_dec80(__xdata const char* signif_str, __xdata exp_t exponent); @@ -85,11 +85,14 @@ void log10_decn(void); void exp_decn(void); void exp10_decn(void); void pow_decn(void); +void sqrt_decn(void); void sin_decn(void); void cos_decn(void); void tan_decn(void); void arctan_decn(void); +void arcsin_decn(void); +void arccos_decn(void); void to_degree_decn(void); void to_radian_decn(void); void pi_decn(void); @@ -111,12 +114,16 @@ void decn_to_str_complete(const dec80* x); void build_decn_at(dec80* dest, const char* signif_str, exp_t exponent); #endif +#ifdef DESKTOP #define PRINT_DEC80(n, v) \ printf(n " %d %5d: ", v.exponent < 0, get_exponent(&v)); \ for (int i = 0; i < DEC80_NUM_LSU; i++) { \ printf("%02d ", v.lsu[i]); \ } \ fputc('\n', stdout); +#else +#define PRINT_DEC80(n, v) +#endif #ifdef __cplusplus } diff --git a/src/decn/decn_tests_trig.cpp b/src/decn/decn_tests_trig.cpp index ceba7a7..28236ae 100644 --- a/src/decn/decn_tests_trig.cpp +++ b/src/decn/decn_tests_trig.cpp @@ -59,6 +59,18 @@ static void atan_test( trig_test(arctan_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return atan(x);}, a_str, a_exp, rtol, atol); } +static void asin_test( + const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3) +{ + trig_test(arcsin_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return asin(x);}, a_str, a_exp, rtol, atol); +} + +static void acos_test( + const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3) +{ + trig_test(arccos_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return acos(x);}, a_str, a_exp, rtol, atol); +} + const char * const pi = "3.141592653589793239"; const char * const pi_threequarters = "2.356194490192344929"; @@ -167,3 +179,25 @@ TEST_CASE("arctan") { atan_test("-3.0", 0); atan_test("0", 0, -1); } + +TEST_CASE("arcsin") { + asin_test("0.001", 0, -1); + asin_test("-0.001", 0, -1); + asin_test("0.7", 0); + asin_test("-0.7", 0); + asin_test("0.1", 0); + asin_test("-0.1", 0); + asin_test("0.9", 0); + asin_test("-0.9", 0); +} + +TEST_CASE("arccos") { + acos_test("0.001", 0); + acos_test("-0.001", 0); + acos_test("0.7", 0); + acos_test("-0.7", 0); + acos_test("0.1", 0); + acos_test("-0.1", 0); + acos_test("0.9", 0); + acos_test("-0.9", 0); +}