From e4ad37623b2e5fdfccd76c15d892a498c6dc37da Mon Sep 17 00:00:00 2001
From: Mirko Scholz <srtlg@users.noreply.github.com>
Date: Thu, 10 Sep 2020 15:38:26 +0200
Subject: [PATCH] implemented arctan

---
 src/decn/decn.c              | 56 +++++++++++++++++++++++++-----------
 src/decn/decn.h              |  9 +++++-
 src/decn/decn_tests_trig.cpp | 22 ++++++++++++++
 3 files changed, 70 insertions(+), 17 deletions(-)

diff --git a/src/decn/decn.c b/src/decn/decn.c
index 1c1faf0..5322ea5 100644
--- a/src/decn/decn.c
+++ b/src/decn/decn.c
@@ -1297,19 +1297,39 @@ void project_decn_into_0_2pi(void) {
 	}
 }
 
+// K. Shirriff, "Reversing Sinclair's amazing 1974 calculator hack - half the ROM of the HP-35"
+// http://files.righto.com/calculator/sinclair_scientific_simulator.html
 #define SIN Tmp2Decn
 #define COS Tmp3Decn
 #define THETA Tmp4Decn
-void sincos_decn(void) {
-	project_decn_into_0_2pi();
-
-	copy_decn(&THETA, &AccDecn);
-	copy_decn(&COS, &DECN_1);
-	set_dec80_zero(&SIN);
-	// 0.0 00 05
-	SIN.lsu[2] = 5;
-	negate_decn(&SIN);
-	while (THETA.exponent >= 0) {
+void sincos_decn(const uint8_t sincos_arctan) {
+	const uint8_t is_negative = AccDecn.exponent < 0;
+	if (sincos_arctan) {
+		set_dec80_zero(&THETA);
+		if (is_negative) negate_decn(&AccDecn);
+		copy_decn(&COS, &AccDecn);
+		copy_decn(&SIN, &DECN_1);
+	} else {
+		project_decn_into_0_2pi();
+		copy_decn(&THETA, &AccDecn);
+		copy_decn(&COS, &DECN_1);
+		set_dec80_zero(&SIN);
+		// 0.0 00 05
+		SIN.lsu[2] = 5;
+		negate_decn(&SIN);
+	}
+	do {
+		if (sincos_arctan) {
+			// THETA is in AccDecn from previous iteration
+			if (COS.exponent < 0) {
+				if (is_negative) negate_decn(&AccDecn);
+				break;
+			}
+		} else {
+			if (THETA.exponent < 0) {
+				break;
+			}
+		}
 		// COS = COS - SIN / 10000
 		copy_decn(&AccDecn, &COS);
 		copy_decn(&BDecn, &SIN);
@@ -1329,32 +1349,36 @@ void sincos_decn(void) {
 		shift_right(&BDecn);
 		add_decn();
 		copy_decn(&SIN, &AccDecn);
-		// THETA = THETA - 0.0 00 1
+		// THETA = THETA -/+ 0.0 00 1
 		copy_decn(&AccDecn, &THETA);
 		set_dec80_zero(&BDecn);
 		BDecn.lsu[2] = 10;
-		negate_decn(&BDecn);
+		if (!sincos_arctan) negate_decn(&BDecn);
 		add_decn();
 		copy_decn(&THETA, &AccDecn);
-	}
+	} while (1);
 }
 
 void sin_decn(void) {
-	sincos_decn();
+	sincos_decn(0);
 	copy_decn(&AccDecn, &SIN);
 }
 
 void cos_decn(void) {
-	sincos_decn();
+	sincos_decn(0);
 	copy_decn(&AccDecn, &COS);
 }
 
 void tan_decn(void) {
-	sincos_decn();
+	sincos_decn(0);
 	copy_decn(&AccDecn, &SIN);
 	copy_decn(&BDecn, &COS);
 	div_decn();
 }
+
+void arctan_decn(void) {
+	sincos_decn(1);
+}
 #undef SIN
 #undef COS
 #undef THETA
diff --git a/src/decn/decn.h b/src/decn/decn.h
index 9f45e86..af9803f 100644
--- a/src/decn/decn.h
+++ b/src/decn/decn.h
@@ -86,10 +86,10 @@ void exp_decn(void);
 void exp10_decn(void);
 void pow_decn(void);
 
-void project_decn_into_0_2pi(void);
 void sin_decn(void);
 void cos_decn(void);
 void tan_decn(void);
+void arctan_decn(void);
 void to_degree_decn(void);
 void to_radian_decn(void);
 
@@ -110,6 +110,13 @@ void decn_to_str_complete(const dec80* x);
 void build_decn_at(dec80* dest, const char* signif_str, exp_t exponent);
 #endif
 
+#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);
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/decn/decn_tests_trig.cpp b/src/decn/decn_tests_trig.cpp
index 91b11fa..ceba7a7 100644
--- a/src/decn/decn_tests_trig.cpp
+++ b/src/decn/decn_tests_trig.cpp
@@ -53,6 +53,12 @@ static void tan_test(
 	trig_test(tan_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return tan(x);}, a_str, a_exp, rtol, atol);
 }
 
+static void atan_test(
+	const char* a_str, int a_exp, double rtol=5e-3, double atol=1e-3)
+{
+	trig_test(arctan_decn, [](bmp::mpfr_float x) -> bmp::mpfr_float {return atan(x);}, a_str, a_exp, rtol, atol);
+}
+
 
 const char * const pi = "3.141592653589793239";
 const char * const pi_threequarters = "2.356194490192344929";
@@ -145,3 +151,19 @@ TEST_CASE("tan") {
 	tan_test("2.5", 0);
 	tan_test("3.0", 0);
 }
+
+TEST_CASE("arctan") {
+	atan_test("0.001", 0);
+	atan_test("-0.001", 0);
+	atan_test("0.7", 0);
+	atan_test("-0.7", 0);
+	atan_test("0.1", 0);
+	atan_test("-0.1", 0);
+	atan_test("1.0", 0);
+	atan_test("-1.0", 0);
+	atan_test("2.0", 0);
+	atan_test("-2.0", 0);
+	atan_test("3.0", 0);
+	atan_test("-3.0", 0);
+	atan_test("0", 0, -1);
+}