洲阁筛 & min_25筛学习笔记

洲阁筛

给定一个积性函数$F(n)$，求$sum_{i = 1}^{n}F(n)$。并且$F(n)$满足在素数和素数次幂的时候易于计算。

显然有：

$sum_{i = 1}^{n} F(n) = sum_{i = 1}^{sqrt{n}}F(i) left(sum_{sqrt{n} < pleqslant n/i, p is a prime} F(p) ight) + sum_{i = 1, i has no prime factor greater than sqrt{n}}^{n} F(i)$。

Part 1

第一部分是要求出$sum_{sqrt{n}<p leqslant n/i, p is a prime} F(p)$

假设小于等于$sqrt{n}$的素数从小到大分别是$P_1, P_2, cdots, P_{m}$。

设$f_{i, j}$表示$[1, j]$中和前$i$个素数互质的数的函数值和。

当$i = 0$时候想办法计算。

考虑当$i > 0$的时候的转移。考虑从$f_{i - 1, j}$中需要删掉最小质因子等于$P_i$的数的贡献。

所以有转移式子$f_{i, j} = f_{i - 1, j} - F(P_i) f_{i - 1, left lfloor j / P_i ight floor}$

如果$F(n)$本身不是完全积性函数，需要把它拆成若干积性函数的和。注意这一部分只用保证素数处取值相同即可。

注意到第一维约有$Oleft (frac{sqrt{n}}{log n} ight)$种取值，第二维有$Oleft (sqrt{n} ight)$种取值。暴力计算的复杂度是$Oleft (frac{n}{log n} ight )$

优化

考虑当$j < P_{i+ 1}$的时候，$f_{i, j} = 1$。考虑把条件放缩成$j < P_i$。

因此，当$P_i leqslant j < P_{i}^2$的时候$f_{i, j} = f_{i - 1, j} - F(P_i)$。

所以做法是：

• 维护一个$F(P_i)$的前缀和
• 考虑转移时的$f_{i - 1, left lfloor j / P_i ight floor}$
  • 当$j geqslant P_{i}^2$的时候暴力计算
  • 当$P_i leqslant j < P_{i}^2$的时候，需要这样的状态的时候减去一段和就行了。
  • 当$j < P_i$的时候，由于$f_{i, 0} = 0$，所以这一部分的值不会改变。

所以对于每个$j$，有效的状态只有$Oleft (frac{sqrt{j}}{log j} ight )$。

时间复杂度为：$sum_{i = 1}^{sqrt{n}} O left (frac{sqrt{i}}{log i} ight ) + sum_{i = 1}^{sqrt{n}} Oleft ( frac{sqrt{n/i}}{log {n/i}} ight)$

由于第一部分显然小于第二部分，所以只用考虑第二部分的和。

所以有：$sum_{i = 1}^{sqrt{n}}O left (frac{sqrt{n/i}}{log n - log i} ight)$

由于$log i leqslant log sqrt{n} = frac{1}{2}log n$

所以：$sum_{i = 1}^{sqrt{n}}O left (frac{sqrt{n/i}}{frac{1}{2}log n} ight) = sum_{i = 1}^{sqrt{n}}Oleft (frac{sqrt{n / i}}{log n} ight) = Oleft ( frac{n^{3/4}}{log n} ight)$

Part 2

这一部分是要求出$sum_{i = 1, i has no prime factor greater than sqrt{n}}^{n} F(i)$

假设小于等于$sqrt{n}$的素数从小到大分别是$P_1, P_2, cdots, P_{m}$。

设$g_{i, j}$表示在$[1, j]$中，质因子只包含$P_{i}, P_{i + 1}, cdots, P_{m}$的数的函数值的和。

转移考虑枚举当前质数$P_i$的指数。所以有：$g_{i, j} = g_{i + 1, j} + sum_{e = 1}^{log_{P_i} j} F(P_{i}^{e}) g_{i + 1, lfloor j / P_i^e floor}​$

对于时间复杂度，考虑对于每一个$j$的转移的枚举量，可以得到大概是$O(frac{sqrt{j}}{log j})$。(不会积分，告辞)

和Part 1一样，得到暴力计算这一部分的复杂度是$O(frac{n}{log n})$。

优化

因为质因子是从大到小枚举的，所以当$j < P_{i}$的时候，$g_{i, j} = 1$。

所以当满足$P_i leqslant j < P_{i}^2$，满足$g_{i, j} = g_{i + 1, j} + F(P_i)$。

