不排序直接升序枚举一个数的因子？
Enumerate factors of a number directly in ascending order without sorting?
有没有一种高效的算法来枚举一个数字n的因数，按升序排列，不排序？我所说的“高效”是指：
该算法从 n 的素数幂分解开始，避免了对除数的强力搜索。
算法的运行时复杂度为 O(d log₂ d) 或更好，其中 d 是 n 的除数。
算法的空间复杂度为O(d).
该算法避免了排序操作。也就是说，因子是按顺序 生成的，而不是乱序生成然后再排序的。尽管使用简单的递归方法枚举然后排序是 O(d log₂ d)，但所有涉及的内存访问的成本非常低在排序中。
一个简单的例子是 n = 360 = 2³ × 3² × 5，它有 d=24 个因子：{ 1, 2 , 3, 4, 5, 6, 8, 9, 10, 12, 15, 18, 20, 24, 30, 36, 40, 45, 60, 72, 90, 120, 180, 360 }.
更严重的例子是n = 278282512406132373381723386382308832000 = 2⁸ × 3⁴ × 5³ × 7² × 11² × 13² × 17 × 19 × 23 × 29 × 31 × 37 × 41 × 43 × 47 × 53 × 59 × 61 × 67 × 71 × 73 × 79，其中有 d=318504960 个因子（显然太多了，无法在这里列出！）。顺便说一下，这个数是最多 2^128 的任何数的最大因数。
我可以发誓几周前我看到了这种算法的描述，带有示例代码，但现在我似乎无法在任何地方找到它。它使用了一些魔术，在输出列表中为每个素因子维护一个祖先索引列表。 （更新：我将因子生成与操作类似的汉明数混淆了。）
更新
我最终在运行时使用了 O(d) 的解决方案，创建 O(d 的内存开销极低) 就地输出，并且比我知道的任何其他方法都快得多。我已经发布了这个解决方案作为答案，带有 C 源代码。它是另一位贡献者 Will Ness 在 Haskell 中介绍的漂亮算法的高度优化简化版本。我选择 Will 的回答作为接受的答案，因为它提供了一个非常优雅的解决方案，符合最初陈述的所有要求。
简而言之：重复从堆中拉出 next-smallest 因子，然后推回该因子的每个倍数，该倍数仍然是 n 的因子。使用技巧来避免出现重复项，以便堆大小永远不会超过 d。时间复杂度为 O(kd log d)，其中 k 是不同质因数的数量。
我们使用的关键 属性 是如果 x 和 y 都是 n 的因数且 y = x*p 对于某个因数 p >= 2 —— 即如果 x 的质因数是 y 的素因子的真子集——那么 x < y。这意味着在将任何此类 y 插入堆之前输出 x 总是安全的。堆仅有效地用于比较两个都不是另一个子多重集的因子。
第一次尝试：重复会减慢速度
首先描述一个产生正确答案但也产生许多重复的算法会很有帮助：
- 设置 prev = NULL。
- 将 1 插入堆 H.
- 从H中取出堆顶t，如果堆为空，则停止。
- 如果 t == prev 然后转到 3。[编辑：固定]
- 输出 t.
- 设置 prev = t.
- 对于 n 的每个不同的质因数 p：
- 如果 n % (t*p) == 0（即如果 t*p 仍然是 n 的一个因子），将 t*p 推到 H 上。
- 转到 3.
上述算法唯一的问题是它可以多次生成相同的因子。例如，如果 n = 30，则因子 15 将生成为因子 5 的 child（通过乘以质因子 3），以及因子 3 的 child（乘以 5)。解决这个问题的一种方法是注意任何重复项到达堆顶时必须在连续块中读出，因此您可以简单地检查堆顶是否等于 just-extracted 值，并且如果是这样，继续提取和丢弃它。但可能有更好的方法：
从源头上消除重复
因子 x 有多少种生成方式？首先考虑 x 不包含重数 > 1 的素因子的情况。在这种情况下，如果它包含 m 个不同的素因子，则有 m-1 "parent" 个因子将其生成为 "child" 在前面的算法中——这些 parent 中的每一个都由 m 个质因数中的 m-1 个子集组成，其余质因数是添加到 child 的质因数。 （如果 x 有一个重数 > 1 的素因子，那么实际上有 m parents。）如果我们有办法确定这些 parents 中的一个恰好是 "chosen one" 实际上将 x 生成为 child，并且此规则导致可以在 parent 弹出时应用于每个 parent 的测试，那么我们就可以避免创建任何重复项。
我们可以使用以下规则：对于任何给定的 x，选择缺少 最大 x 的 m 个因子的潜在 parent y。这形成了一个简单的规则：A parent y 产生一个 child x 当且仅当 x = y*p 对于某个 p 大于或等于 y 中已有的任何质因数。这很容易测试：只需按降序遍历质因数，为每个质因数生成 children，直到我们找到一个已经整除 y 的质因数。在前面的示例中，parent 3 将产生 15，但 parent 5 不会（因为 3 < 5）——因此 15 确实只产生一次。对于 n = 30，完整的树如下所示：
1
/|\
2 3 5
/| \
6 10 15
|
30
注意每个因子只生成一次。
新的duplicate-free算法如下：
- 将 1 插入堆 H.
- 从H中取出堆顶t，如果堆为空，则停止。
- 输出 t.
- 对于 n 的每个不同的素因子 p 按降序排列：
- 如果 n % (t*p) == 0（即如果 t*p 仍然是 n 的一个因子），将 t*p 推到 H 上。
- 如果 t % p == 0（即如果 t 已经包含 p 作为一个因子）则停止。
- 转到 2.
我们可以合并多个流，这样一来就没有重复项了。
从列表 [1] 开始，对于每个唯一的素因子 p，我们将列表乘以 p 迭代 k 次（其中 k是 p) 的多重性以获得 k 个新列表，并将它们与传入列表合并在一起。
在Haskell,
ordfactors n = -- <----------------------<---<---<-----\
foldr g [1] -- g (19,1) (g (7,1) (g (3,2) (g (2,3) [1])))
. reverse -- [(19,1),(7,1),(3,2),(2,3)] ^
. primePowers -- [(2,3),(3,2),(7,1),(19,1)] |
$ n -- 9576 --->--->--->----/
where
g (p,k) xs = mergeAsTree
. take (k+1) -- take 3 [a,b,c,d..] = [a,b,c]
. iterate (map (p*)) -- iterate f x = [x, y, z,...]
$ xs -- where { y=f x; z=f y; ... }
{- g (2,3) [1] = let a0 = [1]
a1 = map (2*) a0 -- [2]
a2 = map (2*) a1 -- [4]
a3 = map (2*) a2 -- [8]
in mergeAsTree [ a0, a1, a2, a3 ] -- xs2 = [1,2,4,8]
g (3,2) xs2 = let b0 = xs2 -- [1,2,4,8]
b1 = map (3*) b0 -- [3,6,12,24]
b2 = map (3*) b1 -- [9,18,36,72]
in mergeAsTree [ b0, b1, b2 ] -- xs3 = [1,2,3,4,6,8,9,12,...]
g (7,1) xs3 = mergeAsTree [ xs3, map (7*) xs3 ] -- xs4 = [1,2,3,4,6,7,8,9,...]
g (19,1) xs4 = mergeAsTree [ xs4, map (19*) xs4 ]
-}
mergeAsTree xs = foldi merge [] xs -- [a,b] --> merge a b
-- [a,b,c,d] --> merge (merge a b) (merge c d)
foldi f z [] = z
foldi f z (x:xs) = f x (foldi f z (pairs f xs))
pairs f (x:y:t) = f x y : pairs f t
pairs _ t = t
foldi 以线性方式排列二进制 merges (which presuppose that there's no duplicates in the streams being merged, for streamlined operation) in a tree-like fashion for speed; whereas foldr。
primePowers n | n > 1 = -- 9576 => [(2,3),(3,2),(7,1),(19,1)]
map (\xs -> (head xs,length xs)) -- ^
. group -- [[2,2,2],[3,3],[7],[19]] |
. factorize -- [2,2,2,3,3,7,19] |
$ n -- 9576 --->--->--->---/
group 是一个标准函数，它将列表中的相邻等值组合在一起，而 factorize 是一个众所周知的试分算法，它以非递减顺序产生一个数的质因数:
factorize n | n > 1 = go n (2:[3,5..]) -- 9576 = 2*2*2*3*3*7*19
where -- produce prime factors only
go n (d:t)
| d*d > n = [n]
| r == 0 = d : go q (d:t)
| otherwise = go n t
where
(q,r) = quotRem n d
factorize _ = []
(.) 是一个函数组合运算符，将其右侧的参数函数的输出链接为左侧函数的输入。它（和 $）可以读作“of”。
[我发布这个答案只是为了形式上的完整性。我已经选择了别人的答案作为接受的答案。]
算法概述。 在寻找生成 in-memory 因子列表（在我的例子中是 64 位无符号值）的最快方法时，我选择了一种实现 two-dimensional 桶排序的混合算法，它利用了排序键的内部知识（即，它们只是整数，因此可以对其进行计算）。具体方法更接近于“ProxMapSort”，但具有两个级别的键（主要，次要）而不是只有一个。 主键 只是该值的以 2 为底的对数。 minor key 是在第二层桶中产生合理传播所需的值的最高有效数字的最小数量。因子首先生成为未分类因子的临时工作数组。接下来，分析它们的分布并分配和填充桶索引数组。最后，使用插入排序将因子直接存储到最终排序的数组中。绝大多数桶只有 1、2 或 3 个值。示例在源代码中给出，附在这个答案的底部。
空间复杂度。 内存利用率大约是 Quicksort-based 解决方案的 4 倍，但这实际上是微不足道的，因为在最坏情况下使用过的最大内存（对于 64 位输入）是 5.5 MB，其中 4.0 MB 只保留了几毫秒。
运行时复杂度。 性能远好于 hand-coded Quicksort-based 解决方案：对于具有非平凡数量的因子的数字，它不规则地约为 2.5快 x 倍。在我的 3.4 GHz 上。 Intel i7，它在0.0052秒内按排序顺序产生184,320个因子18,401,055,938,125,660,800，或每个因子约96个时钟周期，或每秒约3500万个因子。
图表。 内存和运行时性能分析了 47,616 个等价的规范代表 类 最多 2⁶⁴–1 的素数签名。这些是 64 位搜索中的 so-called“高度可分解的数字”space。
对于非平凡的因子计数，总运行时间比 Quicksort-based 解决方案好 ~2.5 倍，如下图所示：
每秒产生的排序因子数本质上是上面的倒数。在大约 2000 个因子的最佳点之后，per-factor 的性能下降，但下降幅度不大。性能受 L1、L2 和 L3 缓存大小以及被分解数的唯一素数的影响，大致随输入值的对数增加。
在此对数-对数图中，内存使用峰值是一条直线，因为它与因子数量的以 2 为底的对数成正比。请注意，内存使用峰值仅在很短的时间内出现； short-lived 工作数组在几毫秒内被丢弃。丢弃临时数组后，剩下的是最终的因子列表，这与 Quicksort-based 解决方案中看到的最小用法相同。
源代码。 下面附上一个 proof-of-concept C 编程语言（特别是 C11）程序。它已经在 x86-64 上使用 Clang/LLVM 进行了测试，尽管它在 GCC 上也应该可以正常工作。
/*==============================================================================
DESCRIPTION
This is a small proof-of-concept program to test the idea of "sorting"
factors using a form of bucket sort. The method is essentially a 2D version
of ProxMapSort that has tuned for vast, nonlinear distributions using two
keys (major, minor) rather than one. The major key is simply the floor of
the base-2 logarithm of the value, and the minor key is derived from the most
significant bits of the value.
INPUT
Input is given on the command line, either as a single argument giving the
number to be factored or an even number of arguments giving the 2-tuples that
comprise the prime-power factorization of the desired number. For example,
the number
75600 = 2^4 x 3^3 x 5^2 x 7
can be given by the following list of arguments:
2 4 3 3 5 2 7 1
Note: If a single number is given, it will require factoring to produce its
prime-power factorization. Since this is just a small test program, a very
crude factoring method is used that is extremely fast for small prime factors
but extremely slow for large prime factors. This is actually fine, because
the largest factor lists occur with small prime factors anyway, and it is the
production of large factor lists at which this program aims to be proficient.
It is simply not interesting to be fast at producing the factor list of a
number like 17293823921105882610 = 2 x 3 x 5 x 576460797370196087, because
it has only 32 factors. Numbers with tens or hundreds of thousands of
factors are much more interesting.
OUTPUT
Results are written to standard output. A list of factors in ascending order
is produced, followed by runtime (in microseconds) required to generate the
list (not including time to print it).
STATISTICS
Bucket size statistics for the 47616 canonical representatives of the prime
signature equivalence classes of 64-bit numbers:
==============================================================
Bucket size Total count of factored Total count of
b numbers needing size b buckets of size b
--------------------------------------------------------------
1 47616 (100.0%) 514306458 (76.2%)
2 47427 (99.6%) 142959971 (21.2%)
3 43956 (92.3%) 16679329 (2.5%)
4 27998 (58.8%) 995458 (0.1%)
5 6536 (13.7%) 33427 (<0.1%)
6 400 (0.8%) 729 (<0.1%)
7 12 (<0.1%) 18 (<0.1%)
--------------------------------------------------------------
~ 47616 (100.0%) 674974643 (100.0%)
--------------------------------------------------------------
Thus, no 64-bit number (of the input set) ever requires more than 7 buckets,
and the larger the bucket size the less frequent it is. This is highly
desirable. Note that although most numbers need at least 1 bucket of size 5,
the vast majority of buckets (99.9%) are of size 1, 2, or 3, meaning that
insertions are extremely efficient. Therefore, the use of insertion sort
for the buckets is clearly the right choice and is arguably optimal for
performance.
AUTHOR
Todd Lehman
2015/05/08
*/
#include <inttypes.h>
#include <limits.h>
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <assert.h>
typedef unsigned int uint;
typedef uint8_t uint8;
typedef uint16_t uint16;
typedef uint32_t uint32;
typedef uint64_t uint64;
#define ARRAY_CAPACITY(x) (sizeof(x) / sizeof((x)[0]))
//-----------------------------------------------------------------------------
// This structure is sufficient to represent the prime-power factorization of
// all 64-bit values. The field names ω and Ω are dervied from the standard
// number theory functions ω(n) and Ω(n), which count the number of unique and
// non-unique prime factors of n, respectively. The field name d is derived
// from the standard number theory function d(n), which counts the number of
// divisors of n, including 1 and n.
//
// The maximum possible value here of ω is 15, which occurs for example at
// n = 7378677391061896920 = 2^3 x 3^2 x 5 x 7 x 11 x 13 x 17 x 19 x 23 x 29
// 31 x 37 x 41 x 43 x 47, which has 15 unique prime factors.
//
// The maximum possible value of Ω here is 63, which occurs for example at
// n = 2^63 and n = 2^62 x 3, both of which have 63 non-unique prime factors.
//
// The maximum possible value of d here is 184320, which occurs at
// n = 18401055938125660800 = 2^7 x 3^4 x 5^2 x 7^2 x 11 x 13 x 17 x 19 x 23 x
// 29 x 31 x 37 x 41.
//
// Maximum possible exponents when exponents are sorted in decreasing order:
//
// Index Maximum Bits Example of n
// ----- ------- ---- --------------------------------------------
// 0 63 6 (2)^63
// 1 24 5 (2*3)^24
// 2 13 4 (2*3*5)^13
// 3 8 4 (2*3*5*7)^8
// 4 5 3 (2*3*5*7*11)^5
// 5 4 3 (2*3*5*7*11*13)^4
// 6 3 2 (2*3*5*7*11*13*17)^3
// 7 2 2 (2*3*5*7*11*13*17*19)^2
// 8 2 2 (2*3*5*7*11*13*17*19*23)^2
// 9 1 1 (2*3*5*7*11*13*17*19*23*29)^1
// 10 1 1 (2*3*5*7*11*13*17*19*23*29*31)^1
// 11 1 1 (2*3*5*7*11*13*17*19*23*29*31*37)^1
// 12 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41)^1
// 13 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41*43)^1
// 14 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41*43*47)^1
// ----- ------- ---- --------------------------------------------
// 15 63 37
//
#pragma pack(push, 8)
typedef struct
{
uint8 e[16]; // Exponents.
uint64 p[16]; // Primes in increasing order.
uint8 ω; // Count of prime factors without multiplicity.
uint8 Ω; // Count of prime factors with multiplicity.
uint32 d; // Count of factors of n, including 1 and n.
uint64 n; // Value of n on which all other fields of this struct depend.
}
PrimePowerFactorization; // 176 bytes with 8-byte packing
#pragma pack(pop)
#define MAX_ω 15
#define MAX_Ω 63
//-----------------------------------------------------------------------------
// Fatal error: print error message and abort.
void fatal_error(const char *format, ...)
{
va_list args;
va_start(args, format);
vfprintf(stderr, format, args);
exit(1);
}
//-----------------------------------------------------------------------------
// Compute 64-bit 2-adic integer inverse.
uint64 uint64_inv(const uint64 x)
{
assert(x != 0);
uint64 y = 1;
for (uint i = 0; i < 6; i++) // 6 = log2(log2(2**64)) = log2(64)
y = y * (2 - (x * y));
return y;
}
//------------------------------------------------------------------------------
// Compute 2 to arbitrary power. This is just a portable and abstract way to
// write a left-shift operation. Note that the use of the UINT64_C macro here
// is actually required, because the result of 1U<<x is not guaranteed to be a
// 64-bit result; on a 32-bit compiler, 1U<<32 is 0 or is undefined.
static inline
uint64 uint64_pow2(x)
{
return UINT64_C(1) << x;
}
//------------------------------------------------------------------------------
// Deduce native word size (int, long, or long long) for 64-bit integers.
// This is needed for abstracting certain compiler-specific intrinsic functions.
#if UINT_MAX == 0xFFFFFFFFFFFFFFFFU
#define UINT64_IS_U
#elif ULONG_MAX == 0xFFFFFFFFFFFFFFFFUL
#define UINT64_IS_UL
#elif ULLONG_MAX == 0xFFFFFFFFFFFFFFFFULL
#define UINT64_IS_ULL
#else
//error "Unable to deduce native word size of 64-bit integers."
#endif
//------------------------------------------------------------------------------
// Define abstracted intrinsic function for counting leading zeros. Note that
// the value is well-defined for nonzero input but is compiler-specific for
// input of zero.
#if defined(UINT64_IS_U) && __has_builtin(__builtin_clz)
#define UINT64_CLZ(x) __builtin_clz(x)
#elif defined(UINT64_IS_UL) && __has_builtin(__builtin_clzl)
#define UINT64_CLZ(x) __builtin_clzl(x)
#elif defined(UINT64_IS_ULL) && __has_builtin(__builtin_clzll)
#define UINT64_CLZ(x) __builtin_clzll(x)
#else
#undef UINT64_CLZ
#endif
//------------------------------------------------------------------------------
// Compute floor of base-2 logarithm y = log_2(x), where x > 0. Uses fast
// intrinsic function if available; otherwise resorts to hand-rolled method.
static inline
uint uint64_log2(uint64 x)
{
assert(x > 0);
#if defined(UINT64_CLZ)
return 63 - UINT64_CLZ(x);
#else
#define S(k) if ((x >> k) != 0) { y += k; x >>= k; }
uint y = 0; S(32); S(16); S(8); S(4); S(2); S(1); return y;
#undef S
#endif
}
//------------------------------------------------------------------------------
// Compute major key, given a nonzero number. The major key is simply the
// floor of the base-2 logarithm of the number.
static inline
uint major_key(const uint64 n)
{
assert(n > 0);
uint k1 = uint64_log2(n);
return k1;
}
//------------------------------------------------------------------------------
// Compute minor key, given a nonzero number, its major key, k1, and the
// bit-size b of major bucket k1. The minor key, k2, is is computed by first
// removing the most significant 1-bit from the number, because it adds no
// information, and then extracting the desired number of most significant bits
// from the remainder. For example, given the number n=1463 and a major bucket
// size of b=6 bits, the keys are computed as follows:
//
// Step 0: Given number n = 0b10110110111 = 1463
//
// Step 1: Compute major key: k1 = floor(log_2(n)) = 10
//
// Step 2: Remove high-order 1-bit: n' = 0b0110110111 = 439
//
// Step 3: Compute minor key: k2 = n' >> (k1 - b)
// = 0b0110110111 >> (10 - 6)
// = 0b0110110111 >> 4
// = 0b011011
// = 27
static inline
uint minor_key(const uint64 n, const uint k1, const uint b)
{
assert(n > 0); assert(k1 >= 0); assert(b > 0);
const uint k2 = (uint)((n ^ uint64_pow2(k1)) >> (k1 - b));
return k2;
}
//------------------------------------------------------------------------------
// Raw unsorted factor.
#pragma push(pack, 4)
typedef struct
{
uint64 n; // Value of factor.
uint32 k1; // Major key.
uint32 k2; // Minor key.
}
UnsortedFactor;
#pragma pop(pack)
//------------------------------------------------------------------------------
// Compute sorted list of factors, given a prime-power factorization.
static uint64 memory_usage;
uint64 *compute_factors(const PrimePowerFactorization ppf)
{
memory_usage = 0;
if (ppf.n == 0)
return NULL;
uint64 *sorted_factors = calloc(ppf.d, sizeof(*sorted_factors));
if (!sorted_factors)
fatal_error("Failed to allocate array of %"PRIu32" factors.", ppf.d);
memory_usage += ppf.d * sizeof(*sorted_factors);
UnsortedFactor *unsorted_factors = malloc(ppf.d * sizeof(*unsorted_factors));
if (!unsorted_factors)
fatal_error("Failed to allocate array of %"PRIu32" factors.", ppf.d);
memory_usage += ppf.d * sizeof(*unsorted_factors);
// These arrays are indexed by the major key of a number.
uint32 major_counts[64]; // Counts of factors in major buckets.
uint32 major_spans[64]; // Counts rounded up to power of 2.
uint32 major_bits[64]; // Base-2 logarithm of bucket size.
uint32 major_indexes[64]; // Indexes into minor array.
memset(major_counts, 0, sizeof(major_counts));
memset(major_spans, 0, sizeof(major_spans));
memset(major_bits, 0, sizeof(major_bits));
memset(major_indexes, 0, sizeof(major_indexes));
// --- Step 1: Produce unsorted list of factors from prime-power
// factorization. At the same time, count groups of factors by their
// major keys.
{
// This array is for counting in the multi-radix number system dictated by
// the exponents of the prime-power factorization. An invariant is that
// e[i] <= ppf.e[i] for all i (0 < i <ppf.ω).
uint8 e[MAX_ω];
for (uint i = 0; i < ppf.ω; i++)
e[i] = 0;
// Initialize inverse-prime-powers. This array allows for division by
// p[i]**e[i] extremely quickly in the main loop below. Note that 2-adic
// inverses are not defined for even numbers (of which 2 is the only prime),
// so powers of 2 must be handled specially.
uint64 pe_inv[MAX_ω];
for (uint i = 0; i < ppf.ω; i++)
{
uint64 pe = 1; for (uint j = 1; j <= ppf.e[i]; j++) pe *= ppf.p[i];
pe_inv[i] = uint64_inv(pe);
}
uint64 n = 1; // Current factor accumulator.
for (uint k = 0; k < ppf.d; k++) // k indexes into unsorted_factors[].
{
//printf("unsorted_factors[%u] = %"PRIu64" j = %u\n", k, n, j);
assert(ppf.n % n == 0);
unsorted_factors[k].n = n;
uint k1 = major_key(n);
assert(k1 < ARRAY_CAPACITY(major_counts));
unsorted_factors[k].k1 = k1;
major_counts[k1] += 1;
// Increment the remainder of the multi-radix number e[].
for (uint i = 0; i < ppf.ω; i++)
{
if (e[i] == ppf.e[i]) // Carrying is occurring.
{
if (ppf.p[i] == 2)
n >>= ppf.e[i]; // Divide n by 2 ** ppf.e[i].
else
n *= pe_inv[i]; // Divide n by ppf.p[i] ** ppf.e[i].
e[i] = 0;
}
else // Carrying is not occurring.
{
n *= ppf.p[i];
e[i] += 1;
break;
}
}
}
assert(n == 1); // n always cycles back to 1, not to ppf.n.
assert(unsorted_factors[ppf.d-1].n == ppf.n);
}
// --- Step 2: Define the major bits array, the major spans array, the major
// index array, and count the total spans.
uint32 total_spans = 0;
{
uint32 k = 0;
for (uint k1 = 0; k1 < ARRAY_CAPACITY(major_counts); k1++)
{
uint32 count = major_counts[k1];
uint32 bits = (count <= 1)? count : uint64_log2(count - 1) + 1;
major_bits[k1] = bits;
major_spans[k1] = (count > 0)? (UINT32_C(1) << bits) : 0;
major_indexes[k1] = k;
k += major_spans[k1];
}
total_spans = k;
}
// --- Step 3: Allocate and populate the minor counts array. Note that it
// must be initialized to zero.
uint32 *minor_counts = calloc(total_spans, sizeof(*minor_counts));
if (!minor_counts)
fatal_error("Failed to allocate array of %"PRIu32" counts.", total_spans);
memory_usage += total_spans * sizeof(*minor_counts);
for (uint k = 0; k < ppf.d; k++)
{
const uint64 n = unsorted_factors[k].n;
const uint k1 = unsorted_factors[k].k1;
const uint k2 = minor_key(n, k1, major_bits[k1]);
assert(k2 < major_spans[k1]);
unsorted_factors[k].k2 = k2;
minor_counts[major_indexes[k1] + k2] += 1;
}
// --- Step 4: Define the minor indexes array.
//
// NOTE: Instead of allocating a separate array, the earlier-allocated array
// of minor indexes is simply repurposed here using an alias.
uint32 *minor_indexes = minor_counts; // Alias the array for repurposing.
{
uint32 k = 0;
for (uint i = 0; i < total_spans; i++)
{
uint32 count = minor_counts[i]; // This array is the same array...
minor_indexes[i] = k; // ...as this array.
k += count;
}
}
// --- Step 5: Populate the sorted factors array. Note that the array must
// be initialized to zero earlier because values of zero are used
// as sentinels in the bucket lists.
for (uint32 i = 0; i < ppf.d; i++)
{
uint64 n = unsorted_factors[i].n;
const uint k1 = unsorted_factors[i].k1;
const uint k2 = unsorted_factors[i].k2;
// Insert factor into bucket using insertion sort (which happens to be
// extremely fast because we know the bucket sizes are always very small).
uint32 k;
for (k = minor_indexes[major_indexes[k1] + k2];
sorted_factors[k] != 0;
k++)
{
assert(k < ppf.d);
if (sorted_factors[k] > n)
{ uint64 t = sorted_factors[k]; sorted_factors[k] = n; n = t; }
}
sorted_factors[k] = n;
}
// --- Step 6: Validate array of sorted factors.
{
for (uint32 k = 1; k < ppf.d; k++)
{
if (sorted_factors[k] == 0)
fatal_error("Produced a factor of 0 at index %"PRIu32".", k);
if (ppf.n % sorted_factors[k] != 0)
fatal_error("Produced non-factor %"PRIu64" at index %"PRIu32".",
sorted_factors[k], k);
if (sorted_factors[k-1] == sorted_factors[k])
fatal_error("Duplicate factor %"PRIu64" at index %"PRIu32".",
sorted_factors[k], k);
if (sorted_factors[k-1] > sorted_factors[k])
fatal_error("Out-of-order factors %"PRIu64" and %"PRIu64" "
"at indexes %"PRIu32" and %"PRIu32".",
sorted_factors[k-1], sorted_factors[k], k-1, k);
}
}
free(minor_counts);
free(unsorted_factors);
return sorted_factors;
}
//------------------------------------------------------------------------------
// Compute prime-power factorization of a 64-bit value. Note that this function
// is designed to be fast *only* for numbers with very simple factorizations,
// e.g., those that produce large factor lists. Do not attempt to factor
// large semiprimes with this function. (The author does know how to factor
// large numbers efficiently; however, efficient factorization is beyond the
// scope of this small test program.)
PrimePowerFactorization compute_ppf(const uint64 n)
{
PrimePowerFactorization ppf;
if (n == 0)
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 };
}
else if (n == 1)
{
ppf = (PrimePowerFactorization){ .p = { 1 }, .e = { 1 },
.ω = 1, .Ω = 1, .d = 1, .n = 1 };
}
else
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = n };
uint64 m = n;
for (uint64 p = 2; p * p <= m; p += 1 + (p > 2))
{
if (m % p == 0)
{
assert(ppf.ω <= MAX_ω);
ppf.p[ppf.ω] = p;
ppf.e[ppf.ω] = 0;
while (m % p == 0)
{ m /= p; ppf.e[ppf.ω] += 1; }
ppf.d *= (1 + ppf.e[ppf.ω]);
ppf.Ω += ppf.e[ppf.ω];
ppf.ω += 1;
}
}
if (m > 1)
{
assert(ppf.ω <= MAX_ω);
ppf.p[ppf.ω] = m;
ppf.e[ppf.ω] = 1;
ppf.d *= 2;
ppf.Ω += 1;
ppf.ω += 1;
}
}
return ppf;
}
//------------------------------------------------------------------------------
// Parse prime-power factorization from a list of ASCII-encoded base-10 strings.
// The values are assumed to be 2-tuples (p,e) of prime p and exponent e.
// Primes must not exceed 2^64 - 1. Exponents must not exceed 2^8 - 1. The
// constructed value must not exceed 2^64 - 1.
PrimePowerFactorization parse_ppf(const uint pairs, const char *const values[])
{
assert(pairs <= MAX_ω);
PrimePowerFactorization ppf;
if (pairs == 0)
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 };
}
else
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = 1 };
for (uint i = 0; i < pairs; i++)
{
ppf.p[i] = (uint64)strtoumax(values[(i*2)+0], NULL, 10);
ppf.e[i] = (uint8)strtoumax(values[(i*2)+1], NULL, 10);
// Validate prime value.
if (ppf.p[i] < 2) // (Ideally this would actually do a primality test.)
fatal_error("Factor %"PRIu64" is invalid.", ppf.p[i]);
// Accumulate count of unique prime factors.
if (ppf.ω > UINT8_MAX - 1)
fatal_error("Small-omega overflow at factor %"PRIu64"^%"PRIu8".",
ppf.p[i], ppf.e[i]);
ppf.ω += 1;
// Accumulate count of total prime factors.
if (ppf.Ω > UINT8_MAX - ppf.e[i])
fatal_error("Big-omega wverflow at factor %"PRIu64"^%"PRIu8".",
ppf.p[i], ppf.e[i]);
ppf.Ω += ppf.e[i];
// Accumulate total divisor count.
if (ppf.d > UINT32_MAX / (1 + ppf.e[i]))
fatal_error("Divisor count overflow at factor %"PRIu64"^%"PRIu8".",
ppf.p[i], ppf.e[i]);
ppf.d *= (1 + ppf.e[i]);
// Accumulate value.
for (uint8 k = 1; k <= ppf.e[i]; k++)
{
if (ppf.n > UINT64_MAX / ppf.p[i])
fatal_error("Value overflow at factor %"PRIu64".", ppf.p[i]);
ppf.n *= ppf.p[i];
}
}
}
return ppf;
}
//------------------------------------------------------------------------------
// Main control. Parse command line and produce list of factors.
int main(const int argc, const char *const argv[])
{
PrimePowerFactorization ppf;
uint values = (uint)argc - 1; // argc is always guaranteed to be at least 1.
if (values == 1)
{
ppf = compute_ppf((uint64)strtoumax(argv[1], NULL, 10));
}
else
{
if (values % 2 != 0)
fatal_error("Odd number of arguments (%u) given.", values);
uint pairs = values / 2;
ppf = parse_ppf(pairs, &argv[1]);
}
// Run for (as close as possible to) a fixed amount of time, tallying the
// elapsed CPU time.
uint64 iterations = 0;
double cpu_time = 0.0;
const double cpu_time_limit = 0.05;
while (cpu_time < cpu_time_limit)
{
clock_t clock_start = clock();
uint64 *factors = compute_factors(ppf);
clock_t clock_end = clock();
cpu_time += (double)(clock_end - clock_start) / (double)CLOCKS_PER_SEC;
if (++iterations == 1)
{
for (uint32 i = 0; i < ppf.d; i++)
printf("%"PRIu64"\n", factors[i]);
}
if (factors) free(factors);
}
// Print the average amount of CPU time required for each iteration.
uint mem_scale = (memory_usage >= 1e9)? 9:
(memory_usage >= 1e6)? 6:
(memory_usage >= 1e3)? 3:
0;
char *mem_units = (mem_scale == 9)? "GB":
(mem_scale == 6)? "MB":
(mem_scale == 3)? "KB":
"B";
printf("%"PRIu64" %"PRIu32" factors %.6f ms %.3f ns/factor %.3f %s\n",
ppf.n,
ppf.d,
cpu_time/iterations * 1e3,
cpu_time/iterations * 1e9 / (double)(ppf.d? ppf.d : 1),
(double)memory_usage / pow(10, mem_scale),
mem_units);
return 0;
}
这个答案给出了一个 C 实现，我认为它是最快和最节省内存的。
算法概述。该算法基于 Will Ness 在 中引入的自下而上合并方法，但进一步简化了列表合并实际上并不存在于内存中的任何地方。每个列表的头元素都经过整理并保存在一个小数组中，而列表的所有其他元素都是根据需要即时构建的。这种“幻影列表”的使用——运行 代码的想象力的虚构——极大地减少了内存占用，以及内存访问量，包括读取和写入，并且还改善了空间局部性，这反过来又显着提高了算法的速度。每个级别的因素都按顺序直接写入输出数组中的最终位置。
基本思想是使用数学归纳法对质数幂因式分解产生因子。例如：
- 要生成 360 的因数，需要计算 72 的因数，然后乘以 5 的相关幂，在本例中为 {1,5}。
- 为了产生 72 的因数，先计算 8 的因数，然后乘以 3 的相关幂，在本例中为 {1,3,9}。
- 为了产生 8 的因数，基本情况 1 乘以 2 的相关幂，在本例中为 {1,2,4,8}。
因此，在每个归纳步骤中，都会在现有因子集和素数幂集之间计算笛卡尔积，并通过乘法将结果还原为整数。
下面是数字 360 的插图。从左到右显示的是存储单元；一排代表一个时间步长。时间被表示为垂直向下流动。
空间复杂度。产生因子的临时数据结构极小：仅O(log₂(n))space 被使用，其中 n 是正在生成其因子的数字。例如，在该算法的 128 位 C 实现中，诸如 333,939,014,887,358,848,058,068,063,658,770,598,400（其以 2 为底的对数≈127.97）之类的数字分配 5.1 GB 来存储其 318,504,960 = 0 [但仅使用 4 = 136 ]bytes 的额外开销来生成该列表。任何 128 位数字最多需要大约 5 KB 的开销。
运行时复杂度。运行时完全取决于素数幂分解的指数（例如，素数签名）。为了获得最佳结果，应首先处理最大的指数，最后处理最小的指数，以便运行时在最后阶段由低复杂度合并主导，这在实践中通常是二进制合并。渐近运行时间为 O(c(n) d(n))，其中 d(n) 是 n 的除数，其中 c(n) 是一些依赖于 n 的素数签名的常数。也就是说，每个与素数签名相关联的等价 class 都有不同的常数。由小指数主导的素数签名具有更小的常数；由大指数支配的素数签名具有更大的常数。因此，运行时复杂度由素数签名聚集。
图表。 运行时性能在 3.4 GHz 上进行了分析。英特尔酷睿 i7 对 n 的 66,591 个值进行采样，具有 d(n) 个独特 d(n) 最多1.6亿。此图分析的 n 的最大值为 14,550,525,518,294,259,162,294,162,737,840,640,000，在 2.35 秒内产生 159,744,000 个因子。
每秒产生的排序因子数本质上是上面的倒数。按素数签名的聚类在数据中很明显。性能还受 L1、L2 和 L3 缓存大小以及物理 RAM 限制的影响。
源代码。 下面附上一个 C 编程语言（特别是 C11）的工作程序。它已经在 x86-64 上使用 Clang/LLVM 进行了测试，尽管它在 GCC 上也应该可以正常工作。
/*==============================================================================
DESCRIPTION
This is a small proof-of-concept program to test the idea of generating the
factors of a number in ascending order using an ultra-efficient sortless
method.
INPUT
Input is given on the command line, either as a single argument giving the
number to be factored or an even number of arguments giving the 2-tuples that
comprise the prime-power factorization of the desired number. For example,
the number
75600 = 2^4 x 3^3 x 5^2 x 7
can be given by the following list of arguments:
2 4 3 3 5 2 7 1
Note: If a single number is given, it will require factoring to produce its
prime-power factorization. Since this is just a small test program, a very
crude factoring method is used that is extremely fast for small prime factors
but extremely slow for large prime factors. This is actually fine, because
the largest factor lists occur with small prime factors anyway, and it is the
production of large factor lists at which this program aims to be proficient.
It is simply not interesting to be fast at producing the factor list of a
number like 17293823921105882610 = 2 x 3 x 5 x 576460797370196087, because
it has only 32 factors. Numbers with tens or hundreds of thousands of
factors are much more interesting.
OUTPUT
Results are written to standard output. A list of factors in ascending order
is produced, followed by runtime required to generate the list (not including
time to print it).
AUTHOR
Todd Lehman
2015/05/10
*/
//-----------------------------------------------------------------------------
#include <inttypes.h>
#include <limits.h>
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#include <math.h>
#include <assert.h>
//-----------------------------------------------------------------------------
typedef unsigned int uint;
typedef uint8_t uint8;
typedef uint16_t uint16;
typedef uint32_t uint32;
typedef uint64_t uint64;
typedef __uint128_t uint128;
#define UINT128_MAX (uint128)(-1)
#define UINT128_MAX_STRLEN 39
//-----------------------------------------------------------------------------
#define ARRAY_CAPACITY(x) (sizeof(x) / sizeof((x)[0]))
//-----------------------------------------------------------------------------
// This structure encode a single prime-power pair for the prime-power
// factorization of numbers, for example 3 to the 4th power.
#pragma pack(push, 8)
typedef struct
{
uint128 p; // Prime.
uint8 e; // Power (exponent).
}
PrimePower; // 24 bytes using 8-byte packing
#pragma pack(pop)
//-----------------------------------------------------------------------------
// Prime-power factorization structure.
//
// This structure is sufficient to represent the prime-power factorization of
// all 128-bit values. The field names ω and Ω are dervied from the standard
// number theory functions ω(n) and Ω(n), which count the number of unique and
// non-unique prime factors of n, respectively. The field name d is derived
// from the standard number theory function d(n), which counts the number of
// divisors of n, including 1 and n.
//
// The maximum possible value here of ω is 26, which occurs at
// n = 232862364358497360900063316880507363070 = 2 x 3 x 5 x 7 x 11 x 13 x 17 x
// 19 x 23 x 29 x 31 x 37 x 41 x 43 x 47 x 53 x 59 x 61 x 67 x 71 x 73 x 79 x
// 83 x 89 x 97 x 101, which has 26 unique prime factors.
//
// The maximum possible value of Ω here is 127, which occurs at n = 2^127 and
// n = 2^126 x 3, both of which have 127 non-unique prime factors.
//
// The maximum possible value of d here is 318504960, which occurs at
// n = 333939014887358848058068063658770598400 = 2^9 x 3^5 x 5^2 x 7^2 x 11^2 x
// 13^2 x 17 x 19 x 23 x 29 x 31 x 37 x 41 x 43 x 47 x 53 x 59 x 61 x 67 x 71 x
// 73 x 79.
//
#pragma pack(push, 8)
typedef struct
{
PrimePower f[32]; // Primes and their exponents.
uint8 ω; // Count of prime factors without multiplicity.
uint8 Ω; // Count of prime factors with multiplicity.
uint32 d; // Count of factors of n, including 1 and n.
uint128 n; // Value of n on which all other fields depend.
}
PrimePowerFactorization; // 656 bytes using 8-byte packing
#pragma pack(pop)
#define MAX_ω 26
#define MAX_Ω 127
//-----------------------------------------------------------------------------
// Fatal error: print error message and abort.
void fatal_error(const char *format, ...)
{
va_list args;
va_start(args, format);
vfprintf(stderr, format, args);
exit(1);
}
//------------------------------------------------------------------------------
uint128 uint128_from_string(const char *const str)
{
assert(str != NULL);
uint128 n = 0;
for (int i = 0; isdigit(str[i]); i++)
n = (n * 10) + (uint)(str[i] - '0');
return n;
}
//------------------------------------------------------------------------------
void uint128_to_string(const uint128 n,
char *const strbuf, const uint strbuflen)
{
assert(strbuf != NULL);
assert(strbuflen >= UINT128_MAX_STRLEN + 1);
// Extract digits into string buffer in reverse order.
uint128 a = n;
char *s = strbuf;
do { *(s++) = '0' + (uint)(a % 10); a /= 10; } while (a != 0);
*s = '[=10=]';
// Reverse the order of the digits.
uint l = strlen(strbuf);
for (uint i = 0; i < l/2; i++)
{ char t = strbuf[i]; strbuf[i] = strbuf[l-1-i]; strbuf[l-1-i] = t; }
// Verify result.
assert(uint128_from_string(strbuf) == n);
}
//------------------------------------------------------------------------------
char *uint128_to_static_string(const uint128 n, const uint i)
{
static char str[2][UINT128_MAX_STRLEN + 1];
assert(i < ARRAY_CAPACITY(str));
uint128_to_string(n, str[i], ARRAY_CAPACITY(str[i]));
return str[i];
}
//------------------------------------------------------------------------------
// Compute sorted list of factors, given a prime-power factorization.
uint128 *compute_factors(const PrimePowerFactorization ppf)
{
const uint128 n = ppf.n;
const uint d = (uint)ppf.d;
const uint ω = (uint)ppf.ω;
if (n == 0)
return NULL;
uint128 *factors = malloc((d + 1) * sizeof(*factors));
if (!factors)
fatal_error("Failed to allocate array of %u factors.", d);
uint128 *const factors_end = &factors[d];
// --- Seed the factors[] array.
factors_end[0] = 0; // Dummy value to simplify looping in bottleneck code.
factors_end[-1] = 1; // Seed value.
if (n == 1)
return factors;
// --- Iterate over all prime factors.
uint range = 1;
for (uint i = 0; i < ω; i++)
{
const uint128 p = ppf.f[i].p;
const uint e = ppf.f[i].e;
// --- Initialize phantom input lists and output list.
assert(e < 128);
assert(range < d);
uint128 *restrict in[128];
uint128 pe[128], f[128];
for (uint j = 0; j <= e; j++)
{
in[j] = &factors[d - range];
pe[j] = (j == 0)? 1 : pe[j-1] * p;
f[j] = pe[j];
}
uint active_list_count = 1 + e;
range *= 1 + e;
uint128 *restrict out = &factors[d - range];
// --- Merge phantom input lists to output list, until all input lists are
// extinguished.
while (active_list_count > 0)
{
if (active_list_count == 1)
{
assert(out == in[0]);
while (out != factors_end)
*(out++) *= pe[0];
in[0] = out;
}
else if (active_list_count == 2)
{
// This section of the code is the bottleneck of the entire factor-
// producing algorithm. Other portions need to be fast, but this
// *really* needs to be fast; therefore, it has been highly optimized.
// In fact, it is by far most frequently the case here that pe[0] is 1,
// so further optimization is warranted in this case.
uint128 f0 = f[0], f1 = f[1];
uint128 *in0 = in[0], *in1 = in[1];
const uint128 pe0 = pe[0], pe1 = pe[1];
if (pe[0] == 1)
{
while (true)
{
if (f0 < f1)
{ *(out++) = f0; f0 = *(++in0);
if (in0 == factors_end) break; }
else
{ *(out++) = f1; f1 = *(++in1) * pe1; }
}
}
else
{
while (true)
{
if (f0 < f1)
{ *(out++) = f0; f0 = *(++in0) * pe0;
if (in0 == factors_end) break; }
else
{ *(out++) = f1; f1 = *(++in1) * pe1; }
}
}
f[0] = f0; f[1] = f1;
in[0] = in0; in[1] = in1;
}
else if (active_list_count == 3)
{
uint128 f0 = f[0], f1 = f[1], f2 = f[2];
uint128 *in0 = in[0], *in1 = in[1], *in2 = in[2];
const uint128 pe0 = pe[0], pe1 = pe[1], pe2 = pe[2];
while (true)
{
if (f0 < f1)
{
if (f0 < f2)
{ *(out++) = f0; f0 = *(++in0) * pe0;
if (in0 == factors_end) break; }
else
{ *(out++) = f2; f2 = *(++in2) * pe2; }
}
else
{
if (f1 < f2)
{ *(out++) = f1; f1 = *(++in1) * pe1; }
else
{ *(out++) = f2; f2 = *(++in2) * pe2; }
}
}
f[0] = f0; f[1] = f1, f[2] = f2;
in[0] = in0; in[1] = in1, in[2] = in2;
}
else if (active_list_count >= 3)
{
while (true)
{
// Chose the smallest multiplier.
uint k_min = 0;
uint128 f_min = f[0];
for (uint k = 0; k < active_list_count; k++)
if (f[k] < f_min)
{ f_min = f[k]; k_min = k; }
// Write the output factor, advance the input pointer, and
// produce a new factor in the array f[] of list heads.
*(out++) = f_min;
f[k_min] = *(++in[k_min]) * pe[k_min];
if (in[k_min] == factors_end)
{ assert(k_min == 0); break; }
}
}
// --- Remove the newly emptied phantom input list. Note that this is
// guaranteed *always* to be the first remaining non-empty list.
assert(in[0] == factors_end);
for (uint j = 1; j < active_list_count; j++)
{
in[j-1] = in[j];
pe[j-1] = pe[j];
f[j-1] = f[j];
}
active_list_count -= 1;
}
assert(out == factors_end);
}
// --- Validate array of sorted factors.
#ifndef NDEBUG
{
for (uint k = 0; k < d; k++)
{
if (factors[k] == 0)
fatal_error("Produced a factor of 0 at index %u.", k);
if (n % factors[k] != 0)
fatal_error("Produced non-factor %s at index %u.",
uint128_to_static_string(factors[k], 0), k);
if ((k > 0) && (factors[k-1] == factors[k]))
fatal_error("Duplicate factor %s at index %u.",
uint128_to_static_string(factors[k], 0), k);
if ((k > 0) && (factors[k-1] > factors[k]))
fatal_error("Out-of-order factors %s and %s at indexes %u and %u.",
uint128_to_static_string(factors[k-1], 0),
uint128_to_static_string(factors[k], 1),
k-1, k);
}
}
#endif
return factors;
}
//------------------------------------------------------------------------------
// Print prime-power factorization of a number.
void print_ppf(const PrimePowerFactorization ppf)
{
printf("%s = ", uint128_to_static_string(ppf.n, 0));
if (ppf.n == 0)
{
printf("0");
}
else
{
for (uint i = 0; i < ppf.ω; i++)
{
if (i > 0)
printf(" x ");
printf("%s", uint128_to_static_string(ppf.f[i].p, 0));
if (ppf.f[i].e > 1)
printf("^%"PRIu8"", ppf.f[i].e);
}
}
printf("\n");
}
//------------------------------------------------------------------------------
int compare_powers_ascending(const void *const pf1,
const void *const pf2)
{
const PrimePower f1 = *((const PrimePower *)pf1);
const PrimePower f2 = *((const PrimePower *)pf2);
return (f1.e < f2.e)? -1:
(f1.e > f2.e)? +1:
0; // Not an error; duplicate exponents are common.
}
//------------------------------------------------------------------------------
int compare_powers_descending(const void *const pf1,
const void *const pf2)
{
const PrimePower f1 = *((const PrimePower *)pf1);
const PrimePower f2 = *((const PrimePower *)pf2);
return (f1.e < f2.e)? +1:
(f1.e > f2.e)? -1:
0; // Not an error; duplicate exponents are common.
}
//------------------------------------------------------------------------------
int compare_primes_ascending(const void *const pf1,
const void *const pf2)
{
const PrimePower f1 = *((const PrimePower *)pf1);
const PrimePower f2 = *((const PrimePower *)pf2);
return (f1.p < f2.p)? -1:
(f1.p > f2.p)? +1:
0; // Error; duplicate primes must never occur.
}
//------------------------------------------------------------------------------
int compare_primes_descending(const void *const pf1,
const void *const pf2)
{
const PrimePower f1 = *((const PrimePower *)pf1);
const PrimePower f2 = *((const PrimePower *)pf2);
return (f1.p < f2.p)? +1:
(f1.p > f2.p)? -1:
0; // Error; duplicate primes must never occur.
}
//------------------------------------------------------------------------------
// Sort prime-power factorization.
void sort_ppf(PrimePowerFactorization *const ppf,
const bool primes_major, // Best false
const bool primes_ascending, // Best false
const bool powers_ascending) // Best false
{
int (*compare_primes)(const void *, const void *) =
primes_ascending? compare_primes_ascending : compare_primes_descending;
int (*compare_powers)(const void *, const void *) =
powers_ascending? compare_powers_ascending : compare_powers_descending;
if (primes_major)
{
mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_powers);
mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_primes);
}
else
{
mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_primes);
mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_powers);
}
}
//------------------------------------------------------------------------------
// Compute prime-power factorization of a 128-bit value. Note that this
// function is designed to be fast *only* for numbers with very simple
// factorizations, e.g., those that produce large factor lists. Do not attempt
// to factor large semiprimes with this function. (The author does know how to
// factor large numbers efficiently; however, efficient factorization is beyond
// the scope of this small test program.)
PrimePowerFactorization compute_ppf(const uint128 n)
{
PrimePowerFactorization ppf;
if (n == 0)
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 };
}
else if (n == 1)
{
ppf = (PrimePowerFactorization){ .f[0] = { .p = 1, .e = 1 },
.ω = 1, .Ω = 1, .d = 1, .n = 1 };
}
else
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = n };
uint128 m = n;
for (uint128 p = 2; p * p <= m; p += 1 + (p > 2))
{
if (m % p == 0)
{
assert(ppf.ω <= MAX_ω);
ppf.f[ppf.ω].p = p;
ppf.f[ppf.ω].e = 0;
while (m % p == 0)
{ m /= p; ppf.f[ppf.ω].e += 1; }
ppf.d *= (1 + ppf.f[ppf.ω].e);
ppf.Ω += ppf.f[ppf.ω].e;
ppf.ω += 1;
}
}
if (m > 1)
{
assert(ppf.ω <= MAX_ω);
ppf.f[ppf.ω].p = m;
ppf.f[ppf.ω].e = 1;
ppf.d *= 2;
ppf.Ω += 1;
ppf.ω += 1;
}
}
return ppf;
}
//------------------------------------------------------------------------------
// Parse prime-power factorization from a list of ASCII-encoded base-10 strings.
// The values are assumed to be 2-tuples (p,e) of prime p and exponent e.
// Primes must not exceed 2^128 - 1 and must not be repeated. Exponents must
// not exceed 2^8 - 1, but can of course be repeated. The constructed value
// must not exceed 2^128 - 1.
PrimePowerFactorization parse_ppf(const uint pairs, const char *const values[])
{
assert(pairs <= MAX_ω);
PrimePowerFactorization ppf;
if (pairs == 0)
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 };
}
else
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = 1 };
for (uint i = 0; i < pairs; i++)
{
ppf.f[i].p = uint128_from_string(values[(i*2)+0]);
ppf.f[i].e = (uint8)strtoumax(values[(i*2)+1], NULL, 10);
// Validate prime value.
if (ppf.f[i].p < 2) // (Ideally this would actually do a primality test.)
fatal_error("Factor %s is invalid.",
uint128_to_static_string(ppf.f[i].p, 0));
// Accumulate count of unique prime factors.
if (ppf.ω > UINT8_MAX - 1)
fatal_error("Small-omega overflow at factor %s^%"PRIu8".",
uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e);
ppf.ω += 1;
// Accumulate count of total prime factors.
if (ppf.Ω > UINT8_MAX - ppf.f[i].e)
fatal_error("Big-omega wverflow at factor %s^%"PRIu8".",
uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e);
ppf.Ω += ppf.f[i].e;
// Accumulate total divisor count.
if (ppf.d > UINT32_MAX / (1 + ppf.f[i].e))
fatal_error("Divisor count overflow at factor %s^%"PRIu8".",
uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e);
ppf.d *= (1 + ppf.f[i].e);
// Accumulate value.
for (uint8 k = 1; k <= ppf.f[i].e; k++)
{
if (ppf.n > UINT128_MAX / ppf.f[i].p)
fatal_error("Value overflow at factor %s.",
uint128_to_static_string(ppf.f[i].p, 0));
ppf.n *= ppf.f[i].p;
}
}
}
return ppf;
}
//------------------------------------------------------------------------------
// Main control. Parse command line and produce list of factors.
int main(const int argc, const char *const argv[])
{
bool primes_major = false;
bool primes_ascending = false;
bool powers_ascending = false;
PrimePowerFactorization ppf;
// --- Parse prime-power sort specifier (if present).
uint value_base = 1;
uint value_count = (uint)argc - 1;
if ((argc > 1) && (argv[1][0] == '-'))
{
static const struct
{
char *str; bool primes_major, primes_ascending, powers_ascending;
}
sort_options[] =
{
// Sorting criteria:
// ----------------------------------------
{ "ep", 0,0,0 }, // Exponents descending, primes descending
{ "Ep", 0,0,1 }, // Exponents ascending, primes descending
{ "eP", 0,1,0 }, // Exponents descending, primes ascending
{ "EP", 0,1,1 }, // Exponents ascending, primes ascending
{ "p", 1,0,0 }, // Primes descending (exponents irrelevant)
{ "P", 1,1,0 }, // Primes ascending (exponents irrelevant)
};
bool valid = false;
for (uint i = 0; i < ARRAY_CAPACITY(sort_options); i++)
{
if (strcmp(&argv[1][1], sort_options[i].str) == 0)
{
primes_major = sort_options[i].primes_major;
primes_ascending = sort_options[i].primes_ascending;
powers_ascending = sort_options[i].powers_ascending;
valid = true;
break;
}
}
if (!valid)
fatal_error("Bad sort specifier: \"%s\"", argv[1]);
value_base += 1;
value_count -= 1;
}
// --- Prime factorization from either a number or a raw prime factorization.
if (value_count == 1)
{
uint128 n = uint128_from_string(argv[value_base]);
ppf = compute_ppf(n);
}
else
{
if (value_count % 2 != 0)
fatal_error("Odd number of arguments (%u) given.", value_count);
uint pairs = value_count / 2;
ppf = parse_ppf(pairs, &argv[value_base]);
}
// --- Sort prime factorization by either the default or the user-overridden
// configuration.
sort_ppf(&ppf, primes_major, primes_ascending, powers_ascending);
print_ppf(ppf);
// --- Run for (as close as possible to) a fixed amount of time, tallying the
// elapsed CPU time.
uint128 iterations = 0;
double cpu_time = 0.0;
const double cpu_time_limit = 0.10;
uint128 memory_usage = 0;
while (cpu_time < cpu_time_limit)
{
clock_t clock_start = clock();
uint128 *factors = compute_factors(ppf);
clock_t clock_end = clock();
cpu_time += (double)(clock_end - clock_start) / (double)CLOCKS_PER_SEC;
memory_usage = sizeof(*factors) * ppf.d;
if (++iterations == 0) //1)
{
for (uint32 i = 0; i < ppf.d; i++)
printf("%s\n", uint128_to_static_string(factors[i], 0));
}
if (factors) free(factors);
}
// --- Print the average amount of CPU time required for each iteration.
uint memory_scale = (memory_usage >= 1e9)? 9:
(memory_usage >= 1e6)? 6:
(memory_usage >= 1e3)? 3:
0;
char *memory_units = (memory_scale == 9)? "GB":
(memory_scale == 6)? "MB":
(memory_scale == 3)? "KB":
"B";
printf("%s %"PRIu32" factors %.6f ms %.3f ns/factor %.3f %s\n",
uint128_to_static_string(ppf.n, 0),
ppf.d,
cpu_time/iterations * 1e3,
cpu_time/iterations * 1e9 / (double)(ppf.d? ppf.d : 1),
(double)memory_usage / pow(10, memory_scale),
memory_units);
return 0;
}
有没有一种高效的算法来枚举一个数字n的因数，按升序排列，不排序？我所说的“高效”是指：
该算法从 n 的素数幂分解开始，避免了对除数的强力搜索。
算法的运行时复杂度为 O(d log₂ d) 或更好，其中 d 是 n 的除数。
算法的空间复杂度为O(d).
该算法避免了排序操作。也就是说，因子是按顺序 生成的，而不是乱序生成然后再排序的。尽管使用简单的递归方法枚举然后排序是 O(d log₂ d)，但所有涉及的内存访问的成本非常低在排序中。
一个简单的例子是 n = 360 = 2³ × 3² × 5，它有 d=24 个因子：{ 1, 2 , 3, 4, 5, 6, 8, 9, 10, 12, 15, 18, 20, 24, 30, 36, 40, 45, 60, 72, 90, 120, 180, 360 }.
更严重的例子是n = 278282512406132373381723386382308832000 = 2⁸ × 3⁴ × 5³ × 7² × 11² × 13² × 17 × 19 × 23 × 29 × 31 × 37 × 41 × 43 × 47 × 53 × 59 × 61 × 67 × 71 × 73 × 79，其中有 d=318504960 个因子（显然太多了，无法在这里列出！）。顺便说一下，这个数是最多 2^128 的任何数的最大因数。
我可以发誓几周前我看到了这种算法的描述，带有示例代码，但现在我似乎无法在任何地方找到它。它使用了一些魔术，在输出列表中为每个素因子维护一个祖先索引列表。 （更新：我将因子生成与操作类似的汉明数混淆了。）
更新
我最终在运行时使用了 O(d) 的解决方案，创建 O(d 的内存开销极低) 就地输出，并且比我知道的任何其他方法都快得多。我已经发布了这个解决方案作为答案，带有 C 源代码。它是另一位贡献者 Will Ness 在 Haskell 中介绍的漂亮算法的高度优化简化版本。我选择 Will 的回答作为接受的答案，因为它提供了一个非常优雅的解决方案，符合最初陈述的所有要求。
简而言之：重复从堆中拉出 next-smallest 因子，然后推回该因子的每个倍数，该倍数仍然是 n 的因子。使用技巧来避免出现重复项，以便堆大小永远不会超过 d。时间复杂度为 O(kd log d)，其中 k 是不同质因数的数量。
我们使用的关键 属性 是如果 x 和 y 都是 n 的因数且 y = x*p 对于某个因数 p >= 2 —— 即如果 x 的质因数是 y 的素因子的真子集——那么 x < y。这意味着在将任何此类 y 插入堆之前输出 x 总是安全的。堆仅有效地用于比较两个都不是另一个子多重集的因子。
第一次尝试：重复会减慢速度
首先描述一个产生正确答案但也产生许多重复的算法会很有帮助：
- 设置 prev = NULL。
- 将 1 插入堆 H.
- 从H中取出堆顶t，如果堆为空，则停止。
- 如果 t == prev 然后转到 3。[编辑：固定]
- 输出 t.
- 设置 prev = t.
- 对于 n 的每个不同的质因数 p：
- 如果 n % (t*p) == 0（即如果 t*p 仍然是 n 的一个因子），将 t*p 推到 H 上。
- 转到 3.
上述算法唯一的问题是它可以多次生成相同的因子。例如，如果 n = 30，则因子 15 将生成为因子 5 的 child（通过乘以质因子 3），以及因子 3 的 child（乘以 5)。解决这个问题的一种方法是注意任何重复项到达堆顶时必须在连续块中读出，因此您可以简单地检查堆顶是否等于 just-extracted 值，并且如果是这样，继续提取和丢弃它。但可能有更好的方法：
从源头上消除重复
因子 x 有多少种生成方式？首先考虑 x 不包含重数 > 1 的素因子的情况。在这种情况下，如果它包含 m 个不同的素因子，则有 m-1 "parent" 个因子将其生成为 "child" 在前面的算法中——这些 parent 中的每一个都由 m 个质因数中的 m-1 个子集组成，其余质因数是添加到 child 的质因数。 （如果 x 有一个重数 > 1 的素因子，那么实际上有 m parents。）如果我们有办法确定这些 parents 中的一个恰好是 "chosen one" 实际上将 x 生成为 child，并且此规则导致可以在 parent 弹出时应用于每个 parent 的测试，那么我们就可以避免创建任何重复项。
我们可以使用以下规则：对于任何给定的 x，选择缺少 最大 x 的 m 个因子的潜在 parent y。这形成了一个简单的规则：A parent y 产生一个 child x 当且仅当 x = y*p 对于某个 p 大于或等于 y 中已有的任何质因数。这很容易测试：只需按降序遍历质因数，为每个质因数生成 children，直到我们找到一个已经整除 y 的质因数。在前面的示例中，parent 3 将产生 15，但 parent 5 不会（因为 3 < 5）——因此 15 确实只产生一次。对于 n = 30，完整的树如下所示：
1
/|\
2 3 5
/| \
6 10 15
|
30
注意每个因子只生成一次。
新的duplicate-free算法如下：
- 将 1 插入堆 H.
- 从H中取出堆顶t，如果堆为空，则停止。
- 输出 t.
- 对于 n 的每个不同的素因子 p 按降序排列：
- 如果 n % (t*p) == 0（即如果 t*p 仍然是 n 的一个因子），将 t*p 推到 H 上。
- 如果 t % p == 0（即如果 t 已经包含 p 作为一个因子）则停止。
- 转到 2.
我们可以合并多个流，这样一来就没有重复项了。
从列表 [1] 开始，对于每个唯一的素因子 p，我们将列表乘以 p 迭代 k 次（其中 k是 p) 的多重性以获得 k 个新列表，并将它们与传入列表合并在一起。
在Haskell,
ordfactors n = -- <----------------------<---<---<-----\
foldr g [1] -- g (19,1) (g (7,1) (g (3,2) (g (2,3) [1])))
. reverse -- [(19,1),(7,1),(3,2),(2,3)] ^
. primePowers -- [(2,3),(3,2),(7,1),(19,1)] |
$ n -- 9576 --->--->--->----/
where
g (p,k) xs = mergeAsTree
. take (k+1) -- take 3 [a,b,c,d..] = [a,b,c]
. iterate (map (p*)) -- iterate f x = [x, y, z,...]
$ xs -- where { y=f x; z=f y; ... }
{- g (2,3) [1] = let a0 = [1]
a1 = map (2*) a0 -- [2]
a2 = map (2*) a1 -- [4]
a3 = map (2*) a2 -- [8]
in mergeAsTree [ a0, a1, a2, a3 ] -- xs2 = [1,2,4,8]
g (3,2) xs2 = let b0 = xs2 -- [1,2,4,8]
b1 = map (3*) b0 -- [3,6,12,24]
b2 = map (3*) b1 -- [9,18,36,72]
in mergeAsTree [ b0, b1, b2 ] -- xs3 = [1,2,3,4,6,8,9,12,...]
g (7,1) xs3 = mergeAsTree [ xs3, map (7*) xs3 ] -- xs4 = [1,2,3,4,6,7,8,9,...]
g (19,1) xs4 = mergeAsTree [ xs4, map (19*) xs4 ]
-}
mergeAsTree xs = foldi merge [] xs -- [a,b] --> merge a b
-- [a,b,c,d] --> merge (merge a b) (merge c d)
foldi f z [] = z
foldi f z (x:xs) = f x (foldi f z (pairs f xs))
pairs f (x:y:t) = f x y : pairs f t
pairs _ t = t
foldi 以线性方式排列二进制 merges (which presuppose that there's no duplicates in the streams being merged, for streamlined operation) in a tree-like fashion for speed; whereas foldr。
primePowers n | n > 1 = -- 9576 => [(2,3),(3,2),(7,1),(19,1)]
map (\xs -> (head xs,length xs)) -- ^
. group -- [[2,2,2],[3,3],[7],[19]] |
. factorize -- [2,2,2,3,3,7,19] |
$ n -- 9576 --->--->--->---/
group 是一个标准函数，它将列表中的相邻等值组合在一起，而 factorize 是一个众所周知的试分算法，它以非递减顺序产生一个数的质因数:
factorize n | n > 1 = go n (2:[3,5..]) -- 9576 = 2*2*2*3*3*7*19
where -- produce prime factors only
go n (d:t)
| d*d > n = [n]
| r == 0 = d : go q (d:t)
| otherwise = go n t
where
(q,r) = quotRem n d
factorize _ = []
(.) 是一个函数组合运算符，将其右侧的参数函数的输出链接为左侧函数的输入。它（和 $）可以读作“of”。
[我发布这个答案只是为了形式上的完整性。我已经选择了别人的答案作为接受的答案。]
算法概述。 在寻找生成 in-memory 因子列表（在我的例子中是 64 位无符号值）的最快方法时，我选择了一种实现 two-dimensional 桶排序的混合算法，它利用了排序键的内部知识（即，它们只是整数，因此可以对其进行计算）。具体方法更接近于“ProxMapSort”，但具有两个级别的键（主要，次要）而不是只有一个。 主键 只是该值的以 2 为底的对数。 minor key 是在第二层桶中产生合理传播所需的值的最高有效数字的最小数量。因子首先生成为未分类因子的临时工作数组。接下来，分析它们的分布并分配和填充桶索引数组。最后，使用插入排序将因子直接存储到最终排序的数组中。绝大多数桶只有 1、2 或 3 个值。示例在源代码中给出，附在这个答案的底部。
空间复杂度。 内存利用率大约是 Quicksort-based 解决方案的 4 倍，但这实际上是微不足道的，因为在最坏情况下使用过的最大内存（对于 64 位输入）是 5.5 MB，其中 4.0 MB 只保留了几毫秒。
运行时复杂度。 性能远好于 hand-coded Quicksort-based 解决方案：对于具有非平凡数量的因子的数字，它不规则地约为 2.5快 x 倍。在我的 3.4 GHz 上。 Intel i7，它在0.0052秒内按排序顺序产生184,320个因子18,401,055,938,125,660,800，或每个因子约96个时钟周期，或每秒约3500万个因子。
图表。 内存和运行时性能分析了 47,616 个等价的规范代表 类 最多 2⁶⁴–1 的素数签名。这些是 64 位搜索中的 so-called“高度可分解的数字”space。
对于非平凡的因子计数，总运行时间比 Quicksort-based 解决方案好 ~2.5 倍，如下图所示：
每秒产生的排序因子数本质上是上面的倒数。在大约 2000 个因子的最佳点之后，per-factor 的性能下降，但下降幅度不大。性能受 L1、L2 和 L3 缓存大小以及被分解数的唯一素数的影响，大致随输入值的对数增加。
在此对数-对数图中，内存使用峰值是一条直线，因为它与因子数量的以 2 为底的对数成正比。请注意，内存使用峰值仅在很短的时间内出现； short-lived 工作数组在几毫秒内被丢弃。丢弃临时数组后，剩下的是最终的因子列表，这与 Quicksort-based 解决方案中看到的最小用法相同。
源代码。 下面附上一个 proof-of-concept C 编程语言（特别是 C11）程序。它已经在 x86-64 上使用 Clang/LLVM 进行了测试，尽管它在 GCC 上也应该可以正常工作。
/*==============================================================================
DESCRIPTION
This is a small proof-of-concept program to test the idea of "sorting"
factors using a form of bucket sort. The method is essentially a 2D version
of ProxMapSort that has tuned for vast, nonlinear distributions using two
keys (major, minor) rather than one. The major key is simply the floor of
the base-2 logarithm of the value, and the minor key is derived from the most
significant bits of the value.
INPUT
Input is given on the command line, either as a single argument giving the
number to be factored or an even number of arguments giving the 2-tuples that
comprise the prime-power factorization of the desired number. For example,
the number
75600 = 2^4 x 3^3 x 5^2 x 7
can be given by the following list of arguments:
2 4 3 3 5 2 7 1
Note: If a single number is given, it will require factoring to produce its
prime-power factorization. Since this is just a small test program, a very
crude factoring method is used that is extremely fast for small prime factors
but extremely slow for large prime factors. This is actually fine, because
the largest factor lists occur with small prime factors anyway, and it is the
production of large factor lists at which this program aims to be proficient.
It is simply not interesting to be fast at producing the factor list of a
number like 17293823921105882610 = 2 x 3 x 5 x 576460797370196087, because
it has only 32 factors. Numbers with tens or hundreds of thousands of
factors are much more interesting.
OUTPUT
Results are written to standard output. A list of factors in ascending order
is produced, followed by runtime (in microseconds) required to generate the
list (not including time to print it).
STATISTICS
Bucket size statistics for the 47616 canonical representatives of the prime
signature equivalence classes of 64-bit numbers:
==============================================================
Bucket size Total count of factored Total count of
b numbers needing size b buckets of size b
--------------------------------------------------------------
1 47616 (100.0%) 514306458 (76.2%)
2 47427 (99.6%) 142959971 (21.2%)
3 43956 (92.3%) 16679329 (2.5%)
4 27998 (58.8%) 995458 (0.1%)
5 6536 (13.7%) 33427 (<0.1%)
6 400 (0.8%) 729 (<0.1%)
7 12 (<0.1%) 18 (<0.1%)
--------------------------------------------------------------
~ 47616 (100.0%) 674974643 (100.0%)
--------------------------------------------------------------
Thus, no 64-bit number (of the input set) ever requires more than 7 buckets,
and the larger the bucket size the less frequent it is. This is highly
desirable. Note that although most numbers need at least 1 bucket of size 5,
the vast majority of buckets (99.9%) are of size 1, 2, or 3, meaning that
insertions are extremely efficient. Therefore, the use of insertion sort
for the buckets is clearly the right choice and is arguably optimal for
performance.
AUTHOR
Todd Lehman
2015/05/08
*/
#include <inttypes.h>
#include <limits.h>
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <assert.h>
typedef unsigned int uint;
typedef uint8_t uint8;
typedef uint16_t uint16;
typedef uint32_t uint32;
typedef uint64_t uint64;
#define ARRAY_CAPACITY(x) (sizeof(x) / sizeof((x)[0]))
//-----------------------------------------------------------------------------
// This structure is sufficient to represent the prime-power factorization of
// all 64-bit values. The field names ω and Ω are dervied from the standard
// number theory functions ω(n) and Ω(n), which count the number of unique and
// non-unique prime factors of n, respectively. The field name d is derived
// from the standard number theory function d(n), which counts the number of
// divisors of n, including 1 and n.
//
// The maximum possible value here of ω is 15, which occurs for example at
// n = 7378677391061896920 = 2^3 x 3^2 x 5 x 7 x 11 x 13 x 17 x 19 x 23 x 29
// 31 x 37 x 41 x 43 x 47, which has 15 unique prime factors.
//
// The maximum possible value of Ω here is 63, which occurs for example at
// n = 2^63 and n = 2^62 x 3, both of which have 63 non-unique prime factors.
//
// The maximum possible value of d here is 184320, which occurs at
// n = 18401055938125660800 = 2^7 x 3^4 x 5^2 x 7^2 x 11 x 13 x 17 x 19 x 23 x
// 29 x 31 x 37 x 41.
//
// Maximum possible exponents when exponents are sorted in decreasing order:
//
// Index Maximum Bits Example of n
// ----- ------- ---- --------------------------------------------
// 0 63 6 (2)^63
// 1 24 5 (2*3)^24
// 2 13 4 (2*3*5)^13
// 3 8 4 (2*3*5*7)^8
// 4 5 3 (2*3*5*7*11)^5
// 5 4 3 (2*3*5*7*11*13)^4
// 6 3 2 (2*3*5*7*11*13*17)^3
// 7 2 2 (2*3*5*7*11*13*17*19)^2
// 8 2 2 (2*3*5*7*11*13*17*19*23)^2
// 9 1 1 (2*3*5*7*11*13*17*19*23*29)^1
// 10 1 1 (2*3*5*7*11*13*17*19*23*29*31)^1
// 11 1 1 (2*3*5*7*11*13*17*19*23*29*31*37)^1
// 12 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41)^1
// 13 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41*43)^1
// 14 1 1 (2*3*5*7*11*13*17*19*23*29*31*37*41*43*47)^1
// ----- ------- ---- --------------------------------------------
// 15 63 37
//
#pragma pack(push, 8)
typedef struct
{
uint8 e[16]; // Exponents.
uint64 p[16]; // Primes in increasing order.
uint8 ω; // Count of prime factors without multiplicity.
uint8 Ω; // Count of prime factors with multiplicity.
uint32 d; // Count of factors of n, including 1 and n.
uint64 n; // Value of n on which all other fields of this struct depend.
}
PrimePowerFactorization; // 176 bytes with 8-byte packing
#pragma pack(pop)
#define MAX_ω 15
#define MAX_Ω 63
//-----------------------------------------------------------------------------
// Fatal error: print error message and abort.
void fatal_error(const char *format, ...)
{
va_list args;
va_start(args, format);
vfprintf(stderr, format, args);
exit(1);
}
//-----------------------------------------------------------------------------
// Compute 64-bit 2-adic integer inverse.
uint64 uint64_inv(const uint64 x)
{
assert(x != 0);
uint64 y = 1;
for (uint i = 0; i < 6; i++) // 6 = log2(log2(2**64)) = log2(64)
y = y * (2 - (x * y));
return y;
}
//------------------------------------------------------------------------------
// Compute 2 to arbitrary power. This is just a portable and abstract way to
// write a left-shift operation. Note that the use of the UINT64_C macro here
// is actually required, because the result of 1U<<x is not guaranteed to be a
// 64-bit result; on a 32-bit compiler, 1U<<32 is 0 or is undefined.
static inline
uint64 uint64_pow2(x)
{
return UINT64_C(1) << x;
}
//------------------------------------------------------------------------------
// Deduce native word size (int, long, or long long) for 64-bit integers.
// This is needed for abstracting certain compiler-specific intrinsic functions.
#if UINT_MAX == 0xFFFFFFFFFFFFFFFFU
#define UINT64_IS_U
#elif ULONG_MAX == 0xFFFFFFFFFFFFFFFFUL
#define UINT64_IS_UL
#elif ULLONG_MAX == 0xFFFFFFFFFFFFFFFFULL
#define UINT64_IS_ULL
#else
//error "Unable to deduce native word size of 64-bit integers."
#endif
//------------------------------------------------------------------------------
// Define abstracted intrinsic function for counting leading zeros. Note that
// the value is well-defined for nonzero input but is compiler-specific for
// input of zero.
#if defined(UINT64_IS_U) && __has_builtin(__builtin_clz)
#define UINT64_CLZ(x) __builtin_clz(x)
#elif defined(UINT64_IS_UL) && __has_builtin(__builtin_clzl)
#define UINT64_CLZ(x) __builtin_clzl(x)
#elif defined(UINT64_IS_ULL) && __has_builtin(__builtin_clzll)
#define UINT64_CLZ(x) __builtin_clzll(x)
#else
#undef UINT64_CLZ
#endif
//------------------------------------------------------------------------------
// Compute floor of base-2 logarithm y = log_2(x), where x > 0. Uses fast
// intrinsic function if available; otherwise resorts to hand-rolled method.
static inline
uint uint64_log2(uint64 x)
{
assert(x > 0);
#if defined(UINT64_CLZ)
return 63 - UINT64_CLZ(x);
#else
#define S(k) if ((x >> k) != 0) { y += k; x >>= k; }
uint y = 0; S(32); S(16); S(8); S(4); S(2); S(1); return y;
#undef S
#endif
}
//------------------------------------------------------------------------------
// Compute major key, given a nonzero number. The major key is simply the
// floor of the base-2 logarithm of the number.
static inline
uint major_key(const uint64 n)
{
assert(n > 0);
uint k1 = uint64_log2(n);
return k1;
}
//------------------------------------------------------------------------------
// Compute minor key, given a nonzero number, its major key, k1, and the
// bit-size b of major bucket k1. The minor key, k2, is is computed by first
// removing the most significant 1-bit from the number, because it adds no
// information, and then extracting the desired number of most significant bits
// from the remainder. For example, given the number n=1463 and a major bucket
// size of b=6 bits, the keys are computed as follows:
//
// Step 0: Given number n = 0b10110110111 = 1463
//
// Step 1: Compute major key: k1 = floor(log_2(n)) = 10
//
// Step 2: Remove high-order 1-bit: n' = 0b0110110111 = 439
//
// Step 3: Compute minor key: k2 = n' >> (k1 - b)
// = 0b0110110111 >> (10 - 6)
// = 0b0110110111 >> 4
// = 0b011011
// = 27
static inline
uint minor_key(const uint64 n, const uint k1, const uint b)
{
assert(n > 0); assert(k1 >= 0); assert(b > 0);
const uint k2 = (uint)((n ^ uint64_pow2(k1)) >> (k1 - b));
return k2;
}
//------------------------------------------------------------------------------
// Raw unsorted factor.
#pragma push(pack, 4)
typedef struct
{
uint64 n; // Value of factor.
uint32 k1; // Major key.
uint32 k2; // Minor key.
}
UnsortedFactor;
#pragma pop(pack)
//------------------------------------------------------------------------------
// Compute sorted list of factors, given a prime-power factorization.
static uint64 memory_usage;
uint64 *compute_factors(const PrimePowerFactorization ppf)
{
memory_usage = 0;
if (ppf.n == 0)
return NULL;
uint64 *sorted_factors = calloc(ppf.d, sizeof(*sorted_factors));
if (!sorted_factors)
fatal_error("Failed to allocate array of %"PRIu32" factors.", ppf.d);
memory_usage += ppf.d * sizeof(*sorted_factors);
UnsortedFactor *unsorted_factors = malloc(ppf.d * sizeof(*unsorted_factors));
if (!unsorted_factors)
fatal_error("Failed to allocate array of %"PRIu32" factors.", ppf.d);
memory_usage += ppf.d * sizeof(*unsorted_factors);
// These arrays are indexed by the major key of a number.
uint32 major_counts[64]; // Counts of factors in major buckets.
uint32 major_spans[64]; // Counts rounded up to power of 2.
uint32 major_bits[64]; // Base-2 logarithm of bucket size.
uint32 major_indexes[64]; // Indexes into minor array.
memset(major_counts, 0, sizeof(major_counts));
memset(major_spans, 0, sizeof(major_spans));
memset(major_bits, 0, sizeof(major_bits));
memset(major_indexes, 0, sizeof(major_indexes));
// --- Step 1: Produce unsorted list of factors from prime-power
// factorization. At the same time, count groups of factors by their
// major keys.
{
// This array is for counting in the multi-radix number system dictated by
// the exponents of the prime-power factorization. An invariant is that
// e[i] <= ppf.e[i] for all i (0 < i <ppf.ω).
uint8 e[MAX_ω];
for (uint i = 0; i < ppf.ω; i++)
e[i] = 0;
// Initialize inverse-prime-powers. This array allows for division by
// p[i]**e[i] extremely quickly in the main loop below. Note that 2-adic
// inverses are not defined for even numbers (of which 2 is the only prime),
// so powers of 2 must be handled specially.
uint64 pe_inv[MAX_ω];
for (uint i = 0; i < ppf.ω; i++)
{
uint64 pe = 1; for (uint j = 1; j <= ppf.e[i]; j++) pe *= ppf.p[i];
pe_inv[i] = uint64_inv(pe);
}
uint64 n = 1; // Current factor accumulator.
for (uint k = 0; k < ppf.d; k++) // k indexes into unsorted_factors[].
{
//printf("unsorted_factors[%u] = %"PRIu64" j = %u\n", k, n, j);
assert(ppf.n % n == 0);
unsorted_factors[k].n = n;
uint k1 = major_key(n);
assert(k1 < ARRAY_CAPACITY(major_counts));
unsorted_factors[k].k1 = k1;
major_counts[k1] += 1;
// Increment the remainder of the multi-radix number e[].
for (uint i = 0; i < ppf.ω; i++)
{
if (e[i] == ppf.e[i]) // Carrying is occurring.
{
if (ppf.p[i] == 2)
n >>= ppf.e[i]; // Divide n by 2 ** ppf.e[i].
else
n *= pe_inv[i]; // Divide n by ppf.p[i] ** ppf.e[i].
e[i] = 0;
}
else // Carrying is not occurring.
{
n *= ppf.p[i];
e[i] += 1;
break;
}
}
}
assert(n == 1); // n always cycles back to 1, not to ppf.n.
assert(unsorted_factors[ppf.d-1].n == ppf.n);
}
// --- Step 2: Define the major bits array, the major spans array, the major
// index array, and count the total spans.
uint32 total_spans = 0;
{
uint32 k = 0;
for (uint k1 = 0; k1 < ARRAY_CAPACITY(major_counts); k1++)
{
uint32 count = major_counts[k1];
uint32 bits = (count <= 1)? count : uint64_log2(count - 1) + 1;
major_bits[k1] = bits;
major_spans[k1] = (count > 0)? (UINT32_C(1) << bits) : 0;
major_indexes[k1] = k;
k += major_spans[k1];
}
total_spans = k;
}
// --- Step 3: Allocate and populate the minor counts array. Note that it
// must be initialized to zero.
uint32 *minor_counts = calloc(total_spans, sizeof(*minor_counts));
if (!minor_counts)
fatal_error("Failed to allocate array of %"PRIu32" counts.", total_spans);
memory_usage += total_spans * sizeof(*minor_counts);
for (uint k = 0; k < ppf.d; k++)
{
const uint64 n = unsorted_factors[k].n;
const uint k1 = unsorted_factors[k].k1;
const uint k2 = minor_key(n, k1, major_bits[k1]);
assert(k2 < major_spans[k1]);
unsorted_factors[k].k2 = k2;
minor_counts[major_indexes[k1] + k2] += 1;
}
// --- Step 4: Define the minor indexes array.
//
// NOTE: Instead of allocating a separate array, the earlier-allocated array
// of minor indexes is simply repurposed here using an alias.
uint32 *minor_indexes = minor_counts; // Alias the array for repurposing.
{
uint32 k = 0;
for (uint i = 0; i < total_spans; i++)
{
uint32 count = minor_counts[i]; // This array is the same array...
minor_indexes[i] = k; // ...as this array.
k += count;
}
}
// --- Step 5: Populate the sorted factors array. Note that the array must
// be initialized to zero earlier because values of zero are used
// as sentinels in the bucket lists.
for (uint32 i = 0; i < ppf.d; i++)
{
uint64 n = unsorted_factors[i].n;
const uint k1 = unsorted_factors[i].k1;
const uint k2 = unsorted_factors[i].k2;
// Insert factor into bucket using insertion sort (which happens to be
// extremely fast because we know the bucket sizes are always very small).
uint32 k;
for (k = minor_indexes[major_indexes[k1] + k2];
sorted_factors[k] != 0;
k++)
{
assert(k < ppf.d);
if (sorted_factors[k] > n)
{ uint64 t = sorted_factors[k]; sorted_factors[k] = n; n = t; }
}
sorted_factors[k] = n;
}
// --- Step 6: Validate array of sorted factors.
{
for (uint32 k = 1; k < ppf.d; k++)
{
if (sorted_factors[k] == 0)
fatal_error("Produced a factor of 0 at index %"PRIu32".", k);
if (ppf.n % sorted_factors[k] != 0)
fatal_error("Produced non-factor %"PRIu64" at index %"PRIu32".",
sorted_factors[k], k);
if (sorted_factors[k-1] == sorted_factors[k])
fatal_error("Duplicate factor %"PRIu64" at index %"PRIu32".",
sorted_factors[k], k);
if (sorted_factors[k-1] > sorted_factors[k])
fatal_error("Out-of-order factors %"PRIu64" and %"PRIu64" "
"at indexes %"PRIu32" and %"PRIu32".",
sorted_factors[k-1], sorted_factors[k], k-1, k);
}
}
free(minor_counts);
free(unsorted_factors);
return sorted_factors;
}
//------------------------------------------------------------------------------
// Compute prime-power factorization of a 64-bit value. Note that this function
// is designed to be fast *only* for numbers with very simple factorizations,
// e.g., those that produce large factor lists. Do not attempt to factor
// large semiprimes with this function. (The author does know how to factor
// large numbers efficiently; however, efficient factorization is beyond the
// scope of this small test program.)
PrimePowerFactorization compute_ppf(const uint64 n)
{
PrimePowerFactorization ppf;
if (n == 0)
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 };
}
else if (n == 1)
{
ppf = (PrimePowerFactorization){ .p = { 1 }, .e = { 1 },
.ω = 1, .Ω = 1, .d = 1, .n = 1 };
}
else
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = n };
uint64 m = n;
for (uint64 p = 2; p * p <= m; p += 1 + (p > 2))
{
if (m % p == 0)
{
assert(ppf.ω <= MAX_ω);
ppf.p[ppf.ω] = p;
ppf.e[ppf.ω] = 0;
while (m % p == 0)
{ m /= p; ppf.e[ppf.ω] += 1; }
ppf.d *= (1 + ppf.e[ppf.ω]);
ppf.Ω += ppf.e[ppf.ω];
ppf.ω += 1;
}
}
if (m > 1)
{
assert(ppf.ω <= MAX_ω);
ppf.p[ppf.ω] = m;
ppf.e[ppf.ω] = 1;
ppf.d *= 2;
ppf.Ω += 1;
ppf.ω += 1;
}
}
return ppf;
}
//------------------------------------------------------------------------------
// Parse prime-power factorization from a list of ASCII-encoded base-10 strings.
// The values are assumed to be 2-tuples (p,e) of prime p and exponent e.
// Primes must not exceed 2^64 - 1. Exponents must not exceed 2^8 - 1. The
// constructed value must not exceed 2^64 - 1.
PrimePowerFactorization parse_ppf(const uint pairs, const char *const values[])
{
assert(pairs <= MAX_ω);
PrimePowerFactorization ppf;
if (pairs == 0)
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 };
}
else
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = 1 };
for (uint i = 0; i < pairs; i++)
{
ppf.p[i] = (uint64)strtoumax(values[(i*2)+0], NULL, 10);
ppf.e[i] = (uint8)strtoumax(values[(i*2)+1], NULL, 10);
// Validate prime value.
if (ppf.p[i] < 2) // (Ideally this would actually do a primality test.)
fatal_error("Factor %"PRIu64" is invalid.", ppf.p[i]);
// Accumulate count of unique prime factors.
if (ppf.ω > UINT8_MAX - 1)
fatal_error("Small-omega overflow at factor %"PRIu64"^%"PRIu8".",
ppf.p[i], ppf.e[i]);
ppf.ω += 1;
// Accumulate count of total prime factors.
if (ppf.Ω > UINT8_MAX - ppf.e[i])
fatal_error("Big-omega wverflow at factor %"PRIu64"^%"PRIu8".",
ppf.p[i], ppf.e[i]);
ppf.Ω += ppf.e[i];
// Accumulate total divisor count.
if (ppf.d > UINT32_MAX / (1 + ppf.e[i]))
fatal_error("Divisor count overflow at factor %"PRIu64"^%"PRIu8".",
ppf.p[i], ppf.e[i]);
ppf.d *= (1 + ppf.e[i]);
// Accumulate value.
for (uint8 k = 1; k <= ppf.e[i]; k++)
{
if (ppf.n > UINT64_MAX / ppf.p[i])
fatal_error("Value overflow at factor %"PRIu64".", ppf.p[i]);
ppf.n *= ppf.p[i];
}
}
}
return ppf;
}
//------------------------------------------------------------------------------
// Main control. Parse command line and produce list of factors.
int main(const int argc, const char *const argv[])
{
PrimePowerFactorization ppf;
uint values = (uint)argc - 1; // argc is always guaranteed to be at least 1.
if (values == 1)
{
ppf = compute_ppf((uint64)strtoumax(argv[1], NULL, 10));
}
else
{
if (values % 2 != 0)
fatal_error("Odd number of arguments (%u) given.", values);
uint pairs = values / 2;
ppf = parse_ppf(pairs, &argv[1]);
}
// Run for (as close as possible to) a fixed amount of time, tallying the
// elapsed CPU time.
uint64 iterations = 0;
double cpu_time = 0.0;
const double cpu_time_limit = 0.05;
while (cpu_time < cpu_time_limit)
{
clock_t clock_start = clock();
uint64 *factors = compute_factors(ppf);
clock_t clock_end = clock();
cpu_time += (double)(clock_end - clock_start) / (double)CLOCKS_PER_SEC;
if (++iterations == 1)
{
for (uint32 i = 0; i < ppf.d; i++)
printf("%"PRIu64"\n", factors[i]);
}
if (factors) free(factors);
}
// Print the average amount of CPU time required for each iteration.
uint mem_scale = (memory_usage >= 1e9)? 9:
(memory_usage >= 1e6)? 6:
(memory_usage >= 1e3)? 3:
0;
char *mem_units = (mem_scale == 9)? "GB":
(mem_scale == 6)? "MB":
(mem_scale == 3)? "KB":
"B";
printf("%"PRIu64" %"PRIu32" factors %.6f ms %.3f ns/factor %.3f %s\n",
ppf.n,
ppf.d,
cpu_time/iterations * 1e3,
cpu_time/iterations * 1e9 / (double)(ppf.d? ppf.d : 1),
(double)memory_usage / pow(10, mem_scale),
mem_units);
return 0;
}
这个答案给出了一个 C 实现，我认为它是最快和最节省内存的。
算法概述。该算法基于 Will Ness 在
基本思想是使用数学归纳法对质数幂因式分解产生因子。例如：
- 要生成 360 的因数，需要计算 72 的因数，然后乘以 5 的相关幂，在本例中为 {1,5}。
- 为了产生 72 的因数，先计算 8 的因数，然后乘以 3 的相关幂，在本例中为 {1,3,9}。
- 为了产生 8 的因数，基本情况 1 乘以 2 的相关幂，在本例中为 {1,2,4,8}。
因此，在每个归纳步骤中，都会在现有因子集和素数幂集之间计算笛卡尔积，并通过乘法将结果还原为整数。
下面是数字 360 的插图。从左到右显示的是存储单元；一排代表一个时间步长。时间被表示为垂直向下流动。
空间复杂度。产生因子的临时数据结构极小：仅O(log₂(n))space 被使用，其中 n 是正在生成其因子的数字。例如，在该算法的 128 位 C 实现中，诸如 333,939,014,887,358,848,058,068,063,658,770,598,400（其以 2 为底的对数≈127.97）之类的数字分配 5.1 GB 来存储其 318,504,960 = 0 [但仅使用 4 = 136 ]bytes 的额外开销来生成该列表。任何 128 位数字最多需要大约 5 KB 的开销。
运行时复杂度。运行时完全取决于素数幂分解的指数（例如，素数签名）。为了获得最佳结果，应首先处理最大的指数，最后处理最小的指数，以便运行时在最后阶段由低复杂度合并主导，这在实践中通常是二进制合并。渐近运行时间为 O(c(n) d(n))，其中 d(n) 是 n 的除数，其中 c(n) 是一些依赖于 n 的素数签名的常数。也就是说，每个与素数签名相关联的等价 class 都有不同的常数。由小指数主导的素数签名具有更小的常数；由大指数支配的素数签名具有更大的常数。因此，运行时复杂度由素数签名聚集。
图表。 运行时性能在 3.4 GHz 上进行了分析。英特尔酷睿 i7 对 n 的 66,591 个值进行采样，具有 d(n) 个独特 d(n) 最多1.6亿。此图分析的 n 的最大值为 14,550,525,518,294,259,162,294,162,737,840,640,000，在 2.35 秒内产生 159,744,000 个因子。
每秒产生的排序因子数本质上是上面的倒数。按素数签名的聚类在数据中很明显。性能还受 L1、L2 和 L3 缓存大小以及物理 RAM 限制的影响。
源代码。 下面附上一个 C 编程语言（特别是 C11）的工作程序。它已经在 x86-64 上使用 Clang/LLVM 进行了测试，尽管它在 GCC 上也应该可以正常工作。
/*==============================================================================
DESCRIPTION
This is a small proof-of-concept program to test the idea of generating the
factors of a number in ascending order using an ultra-efficient sortless
method.
INPUT
Input is given on the command line, either as a single argument giving the
number to be factored or an even number of arguments giving the 2-tuples that
comprise the prime-power factorization of the desired number. For example,
the number
75600 = 2^4 x 3^3 x 5^2 x 7
can be given by the following list of arguments:
2 4 3 3 5 2 7 1
Note: If a single number is given, it will require factoring to produce its
prime-power factorization. Since this is just a small test program, a very
crude factoring method is used that is extremely fast for small prime factors
but extremely slow for large prime factors. This is actually fine, because
the largest factor lists occur with small prime factors anyway, and it is the
production of large factor lists at which this program aims to be proficient.
It is simply not interesting to be fast at producing the factor list of a
number like 17293823921105882610 = 2 x 3 x 5 x 576460797370196087, because
it has only 32 factors. Numbers with tens or hundreds of thousands of
factors are much more interesting.
OUTPUT
Results are written to standard output. A list of factors in ascending order
is produced, followed by runtime required to generate the list (not including
time to print it).
AUTHOR
Todd Lehman
2015/05/10
*/
//-----------------------------------------------------------------------------
#include <inttypes.h>
#include <limits.h>
#include <stdbool.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#include <math.h>
#include <assert.h>
//-----------------------------------------------------------------------------
typedef unsigned int uint;
typedef uint8_t uint8;
typedef uint16_t uint16;
typedef uint32_t uint32;
typedef uint64_t uint64;
typedef __uint128_t uint128;
#define UINT128_MAX (uint128)(-1)
#define UINT128_MAX_STRLEN 39
//-----------------------------------------------------------------------------
#define ARRAY_CAPACITY(x) (sizeof(x) / sizeof((x)[0]))
//-----------------------------------------------------------------------------
// This structure encode a single prime-power pair for the prime-power
// factorization of numbers, for example 3 to the 4th power.
#pragma pack(push, 8)
typedef struct
{
uint128 p; // Prime.
uint8 e; // Power (exponent).
}
PrimePower; // 24 bytes using 8-byte packing
#pragma pack(pop)
//-----------------------------------------------------------------------------
// Prime-power factorization structure.
//
// This structure is sufficient to represent the prime-power factorization of
// all 128-bit values. The field names ω and Ω are dervied from the standard
// number theory functions ω(n) and Ω(n), which count the number of unique and
// non-unique prime factors of n, respectively. The field name d is derived
// from the standard number theory function d(n), which counts the number of
// divisors of n, including 1 and n.
//
// The maximum possible value here of ω is 26, which occurs at
// n = 232862364358497360900063316880507363070 = 2 x 3 x 5 x 7 x 11 x 13 x 17 x
// 19 x 23 x 29 x 31 x 37 x 41 x 43 x 47 x 53 x 59 x 61 x 67 x 71 x 73 x 79 x
// 83 x 89 x 97 x 101, which has 26 unique prime factors.
//
// The maximum possible value of Ω here is 127, which occurs at n = 2^127 and
// n = 2^126 x 3, both of which have 127 non-unique prime factors.
//
// The maximum possible value of d here is 318504960, which occurs at
// n = 333939014887358848058068063658770598400 = 2^9 x 3^5 x 5^2 x 7^2 x 11^2 x
// 13^2 x 17 x 19 x 23 x 29 x 31 x 37 x 41 x 43 x 47 x 53 x 59 x 61 x 67 x 71 x
// 73 x 79.
//
#pragma pack(push, 8)
typedef struct
{
PrimePower f[32]; // Primes and their exponents.
uint8 ω; // Count of prime factors without multiplicity.
uint8 Ω; // Count of prime factors with multiplicity.
uint32 d; // Count of factors of n, including 1 and n.
uint128 n; // Value of n on which all other fields depend.
}
PrimePowerFactorization; // 656 bytes using 8-byte packing
#pragma pack(pop)
#define MAX_ω 26
#define MAX_Ω 127
//-----------------------------------------------------------------------------
// Fatal error: print error message and abort.
void fatal_error(const char *format, ...)
{
va_list args;
va_start(args, format);
vfprintf(stderr, format, args);
exit(1);
}
//------------------------------------------------------------------------------
uint128 uint128_from_string(const char *const str)
{
assert(str != NULL);
uint128 n = 0;
for (int i = 0; isdigit(str[i]); i++)
n = (n * 10) + (uint)(str[i] - '0');
return n;
}
//------------------------------------------------------------------------------
void uint128_to_string(const uint128 n,
char *const strbuf, const uint strbuflen)
{
assert(strbuf != NULL);
assert(strbuflen >= UINT128_MAX_STRLEN + 1);
// Extract digits into string buffer in reverse order.
uint128 a = n;
char *s = strbuf;
do { *(s++) = '0' + (uint)(a % 10); a /= 10; } while (a != 0);
*s = '[=10=]';
// Reverse the order of the digits.
uint l = strlen(strbuf);
for (uint i = 0; i < l/2; i++)
{ char t = strbuf[i]; strbuf[i] = strbuf[l-1-i]; strbuf[l-1-i] = t; }
// Verify result.
assert(uint128_from_string(strbuf) == n);
}
//------------------------------------------------------------------------------
char *uint128_to_static_string(const uint128 n, const uint i)
{
static char str[2][UINT128_MAX_STRLEN + 1];
assert(i < ARRAY_CAPACITY(str));
uint128_to_string(n, str[i], ARRAY_CAPACITY(str[i]));
return str[i];
}
//------------------------------------------------------------------------------
// Compute sorted list of factors, given a prime-power factorization.
uint128 *compute_factors(const PrimePowerFactorization ppf)
{
const uint128 n = ppf.n;
const uint d = (uint)ppf.d;
const uint ω = (uint)ppf.ω;
if (n == 0)
return NULL;
uint128 *factors = malloc((d + 1) * sizeof(*factors));
if (!factors)
fatal_error("Failed to allocate array of %u factors.", d);
uint128 *const factors_end = &factors[d];
// --- Seed the factors[] array.
factors_end[0] = 0; // Dummy value to simplify looping in bottleneck code.
factors_end[-1] = 1; // Seed value.
if (n == 1)
return factors;
// --- Iterate over all prime factors.
uint range = 1;
for (uint i = 0; i < ω; i++)
{
const uint128 p = ppf.f[i].p;
const uint e = ppf.f[i].e;
// --- Initialize phantom input lists and output list.
assert(e < 128);
assert(range < d);
uint128 *restrict in[128];
uint128 pe[128], f[128];
for (uint j = 0; j <= e; j++)
{
in[j] = &factors[d - range];
pe[j] = (j == 0)? 1 : pe[j-1] * p;
f[j] = pe[j];
}
uint active_list_count = 1 + e;
range *= 1 + e;
uint128 *restrict out = &factors[d - range];
// --- Merge phantom input lists to output list, until all input lists are
// extinguished.
while (active_list_count > 0)
{
if (active_list_count == 1)
{
assert(out == in[0]);
while (out != factors_end)
*(out++) *= pe[0];
in[0] = out;
}
else if (active_list_count == 2)
{
// This section of the code is the bottleneck of the entire factor-
// producing algorithm. Other portions need to be fast, but this
// *really* needs to be fast; therefore, it has been highly optimized.
// In fact, it is by far most frequently the case here that pe[0] is 1,
// so further optimization is warranted in this case.
uint128 f0 = f[0], f1 = f[1];
uint128 *in0 = in[0], *in1 = in[1];
const uint128 pe0 = pe[0], pe1 = pe[1];
if (pe[0] == 1)
{
while (true)
{
if (f0 < f1)
{ *(out++) = f0; f0 = *(++in0);
if (in0 == factors_end) break; }
else
{ *(out++) = f1; f1 = *(++in1) * pe1; }
}
}
else
{
while (true)
{
if (f0 < f1)
{ *(out++) = f0; f0 = *(++in0) * pe0;
if (in0 == factors_end) break; }
else
{ *(out++) = f1; f1 = *(++in1) * pe1; }
}
}
f[0] = f0; f[1] = f1;
in[0] = in0; in[1] = in1;
}
else if (active_list_count == 3)
{
uint128 f0 = f[0], f1 = f[1], f2 = f[2];
uint128 *in0 = in[0], *in1 = in[1], *in2 = in[2];
const uint128 pe0 = pe[0], pe1 = pe[1], pe2 = pe[2];
while (true)
{
if (f0 < f1)
{
if (f0 < f2)
{ *(out++) = f0; f0 = *(++in0) * pe0;
if (in0 == factors_end) break; }
else
{ *(out++) = f2; f2 = *(++in2) * pe2; }
}
else
{
if (f1 < f2)
{ *(out++) = f1; f1 = *(++in1) * pe1; }
else
{ *(out++) = f2; f2 = *(++in2) * pe2; }
}
}
f[0] = f0; f[1] = f1, f[2] = f2;
in[0] = in0; in[1] = in1, in[2] = in2;
}
else if (active_list_count >= 3)
{
while (true)
{
// Chose the smallest multiplier.
uint k_min = 0;
uint128 f_min = f[0];
for (uint k = 0; k < active_list_count; k++)
if (f[k] < f_min)
{ f_min = f[k]; k_min = k; }
// Write the output factor, advance the input pointer, and
// produce a new factor in the array f[] of list heads.
*(out++) = f_min;
f[k_min] = *(++in[k_min]) * pe[k_min];
if (in[k_min] == factors_end)
{ assert(k_min == 0); break; }
}
}
// --- Remove the newly emptied phantom input list. Note that this is
// guaranteed *always* to be the first remaining non-empty list.
assert(in[0] == factors_end);
for (uint j = 1; j < active_list_count; j++)
{
in[j-1] = in[j];
pe[j-1] = pe[j];
f[j-1] = f[j];
}
active_list_count -= 1;
}
assert(out == factors_end);
}
// --- Validate array of sorted factors.
#ifndef NDEBUG
{
for (uint k = 0; k < d; k++)
{
if (factors[k] == 0)
fatal_error("Produced a factor of 0 at index %u.", k);
if (n % factors[k] != 0)
fatal_error("Produced non-factor %s at index %u.",
uint128_to_static_string(factors[k], 0), k);
if ((k > 0) && (factors[k-1] == factors[k]))
fatal_error("Duplicate factor %s at index %u.",
uint128_to_static_string(factors[k], 0), k);
if ((k > 0) && (factors[k-1] > factors[k]))
fatal_error("Out-of-order factors %s and %s at indexes %u and %u.",
uint128_to_static_string(factors[k-1], 0),
uint128_to_static_string(factors[k], 1),
k-1, k);
}
}
#endif
return factors;
}
//------------------------------------------------------------------------------
// Print prime-power factorization of a number.
void print_ppf(const PrimePowerFactorization ppf)
{
printf("%s = ", uint128_to_static_string(ppf.n, 0));
if (ppf.n == 0)
{
printf("0");
}
else
{
for (uint i = 0; i < ppf.ω; i++)
{
if (i > 0)
printf(" x ");
printf("%s", uint128_to_static_string(ppf.f[i].p, 0));
if (ppf.f[i].e > 1)
printf("^%"PRIu8"", ppf.f[i].e);
}
}
printf("\n");
}
//------------------------------------------------------------------------------
int compare_powers_ascending(const void *const pf1,
const void *const pf2)
{
const PrimePower f1 = *((const PrimePower *)pf1);
const PrimePower f2 = *((const PrimePower *)pf2);
return (f1.e < f2.e)? -1:
(f1.e > f2.e)? +1:
0; // Not an error; duplicate exponents are common.
}
//------------------------------------------------------------------------------
int compare_powers_descending(const void *const pf1,
const void *const pf2)
{
const PrimePower f1 = *((const PrimePower *)pf1);
const PrimePower f2 = *((const PrimePower *)pf2);
return (f1.e < f2.e)? +1:
(f1.e > f2.e)? -1:
0; // Not an error; duplicate exponents are common.
}
//------------------------------------------------------------------------------
int compare_primes_ascending(const void *const pf1,
const void *const pf2)
{
const PrimePower f1 = *((const PrimePower *)pf1);
const PrimePower f2 = *((const PrimePower *)pf2);
return (f1.p < f2.p)? -1:
(f1.p > f2.p)? +1:
0; // Error; duplicate primes must never occur.
}
//------------------------------------------------------------------------------
int compare_primes_descending(const void *const pf1,
const void *const pf2)
{
const PrimePower f1 = *((const PrimePower *)pf1);
const PrimePower f2 = *((const PrimePower *)pf2);
return (f1.p < f2.p)? +1:
(f1.p > f2.p)? -1:
0; // Error; duplicate primes must never occur.
}
//------------------------------------------------------------------------------
// Sort prime-power factorization.
void sort_ppf(PrimePowerFactorization *const ppf,
const bool primes_major, // Best false
const bool primes_ascending, // Best false
const bool powers_ascending) // Best false
{
int (*compare_primes)(const void *, const void *) =
primes_ascending? compare_primes_ascending : compare_primes_descending;
int (*compare_powers)(const void *, const void *) =
powers_ascending? compare_powers_ascending : compare_powers_descending;
if (primes_major)
{
mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_powers);
mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_primes);
}
else
{
mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_primes);
mergesort(ppf->f, ppf->ω, sizeof(ppf->f[0]), compare_powers);
}
}
//------------------------------------------------------------------------------
// Compute prime-power factorization of a 128-bit value. Note that this
// function is designed to be fast *only* for numbers with very simple
// factorizations, e.g., those that produce large factor lists. Do not attempt
// to factor large semiprimes with this function. (The author does know how to
// factor large numbers efficiently; however, efficient factorization is beyond
// the scope of this small test program.)
PrimePowerFactorization compute_ppf(const uint128 n)
{
PrimePowerFactorization ppf;
if (n == 0)
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 };
}
else if (n == 1)
{
ppf = (PrimePowerFactorization){ .f[0] = { .p = 1, .e = 1 },
.ω = 1, .Ω = 1, .d = 1, .n = 1 };
}
else
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = n };
uint128 m = n;
for (uint128 p = 2; p * p <= m; p += 1 + (p > 2))
{
if (m % p == 0)
{
assert(ppf.ω <= MAX_ω);
ppf.f[ppf.ω].p = p;
ppf.f[ppf.ω].e = 0;
while (m % p == 0)
{ m /= p; ppf.f[ppf.ω].e += 1; }
ppf.d *= (1 + ppf.f[ppf.ω].e);
ppf.Ω += ppf.f[ppf.ω].e;
ppf.ω += 1;
}
}
if (m > 1)
{
assert(ppf.ω <= MAX_ω);
ppf.f[ppf.ω].p = m;
ppf.f[ppf.ω].e = 1;
ppf.d *= 2;
ppf.Ω += 1;
ppf.ω += 1;
}
}
return ppf;
}
//------------------------------------------------------------------------------
// Parse prime-power factorization from a list of ASCII-encoded base-10 strings.
// The values are assumed to be 2-tuples (p,e) of prime p and exponent e.
// Primes must not exceed 2^128 - 1 and must not be repeated. Exponents must
// not exceed 2^8 - 1, but can of course be repeated. The constructed value
// must not exceed 2^128 - 1.
PrimePowerFactorization parse_ppf(const uint pairs, const char *const values[])
{
assert(pairs <= MAX_ω);
PrimePowerFactorization ppf;
if (pairs == 0)
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 0, .n = 0 };
}
else
{
ppf = (PrimePowerFactorization){ .ω = 0, .Ω = 0, .d = 1, .n = 1 };
for (uint i = 0; i < pairs; i++)
{
ppf.f[i].p = uint128_from_string(values[(i*2)+0]);
ppf.f[i].e = (uint8)strtoumax(values[(i*2)+1], NULL, 10);
// Validate prime value.
if (ppf.f[i].p < 2) // (Ideally this would actually do a primality test.)
fatal_error("Factor %s is invalid.",
uint128_to_static_string(ppf.f[i].p, 0));
// Accumulate count of unique prime factors.
if (ppf.ω > UINT8_MAX - 1)
fatal_error("Small-omega overflow at factor %s^%"PRIu8".",
uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e);
ppf.ω += 1;
// Accumulate count of total prime factors.
if (ppf.Ω > UINT8_MAX - ppf.f[i].e)
fatal_error("Big-omega wverflow at factor %s^%"PRIu8".",
uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e);
ppf.Ω += ppf.f[i].e;
// Accumulate total divisor count.
if (ppf.d > UINT32_MAX / (1 + ppf.f[i].e))
fatal_error("Divisor count overflow at factor %s^%"PRIu8".",
uint128_to_static_string(ppf.f[i].p, 0), ppf.f[i].e);
ppf.d *= (1 + ppf.f[i].e);
// Accumulate value.
for (uint8 k = 1; k <= ppf.f[i].e; k++)
{
if (ppf.n > UINT128_MAX / ppf.f[i].p)
fatal_error("Value overflow at factor %s.",
uint128_to_static_string(ppf.f[i].p, 0));
ppf.n *= ppf.f[i].p;
}
}
}
return ppf;
}
//------------------------------------------------------------------------------
// Main control. Parse command line and produce list of factors.
int main(const int argc, const char *const argv[])
{
bool primes_major = false;
bool primes_ascending = false;
bool powers_ascending = false;
PrimePowerFactorization ppf;
// --- Parse prime-power sort specifier (if present).
uint value_base = 1;
uint value_count = (uint)argc - 1;
if ((argc > 1) && (argv[1][0] == '-'))
{
static const struct
{
char *str; bool primes_major, primes_ascending, powers_ascending;
}
sort_options[] =
{
// Sorting criteria:
// ----------------------------------------
{ "ep", 0,0,0 }, // Exponents descending, primes descending
{ "Ep", 0,0,1 }, // Exponents ascending, primes descending
{ "eP", 0,1,0 }, // Exponents descending, primes ascending
{ "EP", 0,1,1 }, // Exponents ascending, primes ascending
{ "p", 1,0,0 }, // Primes descending (exponents irrelevant)
{ "P", 1,1,0 }, // Primes ascending (exponents irrelevant)
};
bool valid = false;
for (uint i = 0; i < ARRAY_CAPACITY(sort_options); i++)
{
if (strcmp(&argv[1][1], sort_options[i].str) == 0)
{
primes_major = sort_options[i].primes_major;
primes_ascending = sort_options[i].primes_ascending;
powers_ascending = sort_options[i].powers_ascending;
valid = true;
break;
}
}
if (!valid)
fatal_error("Bad sort specifier: \"%s\"", argv[1]);
value_base += 1;
value_count -= 1;
}
// --- Prime factorization from either a number or a raw prime factorization.
if (value_count == 1)
{
uint128 n = uint128_from_string(argv[value_base]);
ppf = compute_ppf(n);
}
else
{
if (value_count % 2 != 0)
fatal_error("Odd number of arguments (%u) given.", value_count);
uint pairs = value_count / 2;
ppf = parse_ppf(pairs, &argv[value_base]);
}
// --- Sort prime factorization by either the default or the user-overridden
// configuration.
sort_ppf(&ppf, primes_major, primes_ascending, powers_ascending);
print_ppf(ppf);
// --- Run for (as close as possible to) a fixed amount of time, tallying the
// elapsed CPU time.
uint128 iterations = 0;
double cpu_time = 0.0;
const double cpu_time_limit = 0.10;
uint128 memory_usage = 0;
while (cpu_time < cpu_time_limit)
{
clock_t clock_start = clock();
uint128 *factors = compute_factors(ppf);
clock_t clock_end = clock();
cpu_time += (double)(clock_end - clock_start) / (double)CLOCKS_PER_SEC;
memory_usage = sizeof(*factors) * ppf.d;
if (++iterations == 0) //1)
{
for (uint32 i = 0; i < ppf.d; i++)
printf("%s\n", uint128_to_static_string(factors[i], 0));
}
if (factors) free(factors);
}
// --- Print the average amount of CPU time required for each iteration.
uint memory_scale = (memory_usage >= 1e9)? 9:
(memory_usage >= 1e6)? 6:
(memory_usage >= 1e3)? 3:
0;
char *memory_units = (memory_scale == 9)? "GB":
(memory_scale == 6)? "MB":
(memory_scale == 3)? "KB":
"B";
printf("%s %"PRIu32" factors %.6f ms %.3f ns/factor %.3f %s\n",
uint128_to_static_string(ppf.n, 0),
ppf.d,
cpu_time/iterations * 1e3,
cpu_time/iterations * 1e9 / (double)(ppf.d? ppf.d : 1),
(double)memory_usage / pow(10, memory_scale),
memory_units);
return 0;
}