文本算法 C++

整数分解：优化小因子检查

Patrice T

5.00/5 (1投票)

2019年5月30日

CPOL

5分钟阅读

13461

236

如何优化小因子检查

下载

背景

本文是一系列文章的一部分。每篇文章都重点介绍整数分解的一个方面。

整数分解：试除法算法[^]：关注试除法算法的变体及其各自的效率。
整数分解：可怕的素数列表[^]：关注通过压缩处理大型素数列表的方法。
整数分解：优化小因子检查[^]：关注一种比除法更快地检查小因子的方法。
整数分解：逆向乘法[^]：关注优化试除法算法的另一种方法。

xPrompt 下载和 prg 文件

xPrompt.zip 包含 Windows 版 xHarbour 脚本引擎的二进制文件。这是运行 prg 代码所需的应用程序。

用法
- 将 prg 文件复制到 xPrompt 的同一个目录
- 启动 xPromt
- 输入“do FileName.prg”来运行代码

xHarbour 是免费软件，可以从 xHarbour.org 下载 [^]

xHarbour 是 xBase 系列语言的一部分：xBase - Wikipedia [^]

引言

除法是一种很好的通用算法，用于检查一个素数是否是一个整数的因子，但有时我们可以做得更快。

小素数：比除法更快

在进行整数分解时，通常会在开始主算法之前有一个检查小因子的第一阶段。有时，检查其中一些因子比使用除法更快。

除法是一种有效的通用算法，但对于特定的除数存在更快的算法。这些算法采用**心算**技术。

一些心算技巧可以让你轻松地检查一个数是否能被特定的因子整除。这些技巧同样可以在计算机上用于比除法更快地检查一些因子，因为这些技巧也适用于基数为 2。

检查基数

人类计数

我们以 10 为基数计数。任何人都能够立即判断一个数是否是 10 的倍数，而无需进行除法，因为只看个位数就足够了。而且，由于 2 和 5 是 10 的因子，所以每个人也都能推断出是否是 2 或 5 的倍数。

2 is a factor if unit is 0, 2, 4, 6, 8.
5 is a factor if unit is 0, 5.
Ex: 3141592653
unit is 3, so 2 and 5 are not factor of 3141592653.

成本是常数时间。O(n) = 1。

计算机计数

计算机以 2 为基数计数。个位数可以用来判断 2 是否是该数的因子。快速操作是**按位与**。

2 is a factor if unit is 0.
Ex: 3141592653 => 0b10111011010000001110011001001101
unit is 1, so 2 is not a factor of 3141592653.

成本是常数时间。O(n) = 1。

    if (N & 1) = 0
        return(2)
    endif

	if ((Cand & 1) == 0) {
		return 2;
	}

检查其他因子

人类计数

任何学会手工进行加减乘除的人也学会了如何通过使用**九余法**来验证这些运算。该方法可用于判断一个数是否以 3 作为因子，并且稍加调整后也可以判断 11 等因子。这一切之所以可能，是因为 3 是 10-1 的因子，而 11 是 100-1 的因子。

九余法是通过将一个数的各位数字相加，然后重复此过程直到结果为一位数。

因子 3

该技术是检查 9，即 10-1，而 3 只是 9 的一个因子。

3 is a factor if result is 3, 6 or 9.
Ex: 3141592653 => 3+1+4+1+5+9+2+6+5+3 => 39 => 3+9 => 12 => 1+2 => 3
Result is 3, so 3 is a factor of 3141592653.

n 的成本为 log(n)。O(n) = log(n)。

因子 11

该技术是检查 11，即 10+1。我们使用两位一组的方法来适应。

for 11, we use the fact that 99 (100-1) is a multiple of 11.
11 is a factor if result is 0, 11, 22, 33 ... 88, 99 or if difference of the 2 digits is 0.
3141592653 => 31+41+59+26+53 => 210 => 2+10 => 12
Result is 12, so 11 is not a factor of 3141592653 because 1<>2.

n 的成本为 log(n)。O(n) = log(n)。

因子 101

该技术是检查 101，即 10^2+1。我们使用四位一组的方法来适应。

