演算法?? - Prime
使用加法，?底***?字
例如5 = 1 + 1 + 1 + 1 + 1。
不可***的?元是1。不稀奇。
使用乘法，?底***?字
加法?成乘法之後，事情?化多端，例如5 = 5，例如6 = 2 × 3。可以??，2、3、5、7、11等，是不可***的?元。至於4、6、8、9、10等，可以再***。
不可***的?元叫做「??」！非常稀奇！
接下?要介?的演算法有：?小到大列出??（建立??表）、判?一??是不是??（????）、使用乘法?得?定的?（?因?***）。
Prime Generation: Sieve of Eratosthenes
Sieve of Eratosthenes
?是一??作??表的演算法。??「?法」。
列出所有正整?。?2?始，?掉2的倍?。找下一?未被?掉的?字，找到3，?掉3的倍?。找下一?未被?掉的?字，找到5，?掉5的倍?。如此不?下去，就能?掉所有合?，找到所有??。
??有?限多?，?明省略。我??法找到所有??，通常是?先?立一???，只找到???的所有??。
bool prime[];
void eratosthenes()
for (int i=0; i<; i++) // 初始化
prime[i] =
prime[0] = // 0 和 1 不是??
prime[1] =
// 找下一?未被?掉的?字
for (int i=2; i<; i++)
if (prime[i])
// ?掉??i的倍?，??倍?始。保留原本??。
for (int j=i+i; j<; j+=i)
prime[j] =
欲?掉??i的倍?之?，早已?掉1倍到i-1倍了，直接?i倍?始。
bool prime[];
void eratosthenes()
for (int i=0; i<; i++)
prime[i] =
prime[0] =
prime[1] =
for (int i=2; i<; i++)
if (prime[i])
// ?掉??i的倍?，?i倍?始。小心溢位。
for (int j=i*i; j<; j+=i)
prime[j] =
欲?掉??i的倍?之?，早已?掉「小於i的??、其倍?」倍了，直接?掉「大於等於i的??、其倍?」倍。
乍看下程式?增多而?慢，??上cache miss?少而?快。
bool prime[];
void eratosthenes()
for (int i=0; i<; i++)
prime[i] =
prime[0] =
prime[1] =
for (int i=2; i=i; k--, j-=i)
if (prime[k])
prime[j] =
// ?用小??的prime[k]，避免存取大??的prime[j]，
// 大幅?少cache miss。
一?合?x，必定有一?小於等於sqrt(x)的?因?。所有小於等於sqrt(x)的??，?掉?些??的倍?，就能?掉所有合?了。
bool prime[];
void eratosthenes()
for (int i=0; i<; i++)
prime[i] =
prime[0] =
prime[1] =
// 只需要?掉sqrt()以下的??的倍?。
int sqrt_ = sqrt();
for (int i=2; i=i; k--, j-=i)
if (prime[k])
prime[j] =
?倒true和false，?省初始化??。
// 全域??的?列，每一格?自?初始化?false。
bool sieve[];
void eratosthenes()
int sqrt_ = sqrt();
sieve[0] = sieve[1] =
for (int i=2; i=i; k--, j-=i)
if (!sieve[k])
sieve[j] =
?作??表。?法?束之後，?描一次?列即可。
bool sieve[];
void eratosthenes()
// ?作??表
for (int i=2; i<; i++)
if (!sieve[i])
prime.push_back(i);
UVa 406 516 524 543
使用bitset?取代bool?列
一?int有32?位元，可以?作32??位?使用，?省???空?，?少cache miss。
const int N = ;
int sieve[(N>>5) + 1]; // 加一，以包含??。
GET(int x) { return sieve[x>>5]>>(x&31) & 1; }
void SET(int x) { sieve[x>>5] |= 1<<(x&31); }
void eratosthenes()
// register??字???存放於CPU register，加快存取速度。
register int i, j,
int sqrt_N = sqrt(N);
SET(0); SET(1);
for (i=2; i=i; k--, j-=i)
if (!GET(k))
不?理2的倍?
不?理2的倍?，?省一半???，增?一?速度。
令?列第0格代表?字1、?列第1格代表?字3、?列第2格代表?字5、……，以此?推。
const int N = 1<>6) + 1]; // 不?理2的倍?
GET(int x) { return sieve[x>>5]>>(x&31) & 1; }
void SET(int x) { sieve[x>>5] |= 1<>1; }
I2N(int i) { return (i<>1) // 索引值 = (?字 - 1) / 2
#define I2N(i) (((i)<2 && (x&1) && !GET(N2I(x));
void eratosthenes()
int half_sqrt_N = sqrt(N) / 2;
int half_N = N >> 1;
register int i, j,
// ?3?始。N2I(3) = 1。
for (i=1; i<=half_sqrt_N; i++)
if (!GET(i))
for (j=2*i*(i+1), k=2*i+1; j void eratosthenes()
int half_sqrt_N = sqrt(N) / 2;
register int i, j, k,
for (i=1; i=i; j--, k-=ii)
if (!GET(j))
考????圈索引值j一共有多少?：(N/2 - 2) + (N/3 - 3) + (N/5 - 5) + ... + (N/sqrtN - sqrtN) = O(NloglogN)。
1/1 + 1/2 + 1/3 + ... + 1/N
1/2 + 1/3 + 1/5 + ... + 1/N
= O(loglogN)
1/2 + 1/3 + 1/5 + ... + 1/sqrtN = O(loglogsqrtN) = O(loglogN)
http://en.wikipedia.org/wiki/Divergence_of_the_sum_of_the_reciprocals_of_the_primes
1/1 + 1/2 + 1/3 + ... + 1/n - ln(n)
?近於Euler-Mascheroni constant
1/2 + 1/3 + 1/5 + ... + 1/n - ln(ln(n)) ?近於Meissel-Mertens constant
Prime Generation: Segmented Sieve of Eratosthenes
Segmented Sieve of Eratosthenes
?法?中，以不大於sqrt(N)的??，??大於sqrt(N)的?。有?洞??先、以小搏大的味道。
?法需要大量???空?。?了?省???空?，於是有了分段?理的手法：一、首先求得不大於sqrt(N)的??。二、????切成小段，每段?度sqrt(N)。分???每一段。
空???度?O(N)?成O(sqrt(N))，大幅?少cache miss。
http://www.geeksforgeeks.org/segmented-sieve/
也?可以做?平行演算法的?典?例。
Prime Generation: Linear Sieve Algorithm
一??作??表，一??掉每??的??倍，如此每?合?就只?取一次，????度?到O(N)。
const int N = ;
bool sieve[N];
void linear_sieve()
for (int i=2; i if (!sieve[i]) prime.push_back(i);
for (int j=0; i*prime[j] sieve[i*prime[j]] =
if (i % prime[j] == 0)
Prime Generation: Wheel Factorization
6n±1 Method
?是一?精?版的?法，原理是：只拿2和3?????先??一遍，剩下的?字?用?除法??是不是??。
2和3的最小公倍?是6，我?就把所有?字分?6n、6n+1、6n+2、6n+2、6n+3、6n+4、6n+5六?（n是倍率）。除以六??零的?字?6n，除以六??一的?字?6n+1，以此?推。
可以看出6n、6n+2、6n+3、6n+4?是2或3的倍?，不?於??。因此，只要??6n+1和6n+5是不是??就可以了。（6n+5也可以?成6n-1，意?不?。）
6n-1到6n+1，再到下一?6n-1，再到6n+1，把?些要??的?字由小排到大，可以??之?的差值?是2 4 2 4 2 4 ...不?跳二跳四。?作程式??，就可以直接用加法加二加四，而不必用乘法及加?法?算6n-1、6n+1，如此一?程式的?行效率?好一?。
??的?序是：?字2和3明?的是??，不必??；若是??字5?始??，那?下一?要??的?字就是5+2，再下一?就是5+2+4，再下一?就是5+2+4+2，如此不?下去。
// 以?除法判???
bool isprime(int n)
for (int i=5; i*i<=n; i+=2)
if (n%i == 0)
// 只?查6n±1。?些?字的?隔?2 4 2 4....
void make_prime()
cout << "找到??2";
cout << "找到??3";
for (int i=5, gap=2; i<1000000; i+=gap, gap=6-gap)
if (isprime(i))
cout << "找到??" <<
// 以?除法判???
bool isprime(int n)
for (int i=0; prime[i]*prime[i]<=n; i++)
if (n % prime[i] == 0)
// 只?查6n±1。?些?字的?隔?2 4 2 4....
void make_prime()
cout << "找到??2"; prime.push_back(2);
cout << "找到??3"; prime.push_back(3);
for (int i=5, gap=2; i<1000000; i+=gap, gap=6-gap)
if (isprime(i))
{cout << "找到??" << prime.push_back(i);}
??方法的????度是O(NsqrtN)，空???度是O(1)。事?上6n±1法比?法慢上?多。不?6n±1法不需要?一?超大?列?做?算，?省了很多空?。
Wheel Factorization
6n±1法?行推?，使用?????。
例如取2 3 5的最小公倍?30，建立30n+k系列。又例如取2 3 5 7的最小公倍?210，建立210n+k系列。
取越多??，即可剔除越多合?，事前?省??。
然而，取越多??，也?建立越多??，事中浪???。
?是取???的情?，?者之?必?取得平衡。一般??，大家只取2 3或者2 3 5。
更有甚著，?你?在??要取多少?，普通的?法早就已??行完?了。更有甚者，?你?在??要如何?作?法?，?路上已?有?成的??表了。所?的取得平衡，端看你的?野高度。
Primality Test
Primality Test
????，??一??字是否???。
?????於P??。然而以下介?的都是指???演算法，甚至是不保??果正?的演算法。多?式??、保??果正?的演算法，??行上?搜?AKS Algorithm；但由於方法繁?，??上?人使用。
????，可以直接以?法建立??表，再?判???。然而建立??表需要大量???，因此又?明了其他方法。
Divisibility Primality Test
整除性??法。按照??定?：???法以乘法***、???有任何因?（除了1和自身）。把所有可能的因?拿??除，也就是把介於1到n之?的?拿??除，就能判定??。
此演算法?施了n-2次除法，????度?O(N)。然而前提是：不管n多大，每次除法都是O(1)。
bool divisibiity_test(int n)
// ??所有可能的因?，一一?除（除了1和自身）
for (int d=2; d if (n % d == 0)
因?成?成?（平方根跟自己一?）。?查了小於等於sqrt(n)的因?，就不必再?查大於sqrt(n)的因?。
此演算法?施了sqrt(n)-1次除法，????度?O(sqrtN)。然而前提是：不管n多大，每次除法都是O(1)。
bool divisibiity_test(int n)
// ?查小於等於sqrt(n)的因?（除了1）
int sqrt_n = sqrt(n);
for (int d=2; d<=sqrt_n; ++d)
if (n % d == 0)
因?一定可以??***成?因?，甚至自身就是?因?。?查小於等於sqrt(n)的?因?，效果等同於?查小於等於sqrt(n)的因?。
?因?是??的一部分。上述?容可以?一步改成?查小於等於sqrt(n)的??，一定?涵?到小於等於sqrt(n)的?因?。
若有大量?字需要??，可以先建立??表，??查??。
Square Root Primality Test
??系?的平方根有?特性：
以??n?模?，1的平方根一定等於±1（即是1?n-1）。
以合?n?模?，1的平方根除了等於±1（即是1?n-1），也可能等於其他?字。
平方根??法是?用此特性而想出的方法：
以??n?模?，0到n-1之?，只有1?n-1，平方之後等於1。
以合?n?模?，0到n-1之?，?可能有其它?字，平方之後等於1。
以n?模?，?2到n-2的每一??字，平方之後皆不等於1，就推定n是??。
此演算法的?果不一定正?。通???的?字，可能是合?或??；?法通???的?字，一定是合?。
有些合??被?判成??，例如22就?被判定成??，1?、2?、…、21?模22之後?好都?1。
bool square_root_primality_test(int n)
for (int i=1; i if (i * i % n == 1)
Fermat's Primality Test
??小定理：
若n是??，a小於n，?a? ≡ a (mod n)。
若n是??，a小於n，a不是零，?a??? ≡ 1 (mod n)。
??????法是?用??小定理而想出的方法：
n是??，??小定理一定成立：a??? % n = 1一定成立。
n是合?，??小定理可能成立：a??? % n = 1可能成立。
?a??? % n = 1成立，就推定n是??。
此演算法的?果不一定正?。通???的?字，可能是合?或??；?法通???的?字，一定是合?。
使用各式各?的a??施??，那?判定?果就更??。
bool fermat_primality_test(int n)
// ?便取得一??值?作a，
// 但是必?大於1、小於n，才有意?。
// srand(time(0));
int a = rand() % (n-2) + 2;
// ?算a??? % n的值，存於x中。
int x = 1;
for (int i=1; i x = x * a %
return x == 1;
// return pow(a, n-1, n) == 1;
Primality Test: Miller-Rabin Algorithm
?合Square Root Primality Test?Fermat's Primality Test。
此演算法的?果不一定正?。通???的?字，可能是合?或??；?法通???的?字，一定是合?。
一、?定一?底?a，大於1、小於n，用??行????。
　　待??字n，?是????的模?。
二、令 n-1 = u × 2?，其中t?量大（u?奇?）。
三、?察a^u???字：
　　若等於±1，
　　?表示步?四，所有?字都是1，永不出?平方根??的反例。
　　也表示步?五，最??通?????。
　　推定n是??。
四、依序?察a^(u × 2?)、a^(u × 2?)、…、a^(u × 2???)?些?字，
　　每??字正好是前一??字的平方：
　甲、一旦??有??字的平方等於+1，
　　　?表示?法通?平方根??。
　　　（但是步?五，最??通?????。）
　　　?定n是合?。
　乙、一旦??有??字的平方等於-1，
　　　?表示接下?的?字都是1，永不出?平方根??的反例。
　　　也表示步?五，最??通?????。
　　　推定n是??。
五、?察a^(u × 2?)???字：
　甲、若等於+1，表示通?????，推定n是??。
　乙、若不等於+1，表示?法通?????，?定n是合?。
　回、也就是?，a^(u × 2???)必?等於±1，平方之後才?等於+1。
　　　步?五才能通?????。
　回、由於步?四已??查?a^(u × 2???)是否?±1，
　　　所以步?五可以省略，直接?定n是合?。
bool miller_rabin(int n)
if (n >= 1, t++;
int x = pow(a, u, n); // x = a ^ u %
if (x == 1 || x == n-1)
for (int i=0; i x = mul(x, x, n); // x = x * x %
if (x == 1)
if (x == n-1)
步?二：不?***因?2，最多log(n)步，需?O(logN)。
步?三：??次方，Divide and Conquer，需?O(logN)。
步?四：根?步?二，最多log(n)-1??字，需?O(logN)。
?????度?O(logN)，然而前提是：不管N多大，每次??乘法都是O(1)。
Strong Probable-prime Base（SPRP）
以??合?的角度?看，多取??相?的底?a?施??，判定?果就更??。
事?上已?有?心人士，找出特定的底??合，仔??查了各??字的判定?果，保?判定?果正?。例如底??合{2, 7, 61}就能正?判?2??以下的?字是不是??。
bool isprime(int n)
// ?先判?偶??1，?省一???。
if (n == 2)
if (n >= 1; t++;}
// 推定是??，就?施下一次??；
// ?定是合?，就?上?束。
int sprp[3] = {2, 7, 61}; // 足以涵?int??
for (int k=0; k<3; ++k)
// a?有大於1、小於n-1的限制，
// ?有??意?的?字，?作是通???。
// （a是???，模n後不?等於零，不必特?判?。）
int a = sprp[k] %
if (a == 0 || a == 1 || a == n-1)
long long x = pow(a, u, n);
if (x == 1 || x == n-1)
for (int i = 0; i < t-1; i++)
x = mul(x, x, n);
if (x == 1)
if (x == n-1)
if (x == n-1)
// 通?全部??，?定是??。
UVa 10956 PKU 1811
Integer Factorization
Integer Factorization
把一?正整?***成?因?的?乘?。
n = 2n? × 3n? × 5n? × 7n? × 11n? × …
Factor是「因式」。Factorization是「因式***」。Integer Factorization是「因式***的?象?整?」，?作「整?***」。中??本?作「?因?***」，著眼於***?果而非***?象。
?因?***?於NP??，但是目前?不?定它究竟是P??或是NP-complete??，相?特?。
Fundamental Theorem of Arithmetic（算?基本定理）
凡是正整?都可以藉由?因?***成?一??一?二的式子，不同的n就???不同的(n?, n?, …)，反方向亦同。
根?算?基本定理，凡是?扯到乘法、因?、倍?的???算，都可以改?成比???的?算。
| (n?, n?, …)
| (m?, m?, …)
--------------------------------------------------
| (n? + m?, n? + m?, …)
| (n? - m?, n? - m?, …)
| 大於等於
| (n?, n?, …) ≥ (m?, m?, …)
| n?≥m? and n?≥m? and …
最大公因? | 最小值
| (min(n?, m?), min(n?, m?), …)
最小公倍? | 最大值
| (max(n?, m?), max(n?, m?), …)
| ???必?有零
| min(n?, m?) = 0 and min(n?, m?) = 0 and …
| n?×m? = 0 and n?×m? = 0 and …
算?基本定理?述了另一?世界?，把?字看作是??的?合。??的英文prime有著「原始就有」的意思，便是指??是所有?字的根本。
Trial Division Factorization Method
把所有可能的因?拿??除。
void trial_division(int n)
// ??所有可能的因?，一一?除（除了1和自身）
for (int d=2; d<=n; ++d)
while (n % d == 0)
cout << // 印出?因?
void trial_division(int n)
// ?查小於等於sqrt(n)的因?（除了1）
int sqrt_n = sqrt(n);
for (int d=2; d<=sqrt_n; ++d)
while (n % d == 0)
cout < 1) cout << // n是??
UVa 516 583
Integer Factorization: Pollard's ρ Algorithm
此演算法可以找到n的其中一?因?。
使用一???的???生器f(a???) = a?? + c (mod n)，??各?a??c，?造所有可能的因?，一一?除即可。
以此???生器公式，依序枚?a? a? a? …，模n的情?下，最?必定循?。???通常?成一?ρ的形?，演算法因而得名。ρ?作[ro]，可以?作rho。
?用最大公因?找到因?
因????生器?造的?字a，a恰是n的因?的???小，而a?n有共同因?的???大，所以改用d = gcd(a, n)?找到n的因?d。最大公因?有著?快的演算法，??行速度影?不大。
???生器最?必?出?重覆?字，?生循?。一旦遇到循?，立刻?束枚?，不再?行重覆?算。
另外，此演算法改用abs(ax - ay)，用?字的差值?造所有可能的因?。我不知道如此做的原因。
原?文?用「」??循?。
#define f(x) ((x * x + c) % n) // ???生器
int pollard_rho(int n)
int a0 = rand() % // 最好是+2
int c = rand() %
// 最好不是0和-2
int x = a0, y = a0;
while (true)
// Floyd's Algorithm。
// 一旦出?循?，?回?原本?字n。
x = f(x); y = f(f(y));
if (x == y)
int d = __gcd(abs(x - y), n);
if (d > 1 && d < n)
return -1; // 不??行到?行
#define f(x) ((x * x + c) % n) // ???生器
int pollard_rho(int n)
int a0 = rand() % // 最好是+2
int c = rand() %
// 最好不是0和-2
int x = a0, y = a0;
while (true)
x = f(x); y = f(f(y));
if (x == y)
// 因?出?循?之?，x-y等於0，gcd(0,n)等於n。
// 所以亦得用d==n判?是否出?循?，?少一行程式?。
int d = __gcd(abs(x - y), n);
if (d > 1 && d <= n)
return -1; // 不??行到?行
亦可?用「」??循?，效率?佳。
int pollard_rho(int n)
int a0 = rand() % // 最好是+2
int c = rand() %
// 最好不是0和-2
int x = a0, y = a0;
int i = 1, k = 2;
while (true)
x = (x * x + c) %
if (x == y)
int d = __gcd(abs(x - y), n);
if (d > 1 && d <= n)
if (++i == k) y = x, k <<= 1;
return -1; // 不??行到?行
?者可以思考一下?些??
一、?何???生器不?用f(a???) = a? + c (mod n)??更加??的式子？
二、?何a?至少要是+2？（????，a?最好是+2。）
三、?何c最好不是0和-2？??看????生器公式，代入到x和y之中，?算一下x-y，然後?算一下gcd(abs(x-y), n)。
四、?何x等於y的?候，就要?上?束?圈呢？
五、如果略去abs()，改成gcd(x-y, n)，??果有影??？
似乎只要??各?c，就一定可以??出所有可能的因?。我不知道原因。
vector // 存放n的?因?
int pollard_rho(int n, int c)
int x = 2, y = 2, i = 1, k = 2;
while (true)
x = (mul(x, x, n) + c) %
int d = __gcd(abs(x - y), n);
if (d > 1 && d <= n)
if (++i == k) y = x, k <<= 1;
return -1;
void factorize(int n)
if (n == 1)
// ?先判???。
// 可以使用Miller-Rabin Algorithm with SPRP。
if (isprime(n))
factor.push_back(n);
// 此?只有合?才能?施Pollard's ρ Algorithm，
// 否?永?找不到1?n以外的因?，陷入???圈。
// ??各?c，直到??n的因?d，一定能找到。
for (int c=0; d == ++c)
d = pollard_rho(n, c);
factorize(d);
factorize(n / d);
Integer Factorization: Quadratic Sieve Algorithm
Fermat's Factorization Method
一??字n，***成???x y的乘?。
?用平方差公式，令n = a? - b? = (a+b) (a-b) = x y。?找?好相差n的??平方?a??b?，就能***n。
?作程式??，??整?a，然後推?出b = sqrt(a? - n)，判?b是不是整?。
void fermat(int n)
for (int a = ceil(sqrt(n)); a int bb = a * a -
int b = sqrt(bb);
if (b * b == bb) // 判?b是不是整?
int x = a + b, y = a -
cout << "n =" << x << "*" <<
// ??***x和y，?法再***的?候就找到?因?了。
fermat(x); fermat(y);
另一?更直?的解?方式：
??a?，?算a? - n，判?是不是平方?。如果是平方?，便得到了?好相差n的??平方?a??b?。
a?由小到大依序??。由於a?必?大於等於n，才能相?得到b?，所以由等於n、略大於n的平方??始??。
a? = ceil(sqrt(52)) = 8
a? | a? - n | ? ? |
---|-----|--------|-----|
10 | 100 |
11 | 121 |
12 | 144 |
13 | 169 |
14 | 196 |
52 = 196 - 144 = 14? - 12? = (14+12) (14-12) = 26 × 2
= (a+b) (a-b)
Fermat's Factorization Method推?至??系?
本?是?找?好相差n的??平方?，?在是?找相差kn的??平方?，k是任意整?倍率。放?限制，??更多。
令kn = a? - b? = (a+b) (a-b)。想要***n，可以?算d = gcd(n, a+b)，?a+b?中?取n的因?。或者a-b亦可。
??不好?，d可能等於1或n，?有?到***效果；此?再??其他平方??合即可。
回到一?始。n?kn，可以想成是模n。原???成：?找??平方?，模n同?。
n = a? - b?
n ≡ a? - b? (mod n)
推?成??。模n，就?成差kn。
0 ≡ a? - b? (mod n)
a? ≡ b? (mod n)
移?。式子好看又好用。
?d = gcd(n, a+b) = 1 or n，此?a ≡ ±b (mod n)。
Dixon's Factorization Method
??a?，可是a? - n一直不是平方?，怎??？
各?a? - n相乘，?出平方?，如此便得到了?好相差kn的??平方?a??b?。
a? = ceil(sqrt(26)) = 6
a? | a? - n | ? ? |
---|-----|--------|-----|
16 | 256 |
(6×7×16)? ≡ 230? (mod 26)
{ 16? ≡ 230 (mod 26)
等?左?相乘 6? × 7? × 16? = (6×7×16)? 平方?相乘仍是平方?
等?右?相乘 10×23×230 = 230?
人工?得平方?
得到一?相差 kn 的??平方? a? ? b?！
?了方便?出平方?，?用算?基本定理：a? - n?施?因?***，取次方值，形成向量。平方?的次方值皆?偶?。
?了快速?出平方?，?用?性代?：向量模2，挑其中??向量，XOR等於零。?性?合等於零。解?性方程?。高斯消去法。
??所有?合，?化成高斯消去法，????度?指??成?性，大幅改?！此?精髓！
a? | a? - n | ? ? | factor
| vector (mod 2)
---|-----|--------|-----|----------|-----------------|
| [1 0 1 0 0 ...] | ?
| [0 0 0 0 0 ...] | ?
| [1 0 0 0 0 ...] |
| [0 0 1 0 1 ...] |
16 | 256 |
| [1 0 1 0 0 ...] | ?
24 | 576 |
| 2 5 5 11 | [1 0 0 0 1 ...] |
[ 1 0 1 0 .. 1 .. 1 .. ] [ x?
[ 0 0 0 0 .. 0 .. 0 .. ] [ x?
[ 1 0 0 1 .. 1 .. 0 .. ] [ x?
solve [ 0 0 0 0 .. 0 .. 0 .. ] [
[ 0 0 0 1 .. 0 .. 1 .. ] [ x?? ]
矩?通常很大，於是?用?推法：a? - n?施?因?***，逐步增加?因?；每多一??因?，就解一次?性方程?。
每回合只拿已??底***的a? - n，?解?性方程?；也就是除到剩下1的。
用??????：逐步增加B，每次只拿B-smooth的a? - n，?解?性方程?。
a? | a? - n |
---|-----|--------|--------------|--------------|--------------
10 | 100 |
24 | 576 |
25 | 625 |
最後介?一?加速技巧：剩下的?字，如果不?1，但是?字一?，此?可以融合得到一?新的B-smooth。
有?候剩下的?字很大、仍需?回合才能除?。?用此技巧，得以提早?入一些a? - n，提早找到正?***。
|| round 5
a? | a? - n || ÷11
---|-----|--------||--------------|
14 | 196 |
20 | 400 |
?好都剩17，可以相乘，去掉17，得到一?新的11-smooth。
甚至可以考?消除原本那??。
|| round 5
a? | a? - n || ÷11
---|-----|--------||--------------|
280| *** | ****** ||
1 | 2 5 11 |
Quadratic Sieve Algorithm
整除?因?2的?隔似乎是2。似乎是?法。仔?推?一下：
已知某? a? - n 整除?因? p。?在令?隔?p：
(a+p)? - n
= a? + 2ap + p? - n
= (a? - n) + 2ap + p?
= (a? - n) + p (2a+p)
仍整除?因?p！
在0到p-1之?，找到所有的a，?足a? - n整除?因?p，做??法的起?；跳??隔皆是p。如此一?便?一?漏。
a? - n ≡ 0 (mod p)
a? ≡ n (mod p)
已知n p，解a，找出所有的a。
此形式的??方程式求解，??考「」、「」。
????的?察和推?，一一?除改??法，?省不少??。此即?今世上第二快的?因?***演算法。
Euler's Totient Function
Euler's Totient Function（Euler's φ Function）
?是一?公式。?算1到n的正整??中，跟n互?（最大公因?是一）的?，?共有??。
首先要?n做?因?***：
n = p?a? × p?a? × … × p?a?
where p? … p? are primes
以?因??算Euler's Totient Function。φ?做[fai]，可以?做phi：
φ(n) = n × (1 - 1/p?) × (1 - 1/p?) × … × (1 - 1/p?)
可以?用???序?算，避免溢位：
n ÷ p? × (p?-1) ÷ p? × (p?-1) ÷ … ÷ p? × (p?-1)
如此就不用一?一?去?算最大公因?了，非常有效率！
?因?***?用?除法，?算Euler's Totient Function的????度?O(sqrtN)。?先建立??表，得加速至O(π(sqrtN))，其中π(N)是1到N的????。
φ(p) = p - 1
iff p is prime.
φ(p?) = p? - p???
iff p is prime.
φ(n × m) = φ(n) × φ(m)
iff n and m are relatively prime.
φ(n) = φ(p?a?) × … × φ(p?a?)
iff n = p?a? × … × p?a?
where p? … p? are prime.
未?改良的?法，能求出每??的?因?。?用?法?算Euler's Totient Function，????度?O(NloglogN)。
int phi[10000];
void phi_table()
phi[1] = 1;
for (int i=2; i<10000; ++i)
if (!phi[i])
for (int j=i; j<10000; j+=i)
if (!phi[j]) phi[j] =
phi[j] = phi[j] / i * (i-1);
或者，首先?用?法，?每??求得一??因?；然後?用Euler's Totient Function的性?，?施Dynamic Programming。此?得以套用?性???法，????度?O(N)。
int phi[10000];
void phi_table()
// ?法。?每??求得一??因?。O(NloglogN)。
// 亦可?用?性???法。O(N)。
for (int i=2; i=i; k--, j-=i)
if (!phi[k])
phi[j] = // i是j的其中一??因?
// ?算Euler's Totient Function。O(N)。
phi[1] = 1;
for (int i=2; i<10000; ++i)
if (!phi[i])
phi[i] = i-1;
// i不是??
int j = i / phi[i]; // i***成j*phi[i]
if (j % phi[i] == 0)// j?phi[i]不互?
phi[i] = phi[j] * phi[i];
phi[i] = phi[j] * (phi[i] - 1);
M?bius Function
M?bius Function
用排容原理求一??字的所有因??和。