From 5bcfb8a37128ae5457583016ff5c58717d6178e0 Mon Sep 17 00:00:00 2001 From: Jeff Wang Date: Sat, 30 Jan 2021 18:45:44 -0500 Subject: [PATCH] do all trig calculations in degrees --- src/decn/decn.c | 48 +++-- src/decn/decn_tests_trig.cpp | 342 ++++++++++++++++++++++------------- 2 files changed, 243 insertions(+), 147 deletions(-) diff --git a/src/decn/decn.c b/src/decn/decn.c index 53f5ae1..6e9391f 100644 --- a/src/decn/decn.c +++ b/src/decn/decn.c @@ -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(); } diff --git a/src/decn/decn_tests_trig.cpp b/src/decn/decn_tests_trig.cpp index af1166d..d2243f9 100644 --- a/src/decn/decn_tests_trig.cpp +++ b/src/decn/decn_tests_trig.cpp @@ -22,6 +22,8 @@ namespace bmp = boost::multiprecision; using Catch::Matchers::Equals; +static const bmp::mpfr_float mPI = boost::math::constants::pi(); + 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 distrib(0,99); - std::uniform_int_distribution exp_distrib(-1,0); //restrict range for now + std::uniform_int_distribution exp_distrib(0, 2); //restrict range for now std::uniform_int_distribution 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 distrib(0,99); - std::uniform_int_distribution exp_distrib(-1,0); //restrict range for now + std::uniform_int_distribution exp_distrib(0, 2); //restrict range for now std::uniform_int_distribution 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 distrib(0,99); - std::uniform_int_distribution exp_distrib(-1,0); //restrict range for now + std::uniform_int_distribution exp_distrib(0, 2); //restrict range for now std::uniform_int_distribution 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); }