arithmetic function（number theoretic function）

「算術函數」。數論領域的經典函數。舉例來說：

因數計數函數　　　　divisor counting function σ₀(n)/d(n)/τ(n)
因數總和函數　　　　divisor sum function σ₁(n)/σ(n)
最大公因數分布函數　gcd distribution function
最大公因數總和函數　gcd sum function P(n)

質數計數函數　　　　prime counting function π(n)
互質數計數函數　　　coprime counting function φ(n)
相異質因數計數函數　distinct prime factor counting function ω(n)
質數指示函數　　　　prime indicator function χprime(n)
互質數指示函數　　　Dirichlet character (mod n) χ₀⁡(k)
質因數取捨函數　　　Möbius function μ(n)

正整數冪和函數　　　power sum function sₖ(n)   (k is power)
因數冪和函數　　　　divisor function σₖ(n)     (k is power)
單位根冪和函數　　　sum of powers of roots of unity ηq(n)   (n is power)
單位原根冪和函數　　Ramanujan's sum cq(n)      (n is power)

維基百科整理了一份詳細列表：

一般來說，函數輸入是正整數，函數輸出是整數。算術函數可以視作一串整數數列。

additive function / multiplicative function

算術函數通常是「加性函數」或「乘性函數」，形成遞迴公式。

additive function:       f(a + b) = f(a) × f(b)   when a ⊥ b
multiplicative function: f(a × b) = f(a) × f(b)   when a ⊥ b

加性函數：輸入相加等同輸出相乘，當輸入互質。
乘性函數：輸入相乘等同輸出相乘，當輸入互質。

completely additive function:       f(a + b) = f(a) × f(b)
completely multiplicative function: f(a × b) = f(a) × f(b)

完全加性函數：輸入相加等同輸出相乘。
完全乘性函數：輸入相乘等同輸出相乘。

演算法是篩法，或者說是dynamic programming。

ICPC 7837

scalar multiplication function / exponential function

「倍率函數」與「指數函數」。

scalar multiplication function: f(a + b) = f(a) + f(b)
                                ⇔ f(x) = cx
exponential function:           f(a + b) = f(a) × f(b)
                                ⇔ f(x) = cˣ

倍率函數：輸入相加等同輸出相加。
　　　　　可以改寫成倍率運算。（整數系統，倍率等同乘法。）
指數函數：輸入相加等同輸出相乘。
　　　　　可以改寫成指數運算。

演算法是divide-and-conquer method。

linear function

「加性函數」與「倍性函數」同時成立，稱作「線性函數」。

這些名詞源自不同領域，名稱撞名了，意義不同了。

additive function:               f(a + b) = f(a) + f(b)
homogeneous function of order 1: f(k ⋅ x) = k ⋅ f(x)

加性函數：輸入相加等同輸出相加。
倍性函數：輸入倍率等同輸出倍率。（整數系統，倍率等同乘法。）

演算法是矩陣求解。

list of divisors

n的因數列表。

vector<int> divisors(int n) { vector<int> v; for (int d=1; d<=n; d++) if (n % d == 0) v.push_back(d); return v; }

divisor counting function σ₀(n)

n的因數數量。

int sigma0(int n) { int c = 0; for (int d=1; d<=n; d++) if (n % d == 0) c++; return c; }

divisor sum function σ₁(n)

n的因數總和。

int sigma1(int n) { int sum = 0; for (int d=1; d<=n; d++) if (n % d == 0) sum += d; return sum; }

divisor function σₖ(n)

n的因數k次方和。k=0得到因數數量，k=1得到因數總和。

int sigmak(int n, int k) { int sum = 0; for (int d=1; d<=n; d++) if (n % d == 0) sum += pow(d, k); return sum; }

性質

σ₀(p) = 2                         when p is prime
σ₀(pⁿ) = n + 1                    when p is prime
σ₀(a × b) = σ₀(a) × σ₀(b)         when a and b are coprime
σ₀(n) = σ₀(p₁n₁) × σ₀(p₂n₂) × …   when n = p₁n₁ × p₂n₂ × … is prime factorization

小寫希臘字母σ，唸作/ˈ sɪɡmə/，可以寫作sigma。

計算一項

試除法，時間複雜度O(sqrt(n))。

int sigma0(int n) { int c = 0; int d; for (d=1; d*d<n; d++) if (n % d == 0) c += 2; if (d*d == n) c++; return c; }

數學公式，時間複雜度仍是O(sqrt(n))，但是執行時間更短。

int sigma0(int n) { int c = 1; for (int p=2; p*p<=n; p++) { int expo = 0; while (n % p == 0) { n /= p; expo++; } c *= (expo + 1); } if (n > 1) c *= 2; return c; }

平方根改成立方根，時間複雜度仍是O(sqrt(n))，而且執行時間更久。弄巧成拙。負面教材。

主迴圈O(n^⅓)。最後的質數判定O(sqrt(n)) = O(n^½)。

