simple trigonometric functions as in Sinclair scientific

src/decn/CMakeLists.txt

@@ -21,7 +21,7 @@ target_link_libraries(decn_cover PUBLIC coverage_config)
 add_executable(decn_test decn_test.c ../utils.c)
 target_link_libraries(decn_test decn_cover coverage_config Catch)
 
-add_executable(decn_tests catch_main.cpp decn_tests.cpp ../utils.c)
+add_executable(decn_tests catch_main.cpp decn_tests.cpp decn_tests_trig.cpp ../utils.c)
 target_link_libraries(decn_tests decn_cover coverage_config mpfr Catch)
 
 # decn prototyping

src/decn/decn.c

@@ -84,6 +84,15 @@ const dec80 DECN_LN_10 = {
 	0, {23,  2, 58, 50, 92, 99, 40, 45, 68}
 };
 
+const dec80 DECN_2PI = {
+	0, {62, 83, 18, 53, 7, 17, 95, 86, 48}
+};
+
+// 180/pi = 1rad in degree
+const dec80 DECN_1RAD = {
+	1, {57, 29, 57, 79, 51, 30, 82, 32,  9}
+};
+
 
 void copy_decn(dec80* const dest, const dec80* const src){
 	uint8_t i;
@@ -1254,6 +1263,82 @@ void pow_decn(void) {
 	exp_decn();
 }
 
+void sincos_decn(void) {
+	#define SIN Tmp2Decn
+	#define COS Tmp3Decn
+	#define STP Tmp4Decn
+
+	remove_leading_zeros(&AccDecn);
+
+	// TODO: implement scaling to 0..2pi
+	copy_decn(&BDecn, &DECN_2PI);
+	if (compare_magn() == 1) {
+		set_dec80_NaN(&AccDecn);
+		set_dec80_NaN(&BDecn);
+		return;
+	}
+	set_dec80_zero(&BDecn);
+	if (compare_magn() == -1) {
+		set_dec80_NaN(&AccDecn);
+		set_dec80_NaN(&BDecn);
+		return;
+	}
+
+	copy_decn(&STP, &AccDecn);
+	copy_decn(&COS, &DECN_1);
+	set_dec80_zero(&SIN);
+	// 0.0 00 05
+	SIN.lsu[2] = 5;
+	negate_decn(&SIN);
+	while (STP.exponent >= 0) {
+		// COS = COS - SIN / 10000
+		copy_decn(&AccDecn, &COS);
+		copy_decn(&BDecn, &SIN);
+		shift_right(&BDecn);
+		shift_right(&BDecn);
+		shift_right(&BDecn);
+		shift_right(&BDecn);
+		negate_decn(&BDecn);
+		add_decn();
+		copy_decn(&COS, &AccDecn);
+		// SIN = SIN + COS / 10000
+		copy_decn(&AccDecn, &SIN);
+		copy_decn(&BDecn, &COS);
+		shift_right(&BDecn);
+		shift_right(&BDecn);
+		shift_right(&BDecn);
+		shift_right(&BDecn);
+		add_decn();
+		copy_decn(&SIN, &AccDecn);
+		// STP = STP - 0.0 00 1
+		copy_decn(&AccDecn, &STP);
+		set_dec80_zero(&BDecn);
+		BDecn.lsu[2] = 10;
+		negate_decn(&BDecn);
+		add_decn();
+		copy_decn(&STP, &AccDecn);
+	}
+}
+
+void sin_decn(void) {
+	sincos_decn();
+	copy_decn(&AccDecn, &SIN);
+}
+
+void cos_decn(void) {
+	sincos_decn();
+	copy_decn(&AccDecn, &COS);
+}
+
+void tan_decn(void) {
+	sincos_decn();
+	copy_decn(&AccDecn, &SIN);
+	copy_decn(&BDecn, &COS);
+	div_decn();
+}
+#undef SIN
+#undef COS
+
 
 static void set_str_error(void){
 	Buf[0] = 'E';

src/decn/decn.h

@@ -86,6 +86,10 @@ void exp_decn(void);
 void exp10_decn(void);
 void pow_decn(void);
 
+void sin_decn(void);
+void cos_decn(void);
+void tan_decn(void);
+
 //Buf should hold at least 18 + 4 + 5 + 1 = 28
 #define DECN_BUF_SIZE 28
 extern __xdata char Buf[DECN_BUF_SIZE];

src/decn/decn_tests_trig.cpp

@@ -0,0 +1,138 @@
+
+#include <string>
+#include <boost/multiprecision/mpfr.hpp>
+#include <catch.hpp>
+#include "decn.h"
+namespace bmp = boost::multiprecision;
+using Catch::Matchers::Equals;
+
+static void sin_test(
+	const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3)
+{
+	CAPTURE(a_str); CAPTURE(a_exp);
+	build_dec80(a_str, a_exp);
+	sin_decn();
+	decn_to_str_complete(&AccDecn);
+	CAPTURE(Buf);  // acc / b
+
+	bmp::mpfr_float::default_precision(50);
+	std::string a_full_str(a_str);
+	a_full_str += "e" + std::to_string(a_exp);
+
+	bmp::mpfr_float a_actual(a_full_str);
+	a_actual = sin(a_actual);
+	CAPTURE(a_actual);  // acc / b
+
+	bmp::mpfr_float calculated(Buf);
+	if (rtol >= 0) {
+		bmp::mpfr_float rel_diff = abs((a_actual - calculated) / a_actual);
+		CHECK(rel_diff < rtol);
+	} else {
+		bmp::mpfr_float diff = abs(a_actual - calculated);
+		CHECK(diff < atol);
+	}
+}
+
+static void cos_test(
+	const char* a_str, int a_exp, double rtol=1e-2)
+{
+	CAPTURE(a_str); CAPTURE(a_exp);
+	build_dec80(a_str, a_exp);
+	cos_decn();
+	decn_to_str_complete(&AccDecn);
+	CAPTURE(Buf);  // acc / b
+
+	bmp::mpfr_float::default_precision(50);
+	std::string a_full_str(a_str);
+	a_full_str += "e" + std::to_string(a_exp);
+
+	bmp::mpfr_float a_actual(a_full_str);
+	a_actual = cos(a_actual);
+	CAPTURE(a_actual);  // acc / b
+
+	bmp::mpfr_float calculated(Buf);
+	bmp::mpfr_float rel_diff = abs((a_actual - calculated) / a_actual);
+
+	CHECK(rel_diff < rtol);
+}
+
+static void tan_test(
+	const char* a_str, int a_exp, double rtol=1e-2)
+{
+	CAPTURE(a_str); CAPTURE(a_exp);
+	build_dec80(a_str, a_exp);
+	tan_decn();
+	decn_to_str_complete(&AccDecn);
+	CAPTURE(Buf);  // acc / b
+
+	bmp::mpfr_float::default_precision(50);
+	std::string a_full_str(a_str);
+	a_full_str += "e" + std::to_string(a_exp);
+
+	bmp::mpfr_float a_actual(a_full_str);
+	a_actual = tan(a_actual);
+	CAPTURE(a_actual);  // acc / b
+
+	bmp::mpfr_float calculated(Buf);
+	bmp::mpfr_float rel_diff = abs((a_actual - calculated) / a_actual);
+
+	CHECK(rel_diff < rtol);
+}
+
+TEST_CASE("sin") {
+	sin_test("0.1", 0);
+	sin_test("0.05", 0);
+	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.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);
+}
+
+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);
+}
+
+
+TEST_CASE("tan") {
+	tan_test("0.1", 0);
+	tan_test("0.05", 0);
+	tan_test("0.01", 0);
+	tan_test("0.001", 0);
+	tan_test("0.0001", 0);
+	tan_test("0.00001", 0);
+	tan_test("0.000001", 0, -1);
+	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);
+}