added arcsin and arccos

added a stack for temporaries, otherwise it became to tedious to keep
track which function uses which temporary for later
This commit is contained in:
Mirko Scholz 2020-09-10 18:49:47 +02:00
parent b56523fadf
commit 2903579efd
5 changed files with 135 additions and 35 deletions

View File

@ -180,7 +180,7 @@ ApplicationWindow
["", "", "", ""],
["", "", "", ""],
["R▲", "", "", ""],
["", "", "tan<sup>1</sup>", "►deg"],
["sin<sup>1</sup>", "cos<sup>1</sup>", "tan<sup>1</sup>", "►deg"],
["", "", "", ""]
]

View File

@ -23,6 +23,10 @@
#include "calc.h"
#include "stack_debug.h"
#ifdef DESKTOP
#include <assert.h>
#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
}

View File

@ -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';

View File

@ -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
}

View File

@ -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);
}