[dev] [sbase][PATCH] Add factor(1)

From: Mattias Andrée <maandree_AT_kth.se>
Date: Thu, 25 Feb 2016 21:48:52 +0100

Signed-off-by: Mattias Andrée <maandree_AT_kth.se>
---
 LICENSE  |   1 +
 Makefile |   4 +
 README   |   1 +
 factor.1 |  62 ++++++
 factor.c | 667 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 735 insertions(+)
 create mode 100644 factor.1
 create mode 100644 factor.c
diff --git a/LICENSE b/LICENSE
index 1ff215b..bd75908 100644
--- a/LICENSE
+++ b/LICENSE
_AT_@ -61,3 +61,4 @@ Authors/contributors include:
 © 2015 Wolfgang Corcoran-Mathe <first.lord.of.teal_AT_gmail.com>
 © 2016 Mattias Andrée <maandree_AT_kth.se>
 © 2016 Eivind Uggedal <eivind_AT_uggedal.com>
+© 2016 Joel Mickelin <joel.mickelin_AT_mykolab.ch>
diff --git a/Makefile b/Makefile
index 5fd0b50..50f30cd 100644
--- a/Makefile
+++ b/Makefile
_AT_@ -104,6 +104,7 @@ BIN =\
 	env\
 	expand\
 	expr\
+	factor\
 	false\
 	find\
 	flock\
_AT_@ -189,6 +190,9 @@ $(BIN): $(LIB) $(@:=.o)
 
 $(OBJ): $(HDR) config.mk
 
+sbase: LDFLAGS += -lgmp -lpthread
+factor: LDFLAGS += -lgmp -lpthread
+
 .o:
 	$(CC) $(LDFLAGS) -o $_AT_ $< $(LIB)
 
diff --git a/README b/README
index aa18ede..c5834c9 100644
--- a/README
+++ b/README
_AT_@ -34,6 +34,7 @@ The following tools are implemented:
 =*|o env             .
 #*|o expand          .
 #*|o expr            .
+ *   factor          .
 =*|o false           .
 =    find            .
 =* x flock           .
