7u83/stc_rpncalc
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

do all trig calculations in degrees

Browse Source
This commit is contained in:
Jeff Wang 2021-01-30 18:45:44 -05:00
parent 1e6d786482
commit 5bcfb8a371
2 changed files with 243 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
src/decn/decn.c
View File

@ -85,9 +85,9 @@ const dec80 DECN_LN_10 = {
0, {23, 2, 58, 50, 92, 99, 40, 45, 68}
};
// 2 pi
const dec80 DECN_2PI = {
0, {62, 83, 18, 53, 7, 17, 95, 86, 48}
// pi
const dec80 DECN_PI = {
0, {31, 41, 59, 26, 53, 58, 97, 93, 24}
};
// pi/2
@ -1426,8 +1426,8 @@ void sqrt_decn(void){
#endif //USE_POW_SQRT_IMPL
// see W.E. Egbert, "Personal Calculator Algorithms II: Trigonometric functions"
void project_decn_into_0_2pi(void) {
// normal angle to between 0 and 360 degrees
void normalize_0_360(void) {
const uint8_t is_negative = (AccDecn.exponent < 0);
exp_t exponent;
@ -1436,11 +1436,14 @@ void project_decn_into_0_2pi(void) {
negate_decn(&AccDecn);
}
exponent = get_exponent(&AccDecn);
copy_decn(&BDecn, &DECN_2PI);
//B = 360
set_dec80_zero(&BDecn);
BDecn.lsu[0] = 36;
BDecn.exponent = 2;
if (compare_magn() > 0) {
do {
do {
copy_decn(&BDecn, &DECN_2PI);
//B = 3.6e...
BDecn.exponent = exponent;
if (compare_magn() >= 0) {
negate_decn(&BDecn);
@ -1450,12 +1453,13 @@ void project_decn_into_0_2pi(void) {
}
} while (1);
exponent--;
} while (exponent >= 0);
} while (exponent >= 2);
}
if (is_negative) {
negate_decn(&AccDecn);
copy_decn(&BDecn, &DECN_2PI);
//B = 360
BDecn.exponent = 2;
add_decn();
}
}
@ -1467,13 +1471,14 @@ void project_decn_into_0_2pi(void) {
#define THETA Tmp4Decn
void sincos_decn(const uint8_t sincos_arctan) {
const uint8_t is_negative = AccDecn.exponent < 0;
if (sincos_arctan) {
if (sincos_arctan) { //calculate arctan
set_dec80_zero(&THETA);
if (is_negative) negate_decn(&AccDecn);
copy_decn(&COS, &AccDecn);
set_decn_one(&SIN);
} else {
project_decn_into_0_2pi();
} else { //calculate sin/cos
normalize_0_360();
to_radian_decn();
copy_decn(&THETA, &AccDecn);
set_decn_one(&COS);
set_dec80_zero(&SIN);
@ -1482,13 +1487,13 @@ void sincos_decn(const uint8_t sincos_arctan) {
negate_decn(&SIN);
}
do {
if (sincos_arctan) {
if (sincos_arctan) { //calculate arctan
// THETA is in AccDecn from previous iteration
if (COS.exponent < 0) {
if (is_negative) negate_decn(&AccDecn);
break;
}
} else {
} else { //calculate sin/cos
if (THETA.exponent < 0) {
break;
}
@ -1539,10 +1544,11 @@ void tan_decn(void) {
void arctan_decn(void) {
sincos_decn(1);
to_degree_decn();
}
// see W.E. Egbert, "Personal Calculator Algorithms III: Inverse Trigonometric Functions"
void arcsin_decn(void) {
void arcsin_decn_rad(void) {
st_push_decn(&AccDecn);
copy_decn(&BDecn, &AccDecn);
mult_decn();
@ -1553,15 +1559,20 @@ void arcsin_decn(void) {
recip_decn();
st_pop_decn(&BDecn);
mult_decn();
sincos_decn(1);
}
void arcsin_decn(void) {
arcsin_decn_rad();
to_degree_decn();
}
void arccos_decn(void) {
arcsin_decn();
arcsin_decn_rad();
negate_decn(&AccDecn);
copy_decn(&BDecn, &DECN_PI2);
add_decn();
to_degree_decn();
}
#undef SIN
#undef COS
@ -1579,8 +1590,7 @@ void to_radian_decn(void) {
void pi_decn(void) {
set_dec80_zero(&BDecn);
BDecn.lsu[0] = 5; // 0.5 00 ..
copy_decn(&AccDecn, &DECN_2PI);
copy_decn(&AccDecn, &DECN_PI);
mult_decn();
}

342
src/decn/decn_tests_trig.cpp
View File

@ -22,6 +22,8 @@
namespace bmp = boost::multiprecision;
using Catch::Matchers::Equals;
static const bmp::mpfr_float mPI = boost::math::constants::pi<bmp::mpfr_float>();
static void trig_test(void (*operation)(void), bmp::mpfr_float (*mpfr_operation)(bmp::mpfr_float x),
double rtol, double atol)
@ -50,10 +52,10 @@ static void trig_test(void (*operation)(void), bmp::mpfr_float (*mpfr_operation)
}
static void sin_test(double rtol=5e-3, double atol=1e-3)
static void sin_test(double rtol=6e-3, double atol=1e-3)
{
CAPTURE("sin test");
trig_test(sin_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return sin(x);}, rtol, atol);
trig_test(sin_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return sin(x * mPI / 180);}, rtol, atol);
}
static void sin_test(const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3)
@ -64,10 +66,10 @@ static void sin_test(const char* a_str, int a_exp, double rtol=5e-3, double atol
}
static void cos_test(double rtol=5e-3, double atol=1e-3)
static void cos_test(double rtol=6e-3, double atol=1e-3)
{
CAPTURE("cos test");
trig_test(cos_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return cos(x);}, rtol, atol);
trig_test(cos_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return cos(x * mPI / 180);}, rtol, atol);
}
static void cos_test(const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3)
@ -80,7 +82,7 @@ static void cos_test(const char* a_str, int a_exp, double rtol=5e-3, double atol
static void tan_test(double rtol=5e-3, double atol=1e-3)
{
trig_test(tan_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return tan(x);}, rtol, atol);
trig_test(tan_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return tan(x * mPI / 180);}, rtol, atol);
}
static void tan_test(const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3)
@ -93,7 +95,7 @@ static void tan_test(const char* a_str, int a_exp, double rtol=5e-3, double atol
static void atan_test(double rtol=5e-3, double atol=1e-3)
{
trig_test(arctan_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return atan(x);}, rtol, atol);
trig_test(arctan_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return atan(x)*180/mPI;}, rtol, atol);
}
static void atan_test(const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3)
@ -106,7 +108,7 @@ static void atan_test(const char* a_str, int a_exp, double rtol=5e-3, double ato
static void asin_test(double rtol=5e-3, double atol=1e-3)
{
trig_test(arcsin_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return asin(x);}, rtol, atol);
trig_test(arcsin_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return asin(x)*180/mPI;}, rtol, atol);
}
static void asin_test(const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3)
@ -119,7 +121,7 @@ static void asin_test(const char* a_str, int a_exp, double rtol=5e-3, double ato
static void acos_test(double rtol=5e-3, double atol=1e-3)
{
trig_test(arccos_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return acos(x);}, rtol, atol);
trig_test(arccos_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return acos(x)*180/mPI;}, rtol, atol);
}
static void acos_test(const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3)
@ -129,101 +131,175 @@ static void acos_test(const char* a_str, int a_exp, double rtol=5e-3, double ato
acos_test(rtol, atol);
}
const char * const pi = "3.141592653589793239";
const char * const pi_threequarters = "2.356194490192344929";
const char * const pi_halved = "1.570796326794896619";
const char * const pi_quarter = ".7853981633974483096";
// const char * const pi = "3.141592653589793239";
// const char * const pi_threequarters = "2.356194490192344929";
// const char * const pi_halved = "1.570796326794896619";
// const char * const pi_quarter = ".7853981633974483096";
TEST_CASE("sin") {
sin_test("0.1", 0);
sin_test("0.05", 0, 1e-2);
sin_test("0.01", 0, -1);
sin_test("0.001", 0, -1);
sin_test("0.0001", 0, -1);
sin_test("0.00001", 0, -1);
sin_test("0.000001", 0, -1);
sin_test("0.1", 0, 0.2);
sin_test("0.0", 0, -1);
sin_test("0.2", 0);
sin_test("0.3", 0);
sin_test("0.4", 0);
sin_test("0.9", 0);
sin_test("1.5", 0);
sin_test("2.0", 0);
sin_test("2.5", 0);
sin_test("3.0", 0);
sin_test(pi, 0, -1);
sin_test(pi_quarter, 0);
sin_test(pi_halved, 0);
sin_test(pi_threequarters, 0);
sin_test("1.5", 0, 0.02);
sin_test("2.0", 0, 0.02);
sin_test("2.5", 0, 0.02);
sin_test("3.0", 0, 0.02);
sin_test("10", 0);
sin_test("20", 0);
sin_test("30", 0);
sin_test("40", 0);
sin_test("80", 0);
sin_test("120", 0);
sin_test("160", 0);
sin_test("200", 0);
sin_test("240", 0);
sin_test("280", 0);
sin_test("320", 0);
sin_test("359", 0, 0.02);
sin_test("360", 0, -1, 0.001);
sin_test("361", 0, 0.02);
sin_test("400", 0);
// sin_test(pi, 0, -1);
// sin_test(pi_quarter, 0);
// sin_test(pi_halved, 0);
// sin_test(pi_threequarters, 0);
sin_test("180.0", 0, -1);
sin_test("45.0", 0);
sin_test("90.0", 0);
sin_test("135.0", 0);
sin_test("1000.0", 0);
sin_test("-0.5", 0);
sin_test("-1.5", 0);
sin_test("-2.0", 0);
sin_test("-2.5", 0);
sin_test("-3.0", 0);
sin_test("-0.5", 0, 0.2);
sin_test("-1.5", 0, 0.02);
sin_test("-2.0", 0, 0.02);
sin_test("-2.5", 0, 0.02);
sin_test("-3.0", 0, 0.02);
sin_test("-9.0", 0);
sin_test("-18.0", 0);
sin_test("-27.0", 0);
sin_test("-1000.0", 0);
sin_test("-30", 0);
sin_test("-40", 0);
sin_test("-80", 0);
sin_test("-120", 0);
sin_test("-160", 0);
sin_test("-200", 0);
sin_test("-240", 0);
sin_test("-280", 0);
sin_test("-320", 0);
sin_test("-360", 0, -1, 0.001);
sin_test("-400", 0);
}
TEST_CASE("cos") {
cos_test("0.1", 0);
cos_test("0.05", 0);
cos_test("0.01", 0);
cos_test("0.001", 0);
cos_test("0.0001", 0);
cos_test("0.00001", 0);
cos_test("0.000001", 0);
cos_test("0.0", 0);
cos_test("0.2", 0);
cos_test("0.3", 0);
cos_test("0.4", 0);
cos_test("0.9", 0);
cos_test("1.5", 0);
cos_test("2.0", 0);
cos_test("2.5", 0);
cos_test("3.0", 0);
cos_test(pi, 0);
cos_test(pi_quarter, 0);
cos_test(pi_halved, 0, -1);
cos_test(pi_threequarters, 0);
cos_test("1000.0", 0);
cos_test("-0.5", 0);
cos_test("-1.5", 0);
cos_test("-2.0", 0);
cos_test("-2.5", 0);
cos_test("-3.0", 0);
cos_test("0.1", 0, 0.2);
cos_test("0.0", 0, -1);
cos_test("1.5", 0, 0.02);
cos_test("2.0", 0, 0.02);
cos_test("2.5", 0, 0.02);
cos_test("3.0", 0, 0.02);
cos_test("10", 0);
cos_test("20", 0);
cos_test("30", 0);
cos_test("40", 0);
cos_test("80", 0);
cos_test("120", 0);
cos_test("160", 0);
cos_test("200", 0);
cos_test("240", 0);
cos_test("280", 0, 0.006);
cos_test("320", 0);
cos_test("359", 0, 0.02);
cos_test("360", 0, -1, 0.001);
cos_test("361", 0, 0.02);
cos_test("400", 0);
// cos_test(pi, 0, -1);
// cos_test(pi_quarter, 0);
// cos_test(pi_halved, 0);
// cos_test(pi_threequarters, 0);
cos_test("180.0", 0, -1);
cos_test("45.0", 0);
cos_test("90.0", 0, -1, 0.001);
cos_test("135.0", 0);
cos_test("1000.0", 0, 0.006);
cos_test("-0.5", 0, 0.2);
cos_test("-1.5", 0, 0.02);
cos_test("-2.0", 0, 0.02);
cos_test("-2.5", 0, 0.02);
cos_test("-3.0", 0, 0.02);
cos_test("-9.0", 0);
cos_test("-18.0", 0);
cos_test("-27.0", 0);
cos_test("-1000.0", 0);
cos_test("-30", 0);
cos_test("-40", 0);
cos_test("-80", 0, 0.006);
cos_test("-120", 0);
cos_test("-160", 0);
cos_test("-200", 0);
cos_test("-240", 0);
cos_test("-280", 0);
cos_test("-320", 0);
cos_test("-360", 0, -1, 0.001);
cos_test("-400", 0);
}
TEST_CASE("tan") {
tan_test("0.1", 0);
tan_test("0.05", 0, 1e-2);
tan_test("0.01", 0, 5e-2);
tan_test("0.001", 0, -1);
tan_test("0.0001", 0, -1);
tan_test("0.00001", 0, -1);
tan_test("0.000001", 0, -1);
tan_test("0.1", 0, 0.2);
tan_test("0.0", 0, -1);
tan_test("0.2", 0);
tan_test("0.3", 0);
tan_test("0.4", 0);
tan_test("0.9", 0);
tan_test("1.5", 0);
tan_test("2.0", 0);
tan_test("2.5", 0);
tan_test("3.0", 0);
tan_test("1.5", 0, 0.02);
tan_test("2.0", 0, 0.02);
tan_test("2.5", 0, 0.02);
tan_test("3.0", 0, 0.02);
tan_test("10", 0);
tan_test("20", 0);
tan_test("30", 0);
tan_test("40", 0);
tan_test("80", 0);
tan_test("120", 0);
tan_test("160", 0);
tan_test("200", 0);
tan_test("240", 0);
tan_test("280", 0, 0.006);
tan_test("320", 0);
tan_test("359", 0, 0.02);
tan_test("360", 0, -1, 0.001);
tan_test("361", 0, 0.02);
tan_test("400", 0);
// tan_test(pi, 0, -1);
// tan_test(pi_quarter, 0);
// tan_test(pi_halved, 0);
// tan_test(pi_threequarters, 0);
tan_test("180.0", 0, -1);
tan_test("45.0", 0);
tan_test("90.0", 0, 2);
tan_test("135.0", 0);
tan_test("1000.0", 0, 0.006);
tan_test("-0.5", 0, 0.2);
tan_test("-1.5", 0, 0.02);
tan_test("-2.0", 0, 0.02);
tan_test("-2.5", 0, 0.02);
tan_test("-3.0", 0, 0.02);
tan_test("-9.0", 0);
tan_test("-18.0", 0);
tan_test("-27.0", 0);
tan_test("-1000.0", 0);
tan_test("-30", 0);
tan_test("-40", 0);
tan_test("-80", 0, 0.006);
tan_test("-120", 0);
tan_test("-160", 0);
tan_test("-200", 0);
tan_test("-240", 0);
tan_test("-280", 0);
tan_test("-320", 0);
tan_test("-360", 0, -1, 0.001);
tan_test("-400", 0);
}
TEST_CASE("arctan") {
atan_test("0.001", 0, -1, 2e-3);
atan_test("-0.001", 0, -1, 2e-3);
atan_test("0.001", 0, -1, 0.06);
atan_test("-0.001", 0, -1, 0.06);
atan_test("0.7", 0);
atan_test("-0.7", 0);
atan_test("0.1", 0);
@ -234,12 +310,12 @@ TEST_CASE("arctan") {
atan_test("-2.0", 0);
atan_test("3.0", 0);
atan_test("-3.0", 0);
atan_test("0", 0, -1);
atan_test("0", 0, -1, 0.06);
}
TEST_CASE("arcsin") {
asin_test("0.001", 0, -1);
asin_test("-0.001", 0, -1);
asin_test("0.001", 0, -1, 0.06);
asin_test("-0.001", 0, -1, 0.06);
asin_test("0.7", 0);
asin_test("-0.7", 0);
asin_test("0.1", 0, 1e-2);
@ -265,7 +341,7 @@ static const int NUM_RAND_TRIG_TESTS = 4321; //trig tests are slow
TEST_CASE("sin random"){
std::default_random_engine gen;
std::uniform_int_distribution<int> distrib(0,99);
std::uniform_int_distribution<int> exp_distrib(-1,0); //restrict range for now
std::uniform_int_distribution<int> exp_distrib(0, 2); //restrict range for now
std::uniform_int_distribution<int> sign_distrib(0,1);
for (int j = 0; j < NUM_RAND_TRIG_TESTS; j++){
for (int i = 0; i < DEC80_NUM_LSU; i++){
@ -274,25 +350,27 @@ TEST_CASE("sin random"){
int exp = exp_distrib(gen);
int sign = sign_distrib(gen);
set_exponent(&AccDecn, exp, sign);
remove_leading_zeros(&AccDecn);
int lsu0 = AccDecn.lsu[0];
exp = get_exponent(&AccDecn);
CAPTURE(lsu0);
CAPTURE(exp);
CAPTURE(sign);
if (exp == -1 && lsu0 == 0){
if (exp <= -1){
//very small
sin_test(40);
} else if ((exp == -1 && lsu0 < 10) || (exp == 0 && lsu0 == 0)){
sin_test(4000);
} else if (exp == 0 && lsu0 < 50){
//small
sin_test(0.4);
} else if ((exp == 0 && lsu0 == 31)){
//near pi
sin_test(0.2);
} else if ((exp == 0 && lsu0 == 62)){
//near 2pi
sin_test(0.2);
} else if ((exp == 0 && lsu0 > 62)){
} else if (exp == 2 && lsu0 >= 17 && lsu0 <= 19){
//near 180
sin_test(3);
} else if (exp == 2 && lsu0 >= 35 && lsu0 <= 36){
//near 360
sin_test(12);
} else if (exp == 2 && lsu0 > 50){
//large
sin_test(0.1);
sin_test(35);
} else {
sin_test(0.02);
}
@ -302,7 +380,7 @@ TEST_CASE("sin random"){
TEST_CASE("cos random"){
std::default_random_engine gen;
std::uniform_int_distribution<int> distrib(0,99);
std::uniform_int_distribution<int> exp_distrib(-1,0); //restrict range for now
std::uniform_int_distribution<int> exp_distrib(0, 2); //restrict range for now
std::uniform_int_distribution<int> sign_distrib(0,1);
for (int j = 0; j < NUM_RAND_TRIG_TESTS; j++){
for (int i = 0; i < DEC80_NUM_LSU; i++){
@ -311,20 +389,24 @@ TEST_CASE("cos random"){
int exp = exp_distrib(gen);
int sign = sign_distrib(gen);
set_exponent(&AccDecn, exp, sign);
remove_leading_zeros(&AccDecn);
int lsu0 = AccDecn.lsu[0];
exp = get_exponent(&AccDecn);
CAPTURE(lsu0);
CAPTURE(exp);
CAPTURE(sign);
if (exp == 0 && lsu0 == 15){
//near pi/2
cos_test(0.4);
} else if (exp == 0 && lsu0 == 47){
//near 3/2 * pi
cos_test(0.4);
} else if (exp == 0 && lsu0 == 78){
//near 5/2 * pi
// cos_test(0.4);
cos_test(1.1); //actual rtol is much worse than 0.4, random test happens to hit a bad one
if (exp == 1 && lsu0 >= 89 && lsu0 <= 90){
//very near 90
cos_test(500);
} else if (exp == 1 && lsu0 >= 87 && lsu0 <= 92){
//near 90
cos_test(2);
} else if (exp == 2 && lsu0 >= 26 && lsu0 <= 27){
//near 270
cos_test(500);
} else if (exp == 2 && lsu0 >= 44){
//large
cos_test(20);
} else {
cos_test(0.02);
}
@ -334,7 +416,7 @@ TEST_CASE("cos random"){
TEST_CASE("tan random"){
std::default_random_engine gen;
std::uniform_int_distribution<int> distrib(0,99);
std::uniform_int_distribution<int> exp_distrib(-1,0); //restrict range for now
std::uniform_int_distribution<int> exp_distrib(0, 2); //restrict range for now
std::uniform_int_distribution<int> sign_distrib(0,1);
for (int j = 0; j < NUM_RAND_TRIG_TESTS; j++){
for (int i = 0; i < DEC80_NUM_LSU; i++){
@ -343,35 +425,39 @@ TEST_CASE("tan random"){
int exp = exp_distrib(gen);
int sign = sign_distrib(gen);
set_exponent(&AccDecn, exp, sign);
remove_leading_zeros(&AccDecn);
int lsu0 = AccDecn.lsu[0];
exp = get_exponent(&AccDecn);
CAPTURE(lsu0);
CAPTURE(exp);
CAPTURE(sign);
if (exp == -1 && lsu0 == 0){
if (exp <= -3){
//extremely small
tan_test(5000);
} if (exp <= -1){
//very small
tan_test(40);
} else if ((exp == -1 && lsu0 < 10) || (exp == 0 && lsu0 == 0)){
tan_test(400);
} else if (exp == 0 && lsu0 < 50){
//small
tan_test(0.5);
} else if (exp == 0 && lsu0 == 15){
//near pi/2
tan_test(0.5);
} else if ((exp == 0 && lsu0 == 31)){
//near pi
tan_test(0.2);
} else if (exp == 0 && lsu0 == 47){
//near 3/2 * pi
tan_test(0.5);
} else if ((exp == 0 && lsu0 == 62)){
//near 2pi
tan_test(0.2);
} else if (exp == 0 && lsu0 == 78){
//near 5/2 * pi
// tan_test(0.5);
tan_test(0.6); //actual rtol is much worse than 0.4, random test happens to hit a bad one
} else if ((exp == 0 && lsu0 > 62)){
tan_test(1);
} else if (exp == 1 && lsu0 >= 89 && lsu0 <= 90){
//very near 90
tan_test(5);
} else if (exp == 1 && lsu0 >= 87 && lsu0 <= 92){
//near 90
tan_test(1);
} else if (exp == 2 && lsu0 >= 17 && lsu0 <= 19){
//near 180
tan_test(3);
} else if (exp == 2 && lsu0 >= 26 && lsu0 <= 27){
//near 270
tan_test(5);
} else if (exp == 2 && lsu0 >= 35 && lsu0 <= 37){
//near 360
tan_test(20);
} else if (exp == 2 && lsu0 >= 44){
//large
tan_test(0.1);
tan_test(50);
} else {
tan_test(0.02);
}