和Part 1一样，分成三段，一段暴力转移，一段需要时前缀和，一段不需要维护。

复杂度也是同样的计算方法

/**
 * SPOJ
 * Problem#DIVCNT3
 * Accepted
 * Time: 35190ms
 * Memory: 33792k
 */
#include <bits/stdc++.h>
using namespace std;
typedef bool boolean;

template <typename T>
void pfill(T* pst, const T* ped, T val) {
    for ( ; pst != ped; *(pst++) = val);
}

#define ll long long

typedef class Euler {
    public:
        int num;
        boolean *vis;
        int *pri, *d3, *mp;
        
        Euler(int n) : num(0) {
            d3 = new int[(n + 1)];
            mp = new int[(n + 1)];
            pri = new int[(n + 1)];
            vis = new boolean[(n + 1)];
            pfill(vis, vis + n + 1, false);
            d3[1] = 1;
            for (int i = 2; i <= n; i++) {
                if (!vis[i]) {
                    d3[i] = 4, pri[num++] = i, mp[i] = 1;
                }
                for (int j = 0, _ = n / i, x; j < num && pri[j] <= _; j++) {
                    vis[x = pri[j] * i] = true;
                    if (!(i % pri[j])) {
                        d3[x] = (d3[i / mp[i]] + 3) * d3[mp[i]];
                        mp[x] = mp[i];
                        break;
                    } else {
                        mp[x] = i, d3[x] = d3[i] << 2;
                    }
                } 
            }
        }
        ~Euler() {
            delete[] mp;
            delete[] vis;
            delete[] d3;
            delete[] pri;
        }
} Euler;

const int N = 330000;

Euler *_euler;

int T;
vector<ll> ns;

ll maxn, D;

int cnt_prime;
int *pri, *d3;

inline void init() {
    scanf("%d", &T);
    for (ll x; T--; ) {
        scanf("%lld", &x);
        ns.push_back(x);
        maxn = max(maxn, x);
    }
    D = sqrt(maxn + 0.5) + 1;
    while (D * 1ll * D <= maxn)
        D++;
    _euler = new Euler(D + 23);
    cnt_prime = _euler->num;
    pri = _euler->pri, d3 = _euler->d3;
}

ll n;
int cnt_P[N];
int range0[N], range1[N];
ll f0[N], f1[N], g0[N], g1[N];

void precalc() {
    for (int i = 1; i < D; i++) {
        range0[i] = range0[i - 1];
        while (pri[range0[i]] * 1ll * pri[range0[i]] <= i) {
            range0[i]++;
        }
    }
    for (int i = D - 1; i; i--) {
        ll y = n / i;
        range1[i] = range1[i + 1];
        while (pri[range1[i]] * 1ll * pri[range1[i]] <= y) {
            range1[i]++;
        }
    }
    pfill(cnt_P, cnt_P + D, 0);
    for (int i = 0; pri[i] < D; i++) {
        cnt_P[pri[i]]++;
    }
    for (int i = 1; i < D; i++) {
        cnt_P[i] += cnt_P[i - 1];
    }
}

int nump;
ll get_value_f(int i, ll j) {
    if (j >= D) {
        int rj = n / j;
        return f1[rj] + (max(0, nump - max(range1[rj], i)) << 2);
    }
    return f0[j] + (max(0, cnt_P[j] - max(range0[j], i)) << 2);
}

void calculate_f() {
    for (nump = 0; pri[nump] < D; nump++);
    for (int i = 1; i < D; i++) {
        f0[i] = 1, f1[i] = 1;
    }
    for (int i = nump - 1; ~i; i--) {
        for (int j = 1; j < D && i < range1[j]; j++) {
            ll m = n / j, coef = 1;
            while (m /= pri[i]) {
                f1[j] += (coef += 3) * get_value_f(i + 1, m);
            }
        }
        for (int j = D - 1; j && i < range0[j]; j--) {
            ll m = j, coef = 1;
            while (m /= pri[i]) {
                f0[j] += (coef += 3) * get_value_f(i + 1, m);
            }
        }
    }
    for (int i = 1; i < D; i++) {
        f1[i] = get_value_f(0, n / i);        
    }
    for (int i = 1; i < D; i++) {
        f0[i] = get_value_f(0, i);
    }
}

ll get_value_g(int i, ll j) {
    if (j >= D) {
        int rj = n / j;
        return g1[rj] - max(0, i - range1[rj]);
    }
    return g0[j] - max(0, min(i, cnt_P[j]) - range0[j]);
}

