summaryrefslogtreecommitdiffhomepage
path: root/coreutils/factor.c
blob: 0f838a735da9340d754e379ff4106e0133ca5f16 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
/*
 * Copyright (C) 2017 Denys Vlasenko <vda.linux@googlemail.com>
 *
 * Licensed under GPLv2, see file LICENSE in this source tree.
 */
//config:config FACTOR
//config:	bool "factor"
//config:	default y
//config:	help
//config:	  factor factorizes integers

//applet:IF_FACTOR(APPLET(factor, BB_DIR_USR_BIN, BB_SUID_DROP))

//kbuild:lib-$(CONFIG_FACTOR) += factor.o

//usage:#define factor_trivial_usage
//usage:       "NUMBER..."
//usage:#define factor_full_usage "\n\n"
//usage:       "Print prime factors"

#include "libbb.h"

#if ULLONG_MAX == (UINT_MAX * UINT_MAX + 2 * UINT_MAX)
/* "unsigned" is half as wide as ullong */
typedef unsigned half_t;
#define HALF_FMT ""
#elif ULLONG_MAX == (ULONG_MAX * ULONG_MAX + 2 * ULONG_MAX)
/* long is half as wide as ullong */
typedef unsigned long half_t;
#define HALF_FMT "l"
#else
#error Cant find an integer type which is half as wide as ullong
#endif

static void factorize(const char *numstr)
{
	unsigned long long N, factor2;
	half_t factor;
	unsigned count3;

	/* Coreutils compat */
	numstr = skip_whitespace(numstr);
	if (*numstr == '+')
		numstr++;

	N = bb_strtoull(numstr, NULL, 10);
	if (errno)
		bb_show_usage();

	printf("%llu:", N);

	if (N < 4)
		goto end;
	while (!(N & 1)) {
		printf(" 2");
		N >>= 1;
	}

	count3 = 3;
	factor = 3;
	factor2 = 3 * 3;
	for (;;) {
		unsigned long long nfactor2;

		while ((N % factor) == 0) {
			N = N / factor;
			printf(" %u"HALF_FMT"", factor);
		}
 next_factor:
		/* (f + 2)^2 = f^2 + 4*f + 4 = f^2 + 4*(f+1) */
		nfactor2 = factor2 + 4 * (factor + 1);
		if (nfactor2 < factor2) /* overflow? */
			break;
		factor2 = nfactor2;
		if (factor2 > N)
			break;
		factor += 2;
		/* Rudimentary wheel sieving: skip multiples of 3:
		 * Every third odd number is divisible by three and thus isn't a prime:
		 * 5  7  9  11  13  15  17 19 21 23 25 27 29 31 33 35 37...
		 * ^  ^     ^   ^       ^  ^     ^  _     ^  ^     _  ^ (^primes)
		 */
		--count3;
		if (count3 == 0) {
			count3 = 3;
			goto next_factor;
		}
	}
 end:
	if (N > 1)
		printf(" %llu", N);
	bb_putchar('\n');
}

int factor_main(int argc, char **argv) MAIN_EXTERNALLY_VISIBLE;
int factor_main(int argc UNUSED_PARAM, char **argv)
{
	//// coreutils has undocumented option ---debug (three dashes)
	//getopt32(argv, "");
	//argv += optind;
	argv++;

	if (!*argv)
		//TODO: read from stdin
		bb_show_usage();

	do {
		factorize(*argv);
	} while (*++argv);

	fflush_stdout_and_exit(EXIT_SUCCESS);
}