[gnumeric] Numtheory: speed things up.

[Morten Welinder <mortenw src gnome org>]

commit d4c8c5cfaf395b48695f0cec12d11cccde036ce2
Author: Morten Welinder <terra gnome org>
Date:   Mon Nov 27 10:35:43 2017 -0500

    Numtheory: speed things up.

 NEWS                             |    1 +
 plugins/fn-numtheory/ChangeLog   |    3 +-
 plugins/fn-numtheory/numtheory.c |   65 ++++++++++++++++++-------------------
 3 files changed, 34 insertions(+), 35 deletions(-)
---
diff --git a/NEWS b/NEWS
index 3f09403..5dba03b 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,7 @@ Morten:
        * Improve xlsx export of cell comments.  [#790756]
        * Plug leaks.
        * Fix potential crash on exit.
+       * Speed-up number theory functions.
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.36
diff --git a/plugins/fn-numtheory/ChangeLog b/plugins/fn-numtheory/ChangeLog
index 0f558d8..50c933c 100644
--- a/plugins/fn-numtheory/ChangeLog
+++ b/plugins/fn-numtheory/ChangeLog
@@ -1,7 +1,6 @@
 2017-11-27  Morten Welinder  <terra gnome org>
 
-       * numtheory.c (ithprime): For each new candidate we at most need
-       to consider one extra prime divisor.
+       * numtheory.c (ithprime): switch to explicit sieve.  Much faster.
 
 2017-11-23  Morten Welinder  <terra gnome org>
 
diff --git a/plugins/fn-numtheory/numtheory.c b/plugins/fn-numtheory/numtheory.c
index 2d8258b..b5dbd2a 100644
--- a/plugins/fn-numtheory/numtheory.c
+++ b/plugins/fn-numtheory/numtheory.c
@@ -64,57 +64,56 @@ intpow (int p, int v)
 #define ITHPRIME_LIMIT 10000000
 static guint *prime_table = NULL;
 
-/* Calculate the i-th prime.  Returns TRUE on error.  */
+// Bit-field macros for sieve.  Note that only odd indices are used.
+#define SIEVE_ITEM(u_) s[(u_) >> 4]
+#define SIEVE_BIT(u_) (1u << (((u_) >> 1) & 7))
+
+/* Calculate the i-th prime.  Returns TRUE on too-big-to-handle error.  */
 static gboolean
 ithprime (int i, guint64 *res)
 {
-       static guint computed = 0;
-       static guint allocated = 0;
-
        if (i < 1 || (guint)i > ITHPRIME_LIMIT)
                return TRUE;
 
-       if ((guint)i > computed) {
-               static guint candidate, jlim;
-
-               if ((guint)i > allocated) {
-                       allocated = MAX ((guint)i, 2 * allocated + 100);
-                       allocated = MIN (allocated, ITHPRIME_LIMIT);
-                       prime_table = g_renew (guint, prime_table, allocated);
-                       if (computed == 0) {
-                               prime_table[computed++] = 2;
-                               candidate = prime_table[computed++] = 3;
-                               jlim = 1;
-                               g_assert (candidate < prime_table[jlim] * prime_table[jlim]);
-                       }
-               }
+       if (!prime_table) {
+               // Compute an upper bound of the largest prime we need.
+               // See https://en.wikipedia.org/wiki/Prime_number_theorem
+               guint ub = (guint)
+                       (ITHPRIME_LIMIT *
+                        (log (ITHPRIME_LIMIT) + log (log (ITHPRIME_LIMIT))));
+               guint8 *s = g_new0 (guint8, (ub >> 4) + 1);
+               guint N = 0, c;
+               guint L = (int)sqrt (ub);
 
-               while ((guint)i > computed) {
-                       gboolean prime = TRUE;
-                       guint j;
+               prime_table = g_new (guint, ITHPRIME_LIMIT);
+               prime_table[N++] = 2;
 
-                       candidate += 2;  /* Skip even candidates.  */
+               for (c = 3; N < ITHPRIME_LIMIT; c += 2) {
+                       guint d;
 
-                       // We at most need to consider one extra candidate
-                       if (candidate >= prime_table[jlim] * prime_table[jlim])
-                               jlim++;
+                       if (SIEVE_ITEM (c) & SIEVE_BIT (c))
+                               continue;
 
-                       for (j = 1; j < jlim; j++) {
-                               if (candidate % prime_table[j] == 0) {
-                                       prime = FALSE;
-                                       break;
-                               }
-                       }
+                       prime_table[N++] = c;
+                       if (c > L)
+                               continue; // Prevents square overflow
 
-                       if (prime)
-                               prime_table[computed++] = candidate;
+                       // Tag all odd multiples starting with c^2
+                       for (d = c * c; d <= ub; d += 2 * c)
+                               SIEVE_ITEM (d) |= SIEVE_BIT (d);
                }
+
+               g_free (s);
        }
 
        *res = prime_table[i - 1];
        return FALSE;
 }
 
+#undef SIEVE_ITEM
+#undef SIEVE_BIT
+
+
 /*
  * A function useful for computing multiplicative arithmetic functions.
  * Returns TRUE on error.