[Gmp-commit] /var/hg/gmp: mpn/generic/rootrem.c: Use sqrtrem when k=2, it's f...

[mercurial at gmplib.org]

details:   /var/hg/gmp/rev/6659bc3f7541
changeset: 16829:6659bc3f7541
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Tue Sep 15 08:56:04 2015 +0200
description:
mpn/generic/rootrem.c: Use sqrtrem when k=2, it's faster.

diffstat:

 mpn/generic/rootrem.c |  9 +++++++--
 1 files changed, 7 insertions(+), 2 deletions(-)

diffs (40 lines):

diff -r 7b920a186f15 -r 6659bc3f7541 mpn/generic/rootrem.c
--- a/mpn/generic/rootrem.c	Sun Sep 13 09:28:30 2015 +0200
+++ b/mpn/generic/rootrem.c	Tue Sep 15 08:56:04 2015 +0200
@@ -3,6 +3,7 @@
 
    Contributed by Paul Zimmermann (algorithm) and
    Paul Zimmermann and Torbjorn Granlund (implementation).
+   Marco Bodrato wrote logbased_root to seed the loop. 
 
    THE FUNCTIONS IN THIS FILE ARE INTERNAL, AND HAVE MUTABLE INTERFACES.  IT'S
    ONLY SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT'S ALMOST
@@ -93,6 +94,8 @@
   ASSERT (up[un - 1] != 0);
   ASSERT (k > 1);
 
+  if (UNLIKELY (k == 2))
+    return mpn_sqrtrem (rootp, remp, up, un);
   /* (un-1)/k > 2 <=> un > 3k <=> (un + 2)/3 > k */
   if (remp == NULL && (un + 2) / 3 > k)
     /* Pad {up,un} with k zero limbs.  This will produce an approximate root
@@ -252,8 +255,8 @@
 
   kk = k * xnb;		/* number of truncated bits in the input */
 
-  /* FIXME: Should we skipp all loops when xnb > snb ? */
-  for (uh = k - 1, logk = 2; (uh >>= 1) != 0; ++logk )
+  /* FIXME: Should we skip the next two loops when xnb <= snb ? */
+  for (uh = (k - 1) / 2, logk = 3; (uh >>= 1) != 0; ++logk )
     ;
   /* logk = ceil(log(k)/log(2)) + 1 */
 
@@ -292,6 +295,8 @@
      So the worst case is S'/S=3/2, thus S'^k <= (3/2)^k * S^k. Since S^k
      fits in un limbs, the number of extra limbs needed is bounded by
      ceil(k*log2(3/2)/GMP_NUMB_BITS). */
+  /* THINK: with the use of logbased_root, maybe the constant is
+     258/256 instead of 3/2 ? log2(258/256) < 1/89 < 1/64 */
 #define EXTRA 2 + (mp_size_t) (0.585 * (double) k / (double) GMP_NUMB_BITS)
   TMP_ALLOC_LIMBS_3 (scratch, un + 1, /* used by mpn_div_q */
 		     qp, un + EXTRA,  /* will contain quotient and remainder