From f33e30a8a45d4c1d1cb457f01cafe5dd75e9f483 Mon Sep 17 00:00:00 2001
From: Ben Lynn <benlynn@gmail.com>
Date: Tue, 16 Sep 2008 03:22:28 -0700
Subject: [PATCH] Quadratic algorithm supports sign, as does decimal
 conversion.

Better bihom_test.
---
 bihom.c       | 24 ++++++++++++++++++++++++
 bihom_test.c  | 19 +++++++++++++++++--
 cf.c          |  4 ++++
 cf.h          |  1 +
 famous_test.c | 10 +++++-----
 mobius.c      |  1 +
 test.c        | 19 ++++++++++++++++++-
 7 files changed, 70 insertions(+), 8 deletions(-)

diff --git a/bihom.c b/bihom.c
index 652040d..96c8942 100644
--- a/bihom.c
+++ b/bihom.c
@@ -90,6 +90,30 @@ static void *bihom(cf_t cf) {
     mpz_set(p->r0, p->p0);  mpz_set(p->r1, p->p1);
     mpz_set(p->p0, t0);     mpz_set(p->p1, t1);
   }
+  // Determine sign.
+  while (mpz_sgn(p->p1) != mpz_sgn(p->q1)
+      || mpz_sgn(p->q1) != mpz_sgn(p->r1)
+      || mpz_sgn(p->r1) != mpz_sgn(p->s1)
+      || mpz_sgn(p->p0) != mpz_sgn(p->q0)
+      || mpz_sgn(p->q0) != mpz_sgn(p->r0)
+      || mpz_sgn(p->r0) != mpz_sgn(p->s0)) {
+    move_right();
+    move_down();
+  }
+  if (mpz_sgn(p->p0) < 0) {
+    mpz_neg(p->p0, p->p0);
+    mpz_neg(p->q0, p->q0);
+    mpz_neg(p->r0, p->r0);
+    mpz_neg(p->s0, p->s0);
+    cf_flip_sign(cf);
+  }
+  if (mpz_sgn(p->p1) < 0) {
+    mpz_neg(p->p1, p->p1);
+    mpz_neg(p->q1, p->q1);
+    mpz_neg(p->r1, p->r1);
+    mpz_neg(p->s1, p->s1);
+    cf_flip_sign(cf);
+  }
 
   int recur() {
     if (!mpz_sgn(p->s1)) {
diff --git a/bihom_test.c b/bihom_test.c
index 874beee..10548f1 100644
--- a/bihom_test.c
+++ b/bihom_test.c
@@ -22,7 +22,7 @@ int main() {
 
   cf_t b = cf_new_bihom(e, pi, addarray);
 
-  CF_EXPECT_DEC(b, "58598744820488384738");
+  CF_EXPECT_DEC(b, "5.8598744820488384738");
 
   cf_free(b);
   cf_free(e);
@@ -38,11 +38,26 @@ int main() {
   mpz_set_si(addarray[7], 1);
   b = cf_new_bihom(s1, c1, addarray);
 
-  CF_EXPECT_DEC(b, "09092974268256816953");
+  CF_EXPECT_DEC(b, "0.9092974268256816953");
 
   cf_free(b);
   cf_free(c1);
   cf_free(s1);
+
+  // Check 2 (cos 1)^2 - 1 = cos 2 =
+  s1 = cf_new_cos1();  // TODO: Implement tee, use that instead.
+  c1 = cf_new_cos1();
+  mpz_set_si(addarray[0], 2);
+  mpz_set_si(addarray[3], -1);
+  mpz_set_si(addarray[7], 1);
+  b = cf_new_bihom(s1, c1, addarray);
+
+  CF_EXPECT_DEC(b, "-0.41614683654714238699");
+
+  cf_free(b);
+  cf_free(c1);
+  cf_free(s1);
+
   for (int i = 0; i < 8; i++) {
     mpz_clear(addarray[i]);
   }
diff --git a/cf.c b/cf.c
index 2bad1b5..780f2c5 100644
--- a/cf.c
+++ b/cf.c
@@ -44,6 +44,10 @@ void *cf_data(cf_t cf) {
   return cf->data;
 }
 
+void cf_set_sign(cf_t cf, int sign) {
+  cf->sign = sign;
+}
+
 int cf_sign(cf_t cf) {
   return cf->sign;
 }
diff --git a/cf.h b/cf.h
index 1fd64e3..c319414 100644
--- a/cf.h
+++ b/cf.h
@@ -14,6 +14,7 @@ static inline cf_t cf_new_const(void *(*func)(cf_t)) {
 }
 void cf_free(cf_t cf);
 
+void cf_set_sign(cf_t cf, int sign);
 int cf_sign(cf_t cf);
 int cf_flip_sign(cf_t cf);
 void cf_get(mpz_t z, cf_t cf);
diff --git a/famous_test.c b/famous_test.c
index ec2ac79..10db391 100644
--- a/famous_test.c
+++ b/famous_test.c
@@ -9,11 +9,11 @@
 #include "test.h"
 
 int main() {
-  CF_NEW_EXPECT_DEC(cf_new_e, "2718281828");
-  CF_NEW_EXPECT_DEC(cf_new_pi, "31415926535897932384");
-  CF_NEW_EXPECT_DEC(cf_new_tan1, "15574077246549022305");
+  CF_NEW_EXPECT_DEC(cf_new_e, "2.718281828");
+  CF_NEW_EXPECT_DEC(cf_new_pi, "3.1415926535897932384");
+  CF_NEW_EXPECT_DEC(cf_new_tan1, "1.5574077246549022305");
 
-  CF_NEW_EXPECT_DEC(cf_new_sin1, "08414709848078965066");
-  CF_NEW_EXPECT_DEC(cf_new_cos1, "05403023058681397174");
+  CF_NEW_EXPECT_DEC(cf_new_sin1, "0.8414709848078965066");
+  CF_NEW_EXPECT_DEC(cf_new_cos1, "0.5403023058681397174");
   return 0;
 }
diff --git a/mobius.c b/mobius.c
index 23b8373..8507449 100644
--- a/mobius.c
+++ b/mobius.c
@@ -242,6 +242,7 @@ static void *mobius_decimal(cf_t cf) {
   mpz_t t0, t1, t2; mpz_init(t2); mpz_init(t1); mpz_init(t0);
 
   // Determine the sign.
+  cf_set_sign(cf, cf_sign(input));
   while (mpz_sgn(pq->pold) != mpz_sgn(pq->p)
       || mpz_sgn(pq->qold) != mpz_sgn(pq->q)) {
     cf_get(denom, input);
diff --git a/test.c b/test.c
index 075c714..ec4d43c 100644
--- a/test.c
+++ b/test.c
@@ -10,11 +10,28 @@ void cf_expect_dec(cf_t x, char *result, char *filename, int line) {
   cf_t conv;
   conv = cf_new_cf_to_decimal(x);
   int len = strlen(result);
+
   char s[len + 1];
-  for (int i = 0; i <= len; i++) {
+  cf_get(z, conv);
+  int i = 0;
+  if (cf_sign(conv) < 0) {
+    s[i] = '-';
+    i++;
+  }
+  if (i >= len) goto testit;
+  i += gmp_snprintf(s + i, len - i + 1, "%Zd", z);
+  if (i >= len) goto testit;
+  s[i] = '.';
+  i++;
+  if (i >= len) goto testit;
+
+  while (i < len) {
     cf_get(z, conv);
     s[i] = mpz_get_ui(z) + '0';
+    i++;
   }
+
+testit:
   s[len] = 0;
   if (strcmp(s, result)) {
     fprintf(stderr, "\n%s:%d: bad continued fraction decimal expansion\n",
-- 
2.11.4.GIT