From 5a76cbf4393e241d1b6484d00c4d5e8ccc104c13 Mon Sep 17 00:00:00 2001
From: Marco Bodrato <bodrato@mail.dm.unipi.it>
Date: Mon, 9 Mar 2020 17:01:17 +0100
Subject: mpz/bin_ui.c (mpz_bin_ui): Siplify special cases.

---
 mpz/bin_ui.c | 58 +++++++++++++++++++++++++++++++---------------------------
 1 file changed, 31 insertions(+), 27 deletions(-)

(limited to 'mpz')

diff --git a/mpz/bin_ui.c b/mpz/bin_ui.c
index 143e2988c..e40af0a16 100644
--- a/mpz/bin_ui.c
+++ b/mpz/bin_ui.c
@@ -69,7 +69,7 @@ posmpz_inc_ui (mpz_ptr r, unsigned long in)
 #else
   ASSERT (SIZ (r) > 0);
   MPN_INCR_U (PTR (r), SIZ (r) + 1, in);
-  SIZ (r) += (PTR (r)[SIZ (r)] != 0);
+  SIZ (r) += PTR (r)[SIZ (r)];
 #endif
 }
 
@@ -379,36 +379,40 @@ mpz_bin_ui (mpz_ptr r, mpz_srcptr n, unsigned long int k)
 	  MPZ_NEWALLOC (r, 1)[0] = 1;
 	}
 #if APARTAJ_KALKULOJ > 2
-      else if (k == 2)
+      else if (k > 1)
 	{
-	  mpz_add_ui (ni, ni, 1);
-	  mpz_mul (r, ni, ni);
-	  mpz_add (r, r, ni);
-	  posmpz_rsh1 (r);
-	}
-#endif
-#if APARTAJ_KALKULOJ > 3
-      else if (k > 2)
-	{ /* k = 3, 4 */
-	  mpz_add_ui (ni, ni, 2); /* n+1 */
-	  mpz_mul (r, ni, ni); /* (n+1)^2 */
-	  mpz_sub_ui (r, r, 1); /* (n+1)^2-1 */
-	  if (k == 3)
+	  mpz_add_ui (ni, ni, 1 + (APARTAJ_KALKULOJ > 2 && k > 2));
+	  mpz_mul (r, ni, ni); /* r = (n + (k>2))^2 */
+	  if (APARTAJ_KALKULOJ == 2 || k == 2)
 	    {
-	      mpz_mul (r, r, ni); /* ((n+1)^2-1)(n+1) = n(n+1)(n+2) */
-	      /* mpz_divexact_ui (r, r, 6); /\* 6=3<<1; div_by3 ? *\/ */
-	      mpn_pi1_bdiv_q_1 (PTR(r), PTR(r), SIZ(r), 3, GMP_NUMB_MASK/3*2+1, 1);
-	      MPN_NORMALIZE_NOT_ZERO (PTR(r), SIZ(r));
+	      mpz_add (r, r, ni); /* n^2+n= n(n+1) */
+	      posmpz_rsh1 (r);
 	    }
-	  else /* k = 4 */
-	    {
-	      mpz_add (ni, ni, r); /* (n+1)^2+n */
-	      mpz_mul (r, ni, ni); /* ((n+1)^2+n)^2 */
-	      mpz_sub_ui (r, r, 1); /* ((n+1)^2+n)^2-1 = n(n+1)(n+2)(n+3) */
-	      /* mpz_divexact_ui (r, r, 24); /\* 24=3<<3; div_by3 ? *\/ */
-	      mpn_pi1_bdiv_q_1 (PTR(r), PTR(r), SIZ(r), 3, GMP_NUMB_MASK/3*2+1, 3);
-	      MPN_NORMALIZE_NOT_ZERO (PTR(r), SIZ(r));
+#if APARTAJ_KALKULOJ > 3
+#if APARTAJ_KALKULOJ != 5
+#error Not implemented! 3 < APARTAJ_KALKULOJ != 5
+#endif
+	  else /* k > 2 */
+	    { /* k = 3, 4 */
+	      mpz_sub_ui (r, r, 1); /* (n+1)^2-1 */
+	      if (k == 3)
+		{
+		  mpz_mul (r, r, ni); /* ((n+1)^2-1)(n+1) = n(n+1)(n+2) */
+		  /* mpz_divexact_ui (r, r, 6); /\* 6=3<<1; div_by3 ? *\/ */
+		}
+	      else /* k = 4 */
+		{
+		  mpz_add (ni, ni, r); /* (n+1)^2+n */
+		  mpz_mul (r, ni, ni); /* ((n+1)^2+n)^2 */
+		  /* We should subtract one: ((n+1)^2+n)^2-1 = n(n+1)(n+2)(n+3). */
+		  /* PTR (r) [0] ^= 1; would suffice, but it is not even needed, */
+		  /* because the next division will shift away this bit anyway.  */
+		  /* mpz_divexact_ui (r, r, 24); /\* 24=3<<3; div_by3 ? *\/ */
+		}
+	      mpn_pi1_bdiv_q_1 (PTR(r), PTR(r), SIZ(r), 3, GMP_NUMB_MASK/3*2+1, 1 | (k>>1));
+	      SIZ(r) -= PTR(r) [SIZ(r) - 1] == 0;
 	    }
+#endif
 	}
 #endif
       else
-- 
cgit v1.2.1