Number Theoretic Function

「數論函數」。數論領域的經典函數。諸如因數數量函數σ₀(n)、因數總和函數σ₁(n)、質因數數量函數ω(n)、質數數量函數π(n)、互質數量函數φ(n)、質因數排容函數μ(n)、……。維基百科整理了一份詳細列表：

https://en.wikipedia.org/wiki/Arithmetic_function

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

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

Linear Function

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

只是順便介紹一下。由於是不同領域，加性函數撞名了。

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

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

演算法是矩陣求解。

Scalar Function

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

只是順便介紹一下。由於是不同領域，名稱不一致了。

https://math.stackexchange.com/questions/22069/

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

Coprime Counting Function φ(n)

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))，其中ω(n)是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)。

int phi[10000]; void phi_table() { phi[0] = 0; 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] = j; phi[j] -= phi[j] / i; } } void phi_table() { for (int i=0; i<10000; ++i) phi[i] = i; for (int i=2; i<10000; ++i) if (phi[i] == i) for (int j=i; j<10000; j+=i) phi[j] -= phi[j] / i; }

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

int phi[10000]; bool sieve[10000]; vector<int> prime; void phi_table() { phi[0] = 0; phi[1] = 1; for (int i=2; i<=10000; 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)。

int lpf[10000]; int phi[10000]; vector<int> prime; // O(NloglogN) void least_prime_factor_table() { for (int i=2; i<10000; ++i) if (!lpf[i]) for (int j=i; j<10000; j+=i) if (!lpf[j]) lpf[j] = i; } // O(N) void least_prime_factor_table() { for (int i=2; i<10000; i++) { if (!lpf[i]) { prime.push_back(i); lpf[i] = i; } for (int j=0; i*prime[j]<10000; 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<10000; ++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)

質因數重複為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))，其中ω(n)是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)。

int mu[10000]; bool sieve[10000]; vector<int> prime; void mu_table() { mu[0] = 0; mu[1] = 1; for (int i=2; i<=10000; 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)。

int lpf[10000]; int mu[10000]; void mu_table() { least_prime_factor_table(); mu[0] = 0; mu[1] = 1; for (int i=2; i<10000; 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

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))。預先建立質數表，降為O(π(sqrt(n)))，其中π(n)是1到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; }

數學公式：質因數分解。時間複雜度依舊，但是執行時間更短。

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(n^⅓)。

https://www.geeksforgeeks.org/count-divisors-n-on13/

計算每一項

篩法。時間複雜度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]++; }

Summatory Divisor Function Dₖ(n)

Divisor Function前綴和。前綴和與相鄰差互為反運算！

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

計算一項

時間複雜度O(sqrt(n))。

D0(n) = floor(n/1) + floor(n/2) + ... + floor(n/n)

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

時間複雜度O(sqrt(sqrt(n))) = O(n^¼)。

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

int D0(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; }

luogu P2424 P6028

Coprime Counting Function

f(n)   =   ∑   [gcd(i,n) = 1]   Coprime Counting Function
         1≤i≤n                  Euler's Totient Function

f(n)   =   ∑   [gcd(i,j) = 1]   Coprime Triangular Counting Function
        1≤i<j≤n

f(n,m) =   ∑   [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()); 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()); 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()); int sum = 0; int l, r; for (l=1; 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)   =   ∑   [gcd(i,n) = d]   GCD Distribution Function
           1≤i≤n

f(d;n)   =   ∑   [gcd(i,j) = d]   GCD Triangular Distribution Function
          1≤i<j≤n

f(d;n,m) =   ∑   [gcd(i,j) = d]   GCD Double Distribution Function
           1≤i≤n
           1≤j≤m

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)   =   ∑   gcd(i,n)   GCD Sum Function
         1≤i≤n            Pillai's Arithmetical Function

F(n)   =   ∑   gcd(i,j)   GCD Triangular Sum Function
        1≤i<j≤n

F(n,m) =   ∑   gcd(i,j)   GCD Double Sum Function
         1≤i≤n
         1≤j≤m

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

// 1¹ + 2¹ + ... + n¹ int S1(int n) { return (1 + n) * n / 2; // 上底加下底乘以高除以二 } // ∑ gcd(i,n) = ∑ d φ(n/d) = ∑ (n/d) φ(d) // 1≤i≤n d|x d|x 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; int l, r; for (l=1; 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; int l, r; for (l=1; 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

Prime Counting Function π(n)

https://min-25.hatenablog.com/entry/2018/11/11/172216

luogu P5493