for 101, we use the fact that 9999 (10000-1) is a multiple of 101.
101 is a factor if result is 0, 101, 202, 303 ... 9999. or if difference of pairs of digits is 0.
3141592653 => 31+4159+2653 => 6843
So 101 is not a factor of 3141592653 because 68<>43.

n 的成本为 log(n)。O(n) = log(n)。

级联

当一个检查中的数字位数是另一个检查中数字位数的倍数时，这些因子检查可以组合。

Because 10000 is a power of 100, the result for factor 101 can be used to check factor 11 and 9.
3141592653 => 31+4159+2653 => 6843 => 101 is not a factor

Because 10000 is a power of 100, the result for 10000 can be reused for 100.
6843 => 68+43 => 111 => 1+11 => 12 => 11 is not a factor

Because 100 is a power of 10, the result for 100 can be reused for 10.
12 => 1+2 => 3 => 3 is a factor

当在 101 之后进行检查时，因子 11 和 3 的成本是常数时间。

计算机计数

“九余法”的原理可以应用于基数为 2 的情况。处理器有非常有效的工具可以用来处理数字，这些工具是加/减、位操作和移位。

有趣的小素数

3	2^2-1=3=1*3
5	2^2+1=5=1*5
7	2^3-1=7=1*7
11	2^5+1=33=1*3*11
13	2^6+1=65=1*5*13
17	2^4+1=17=1*17
19	2^9+1=513=1*3*3*3*19
23	2^11-1=2047=1*23*89
29	2^14+1=16385=1*5*29*113
31	2^5-1=31=1*31
37	2^18+1=262145=1*5*13*37*109
41	2^10+1=1025=1*5*5*41
43	2^7+1=129=1*3*43
47	2^23-1=8388607=1*47*178481
53	2^26+1=67108865=1*5*53*157*1613
59	2^29+1=536870913=1*3*59*3033169

级联因子 17、5 和 3

Ex: 3141592653 to binary 0b10111011010000001110011001001101

Reduce by 2 ^ 16 as intermediate optimization
0b10111011010000001110011001001101 => 0b1011101101000000+0b1110011001001101
=> 0b11010000110001101 => 0b1+0b1010000110001101 => 0b1010000110001110

Check 17: Reduce by 2 ^ 8
0b10100001+0b10001110 => 0b1000101110 => 0b10+0b00101110 => 0b00110000 => 0b0011 0b0000
17 is not a factor of 3141592653 because 0b0011<>0b0000

Check 5: Reduce by 2 ^ 4
0b00110000 => 0b0011+0b0000 => 0b0011 => 0b00 0b11
5 is not a factor of 3141592653 because 0b00<>0b11

Check 3: Reduce by 2 ^ 2
0b0011 => 0b00+0b11 => 0b11
3 is a factor of 3141592653

n 的成本为 log(n) + 常数。O(n) = log(n)。

中间优化

在现代处理器中，寄存器是 64 位的，简单的操作需要 1 个时钟周期，无论操作的大小如何，成本都是相同的。当一个人想将一个 64 位的值缩减到一个 8 位的值时，需要 7 次操作。

通过使用中间步骤，我们可以每次将操作的大小减半。一个 64 位的值首先缩减到 32 位，然后到 16 位，最后到 8 位。这需要 3 次操作。

算法

有 4 种情况需要处理

因子是 2 的幂减 1

缩减到 2 的幂。
将结果与因子进行比较。

因子是 2 的幂减 1 的因子

缩减到 2 的幂。
进行模运算。

因子是 2 的幂加 1

缩减到 2 的幂的两倍。
测试两个一半是否相同。

因子是 2 的幂加 1 的因子

缩减到 2 的幂的两倍。
对两个一半的差进行模运算。

以及优化

级联

进行拓扑排序。

中间优化

在需要的地方添加中间缩减。

生成器

由于事情可能会变得复杂，因此自动化代码生成更容易。

要集成的素数列表

