In this article, I show how I adapted a mental calculation technique to check some small factors faster than division.
Download
Background
This article is part of a set of articles. Each article highlight an aspect around Integer Factorization.
The xPrompt Download and prg files
xPrompt.zip contains the binary of the xHarbour Scripting Engine for Windows. This is the app needed to run the prg code.
Usage :
- Copy prg file in same directory than xPrompt
- Launch xPromt
- Type "do FileName.prg" to run the code
xHarbour is freeware and can be dowloaded at xHarbour.org[^]
xHarbour is part of the xBase family languages:
xBase - Wikipedia[
^]
Introduction
Division is a good general purpose algorithm to check if a prime is factor of an integer, but sometimes we can get faster.
Small primes: Faster than division
When dealing with integer factorization, it is common to have a first stage for checking small factors before starting main algorithm. It happen that it is possible to check some of those factors faster than with a division.
Division is an efficient general purpose algorithm, but faster algorithms exist for particular divisors. Those algorithms resort to Mental Calculation techniques.
Some mental calculation techniques allow you to check easily the divisibility of a number for specific factors. The same techniques can be used on computer to check a few factors faster than with divisions because those thechniques also works with base 2.
Checking the Base
Humans Counting
Base 10 is how we count. Everyone can tell instantly if a number is multiple of 10 without resorting to division because the unit digit is all what matters. And everyone can also deduce if multiple of 2 and 5 as they are factors of 10.
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.
Cost is constant time. O(n) = 1.
Computer Counting
Base 2 is how computers count. The unit allow to know if 2 is a factor of a number. The fast operation is a bitwise AND.
2 is a factor if unit is 0.
Ex: 3141592653 => 0b10111011010000001110011001001101
unit is 1, so 2 is not a factor of 3141592653.
Cost is constant time. O(n) = 1.
if (N & 1) = 0
return(2)
endif
if ((Cand & 1) == 0) {
return 2;
}
Checking Other Factors
Humans Counting
Everyone having learned to do additions and multiplications by hand have also learned how to proof check those operations by using the casting out nines method. The method can be used to know if a number have 3 as factor, and also 11 and others with a little teaking. All this is possible because, 3 is a factor of 10 - 1 and 11 is a factor of 100 - 1.
The casting out nines method is about adding the digits of a number and reapply the method until the result is 1 digit.
Factor 3
The techbique is about checking 9 which is base 10 - 1, 3 being only a factor of 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.
Cost for n is log(n). O(n) = log(n).
Factor 11
The technique is about checking 11 which is base 10 + 1. We adapt the method with packets of 2 digits.
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.
Cost for n is log(n). O(n) = log(n).
Factor 101
The technique is about checking 101 which is base 10^2 + 1. We adapt the method with packets of 4 digits.
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.
Cost for n is log(n). O(n) = log(n).
Cascade
Those factor check can be combined when the number of digits in a check is a multiple of the number of digit of another.
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
The cost for factors 11 and 3 are in constant time when done after 101.
Computer Counting
The principle of "casting out nines" can be applied to base 2. Processors have very efficient tools that can be used to butcher numbers, those tools are addition/subtraction, bitwise operation and shifting.
Interesting Small Primes
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
Cascade Factors 17, 5 and 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
Cost for n is log(n) + constant. O(n) = log(n).
Intermediate oprimization
In modern processors, registers are 64 bits and simple operation take 1 clockcycle, no matter the size of an operation, the cost is the same. When one want to reduce a 64 bits value to a 8 bits value, it takes 7 operations.
By using intermediate steps, we can divide the size of operation by 2 every time. A 64 bits value is first reduced to 32 bits, the 16 bits, and finally to 8 bits. That 3 operations.
Algorithm
There are 4 cases to handle:
- Factor is a power of 2 minus 1
Reduce to power of 2.
Compare result with factor.
- Factor is a factor of a power of 2 minus 1
Reduce to power of 2.
do a modulo.
- Factor is a power of 2 plus 1
Reduce to double of power of 2.
test if both half are same.
- Factor is a factor of a power of 2 plus 1
Reduce to double of power of 2.
do a modulo on the difference of both half.
and optimizations:
Generator
As things can get complicated, It was easier to automate the code generation.
List of prime to integrate
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.}}
Topological sort for cascade
for scan=2 to len(list)
ptr= 0
for place= 1 to scan-1
if list[scan,5] % list[place,5] = 0
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
Adding intermediate steps
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
Code generator
? 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]
? " 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
? " 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]
? " tmp = FactReduce( N"+str(list[scan,5],,,.T.)+", Prime);"
? " if (tmp == 0) {"
? " return Prime;"
? " }"
elseif list[scan,3]= -1
? " if ( N"+str(list[scan,5],,,.T.)+" == Prime) {"
? " return Prime;"
? " }"
endif
?
next
? " return Cand;"
? "}"
Resulting code
function SmallPrimes(N)
if (N & 1) = 0
return(2)
endif
N16= Reduce(N, 16)
if (N16 >> 8) = (N16 & 255)
return(257)
endif
N8= Reduce(N16, 8)
if (N8 >> 4) = (N8 & 15)
return(17)
endif
N4= Reduce(N8, 4)
if (N4 >> 2) = (N4 & 3)
return(5)
endif
N2= Reduce(N4, 2)
if N2 = 3
return(3)
endif
N12= Reduce(N, 12)
Tmp= abs((N12 >> 6) - (N12 & 63))
while Tmp >= 13
Tmp= Tmp- 13
enddo
if Tmp = 0
return(13)
endif
N3= Reduce(N12, 3)
if N3 = 7
return(7)
endif
N10= Reduce(N, 10)
Tmp= abs((N10 >> 5) - (N10 & 31))
while Tmp >= 11
Tmp= Tmp- 11
enddo
if Tmp = 0
return(11)
endif
N5= Reduce(N10, 5)
if N5 = 31
return(31)
endif
N7= Reduce(N, 7)
if N7 = 127
return(127)
endif
return 1
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)
inline long long FactModulo(long long cand, long long fact) {
while (cand >= fact) {
if (cand & 1) {
cand = cand - fact;
}
cand = cand >> 1;
}
return cand;
}
inline long long FactFold(long long cand, long long shift, long long Mask) {
do {
cand = (cand >> shift) + (cand & Mask);
} while (cand >> shift);
return cand;
}
long long IsSFactNonDiv(long long Cand) {
long long Mask, Prime, tmp;
int Size;
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;
}
Benchmark
Runtime is so fast that timing is done on repeated execution of routines. Timing done on a i3 3GHz.
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
On last three numbers, my method is 2.5 times faster than division.
Points of Interest
Interest is limited if value is within register size, is value is big integer, the technique can be extended with more gains.
Unfortunately, it works only on a limited set of primes.
External Links:
History
- 20190501: First draft
- 20201230: Second version
I am a professional programmer.
Problem analyse is certainly what I am best at.
My main programming expertise is in the xBase languages (dBase, Clipper, FoxPro, Harbour, xHarbour), then VBA for Excel and advanced Excel WorkBooks.
I also have knowledge on C/C++, d language, HTML, SVG, XML, XSLT, Javascript, PHP, BASICs, Python, COBOL, Assembly.
My personal interests goes to algorithm optimization and puzzles.