summaryrefslogtreecommitdiffstats
path: root/security/nss/lib/freebl/mpi/utils/sieve.c
diff options
context:
space:
mode:
authorMatt A. Tobin <mattatobin@localhost.localdomain>2018-02-02 04:16:08 -0500
committerMatt A. Tobin <mattatobin@localhost.localdomain>2018-02-02 04:16:08 -0500
commit5f8de423f190bbb79a62f804151bc24824fa32d8 (patch)
tree10027f336435511475e392454359edea8e25895d /security/nss/lib/freebl/mpi/utils/sieve.c
parent49ee0794b5d912db1f95dce6eb52d781dc210db5 (diff)
downloadUXP-5f8de423f190bbb79a62f804151bc24824fa32d8.tar
UXP-5f8de423f190bbb79a62f804151bc24824fa32d8.tar.gz
UXP-5f8de423f190bbb79a62f804151bc24824fa32d8.tar.lz
UXP-5f8de423f190bbb79a62f804151bc24824fa32d8.tar.xz
UXP-5f8de423f190bbb79a62f804151bc24824fa32d8.zip
Add m-esr52 at 52.6.0
Diffstat (limited to 'security/nss/lib/freebl/mpi/utils/sieve.c')
-rw-r--r--security/nss/lib/freebl/mpi/utils/sieve.c243
1 files changed, 243 insertions, 0 deletions
diff --git a/security/nss/lib/freebl/mpi/utils/sieve.c b/security/nss/lib/freebl/mpi/utils/sieve.c
new file mode 100644
index 000000000..57768af9e
--- /dev/null
+++ b/security/nss/lib/freebl/mpi/utils/sieve.c
@@ -0,0 +1,243 @@
+/*
+ * sieve.c
+ *
+ * Finds prime numbers using the Sieve of Eratosthenes
+ *
+ * This implementation uses a bitmap to represent all odd integers in a
+ * given range. We iterate over this bitmap, crossing off the
+ * multiples of each prime we find. At the end, all the remaining set
+ * bits correspond to prime integers.
+ *
+ * Here, we make two passes -- once we have generated a sieve-ful of
+ * primes, we copy them out, reset the sieve using the highest
+ * generated prime from the first pass as a base. Then we cross out
+ * all the multiples of all the primes we found the first time through,
+ * and re-sieve. In this way, we get double use of the memory we
+ * allocated for the sieve the first time though. Since we also
+ * implicitly ignore multiples of 2, this amounts to 4 times the
+ * values.
+ *
+ * This could (and probably will) be generalized to re-use the sieve a
+ * few more times.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+
+typedef unsigned char byte;
+
+typedef struct {
+ int size;
+ byte *bits;
+ long base;
+ int next;
+ int nbits;
+} sieve;
+
+void sieve_init(sieve *sp, long base, int nbits);
+void sieve_grow(sieve *sp, int nbits);
+long sieve_next(sieve *sp);
+void sieve_reset(sieve *sp, long base);
+void sieve_cross(sieve *sp, long val);
+void sieve_clear(sieve *sp);
+
+#define S_ISSET(S, B) (((S)->bits[(B) / CHAR_BIT] >> ((B) % CHAR_BIT)) & 1)
+#define S_SET(S, B) ((S)->bits[(B) / CHAR_BIT] |= (1 << ((B) % CHAR_BIT)))
+#define S_CLR(S, B) ((S)->bits[(B) / CHAR_BIT] &= ~(1 << ((B) % CHAR_BIT)))
+#define S_VAL(S, B) ((S)->base + (2 * (B)))
+#define S_BIT(S, V) (((V) - ((S)->base)) / 2)
+
+int
+main(int argc, char *argv[])
+{
+ sieve s;
+ long pr, *p;
+ int c, ix, cur = 0;
+
+ if (argc < 2) {
+ fprintf(stderr, "Usage: %s <width>\n", argv[0]);
+ return 1;
+ }
+
+ c = atoi(argv[1]);
+ if (c < 0)
+ c = -c;
+
+ fprintf(stderr, "%s: sieving to %d positions\n", argv[0], c);
+
+ sieve_init(&s, 3, c);
+
+ c = 0;
+ while ((pr = sieve_next(&s)) > 0) {
+ ++c;
+ }
+
+ p = calloc(c, sizeof(long));
+ if (!p) {
+ fprintf(stderr, "%s: out of memory after first half\n", argv[0]);
+ sieve_clear(&s);
+ exit(1);
+ }
+
+ fprintf(stderr, "%s: half done ... \n", argv[0]);
+
+ for (ix = 0; ix < s.nbits; ix++) {
+ if (S_ISSET(&s, ix)) {
+ p[cur] = S_VAL(&s, ix);
+ printf("%ld\n", p[cur]);
+ ++cur;
+ }
+ }
+
+ sieve_reset(&s, p[cur - 1]);
+ fprintf(stderr, "%s: crossing off %d found primes ... \n", argv[0], cur);
+ for (ix = 0; ix < cur; ix++) {
+ sieve_cross(&s, p[ix]);
+ if (!(ix % 1000))
+ fputc('.', stderr);
+ }
+ fputc('\n', stderr);
+
+ free(p);
+
+ fprintf(stderr, "%s: sieving again from %ld ... \n", argv[0], p[cur - 1]);
+ c = 0;
+ while ((pr = sieve_next(&s)) > 0) {
+ ++c;
+ }
+
+ fprintf(stderr, "%s: done!\n", argv[0]);
+ for (ix = 0; ix < s.nbits; ix++) {
+ if (S_ISSET(&s, ix)) {
+ printf("%ld\n", S_VAL(&s, ix));
+ }
+ }
+
+ sieve_clear(&s);
+
+ return 0;
+}
+
+void
+sieve_init(sieve *sp, long base, int nbits)
+{
+ sp->size = (nbits / CHAR_BIT);
+
+ if (nbits % CHAR_BIT)
+ ++sp->size;
+
+ sp->bits = calloc(sp->size, sizeof(byte));
+ memset(sp->bits, UCHAR_MAX, sp->size);
+ if (!(base & 1))
+ ++base;
+ sp->base = base;
+
+ sp->next = 0;
+ sp->nbits = sp->size * CHAR_BIT;
+}
+
+void
+sieve_grow(sieve *sp, int nbits)
+{
+ int ns = (nbits / CHAR_BIT);
+
+ if (nbits % CHAR_BIT)
+ ++ns;
+
+ if (ns > sp->size) {
+ byte *tmp;
+ int ix;
+
+ tmp = calloc(ns, sizeof(byte));
+ if (tmp == NULL) {
+ fprintf(stderr, "Error: out of memory in sieve_grow\n");
+ return;
+ }
+
+ memcpy(tmp, sp->bits, sp->size);
+ for (ix = sp->size; ix < ns; ix++) {
+ tmp[ix] = UCHAR_MAX;
+ }
+
+ free(sp->bits);
+ sp->bits = tmp;
+ sp->size = ns;
+
+ sp->nbits = sp->size * CHAR_BIT;
+ }
+}
+
+long
+sieve_next(sieve *sp)
+{
+ long out;
+ int ix = 0;
+ long val;
+
+ if (sp->next > sp->nbits)
+ return -1;
+
+ out = S_VAL(sp, sp->next);
+#ifdef DEBUG
+ fprintf(stderr, "Sieving %ld\n", out);
+#endif
+
+ /* Sieve out all multiples of the current prime */
+ val = out;
+ while (ix < sp->nbits) {
+ val += out;
+ ix = S_BIT(sp, val);
+ if ((val & 1) && ix < sp->nbits) { /* && S_ISSET(sp, ix)) { */
+ S_CLR(sp, ix);
+#ifdef DEBUG
+ fprintf(stderr, "Crossing out %ld (bit %d)\n", val, ix);
+#endif
+ }
+ }
+
+ /* Scan ahead to the next prime */
+ ++sp->next;
+ while (sp->next < sp->nbits && !S_ISSET(sp, sp->next))
+ ++sp->next;
+
+ return out;
+}
+
+void
+sieve_cross(sieve *sp, long val)
+{
+ int ix = 0;
+ long cur = val;
+
+ while (cur < sp->base)
+ cur += val;
+
+ ix = S_BIT(sp, cur);
+ while (ix < sp->nbits) {
+ if (cur & 1)
+ S_CLR(sp, ix);
+ cur += val;
+ ix = S_BIT(sp, cur);
+ }
+}
+
+void
+sieve_reset(sieve *sp, long base)
+{
+ memset(sp->bits, UCHAR_MAX, sp->size);
+ sp->base = base;
+ sp->next = 0;
+}
+
+void
+sieve_clear(sieve *sp)
+{
+ if (sp->bits)
+ free(sp->bits);
+
+ sp->bits = NULL;
+}