factor: switch from mp to single when doable

This significantly improves performance when a number exceeds
2**(W_TYPE_SIZE - 1) and is the product of a prime less than
FIRST_OMITTED_PRIME and another prime less than 2**(W_TYPE_SIZE - 1).
On my platform, for example, it doubled the speed of factoring
4999 * (2**128 - 159).
* src/factor.c (mp_size, mp_finish_in_single): New functions.
(mp_factor_using_division, mp_factor_using_pollard_rho):
Finish using single precision when possible.
* tests/factor/factor.pl (lt-5000-times-128-bit): New test.

diff --git a/src/factor.c b/src/factor.c
index e32cf7c3af47d5bf7ce7791f2599796e15ec829c..17af10745caf6a119ad7cf472f79ca968f15d5b3 100644 (file)
--- a/src/factor.c
+++ b/src/factor.c
@@ -838,6 +838,48 @@ factor_using_division (mp_limb_t t1, mp_limb_t t0,
   return make_uuint (t1, t0);
 }
 
+/* Return the number of limbs in positive N.  */
+static mp_size_t
+mp_size (mpz_t n)
+{
+  /* Tell the compiler that N is positive; this can speed up access to N.  */
+  assume (0 < mpz_sgn (n));
+
+  return mpz_size (n);
+}
+
+/* If N is small enough to be factorable by 'factor',
+   add its factors to MP_FACTORS and return true.
+   Otherwise, return false.  */
+static bool
+mp_finish_in_single (mpz_t n, struct mp_factors *mp_factors)
+{
+  if (2 < mp_size (n))
+    return false;
+  mp_limb_t n1 = mpz_getlimbn (n, 1);
+  if (n1 >> (W_TYPE_SIZE - 1))
+    return false;
+  mp_limb_t n0 = mpz_getlimbn (n, 0);
+  mpz_set_ui (n, 1);
+
+  struct factors factors;
+  factor (n1, n0, &factors);
+
+  if (hi_is_set (&factors.plarge))
+    {
+      mpz_t p = MPZ_ROINIT_N (factors.plarge.uu, 2);
+      mp_factor_insert (mp_factors, p, 1);
+    }
+
+  for (int i = factors.nfactors; 0 < i; i--)
+    {
+      mpz_t p = MPZ_ROINIT_N (&factors.p[i - 1], 1);
+      mp_factor_insert (mp_factors, p, factors.e[i - 1]);
+    }
+
+  return true;
+}
+
 static void
 mp_factor_using_division (mpz_t t, struct mp_factors *factors)
 {
@@ -848,13 +890,22 @@ mp_factor_using_division (mpz_t t, struct mp_factors *factors)
     {
       mpz_fdiv_q_2exp (t, t, m);
       mp_factor_insert_ui (factors, 2, m);
+      if (mp_finish_in_single (t, factors))
+        return;
     }
 
   unsigned long int d = 3;
   for (idx_t i = 1; i <= PRIMES_PTAB_ENTRIES;)
     {
       for (m = 0; mpz_divisible_ui_p (t, d); m++)
-        mpz_tdiv_q_ui (t, t, d);
+        {
+          mpz_tdiv_q_ui (t, t, d);
+          if (mp_finish_in_single (t, factors))
+            {
+              mp_factor_insert_ui (factors, d, m + 1);
+              return;
+            }
+        }
       if (m)
         mp_factor_insert_ui (factors, d, m);
       d += primes_diff[i++];
@@ -1725,16 +1776,20 @@ mp_factor_using_pollard_rho (mpz_t n, unsigned long int a,
 
       mpz_divexact (n, n, t);   /* divide by t, before t is overwritten */
 
-      if (!mp_prime_p (t))
-        {
-          devmsg ("[composite factor--restarting pollard-rho] ");
-          mp_factor_using_pollard_rho (t, a + 1, factors);
-        }
-      else
+      if (!mp_finish_in_single (t, factors))
         {
-          mp_factor_insert (factors, t, 1);
+          if (mp_prime_p (t))
+            mp_factor_insert (factors, t, 1);
+          else
+            {
+              devmsg ("[composite factor--restarting pollard-rho] ");
+              mp_factor_using_pollard_rho (t, a + 1, factors);
+            }
         }
 
+      if (mp_finish_in_single (n, factors))
+        break;
+
       if (mp_prime_p (n))
         {
           mp_factor_insert (factors, n, 1);

diff --git a/tests/factor/factor.pl b/tests/factor/factor.pl
index c45a8fe4d1744908b4eeff51dd70a9f42e04753e..782deefe6c0d75360d20ad795b9ac333a68eb045 100755 (executable)
--- a/tests/factor/factor.pl
+++ b/tests/factor/factor.pl
@@ -92,6 +92,9 @@ my @Tests =
      ['bug-with-128-bit-uintmax_t',
       '340282366920938463463374607431768211355',
       {OUT => '5 31 2195370109167344925570158757624311041'}],
+     ['lt-5000-times-128-bit',
+      '1701071552237771378853409662551409288273703',
+      {OUT => '4999 340282366920938463463374607431768211297'}],
      ['h-1', '-h 3000', {OUT => '2^3 3 5^3'}],
      ['h-2', '3000 --exponents', {OUT => '2^3 3 5^3'}],
     );