*	List of: Prime, power of 2, +-1, is multiple
list={ {3, 2, -1, .F.}, {5, 2, +1, .F.}, {7, 3, -1, .F.}, {11, 5, +1, .T.}, ;
	{13, 6, +1, .T.}, {17, 4, +1, .F.}, {19, 9, +1, .T.}, {23, 11, -1, .T.}, ;
	{29, 14, +1, .T.}, {31, 5, -1, .F.}, {37, 18, +1, .T.}, {41, 10, +1, .T.}, ;
	{43, 7, +1, .T.}, {47, 23, -1, .T.}, {53, 26, +1, .T.}, {59, 29, +1, .T.}}

用于级联的拓扑排序

*	tri topologique
*	Optimization by chaining operations
for scan=2 to len(list)
	ptr= 0
	for place= 1 to scan-1
		if list[scan,5] % list[place,5] = 0 // scan multiple de place
			if list[scan,5] = list[place,5]
				tmp= list[scan]
				adel(list, scan, .T.)
				ains(list, place, tmp, .T.)
				ptr=0
				exit
			endif
			if ispower(list[scan,5] / list[place,5], 2)= 1
				tmp= list[scan]
				adel(list, scan, .T.)
				ains(list, place, tmp, .T.)
				ptr=0
				exit
			endif
		endif
		if list[place,5] % list[scan,5] = 0
			if ispower(list[place,5] / list[scan,5], 2)= 1
				ptr= place
			endif
		endif
	next
	if ptr != 0
		tmp= list[scan]
		adel(list, scan, .T.)
		ains(list, ptr+1, tmp, .T.)
	endif
next

添加中间步骤

*	pad list with intermediate opeartions
for scan=len(list) to 2 step -1
	while list[scan-1,5] % list[scan,5] = 0
		tmp= list[scan-1,5] / list[scan,5]
		tmpp= ispower(list[scan-1,5] / list[scan,5], 2)
		if tmp > 2 .and. tmpp=1
			ains(list, scan, {0,0,0,.F.,list[scan,5]*2}, .T.)
		else
			exit
		endif
	enddo
	while list[scan-1,5] < list[scan,5]
		if list[scan,5] < int_size/2
			ains(list, scan, {0,0,0,.F.,list[scan,5]*2}, .T.)
		else
			exit
		endif
	enddo
next
while list[1,5] < int_size/2
	ains(list, 1, {0,0,0,.F.,list[1,5]*2}, .T.)
enddo

代码生成器

*	generate C code
? int_type+" IsSFactNonDiv("+int_type+" Cand) {"
? "  "+int_type+" Mask, Prime, tmp;"
? "  int Size;"
? "  // fact 2"
? "  if ((Cand & 1) == 0) {"
? "    return 2;"
? "  }"
?

for scan=1 to len(list)
	if list[scan,1] != 0
		? "  Prime="+str(list[scan,1],,,.T.)+";"
	endif
	if scan=1
		org= "Cand"
	elseif list[scan,5]*2=list[scan-1,5]
		org= "N"+str(list[scan,5]*2,,,.T.)
	else
		org= "Cand"
	endif
	if (scan=1 .or. list[scan,5] != list[scan-1,5])
		if list[scan,5] < int_size
			? "  Size="+str(list[scan,5],,,.T.)+";"
			? "  Mask= (0ULL-1) >> ("+str(int_size,,,.T.)+"-Size);"
	
			if list[scan,3]=0
				? "  "+int_type+" N"+str(list[scan,5],,,.T.)+"= ("+org+" >> Size) + ("+org+" & Mask);"
			else
				? "  "+int_type+" N"+str(list[scan,5],,,.T.)+"= FactFold("+org+", Size, Mask);"
			endif
		else
			? "  "+int_type+" N"+str(list[scan,5],,,.T.)+"= "+org+";"
		endif
	endif

	if list[scan,3]=  1 .and. list[scan,4] // P2+1 et multiple
		? "  Size="+str(list[scan,2],,,.T.)+";"
		? "  Mask= (0ULL-1) >> ("+str(int_size,,,.T.)+"-Size);"
		? "  tmp = abs(( N"+str(list[scan,5],,,.T.)+" >> Size) - ( N"+str(list[scan,5],,,.T.)+" & Mask));"
		? "  tmp = FactReduce(tmp, Prime);"
		? "  if (tmp == 0) {"
		? "    return Prime;"
		? "  }"
	elseif list[scan,3]= 1 // P2+1
		? "  Size="+str(list[scan,2],,,.T.)+";"
		? "  Mask= (0ULL-1) >> ("+str(int_size,,,.T.)+"-Size);"
		? "  if (( N"+str(list[scan,5],,,.T.)+" >> Size) == ( N"+str(list[scan,5],,,.T.)+" & Mask)) {"
		? "    return Prime;"
		? "  }"
	elseif list[scan,3]= -1 .and. list[scan,4] // P2-1 et multiple
		? "  tmp = FactReduce( N"+str(list[scan,5],,,.T.)+", Prime);"
		? "  if (tmp == 0) {"
		? "    return Prime;"
		? "  }"
	elseif list[scan,3]= -1 // P2-1
		? "  if ( N"+str(list[scan,5],,,.T.)+" == Prime) {"
		? "    return Prime;"
		? "  }"
	endif

	?