diff --git a/factor.1 b/factor.1
new file mode 100644
index 0000000..7386d39
--- /dev/null
+++ b/factor.1
_AT_@ -0,0 +1,62 @@
+.Dd 2016-02-25
+.Dt FACTOR 1
+.Os sbase
+.Sh NAME
+.Nm factor
+.Nd prime factor numbers
+.Sh SYNOPSIS
+.Nm
+.Op Fl c Ar N
+.Op Fl p Ar N
+.Op Fl q
+.Op Ar number ...
+.Sh DESCRIPTION
+.Nm
+prints the prime factors of
+.Ar number .
+If no
+.Ar number
+has been specified, the numbers to factor are
+line by line from stdin.
+.Sh OPTIONS
+.Bl -tag -width Ds
+.It Fl c Ar N
+Select the certainty for prime factors. The probably
+that a composite number will be listed as a prime
+factor will be less than 4^(-\fIN\fP).
+.It Fl p Ar N
+Use at most
+.Ar N
+threads to factor each
+.Ar number .
+This is seldom useful.
+Note, the input numbers are not factored concurrently,
+only there factors.
+.It Fl q
+Mark factors with \(aq?\(aq if it is not certain that
+they are prime.
+.El
+.Sh EXIT STATUS
+.Bl -tag -width Ds
+.It 0
+All numbers were successfully factorized.
+.It > 1
+An error occurred.
+.El
+.Sh NOTES
+The factors are not output in ascending order,
+but in order of discovery.
+.Pp
+The prime test is two-stage. The first stage
+is unaffected by the
+.Fl c
+flag. In this stage it is possible that the
+number is designated as a prime with complete
+certainty. In second stage, which is affected
+by the
+.Fl c
+flag, the number cannot be designated as a
+prime with complete certainty.
+.Sh SEE ALSO
+.Xr bc 1
+
diff --git a/factor.c b/factor.c
new file mode 100644
index 0000000..350dfc1
--- /dev/null
+++ b/factor.c
_AT_@ -0,0 +1,667 @@
+/* See LICENSE file for copyright and license details. */
+#include <stdlib.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <ctype.h>
+#include <errno.h>
+#include <semaphore.h>
+#include <pthread.h>
+
+#include "util.h"
+
+/*
+ * Optimisations that have been excluded.
+ * 
+ * For composite N, find R, B ∈ ℕ : B↑R = N. If found, continue by
+ * multiplying root_order by R and factor B instead of N. This is
+ * done by calculating N↑−P for all primes P : P ≤ log₃ N, if the
+ * result is an integer, replace N with N↑−P, multiply R with
+ * P and try P again. The final N is B. This is useful in some
+ * cases, but not overall.
+ */
+
+static int certainty = 5;
+static int parallel = 1;
+static const char *questioned = "";
+
+#if !defined(USE_GMP) && !defined(USE_TOMMATH)
+# define USE_GMP /* GMP is significantly fast than libtommath */
+#endif
+
+#if defined(USE_GMP)
+# include <gmp.h>
+# define lowest_bit(x)               mpz_scan1(x, 0)
+# define shift_right(r, x, steps)    mpz_tdiv_q_2exp(r, x, steps)
+# define prime_test(x)               mpz_probab_prime_p(x, certainty)
+# define to_string(x)                mpz_get_str(strbuf, 10, x)
+# define div_mod(q, r, n, d)         mpz_tdiv_qr(q, r, n, d)
+# define bit_count(x)                mpz_sizeinbase(x, 2)
+# define gcd(r, a, b)                mpz_gcd(r, a, b)
+# define zabs(r, x)                  mpz_abs(r, x)
+# define zmul(r, a, b)               mpz_mul(r, a, b)
+# define zmul_i(r, a, b)             mpz_mul_ui(r, a, b)
+# define zadd(r, a, b)               mpz_add(r, a, b)
+# define zadd_i(r, a, b)             mpz_add_ui(r, a, b)
+# define zsub(r, a, b)               mpz_sub(r, a, b)
+# define zmod(r, a, b)               mpz_mod(r, a, b)
+# define zsqrt(r, x)                 mpz_sqrt(r, x)
+# define zsqrt_test(r, x)            mpz_root(r, x, 2)
+# define zinit(x)                    mpz_init(x)
+# define zfree(x)                    mpz_clear(x)
+# define zset(r, x)                  mpz_set(r, x)
+# define zset_i(r, x)                mpz_set_ui(r, x)
+# define zcmp(a, b)                  mpz_cmp(a, b)
+# define zcmp_i(a, b)                mpz_cmp_ui(a, b)
+# define zparse(r, s)                mpz_init_set_str(r, s, 10);
+typedef mpz_t bigint_t;
+#elif defined(USE_TOMMATH)
+# include <tommath.h>
+# define lowest_bit(x)               mp_cnt_lsb(x)
+# define shift_right(r, x, steps)    mp_div_2d(x, steps, r, 0)
+static int prime_test(mp_int *x)     { int ret; mp_prime_is_prime(x, certainty, &ret); return ret;}
+# define to_string(x)                (mp_todecimal(x, strbuf), strbuf)
+# define div_mod(q, r, n, d)         mp_div(n, d, q, r)
+# define bit_count(x)                mp_count_bits(x)
+# define gcd(r, a, b)                mp_gcd(a, b, r)
+# define zabs(r, x)                  mp_abs(x, r)
+# define zmul(r, a, b)               mp_mul(a, b, r)
+# define zmul_i(r, a, b)             (zset_i(ctx->tmp, b), zmul(r, a, ctx->tmp))
+# define zadd(r, a, b)               mp_add(a, b, r)
+# define zadd_i(r, a, b)             (zset_i(ctx->tmp, b), zadd(r, a, ctx->tmp))
+# define zsub(r, a, b)               mp_sub(a, b, r)
+# define zmod(r, a, b)               mp_mod(a, b, r)
+# define zsqrt(r, x)                 mp_sqrt(x, r)
+static int zsqrt_test(mp_int *r, mp_int *x)  { int ret; mp_is_square(x, &ret); zsqrt(r, x); return ret; }
+# define zinit(x)                    mp_init(x)
+# define zfree(x)                    mp_clear(x)
+/*# define zset(r, x)                  mp_copy(x, r) /\* TODO Is this really correct? */
+static void zset(mp_int *r, mp_int *x)  { mp_zero(r); zadd(r, r, x); }
+# define zset_i(r, x)                mp_set_int(r, x)
+# define zcmp(a, b)                  mp_cmp(a, b)
+# define zcmp_i(a, b)                (zset_i(ctx->tmp, b), zcmp(a, ctx->tmp))
+# define zparse(r, s)                (zinit(r), mp_read_radix(r, s, 10))
+typedef mp_int bigint_t[1];
+#endif
+
+#define elementsof(x)               (sizeof(x) / sizeof(*x))
+#define is_factorised(x)            (!zcmp_i(x, 1))
+
+enum { NO = 0, MAYBE, YES };
+
+enum { POLLARDS_RHO_INITIAL_SEED = 0 };
+enum { POLLARDS_RHO_SEED_INCREASEMENT = 500 };
+
+struct context {
+	bigint_t div_n, div_q, div_r, div_d;
+	bigint_t *div_stack;
+	size_t div_stack_size;
+	bigint_t factor, d, x, y, conj_a, conj_b;
+#ifdef USE_TOMMATH
+	bigint_t tmp;
+#endif
+};
+
+struct thread_data {
+	bigint_t integer;
+	size_t root_order;
+};
+
+#define LIST_CONSTANTS X(3) X(5) X(7) X(11) X(13) X(17)
+static bigint_t constants[18];
+
+#define _5(x) x, x, x, x, x
+#define _25(x) _5(_5(x))
+#define _50(x) _25(x), _25(x)
+static const long pollards_rho_seeds[] = {
+	[0] = _50(2),
+	[4] = _50(11),
+	[8] = _50(101),
+	[16] = _50(503),
+	[54] = _25(4993),
+	[64] = _5(6029),
+	[71] = _5(6500),
+	[72] = _5(7001),
+	[74] = _5(7559),
+	[78] = _5(8017),
+	[82] = _5(7559),
+	[83] = _5(7000),
+	[84] = _5(6500),
+	[85] = _5(6029),
+	[86] = _5(10711),
+	[89] = _5(17041),
+	[92] = _5(34511)
+};
+
+static char *strbuf;
+static void (*output_primes)(bigint_t, int, size_t);
+static void (*subfactor)(struct context *, bigint_t, size_t, bigint_t, bigint_t);
+static sem_t semaphore;
+static pthread_mutex_t print_mutex;
+static pthread_mutex_t counter_mutex;
+static pthread_cond_t counter_condition;
+static int running = 0;
+#ifdef DEBUG
+static bigint_t result;
+static bigint_t expected;
+#endif
+
+static void
+context_init(struct context *ctx, bigint_t integer)
+{
+	size_t n;
+
+	if (!integer) {
+		ctx->div_stack_size = 0;
+		ctx->div_stack = 0;
+	} else {
+		n = ctx->div_stack_size = bit_count(integer);
+		ctx->div_stack = emalloc(n * sizeof(bigint_t));
+		while (n--)
+			zinit(ctx->div_stack[n]);
+	}
+
+	zinit(ctx->div_n);
+	zinit(ctx->div_q);
+	zinit(ctx->div_r);
+	zinit(ctx->div_d);
+
+	zinit(ctx->factor);
+	zinit(ctx->d);
+	zinit(ctx->x);
+	zinit(ctx->y);
+	zinit(ctx->conj_a);
+	zinit(ctx->conj_b);
+
+#ifdef USE_TOMMATH
+	zinit(ctx->tmp);
+#endif
+}
+
+static void
+context_reinit(struct context *ctx, bigint_t integer)
+{
+	size_t i, n = bit_count(integer);
+
+	if (n > ctx->div_stack_size) {
+		i = ctx->div_stack_size;
+		ctx->div_stack_size = n;
+		ctx->div_stack = erealloc(ctx->div_stack, n * sizeof(bigint_t));
+		while (i < n)
+			zinit(ctx->div_stack[i++]);
+	}
+}
+
+static void
+context_free(struct context *ctx)
+{
+	size_t n;
+
+	for (n = ctx->div_stack_size; n--;)
+		zfree(ctx->div_stack[n]);
+	free(ctx->div_stack);
+
+	zfree(ctx->div_n);
+	zfree(ctx->div_q);
+	zfree(ctx->div_r);
+	zfree(ctx->div_d);
+
+	zfree(ctx->factor);
+	zfree(ctx->d);
+	zfree(ctx->x);
+	zfree(ctx->y);
+	zfree(ctx->conj_a);
+	zfree(ctx->conj_b);
+
+#ifdef USE_TOMMATH
+	zfree(ctx->tmp);
+#endif
+}
+
+static void
+parallel_init(void)
+{
+	if (sem_init(&semaphore, 0, parallel))
+		eprintf("sem_init:");
+	if ((errno = pthread_mutex_init(&print_mutex, 0)))
+		eprintf("pthread_mutex_init:");
+	if ((errno = pthread_mutex_init(&counter_mutex, 0)))
+		eprintf("pthread_mutex_init:");
+	if ((errno = pthread_cond_init(&counter_condition, 0)))
+		eprintf("pthread_cond_init:");
+}
+
+static void
+parallel_term(void)
+{
+	if (sem_destroy(&semaphore))
+		eprintf("sem_destroy:");
+	if ((errno = pthread_mutex_destroy(&print_mutex)))
+		eprintf("pthread_mutex_destroy:");
+	if ((errno = pthread_mutex_destroy(&counter_mutex)))
+		eprintf("pthread_mutex_destroy:");
+	if ((errno = pthread_cond_destroy(&counter_condition)))
+		eprintf("pthread_cond_destroy:");
+}
+
+static void
+output_primes_nonparallel(bigint_t factor, int is_prime, size_t power)
+{
+	const char *fstr = to_string(factor);
+	const char *qstr = is_prime == MAYBE ? questioned : "";
+	while (power--) {
+		printf(" %s%s", fstr, qstr);
+#ifdef DEBUG
+		zmul(result, result, factor);
+#endif
+	}
+}
+
+static void
+output_primes_parallel(bigint_t factor, int is_prime, size_t power)
+{
+	const char *fstr, *qstr = is_prime == MAYBE ? questioned : "";
+	if (pthread_mutex_lock(&print_mutex))
+		eprintf("pthread_mutex_lock:");
+	fstr = to_string(factor);
+	while (power--) {
+		printf(" %s%s", fstr, qstr);
+#ifdef DEBUG
+		zmul(result, result, factor);
+#endif
+	}
+	if (pthread_mutex_unlock(&print_mutex))
+		eprintf("pthread_mutex_unlock:");
+}
+
+static ssize_t
+iterated_division(struct context *ctx, bigint_t remainder, bigint_t numerator,
+		  bigint_t denominator, size_t root_order, int is_prime)
+{
+	/*
+	 * Just like n↑m by squaring, excepted this is iterated division.
+	 */
+
+	const char *dstr = root_order ? to_string(denominator) : 0;
+	const char *nstr = is_prime == MAYBE ? questioned : "";
+	size_t partial_times = 1, times = 0, out, i;
+	bigint_t *n = &ctx->div_n, *q = &ctx->div_q, *r = &ctx->div_r, *d = &ctx->div_d;
+	bigint_t *div_stack = ctx->div_stack;
+
+	zset(*n, numerator);
+	zset(*d, denominator);
+	zset(*div_stack++, denominator);
+
+	for (;;) {
+		zmul(*d, *d, *d);
+		if (zcmp(*d, *n) <= 0) {
+			zset(*div_stack++, *d);
+			partial_times <<= 1;
+		} else {
+			out = root_order * partial_times;
+			for (; partial_times; out >>= 1, partial_times >>= 1) {
+				div_mod(*q, *r, *n, *--div_stack);
+				if (!zcmp_i(*r, 0)) {
+					for (i = 0; i < out; i++)
+						printf(" %s%s", dstr, nstr);
+					times |= partial_times;
+					zset(*n, *q);
+				}
+			}
+			zset(remainder, *n);
+			return times;
+		}
+	}
+}
+
+static void
+subfactor_proper(struct context *ctx, bigint_t factor, bigint_t c, bigint_t next, size_t root_order)
+{
+	size_t bits, cd;
+	bigint_t *d = &ctx->d, *x = &ctx->x, *y = &ctx->y;
+	bigint_t *conj_a = &ctx->conj_a, *conj_b = &ctx->conj_b;
+	int is_prime;
+	long seed;
+
+beginning:
+	bits = bit_count(factor);
+
+	if (bits < elementsof(pollards_rho_seeds))
+		seed = pollards_rho_seeds[bits];
+	else
+		seed = pollards_rho_seeds[elementsof(pollards_rho_seeds) - 1];
+
+	zadd_i(*x, c, seed);
+	zset(*y, *x);
+
+	for (;;) {
+		/*
+		 * There may exist a number b = (A = ⌊√n⌋ + 1)² − n such that B = √b is an integer
+		 * in which case n = (A − B)(A + B)  [n is(!) odd composite]. If not, the only the
+		 * trivial iteration of A (A += 1) seems to be the one not consuming your entire
+		 * CPU and it is also guaranteed to find the factors, but it is just so slow.
+		 */
+		zsqrt(*conj_a, factor);
+		zadd_i(*conj_a, *conj_a, 1);
+		zmul(*conj_b, *conj_a, *conj_a);
+		zsub(*conj_b, *conj_b, factor);
+
+		if (zsqrt_test(*conj_b, *conj_b)) {
+			zadd(next, *conj_a, *conj_b);
+			zsub(factor, *conj_a, *conj_b);
+			subfactor(ctx, next, root_order, 0, 0);
+			subfactor(ctx, factor, root_order, c, next);
+			break;
+		}
+
+		/* Pollard's rho algorithm with Floyd's cycle-finding algorithm and seed.
+		 * A special-purpose integer factorisation algorithm used for factoring
+		 * integers with small factors. */
+		do {
+			zmul(*x, *x, *x), zadd_i(*x, *x, seed);
+			zmul(*y, *y, *y), zadd_i(*y, *y, seed);
+			zmul(*y, *y, *y), zadd_i(*y, *y, seed);
+			zmod(*x, *x, factor);
+			zmod(*y, *y, factor);
+
+			zsub(*d, *x, *y);
+			zabs(*d, *d);
+			gcd(*d, *d, factor);
+		} while (!zcmp_i(*d, 1));
+
+		if (!zcmp(factor, *d)) {
+			if ((is_prime = prime_test(factor))) {
+				output_primes(factor, is_prime, root_order);
+				break;
+			} else {
+				zadd_i(c, c, POLLARDS_RHO_SEED_INCREASEMENT);
+				goto beginning;
+			}
+		}
+
+		if ((is_prime = prime_test(factor))) {
+			iterated_division(ctx, factor, factor, *d, root_order, is_prime);
+			if (is_factorised(factor))
+				break;
+		} else {
+			cd = iterated_division(ctx, factor, factor, *d, 0, 0);
+			zset(next, *d);
+			subfactor(ctx, next, root_order * cd, 0, 0);
+			if (is_factorised(factor))
+				break;
+		}
+
+		if ((is_prime = prime_test(factor))) {
+			output_primes(factor, is_prime, root_order);
+			break;
+		}
+	}
+}
+
+static void
+subfactor_nonparallel(struct context *ctx, bigint_t integer, size_t root_order,
+		      bigint_t reuse_seed, bigint_t reuse_next)
+{
+	if (reuse_seed) {
+		zset_i(reuse_seed, POLLARDS_RHO_INITIAL_SEED);
+		subfactor_proper(ctx, integer, reuse_seed, reuse_next, root_order);
+	} else {
+		bigint_t seed, next;
+		zinit(seed);
+		zset_i(seed, POLLARDS_RHO_INITIAL_SEED);
+		zinit(next);
+		subfactor_proper(ctx, integer, seed, next, root_order);
+		zfree(seed);
+		zfree(next);
+	}
+}
+
+static void *
+subfactor_parallel_threaded(void *data_)
+{
+	struct thread_data *data = data_;
+	struct context ctx;
+
+	if (sem_wait(&semaphore))
+		eprintf("sem_wait:");
+
+	output_primes(data->integer, 2, 1);
+
+	context_init(&ctx, data->integer);
+	subfactor_nonparallel(&ctx, data->integer, data->root_order, 0, 0);
+	context_free(&ctx);
+	zfree(data->integer);
+	free(data);
+
+	if (sem_post(&semaphore))
+		eprintf("sem_post:");
+
+	if ((errno = pthread_mutex_lock(&counter_mutex)))
+		eprintf("pthread_mutex_lock:");
+	if (--running == 0) {
+		if ((errno = pthread_cond_signal(&counter_condition)))
+			eprintf("pthread_cond_signal:");
+	}
+	if ((errno = pthread_mutex_unlock(&counter_mutex)))
+		eprintf("pthread_mutex_unlock:");
+
+	return NULL;
+}
+
+static void
+subfactor_parallel(struct context *ctx, bigint_t integer, size_t root_order,
+		   bigint_t reuse_seed, bigint_t reuse_next)
+{
+	if (reuse_seed) {
+		subfactor_nonparallel(ctx, integer, root_order, reuse_seed, reuse_next);
+	} else {
+		struct thread_data *data = emalloc(sizeof(*data));
+		pthread_t thread;
+		zinit(data->integer);
+		zset(data->integer, integer);
+		data->root_order = root_order;
+
+		if ((errno = pthread_mutex_lock(&counter_mutex)))
+			eprintf("pthread_mutex_lock:");
+		running++;
+
+		if ((errno = pthread_mutex_unlock(&counter_mutex)))
+			eprintf("pthread_mutex_unlock:");
+
+		if ((errno = pthread_create(&thread, 0, subfactor_parallel_threaded, data)))
+			eprintf("pthread_create:");
+		if ((errno = pthread_detach(thread)))
+			eprintf("pthread_detach:");
+	}
+}
+
+static int
+factor(struct context *ctx, char *integer_str, bigint_t reusable_seed, bigint_t reusable_next)
+{
+	bigint_t integer;
+	size_t i, power;
+	int is_prime;
+
+	if (!*integer_str)
+		goto invalid;
+	for (i = 0; integer_str[i]; i++)
+		if (!isdigit(integer_str[i]))
+			goto invalid;
+
+	zparse(integer, integer_str);
+#ifdef DEBUG
+	zset_i(result, 1);
+	zset(expected, integer);
+#endif
+
+	strbuf = integer_str;
+
+	while (*integer_str == '0' && *integer_str != 0) integer_str++;
+	printf("%s:", integer_str);
+
+	/* Behave like GNU factor: print empty set for 0 and 1, and pretend 0 is positive. */
+	if (zcmp_i(integer, 1) <= 0)
+		goto done;
+
+	/* Remove factors of two. */
+	power = lowest_bit(integer);
+	if (power > 0) {
+		shift_right(integer, integer, power);
+		while (power--) {
+			printf(" 2");
+#ifdef DEBUG
+			zmul_i(result, result, 2);
+#endif
+		}
+		if (is_factorised(integer))
+			goto done;
+	}
+
+	context_reinit(ctx, integer);
+
+	/* Remove factors of other tiny primes. */
+#ifdef DEBUG
+# define print_prime(factor)  printf(" "#factor), zmul_i(result, result, factor);
+#else
+# define print_prime(factor)  printf(" "#factor);
+#endif
+#define X(factor)\
+	power = iterated_division(ctx, integer, integer, constants[factor], 0, YES);\
+	if (power > 0) {\
+		while (power--)\
+			print_prime(factor);\
+		if (is_factorised(integer))\
+			goto done;\
+	}
+	LIST_CONSTANTS;
+#undef X
+
+	if ((is_prime = prime_test(integer))) {
+		output_primes(integer, is_prime, 1);
+		goto done;
+	}
+
+	if (parallel < 2) {
+		subfactor(ctx, integer, 1, reusable_seed, reusable_next);
+	} else {
+		if (sem_trywait(&semaphore))
+			eprintf("sem_trywait:");
+
+		running++;
+		subfactor(ctx, integer, 1, reusable_seed, reusable_next);
+
+		if (sem_post(&semaphore))
+			eprintf("sem_post:");
+
+		if ((errno = pthread_mutex_lock(&counter_mutex)))
+			eprintf("pthread_mutex_lock:");
+		if (--running != 0) {
+			if ((errno = pthread_cond_wait(&counter_condition, &counter_mutex)))
+				eprintf("pthread_cond_wait:");
+		}
+		if ((errno = pthread_mutex_unlock(&counter_mutex)))
+			eprintf("pthread_mutex_unlock:");
+	}
+
+#ifdef DEBUG
+	if (zcmp(result, expected))
+		fprintf(stderr, "\033[1;31mIncorrect factorisation of %s\033[m\n", to_string(expected));
+#endif
+
+done:
+	zfree(integer);
+	printf("\n");
+	return 0;
+invalid:
+	weprintf("%s is not a valid non-negative integer\n", integer_str);
+	return 1;
+}
+
+static void
+usage(void)
+{
+	eprintf("usage: %s [-c N] [-p N] [-q] [number ...]\n", argv0);
+}
+
+int
+main(int argc, char *argv[])
+{
+	long temp;
+	int ret = 0;
+	struct context ctx;
+	bigint_t reusable_seed, reusable_next;
+
+	ARGBEGIN {
+	case 'c':
+		temp = strtol(EARGF(usage()), NULL, 10);
+		if (temp < 1)
+			eprintf("value of -c must be a positive integer\n");
+		certainty = temp > INT_MAX ? INT_MAX : (int)temp;
+		break;
+	case 'p':
+		temp = atoi(EARGF(usage()));
+		if (temp < 1)
+			eprintf("value of -p must be a positive integer\n");
+		parallel = temp > INT_MAX ? INT_MAX : (int)temp;
+		parallel = temp > SEM_VALUE_MAX ? SEM_VALUE_MAX : (int)temp;
+		break;
+	case 'q':
+		questioned = "?";
+		break;
+	default:
+		usage();
+	} ARGEND;
+
+#define X(x)  zinit(constants[x]), zset_i(constants[x], x);
+	LIST_CONSTANTS;
+#undef X
+
+	if (parallel > 1) {
+		output_primes = output_primes_parallel;
+		subfactor = subfactor_parallel;
+		parallel_init();
+	} else {
+		output_primes = output_primes_nonparallel;
+		subfactor = subfactor_nonparallel;
+	}
+
+	context_init(&ctx, 0);
+	zinit(reusable_seed);
+	zinit(reusable_next);
+#ifdef DEBUG
+	zinit(result);
+	zinit(expected);
+#endif
+
+	if (*argv) {
+		while (*argv)
+			ret |= factor(&ctx, *argv++, reusable_seed, reusable_next);
+	} else {
+		ssize_t n;
+		size_t size = 0;
+		char *line = 0;
+		while ((n = getline(&line, &size, stdin)) >= 0) {
+			if (n > 0 && line[n - 1] == '\n')
+				n--;
+			line[n] = 0;
+			ret |= *line ? factor(&ctx, line, reusable_seed, reusable_next) : 0;
+		}
+		free(line);
+	}
+
+	context_free(&ctx);
+	zfree(reusable_seed);
+	zfree(reusable_next);
+#ifdef DEBUG
+	zfree(result);
+	zfree(expected);
+#endif
+
+	if (parallel > 1)
+		parallel_term();
+
+#define X(x)  zfree(constants[x]);
+	LIST_CONSTANTS;
+#undef X
+
+	return fshut(stdout, "<stdout>") || ret;
+}
-- 
2.7.1
Received on Thu Feb 25 2016 - 21:48:52 CET

This archive was generated by hypermail 2.3.0 : Thu Feb 25 2016 - 22:00:14 CET