// https://www.geeksforgeeks.org/count-divisors-n-on13/ int sigma0(int n) { int c = 1; for (int p=0; p*p*p<=n; p++) { int expo = 0; while (n % p == 0) { n /= p; expo++; } c *= (expo + 1); } if (isprime(n)) c *= 2; else if (isprimesquare(n)) c *= 3; else if (n != 1) c *= 4; return c; }

預先建立質數表，降為O(π(sqrt(n)))。π(n)是1到n的質數數量，大約是O(n/log(n))。

預先質因數分解，降為O(ω(n))。ω(n)是n的相異質因數數量，大約是O(log(n)/loglog(n))。

計算每一項

篩法。時間複雜度O(NlogN)。

int sigma0[10000]; void divisor_generation() { for (int i=1; i<=10000; i++) // 窮舉因數 for (int j=i; j<=10000; j+=i) // 窮舉整數 // divisors[j].push_back(i); // sigma1[j] += i; sigma0[j]++; }

list of GCDs

n的最大公因數列表。1到n的正整數當中，每個數與n的最大公因數。

vector<int> gcds(int n) { vector<int> v; for (int i=1; i<=n; i++) v.push_back(gcd(i, n)); return v; }

GCD distribution function

n的最大公因數列表當中，某種最大公因數的數量。

int gcd_distribution(int n, int d) { int c = 0; for (int i=1; i<=n; i++) if (gcd(i, n) == d) c++; return c; }

GCD sum function（Pillai's arithmetical function）

n的最大公因數列表當中，每個最大公因數的總和。

int gcd_sum(int n) { int sum = 0; for (int i=1; i<=n; i++) sum += gcd(i, n); return sum; }

計算一項

同除以d，問題變成coprime counting function。

int gcd_distribution(int n, int d) { return coprime_count(n/d); }

UVa 12425

計算每一項

dynamic programming。時間複雜度O(NlogN)。

int table[10000]; void gcd_distribution_table(int n) { table[0] = 0; for (int i=n; i>=1; i--) { table[i] = n / i; for (int j=i*2; j<=n; j+=i) table[i] -= table[j]; } }

coprime counting function φ(n)
（Euler's totient function）

1到n的正整數當中，跟n互質（最大公因數是一）的數，總共有幾個。

int phi(int n) { int c = 0; for (int i=1; i<=n; i++) if (gcd(i, n) == 1) c++; return c; }

性質

φ(p) = p - 1                   when p is prime
φ(pⁿ) = pⁿ - pⁿ⁻¹              when p is prime
φ(a × b) = φ(a) × φ(b)         when a and b are coprime
φ(n) = φ(p₁n₁) × φ(p₂n₂) × …   when n = p₁n₁ × p₂n₂ × … is prime factorization
φ(a × b) = φ(a) × φ(b) × gcd(a,b) ÷ φ(gcd(a,b))

小寫希臘字母φ，唸作/faɪ/，可以寫作phi。

計算一項

數學公式。不用一個一個去計算最大公因數，非常有效率！

首先是質因數分解：

n = p₁n₁ × p₂n₂ × ... × pₖnₖ   where p₁ ... pₖ are primes

以質因數計算coprime counting function：

φ(n) = n × (1 - 1/p₁) × (1 - 1/p₂) × ... × (1 - 1/pₖ)

採用這個順序計算，避免溢位：

n ÷ p₁ × (p₁-1) ÷ p₂ × (p₂-1) ÷ ... ÷ pₖ × (pₖ-1)

試除法，時間複雜度O(sqrt(n))。

預先建立質數表，降為O(π(sqrt(n)))。π(n)是1到n的質數數量，大約是O(n/log(n))。

預先質因數分解，降為O(ω(n))。ω(n)是n的相異質因數數量，大約是O(log(n)/loglog(n))。

// > O(n) int phi(int n) { int c = n; for (int p=2; p<=n; p++) if (n % p == 0 && isprime(p)) c -= c / p; return c; } // O(sqrt(n)) int phi(int n) { int c = n; for (int p=2; p*p<=n; p++) if (n % p == 0) { while (n % p == 0) n /= p; c -= c / p; } if (n > 1) c -= c / n; return c; }

UVa 10299 10179 11064

計算每一項

篩法。時間複雜度O(NloglogN)。

const int N = 10000; int phi[N]; void phi_table() { phi[0] = 0; phi[1] = 1; for (int i=2; i<N; ++i) if (!phi[i]) for (int j=i; j<N; j+=i) { if (!phi[j]) phi[j] = j; phi[j] -= phi[j] / i; } } void phi_table() { for (int i=0; i<N; ++i) phi[i] = i; for (int i=2; i<N; ++i) if (phi[i] == i) for (int j=i; j<N; j+=i) phi[j] -= phi[j] / i; }

線性篩法。時間複雜度O(N)。

const int N = 10000; int phi[N]; bool sieve[N]; vector<int> prime; void phi_table() { phi[0] = 0; phi[1] = 1; for (int i=2; i<=N; i++) { if (!sieve[i]) { prime.push_back(i); phi[i] = i - 1; } for (int p: prime) { if (i * p >= N) break; sieve[i * p] = true; if (i % p == 0) { // φ(a × b) = φ(a) × φ(b) × gcd(a,b) ÷ φ(gcd(a,b)) // gcd(i, p) = p phi[i * p] = phi[i] * p; break; } else // φ(a × b) = φ(a) × φ(b) when a ⊥ b // phi[i * p] = phi[i] * phi[p]; phi[i * p] = phi[i] * (p - 1); } } }

首先運用線性篩法，求得每個數的最小質因數。然後運用乘性函數，實施dynamic programming。時間複雜度O(N)。

const int N = 10000; int lpf[N]; int phi[N]; vector<int> prime; // O(NloglogN) void least_prime_factor_table() { for (int i=2; i<N; ++i) if (!lpf[i]) for (int j=i; j<N; j+=i) if (!lpf[j]) lpf[j] = i; } // O(N) void least_prime_factor_table() { for (int i=2; i<N; i++) { if (!lpf[i]) { prime.push_back(i); lpf[i] = i; } for (int j=0; i*prime[j]<N; j++) { lpf[i*prime[j]] = i; if (i % prime[j] == 0) break; } } } void phi_table() { least_prime_factor_table(); phi[0] = 0; phi[1] = 1; for (int i=2; i<N; ++i) if (lpf[i] == i) // i是質數 phi[i] = i-1; else // i不是質數 { int j = i / lpf[i]; // i = j × lpf[i] if (j % lpf[i] != 0) // j ⊥ lpf[i] // if (lpf[j] != lpf[i]) // j ⊥ lpf[i] // φ(a × b) = φ(a) × φ(b) when a ⊥ b phi[i] = phi[j] * (lpf[i] - 1); else // φ(a × b) = φ(a) × φ(b) × gcd(a,b) ÷ φ(gcd(a,b)) // gcd(j, lpf[i]) = lpf[i] phi[i] = phi[j] * lpf[i]; } }

UVa 10820 10990 11327 11424 11426 12493

Möbius function μ(n)

n進行質因數分解。質因數重複為0、不重複為±1。不重複的情況下，偶數個質因數為+1、奇數個質因數為-1。

性質

μ(p¹) = -1                     when p is prime
μ(p²) = μ(p³) = ... = 0        when p is prime
μ(a × b) = μ(a) × μ(b)         when a and b are coprime
μ(n) = μ(p₁n₁) × μ(p₂n₂) × …   when n = p₁n₁ × p₂n₂ × … is prime factorization

小寫希臘字母μ，唸作/mju:/，可以寫作mu。

計算一項

試除法，時間複雜度O(sqrt(n))。

預先建立質數表，降為O(π(sqrt(n)))。π(n)是1到n的質數數量，大約是O(n/log(n))。

預先質因數分解，降為O(ω(n))。ω(n)是n的相異質因數數量，大約是O(log(n)/loglog(n))。

// > O(n) int mu(int n) { if (n == 0) return 0; int c = 0; for (int p=2; p<=n; p++) if (n % p == 0 && isprime(p)) { if (n % (p * p) == 0) return 0; c++; } return (c % 2) ? -1 : 1; } // O(sqrt(n)) int mu(int n) { if (n == 0) return 0; int c = 0; for (int p=2; p*p<=n; p++) if (n % p == 0) { n /= p; if (n % p == 0) return 0; c++; } if (n > 1) c++; return (c % 2) ? -1 : 1; }

計算每一項

線性篩法。時間複雜度O(N)。

const int N = 10000; int mu[N]; bool sieve[N]; vector<int> prime; void mu_table() { mu[0] = 0; mu[1] = 1; for (int i=2; i<=N; i++) { if (!sieve[i]) { prime.push_back(i); mu[i] = -1; } for (int p: prime) { if (i * p >= N) break; sieve[i * p] = true; if (i % p == 0) { mu[i * p] = 0; break; } else // mu[i * p] = mu[i] * mu[p]; mu[i * p] = -mu[i]; } } }

首先運用線性篩法，求得每個數的最小質因數。然後運用乘性函數，實施dynamic programming。時間複雜度O(N)。

const int N = 10000; int lpf[N]; int mu[N]; void mu_table() { least_prime_factor_table(); mu[0] = 0; mu[1] = 1; for (int i=2; i<N; i++) if (lpf[i / lpf[i]] == lpf[i]) mu[i] = 0; else // mu[i] = mu[i / lpf[i]] * mu[lpf[i]]; mu[i] = -mu[i / lpf[i]]; }

ICPC 2116 luogu P4318

prime counting function π(n)

1到n的正整數當中，總共有幾個質數。

int pi(int n) { int c = 0; for (int i=1; i<=n; i++) if (isprime(i)) c++; return c; }

計算一項

Meissel–Lehmer algorithm。可能是O(n^¾ / log(n))。

Lagarias–Miller–Odlyzko algorithm。可能是O(n^⅔ / log²(n))。

luogu P3912 P7884

計算每一項

線性篩法。時間複雜度O(N)。

summatory function

「累和函數」。加性積分的結果。累積和、前綴和。

算術函數的累和函數，通常擁有快速演算法。

g = ∑ f     f的加性積分是g
f = Δ g     g的加性微分是f

g(n) = sum f(i)          f的加性積分是g
       i≤n

f(n) = g(n) - g(n-1)     g的加性微分是f

Möbius transform

「莫比烏斯變換」。乘性積分的結果。因數項的和。

算術函數的莫比烏斯變換，通常擁有特殊數學公式。

乘性卷積當中，常數函數𝟏的效果是乘性積分。

乘性卷積當中，質因數取捨函數μ的效果是乘性微分。

g = f ∗ 𝟏     f的乘性積分是g
f = g ∗ μ     g的乘性微分是f

g(n) = sum f(d)            f的乘性積分是g
       d|n

f(n) = sum g(d) μ(n/d)     g的乘性微分是f
       d|n

質因數取捨函數μ的乘性積分是脈衝函數ε。

互質數計數函數φ的乘性積分是恆等函數ι。

ε = μ ∗ 𝟏     μ的乘性積分是ε
ι = φ ∗ 𝟏     φ的乘性積分是ι

sum μ(d) = { 1 if n = 1      μ的乘性積分是ε
d|n        { 0 otherwise     為了書寫方便，交換等號兩側。

n = sum φ(d)     φ的乘性積分是ι
    d|n

常數函數𝟏的乘性積分是因數計數函數σ₀。

恆等函數ι的乘性積分是因數總和函數σ₁。

σ₀ = 𝟏 ∗ 𝟏     𝟏的乘性積分是σ₀
σ₁ = ι ∗ 𝟏     ι的乘性積分是σ₁

σ₀(n) = sum 1     𝟏的乘性積分是σ₀
        d|n

σ₁(n) = sum d     ι的乘性積分是σ₁
        d|n

順便介紹一道經典公式φ = μ ∗ ι。可以利用代數來證明。

{ ε = μ ∗ 𝟏                 μ的乘性積分是ε
{ φ ∗ 𝟏 = ι                 φ的乘性積分是ι
ε ∗ φ ∗ 𝟏 = μ ∗ 𝟏 ∗ ι       左側相卷等於右側相卷
(ε ∗ φ) ∗ 𝟏 = (μ ∗ ι) ∗ 𝟏   交換律、結合律
ε ∗ φ = μ ∗ ι               左右兩側同卷𝟏的反元素
φ = μ ∗ ι                   ε是單位元素

名詞解釋

數論領域和分析領域都有Möbius transform這個詞彙，但是意義不同。此處介紹的是數論領域的Möbius transform。

transform作為動詞，是指轉換。transform作為名詞，是指轉換結果。transformation從動詞變名詞，是指轉換這件事情。

類似詞彙還有inverse與inversion、count與counting、sort與sorting。

power sum function sₖ(n)

「冪和函數」。正整數1到n，各自k次方的總和。

計算一項：時間複雜度O(n)。

// 1¹ + 2¹ + ... + n¹ int s1(int n) { int sum = 0; for (int i=1; i<=n; i++) sum += i; return sum; }

計算一項：時間複雜度O(1)。

// 1¹ + 2¹ + ... + n¹ int s1(int n) { return (1 + n) * n / 2; // 上底加下底乘以高除以二 }

list of quotients【尚無正式名稱】

「商數列表」。quo(n,1) quo(n,2) ... quo(n,n)。

正整數n，分別除以1到n，所得到的商數們。

n | list of quotients | formula
--| ------------------| --------------------------------------------
1 | 1                 | quo(1,1)
2 | 2 1               | quo(2,1) quo(2,2)
3 | 3 1 1             | quo(3,1) quo(3,2) quo(3,3)
4 | 4 2 1 1           | quo(4,1) quo(4,2) quo(4,3) quo(4,4)
5 | 5 2 1 1 1         | quo(5,1) quo(5,2) quo(5,3) quo(5,4) quo(5,5)

n的商數列表最多2sqrt(n)種數值。

分成兩種情況：
一、除數小於等於sqrt(n)：除數總共sqrt(n)種，於是商數最多sqrt(n)種。
二、除數大於sqrt(n)：商數小於sqrt(n)，於是商數最多sqrt(n)種。

求得每種數值所在位置，時間複雜度O(sqrt(n))。

中文網路稱作「整除分块」。

// quo(n,1) quo(n,2) ... quo(n,n) // floor(n/1) floor(n/2) ... floor(n/n) void list_of_quotients(int n) { for (int l=1, r; l<=n; l=r+1) { r = n / (n / l); cout << (n / l) << ": [" << l << ',' << r << ']'; } }

quotient sum function Q(n)【OEIS A006218】

「商數總和函數」。商數列表總和。

計算一項：時間複雜度O(n)。

// quo(n,1) + quo(n,2) + ... + quo(n,n) // floor(n/1) + floor(n/2) + ... + floor(n/n) int Q(int n) { int sum = 0; for (int i=1; i<=n; i++) sum += n / i; return sum; }

計算一項：援引方才的演算法，時間複雜度O(sqrt(n))。

int Q(int n) { int sum = 0; for (int l=1, r; l<=n; l=r+1) { r = n / (n / l); sum += (r - l + 1) * (n / l); } return sum; }

計算一項：時間複雜度O(sqrt(n))。

Q(n) = (floor(n/1) + ... + floor(n/r)) × 2 - r²
where r = floor(sqrt(n))

int Q(int n) { int sum = 0; int l, r; for (l=1; l*l<=n; l=r+1) { r = n / (n / l); sum += (r - l + 1) * (n / l); } return sum * 2 - r * r; }

計算一項：凸函數曲線下方格點，時間複雜度可能是O(n^⅓)。

list of remainders【尚無正式名稱】

「餘數列表」。rem(n,1) rem(n,2) ... rem(n,n)。

正整數n，分別除以1到n，所得到的餘數們。

n | list of remainders | formula
--| -------------------| --------------------------------------------
1 | 0                  | rem(1,1)
2 | 0 0                | rem(2,1) rem(2,2)
3 | 0 1 0              | rem(3,1) rem(3,2) rem(3,3)
4 | 0 0 1 0            | rem(4,1) rem(4,2) rem(4,3) rem(4,4)
5 | 0 1 2 1 0          | rem(5,1) rem(5,2) rem(5,3) rem(5,4) rem(5,5)

remainder sum function R(n)【OEIS A004125】

「餘數總和函數」。餘數列表總和。

計算一項：時間複雜度O(n)。

// rem(n,1) + rem(n,2) + ... + rem(n,n) int R(int n) { int sum = 0; for (int i=1; i<=n; i++) sum += n % i; return sum; }

餘數不易計算rem(n,1) rem(n,2) ... rem(n,n)，整除部分容易計算n-rem(n,1) n-rem(n,2) ... n-rem(n,n)。此即商數列表的變種1⋅quo(n,1) 2⋅quo(n,2) ... n⋅quo(n,n)。其總和扣掉n²再變號，就是正確答案。

計算一項：援引方才的演算法，時間複雜度O(sqrt(n))。

// rem(n,1) + rem(n,2) + ... + rem(n,n) int R(int n) { int sum = 0; for (int l=1, r; l<=n; l=r+1) { r = n / (n / l); sum += (r - l + 1) * (n / l) * (l + r) / 2; } return n*n - sum; }

計算一項：凸函數曲線下方格點，時間複雜度可能是O(n^⅓)。

summatory divisor counting function D₀(n)

計算一項：時間複雜度O(n sqrt(n))。

int D0(int n) { int sum = 0; for (int i=1; i<=n; i++) sum += sigma0(i); return sum; }

商數總和函數、因數計數函數累積和，恰好相等。

證明方式：建立因數列表、套用篩法。

計算一項：時間複雜度O(sqrt(n))。

// D0(n) = Q(n) int D0(int n) { int sum = 0; for (int l=1, r; l<=n; l=r+1) { r = n / (n / l); sum += (r - l + 1) * (n / l); } return sum; }

因數計數函數，得以快速計算一項。

int sigma0(int n) { if (n == 1) return D0(1); else return D0(n) - D0(n-1); }

luogu P2424 P6028

summatory divisor sum function D₁(n)

計算一項：時間複雜度O(n sqrt(n))。

int D1(int n) { int sum = 0; for (int i=1; i<=n; i++) sum += sigma1(i); return sum; }

餘數總和函數、因數總和函數累積和，恰好互補，相加是n²。

計算一項：時間複雜度O(sqrt(n))。

int D1(int n) { int sum = 0; for (int l=1, r; l<=n; l=r+1) { r = n / (n / l); sum += (r - l + 1) * (n / l) * (l + r) / 2; } return sum; }

因數總和函數，得以快速計算一項。

int sigma1(int n) { if (n == 1) return D1(1); else return D1(n) - D1(n-1); }

coprime counting function

f(n)   =  sum  [gcd(i,n) = 1]   coprime counting function
         1≤i≤n                  Euler's totient function

f(n)   =  sum  [gcd(i,j) = 1]   coprime triangular counting function
        1≤i<j≤n

f(n,m) =  sum  [gcd(i,j) = 1]   coprime double counting function
         1≤i≤n
         1≤j≤m

預處理O(n)、計算一項O(1)。只能處理n = m的情況。

int coprime_triangular_count(int n) { if (n == 0) return 0; vector<int> phi = ::phi(n); partial_sum(phi.begin(), phi.end(), phi.begin()); return phi[n] - 1; // remove gcd(1,1) } int coprime_double_count(int n) { if (n == 0) return 0; vector<int> phi = ::phi(n); partial_sum(phi.begin(), phi.end(), phi.begin()); return phi[n] * 2 - 1; // remove duplicate gcd(1,1) }

預處理O(n)、計算一項O(sqrt(n))。

int coprime_double_count(int n, int m) { if (n > m) swap(n, m); vector<int> mu = ::mu(n); partial_sum(mu.begin(), mu.end(), mu.begin()); int sum = 0; for (int l=1, r; l<=n; l=r+1) { r = min(n/(n/l), m/(m/l)); sum += (mu[r] - mu[l-1]) * (n/l) * (m/l); } return sum; }

沒有預處理、計算一項O(n log(n))。

int coprime_count(int n) { if (n == 0) return 0; // GCD distribution function vector<int> f(n+1); for (int i=n; i>=1; i--) { f[i] = n / i; for (int j=i*2; j<=n; j+=i) f[i] -= f[j]; } return f[1]; } int coprime_double_count(int n, int m) { if (n == 0 || m == 0) return 0; if (n > m) swap(n, m); // GCD double distribution function vector<int> f(n+1); for (int i=n; i>=1; i--) { f[i] = (long long)(n/i) * (m/i); for (int j=i*2; j<=n; j+=i) f[i] -= f[j]; } return f[1]; }

luogu P2158

GCD distribution function

f(d;n)   =  sum  [gcd(i,n) = d]   GCD distribution function
           1≤i≤n

f(d;n)   =  sum  [gcd(i,j) = d]   GCD triangular distribution function
          1≤i<j≤n

f(d;n,m) =  sum  [gcd(i,j) = d]   GCD double distribution function
           1≤i≤n
           1≤j≤m

int gcd_distribution(int n, int d) { return coprime_count(n/d); } int gcd_double_distribution(int n, int m, int d) { return coprime_double_count(n/d, m/d); }

luogu P4450 P2522 P2568 P2257 P3172

GCD sum function

F(n)   =  sum  gcd(i,n)   GCD sum function
         1≤i≤n            Pillai's arithmetical function

F(n)   =  sum  gcd(i,j)   GCD triangular sum function
        1≤i<j≤n

F(n,m) =  sum  gcd(i,j)   GCD double sum function
         1≤i≤n
         1≤j≤m

拆開成乘性卷積。然後整除分塊。

F(n) = sum gcd(i,n) = sum d φ(n/d) = sum (n/d) φ(d)
      1≤i≤n           d|n            d|n

預處理O(n)、計算一項O(sqrt(n))。

int gcd_sum(int n) { vector<int> phi = ::phi(n); int sum = 0; int d; for (d=1; d*d<n; d++) { sum += d * phi[n/d]; sum += n/d * phi[d]; } if (d*d == n) sum += d * phi[d]; return sum; } int gcd_triangular_sum(int n) { vector<int> phi = ::phi(n); partial_sum(phi.begin(), phi.end(), phi.begin()); int sum = 0; for (int l=1, r; l<=n; l=r+1) { r = n / (n / l); sum += (phi[n/l] - 1) * (s1(r) - s1(l-1)); } return sum; } int gcd_double_sum(int n) { vector<int> phi = ::phi(n); partial_sum(phi.begin(), phi.end(), phi.begin()); int sum = 0; for (int l=1, r; l<=n; l=r+1) { r = n / (n / l); sum += (phi[r] - phi[l-1]) * (n/l) * (n/l); } return sum; }

沒有預處理、計算一項O(n log(n))。

int gcd_sum(int n) { vector<int> f(n+1); int sum = 0; for (int i=n; i>=1; i--) { f[i] = n / i; for (int j=i*2; j<=n; j+=i) f[i] -= f[j]; sum += f[i] * i; } return sum; } int gcd_double_sum(int n, int m) { vector<int> f(n+1); int sum = 0; for (int i=n; i>0; --i) { f[i] = (n/i) * (m/i); for (int j=i*2; j<=n; j+=i) f[i] -= f[j]; sum += f[i] * i; } return sum; } int gcd_triangular_sum(int n) { // remove diagonal gcd(i,i) = i return (gcd_double_sum(n) - s1(n)) / 2; } int gcd_double_sum(int n) { // add diagonal gcd(i,i) = i return gcd_triangular_sum(n) * 2 + s1(n); }

UVa 11417 11424 11426 luogu P1390 P2398 P3768

summatory multiplicative convolution

乘性卷積的累和，擁有特殊演算法，快速計算一項。

sum (a∗b)(i) = sum sum a(d)b(i/d)
i≤n            i≤n d|i

sum (a∗b)(i) = sum   sum  a(d)b(i) = sum a(d)  sum  b(i)
i≤n            d≤n i≤⌊n/d⌋           d≤n     i≤⌊n/d⌋

sum (a∗b)(i) = sum a(d)B(⌊n/d⌋)
i≤n            d≤n

已知數列a與其累和A，已知數列b與其累和B，求得乘性卷積的累和sum (a∗b)。

一、普通做法，先乘性卷積再累和，必須計算每一項，時間複雜度O(n sqrt(n))。

因數成雙成對，窮舉一半足矣。

二、如果a累和與b累和，可以快速計算一項，需時O(f(n))與O(g(n))，那麼運用整除分塊降為O((f(n) + g(n)) sqrt(n))。

運用整除分塊，分成2sqrt(n)塊，每段計算一項。

summatory multiplicative function

乘性函數的累和，擁有特殊演算法，快速計算一項。

sum (a∗b)(i) = sum a(d)B(⌊n/d⌋)
i≤n            d≤n

sum (a∗b)(i) = a(1)B(⌊n/1⌋) + sum a(d)B(⌊n/d⌋)
i≤n                          2≤d≤n

B(n) = (sum (a∗b)(i) - sum a(d)B(⌊n/d⌋)) / a(1)
        i≤n           2≤d≤n

已知數列a與其累和A，已知sum (a∗b)，已知數列b恰是乘性函數，求得其累和B。

一、普通做法，先線性篩法再累和，必須計算每一項，時間複雜度O(n)。

二、尋找適當的數列a。如果a累和與a∗b累和，可以快速計算一項，假定是O(1)，那麼運用整除分塊降為O(n^¾)。

沿用上個小節的推導結果。運用分部積分，拆出第一項，形成b累和的遞迴公式。運用top-down動態規劃，計算遞迴公式。運用整除分塊，處理遞迴公式當中的乘性卷積。事先推導a累和與a∗b累和的數學公式，以便快速計算一項。

快速計算一項，假定是O(f(n))與O(g(n))，求得時間複雜度O((f(n) + g(n)) n^¾)。假定是O(1)，方便求得時間複雜度O(n^¾)。

三、又b是乘性函數，那麼再運用預計算降為O(n^⅔)。

先線性篩法再累和，預先計算B的前m項。時間複雜度成為O(m + n/sqrt(m))。令m = n/sqrt(m)讓時間複雜度達到最低。

中文網路稱作「杜教筛」。

summatory coprime counting function

回到正題。φ的累和，有兩種演算法：

一、乘性卷積的累和（正向）：採用公式φ = μ ∗ ι。

sum φ(i) = sum (μ∗ι)(i) = sum μ(d)I(⌊n/d⌋)
i≤n        i≤n            d≤n

二、乘性函數的累和（反向）：採用公式ι = 𝟏 ∗ φ。

Φ(n) = (sum ι(i) - sum 𝟏(d)Φ(⌊n/d⌋)) / 𝟏(1)
        i≤n       2≤d≤n

反向的優點是不必另外計算μ。

計算一項：時間複雜度O(n^⅔)。

vector<long long> summatory(vector<int>& a) { vector<long long> b(a.size()); b[0] = 0; for (int i=1; i<(int)a.size(); i++) b[i] = b[i-1] + a[i]; return b; } long long sumphi(int n) { static const int N⅔ = 5e6; static vector<long long> table = summatory(phi_table(N⅔)); if (n <= N⅔) return table[n]; static unordered_map<int,long long> memo; if (memo.count(n)) return memo[n]; long long sum = s1(n); for (long long l=2, r; l<=n; l=r+1) { r = n / (n / l); sum -= (r - l + 1) * sumphi(n / l); } return memo[n] = sum; }

luogu P4213

summatory Möbius function

反向。公式ε = 𝟏 ∗ μ。

計算一項：時間複雜度O(n^⅔)。

long long summu(int n) { static const int N⅔ = 5e6; static vector<int> table = summatory(mu_table(N⅔)); if (n <= N⅔) return table[n]; static unordered_map<int,long long> memo; if (memo.count(n)) return memo[n]; long long sum = 1; for (long long l=2, r; l<=n; l=r+1) { r = n / (n / l); sum -= (r - l + 1) * summu(n / l); } return memo[n] = sum; }

summatory prime counting function

luogu P5493

coprime counting function

[totient multiplicative]
φ(mn) = φ(m)φ(n)
與m互質的數a  與n互質的數b 與mn互質的數x
{ x = a (mod m)     m個m個排整齊剛好是一堆column
{ x = b (mod n)     n個n個排整齊剛好是一堆column
窮舉a與b實施中國餘數定理得到x

[totient divisor sum]
sum φ(d) = n
d|n
1/20 2/20 .... 20/20
約分之後，分母只可能是因數
分母是20的有φ(20)個
分母是10的有φ(10)個
分母是5的有φ(5)個
分母是4 2 1 有 φ(4) φ(2) φ(1)個

[totient gcd]
 sum gcd(i,n) = sum φ(n) * n/d
i=1⋯n           d|n
gcd數列累積和   φ(n)數列只取因數項
原理同[euler divisor sum]  n/d是gcd的值

[totient 取捨]
φ(n) = sum μ(d) * n/d = n sum μ(d)/d
       d|n                d|n
順帶一提，隨便一個poset都可以定義取捨，μ(n)是定在因數關係的poset。

[Dirichlet product]
sum f(d) * g(n/d)
d|n

[滴哩卷積的計算]
令𝟏(n)=1，補到原本的等式裡面，弄成卷積的樣子
sum μ(d) * 𝟏(n/d) = { 1 when n = 1    [möbius divisor]
d|n                 { 0 otherwise
sum φ(d) * 𝟏(n/d) = n　　             [totient divisor]
d|n
sum μ(d) * n/d = φ(n)　　　           [totient 取捨]
d|n

μ(n)和𝟏(n)=1的卷積是脈衝函數。μ(n)的倒數是水平線𝟏(n)=1。非常奇葩。
φ(n)和𝟏(n)=1的卷積是45度斜線ι(n)=n。也很奇葩。
μ(n)和ι(n)=n的卷積是φ(n)。

[滴哩卷積的計算]
  ι    ι
μ -> φ -> gcd累積和

設函數ι(n)=n
 sum gcd(i,n) = φ(n) * ι(n)　　[euler gcd]
i=1⋯n
 sum f(gcd(i,n)) = φ(n) * f(n)
i=1⋯n
https://math.stackexchange.com/questions/135351/

如何簡單表示ι(n)的倒數?
http://math.stackexchange.com/questions/423317/
http://oeis.org/A055615

如何簡單表示φ(n)的倒數?
https://math.stackexchange.com/questions/3073401/
http://oeis.org/A023900

sum of powers of roots of unity

單位根，當作自然數來用。

單位根加總，當作生成函數來用。

小寫希臘字母η，唸作/ˈ i:tə/，可以寫作eta。

Ramanujan's sum

[因數指示函數，中獎結果從1改成因數本身]

                nk
η_q(n) = sum w_q      n個複弦波，n是倍率，q是取樣數，k是索引值，通通加起來
        k=1⋯q         首先定參一下取樣數是多少

       = { 0  if q∤n  宛如餘數加法循環，各數(各等分角)剛好出現一遍，加起來是0
         { q  if q|n  取樣通通抽中0度角=1。

[拉馬努金和]

                nk
c_q(n) = sum w_q        針對每一種波，取樣時取不可約分(互質)的那些點，φ(q)個
         k⊥q

         if q is prime  [totient φ(q)=q-1 if q is prime 的推廣]
       = { -1   if q∤n  全部的取樣和 - 可以約分的取樣和 = 不可約分的取樣和
         { q-1  if q|n      0或q         1(只有第n點)         -1或q-1

[乘性積分/乘性微分]

η_q(n) = sum c_d(n)           乘性積分，原理是[totient divisor sum]
         d|q                  每種因數d，用d次單位原根，湊滿q次單位根

sum η_d(n) * μ(q/d) = c_q(n)  乘性微分
d|q
sum d      * μ(q/d) = c_q(n)  照η的公式，等於0的就刪掉，剩下有因數的
d|gcd(q,n)
sum μ(q/d)   * d    = c_q(n)  調整一下寫法
sum μ(d)     * q/d  = c_q(n)  調整一下寫法
　　　　　　　　　　　　　　　　
[是因數的，滴哩卷積/取捨後，得到互質的]

φ(q) = c_q(0)           平坦波，取樣必定都是0度角=1。照定義來，取小於n的互質數
μ(q) = c_q(1)           一倍波。這很奇葩，需要想一下。
c_q(0) = sum c_d(1) * q/d   [totient 取捨]
         d|q
c_q(n) = sum c_d(1) * q/d   前面推導的式子 (n=0是特例，剛好是totient取捨)
        d|gcd(q,n)          這可以DP：比較小的取樣結果，遞推比較大的取樣結果

[傅立葉轉換]

滴哩卷積下，μ(n)和φ(d)和g(n)=1和g(n)=n經過某種傅立葉轉換之後不知道是啥鬼。

循環卷積下，gcd數列dft的第一項(不是第零項)是φ(n)
　　　　　　gcd數列與一倍頻率複弦波的內積是φ(n)
不知道怎麼證

[傅立葉轉換，針對特例gcd]

注意到沒有第零項
                     ij
gcd(i,n) = 1/n sum (w  * sum c_d(j) * n/d )  n是數量，i是索引值
              j=1⋯n      d|n

                 ij
         = sum (w  * sum c_d(j)/d )  外面1/n跟裡面n相消
          j=1⋯n      d|n

注意到j=1的那一項
sum c_d(1) * n/d = c_d(0) = φ(n)
d|n

還可以算複數系統的最大公因數???

list of GCDs的傅立葉轉換是coprime counting function

φ(n) =  sum  { gcd(n, k) ÷ e𝑖⋅2π⋅(k/N)⋅n }
      k=0⋯N-1

《The Fourier transform of functions of the greatest common divisor》

輸入不是一維數列而是二維陣列。嚴格來說這不是傅立葉轉換，也不存在快速傅立葉轉換。

coprime counting function存在快速演算法。這個轉換似乎派不上用場。