next
? "  return Cand;"
? "}"

生成的代码

	// Fast checking of small primes
function SmallPrimes(N)
    //  Check factor 2
    if (N & 1) = 0
        return(2)
    endif

    //  Check factor 257
    N16= Reduce(N, 16)
    if  (N16 >> 8) = (N16 & 255)
        return(257)
    endif
    //  Check factor 17
    N8= Reduce(N16, 8)
    if  (N8 >> 4) = (N8 & 15)
        return(17)
    endif
    //  Check factor 5
    N4= Reduce(N8, 4)
    if  (N4 >> 2) = (N4 & 3)
        return(5)
    endif
    //  Check factor 3
    N2= Reduce(N4, 2)
    if  N2 = 3
        return(3)
    endif

    // Check factor 13
    N12= Reduce(N, 12)
    Tmp= abs((N12 >> 6) - (N12 & 63))
    while Tmp >= 13
      Tmp= Tmp- 13
    enddo
    if  Tmp = 0
        return(13)
    endif
    //  Check factor 7
    N3= Reduce(N12, 3)
    if  N3 = 7
        return(7)
    endif
    
    //  Check factor 11
    N10= Reduce(N, 10)
    Tmp= abs((N10 >> 5) - (N10 & 31))
    while Tmp >= 11
      Tmp= Tmp- 11
    enddo
    if  Tmp = 0
        return(11)
    endif
    //  Check factor 31
    N5= Reduce(N10, 5)
    if  N5 = 31
        return(31)
    endif
    
    //  Check factor 127
    N7= Reduce(N, 7)
    if  N7 = 127
        return(127)
    endif
    
    return 1

    // Binary 'casting out nines'
function Reduce(N, Size)
    local Tmp, Mask
    Tmp= N
    Mask= (1 << Size) -1
    while (Tmp >> Size) <> 0
      tmp= (Tmp >> Size) + (Tmp & Mask)
    enddo
    return(Tmp)

// My modulo routine
inline long long FactModulo(long long cand, long long fact) {
	// my modulo
	while (cand >= fact) {
		if (cand & 1) {
			cand = cand - fact;
		}
		cand = cand >> 1;
	}
	return cand;
}

// the reduction routine
inline long long FactFold(long long cand, long long shift, long long Mask) {
	// fold integer on itself
	do {
		cand = (cand >> shift) + (cand & Mask);
	} while (cand >> shift);
	return cand;
}

