--- 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.1Received 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