void calculate_g() {
    for (int i = 1; i < D; i++) {
        g0[i] = i, g1[i] = n / i;
    }
    int cp = 0;
    for (int i = 0; pri[i] < D; i++, cp++) {
        for (int j = 1; j < D && i < range1[j]; j++) {
            g1[j] -= get_value_g(i, n / j / pri[i]);
        }
        for (int j = D - 1; j && i < range0[j]; j--) {
            g0[j] -= get_value_g(i, j / pri[i]);
        }
    }
    for (int i = 1; i < D; i++) {
        (g1[i] = get_value_g(cp, n / i) - 1) <<= 2;
    }
    for (int i = 1; i < D; i++) {
        (g0[i] = get_value_g(cp, i) - 1) <<= 2;
    }
}

void Solve(ll _n) {
    n = _n;
    D = sqrt(n + 0.5) + 1;
    while (D * 1ll * D <= n)
        D++;
    precalc();
    calculate_g();
    calculate_f();
    ll ans = f1[1];
    for (int i = 1; i < D; i++) {
        ans += d3[i] * g1[i];
//        cerr << d3[i] << " " << g1[i] << '
'; 
    }
    printf("%lld
", ans);
}

inline void solve() {
    for (int i = 0; i < (signed) ns.size(); i++) {
        Solve(ns[i]);
    }
}

int main() {
    init();
    solve();
    return 0;
}

Min25筛

Part 1

假设小于等于$sqrt{n}$的素数从小到大分别是$P_1, P_2, cdots, P_{m}$。

设$g_{i, j} = sum_{k = 2}^{j} [k的最小质因子大于P_i或者k是质数] F(k)$

转移式子：

$$g_{i, j} =
egin{cases}g_{i - 1, j} - F(P_i)(g_{i - 1, lfloor j / P_i floor} - g_{i - 1, P_i - 1}) & (j geqslant P_i^2) \ 
g_{i - 1, j} & (1 leqslant j < P_i^2)
end{cases}
$$

和洲阁筛类似，不过第一种情况还需要抠掉素数的贡献。

这里也要求$F$是一个完全积性函数。

我们考虑求出$1, 2, cdots, sqrt{n}, lfloor n / sqrt{n} floor, cdots, lfloor n / 2 floor, n$处的值。

Part 2

这里考虑计算$S(n, i) = sum_{i = 2}^{n}[n的最小质因子大于等于P_i]F(n)$。

那么显然答案等于$S(n, 1) + F(1)$。

这里简单粗暴一点，如果是素数，那么可以用第一部分计算的答案再减去某个前缀和得到，否则枚举这个合数的最小质因子。不难得到这样一个式子：

$$
S(n, i) = g_{m, n} - sum_{j = 1}^{i - 1}F(P_j) + sum_{j = i wedge P_j^2 leqslant n}sum_{e = 1wedge P_{j}^{e + 1}leqslant n} F(P_{j}^{e}) S(lfloor n/P_{j}^e floor, j + 1) + F(P_{j}^{e + 1})
$$

前一半是计算素数的贡献。

考虑计算合数的贡献，后一半是先枚举最小质因子，因为它是一个合数最小质因子，所以$p^2 leqslant n$

$p^{e + 1} leqslant n$是类似的原因。

时间复杂度被证明为是$O(frac{n}{Poly(log n)})$。

（似乎目前min_25在$10^{13}$范围内吊打其他筛）

/**
 * loj
 * Problem#6053
 * Accepted
 * Time: 332ms 
 * Memory: 2684k
 */
#include <bits/stdc++.h>
using namespace std;
typedef bool boolean;

#define ll long long

void exgcd(int a, int b, int& x, int& y) {
	if (!b) {
		x = 1, y = 0;
	} else {
		exgcd(b, a % b, y, x);
		y -= (a / b) * x;
	}
}

int inv(int a, int n) {
	int x, y;
	exgcd(a, n, x, y);
	return (x < 0) ? (x + n) : (x);
}

const int Mod = 1e9 + 7;

template <const int Mod = :: Mod>
class Z {
	public:
		int v;

		Z() : v(0) {	}
		Z(int x) : v(x){	}
		Z(ll x) : v(x % Mod) {	}

		friend Z operator + (const Z& a, const Z& b) {
			int x;
			return Z(((x = a.v + b.v) >= Mod) ? (x - Mod) : (x));
		}
		friend Z operator - (const Z& a, const Z& b) {
			int x;
			return Z(((x = a.v - b.v) < 0) ? (x + Mod) : (x));
		}
		friend Z operator * (const Z& a, const Z& b) {
			return Z(a.v * 1ll * b.v);
		}
		friend Z operator ~(const Z& a) {
			return inv(a.v, Mod);
		}
		friend Z operator - (const Z& a) {
			return Z(0) - a;
		}
		Z& operator += (Z b) {
			return *this = *this + b;
		}
		Z& operator -= (Z b) {
			return *this = *this - b;
		}
		Z& operator *= (Z b) {
			return *this = *this * b;
		}
		friend boolean operator == (const Z& a, const Z& b) {
			return a.v == b.v;
		} 
};

Z<> qpow(Z<> a, int p) {
	Z<> rt = Z<>(1), pa = a;
	for ( ; p; p >>= 1, pa = pa * pa) {
		if (p & 1) {
			rt = rt * pa;
		}
	}
	return rt;
}

typedef Z<> Zi;

const Zi inv2 ((Mod + 1) >>1);

const int Dmx = 1e5 + 5;

ll n;
int D;
int pnum;
int pri[Dmx], sf[Dmx];

void Euler(int n) {
	static bitset<Dmx + 23> vis;
	for (int i = 2; i <= n; i++) {
		if (!vis.test(i)) {
			pri[pnum++] = i;
		}
		for (int *p = pri, *_ = pri + pnum, x; p != _ && (x = *p * i) <= n; p++) {
			vis.set(x);
			if (!(i % *p)) {
				break;
			}
		}
	}
	for (int i = 1; i <= pnum; i++) {
		sf[i] = sf[i - 1] + (pri[i - 1] ^ 1);
	}
}

template <typename T>
class SieveArray {
	public:
		ll n;
		vector<T> a0;
		vector<T> a1;
		
		void init(ll n, int D) {
			this->n = n;
			a0.assign(D, 0);
			a1.assign(D, 0);
		}
		
		T& operator [] (ll p) {
			return (p >= (signed) a0.size()) ? a1[n / p] : a0[p]; 
		}
};

SieveArray<int> range;
SieveArray<Zi> f0, f1;

Zi S(ll n, int m) {
	if (pri[m] > n) {
		return 0;
	}
	Zi rt = f1[n] - sf[m];
	for (int j = m; j < pnum && 1ll * pri[j] * pri[j] <= n; j++) {
		int p = pri[j];
		ll nn = n, p2 = 1ll * p * p;
		for (int e = 1; p2 <= nn; e++) {
			rt += S(nn /= p , j + 1) * (p ^ e) + (p ^ (e + 1));
		}
	}
	return rt;
}

int main() {
	scanf("%lld", &n);
	D = sqrt(D + 0.5) + 1;
	while (1ll * D * D <= n)
		D++;
	Euler(D + 3);
	range.init(n, D);
	f0.init(n, D);
	f1.init(n, D);
	ll t = 0;
	range[0] = 0;
	for (int i = 1; i < D; i++) {
		while (t * t <= i)
			t++;
		range.a0[i] = t;
	}
	t = 0;
	for (int i = D; --i; ) {
		ll v = n / i;
		while (t * t <= v)
			t++;
		range.a1[i] = t;
	}
	for (int i = 1; i < D; i++) {
		f0.a0[i] = i - 1;
		f1.a0[i] = ((Zi(i) * (i + 1)) * inv2) - 1;
		ll v = n / i;
		f0.a1[i] = v - 1;
		f1.a1[i] = ((Zi(v) * (v + 1)) * inv2) - 1;
	}
	for (int i = 0; i < pnum; i++) {
		int p = pri[i];
		if (p >= D) {
			break;
		}
		ll v;
		for (int j = 1; j < D && p < range.a1[j]; j++) {
			v = n / (1ll * j * p);
			f0.a1[j] -= f0[v] - f0.a0[p - 1];
			f1.a1[j] -= (f1[v] - f1.a0[p - 1]) * p;
		}
		for (int j = D; p < range.a0[--j]; ) {
			f0.a0[j] -= f0.a0[j / p] - f0.a0[p - 1];
			f1.a0[j] -= (f1.a0[j / p] - f1.a0[p - 1]) * p;
		}
	}
	for (int i = 1; i < D; i++) {
		f1.a0[i] -= f0.a0[i];
		f1.a0[i] += (i >= 2) * 2;
		f1.a1[i] -= f0.a1[i];
		f1.a1[i] += (2 * i <= n) * 2;
	}
	Zi ans = S(n, 0) + 1;
	printf("%d
", ans.v);
	return 0;
}