// My factor checker, result of automated code generation
long long IsSFactNonDiv(long long Cand) {
	long long Mask, Prime, tmp;
	int Size;
	// fact 2
	if ((Cand & 1) == 0) {
		return 2;
	}

	Size = 32;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N32 = (Cand >> Size) + (Cand & Mask);

	Size = 16;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N16 = (N32 >> Size) + (N32 & Mask);

	Prime = 17;
	Size = 8;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N8 = FactFold(N16, Size, Mask);
	Size = 4;
	Mask = (0ULL - 1) >> (64 - Size);
	if ((N8 >> Size) == (N8 & Mask)) {
		return Prime;
	}

	Prime = 5;
	Size = 4;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N4 = FactFold(N8, Size, Mask);
	Size = 2;
	Mask = (0ULL - 1) >> (64 - Size);
	if ((N4 >> Size) == (N4 & Mask)) {
		return Prime;
	}

	Prime = 3;
	Size = 2;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N2 = FactFold(N4, Size, Mask);
	if (N2 == Prime) {
		return Prime;
	}

	Size = 48;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N48 = (Cand >> Size) + (Cand & Mask);

	Size = 24;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N24 = (N48 >> Size) + (N48 & Mask);

	Prime = 13;
	Size = 12;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N12 = FactFold(N24, Size, Mask);
	Size = 6;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N12 >> Size) - (N12 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Size = 6;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N6 = (N12 >> Size) + (N12 & Mask);

	Prime = 7;
	Size = 3;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N3 = FactFold(N6, Size, Mask);
	if (N3 == Prime) {
		return Prime;
	}

	Size = 40;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N40 = (Cand >> Size) + (Cand & Mask);

	Prime = 41;
	Size = 20;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N20 = FactFold(N40, Size, Mask);
	Size = 10;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N20 >> Size) - (N20 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Prime = 11;
	Size = 10;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N10 = FactFold(N20, Size, Mask);
	Size = 5;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N10 >> Size) - (N10 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Prime = 31;
	Size = 5;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N5 = FactFold(N10, Size, Mask);
	if (N5 == Prime) {
		return Prime;
	}

	Prime = 37;
	Size = 36;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N36 = FactFold(Cand, Size, Mask);
	Size = 18;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N36 >> Size) - (N36 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Prime = 19;
	Size = 18;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N18 = FactFold(N36, Size, Mask);
	Size = 9;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N18 >> Size) - (N18 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Prime = 23;
	Size = 11;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N11 = FactFold(Cand, Size, Mask);
	tmp = FactModulo(N11, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Size = 56;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N56 = (Cand >> Size) + (Cand & Mask);

	Prime = 29;
	Size = 28;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N28 = FactFold(N56, Size, Mask);
	Size = 14;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N28 >> Size) - (N28 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Prime = 43;
	Size = 14;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N14 = FactFold(N28, Size, Mask);
	Size = 7;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N14 >> Size) - (N14 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Size = 46;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N46 = (Cand >> Size) + (Cand & Mask);

	Prime = 47;
	Size = 23;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N23 = FactFold(N46, Size, Mask);
	tmp = FactModulo(N23, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Prime = 53;
	Size = 52;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N52 = FactFold(Cand, Size, Mask);
	Size = 26;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N52 >> Size) - (N52 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	Prime = 59;
	Size = 58;
	Mask = (0ULL - 1) >> (64 - Size);
	long long N58 = FactFold(Cand, Size, Mask);
	Size = 29;
	Mask = (0ULL - 1) >> (64 - Size);
	tmp = abs((N58 >> Size) - (N58 & Mask));
	tmp = FactModulo(tmp, Prime);
	if (tmp == 0) {
		return Prime;
	}

	return Cand;
}

基准测试

运行速度如此之快，以至于计时是在例程重复执行的基础上进行的。计时是在一台 3GHz 的 i3 处理器上完成的。

Timing
Integer				Factor				Division	My Method
6561				3					9 ms		6 ms
125					5					18.001 ms	4 ms
2401				7					39.002 ms	9 ms
14641				11					41.002 ms	20.001 ms
3601				13					42.002 ms	11 ms
83521				17					52.003 ms	3 ms
49999				49999				147.008 ms	59.003 ms
4611686018427387899	4611686018427387899	143.008 ms	55.003 ms
4611686018427387877	4611686018427387877	143.008 ms	55.003 ms

对于最后三个数字，我的方法比除法快 2.5 倍。

关注点

如果值在寄存器大小内，那么兴趣有限。如果值是大数据，这项技术可以扩展并获得更多收益。

不幸的是，它只适用于有限的素数集。

外部链接

整数因式分解：整数因式分解 - 维基百科[^]
九余法：九余法 - Wikipedia [^]

历史

20190501：初稿
20201230：第二版