constexpr ll cfloor(long double x) {
    ll ret = (ll)x;
    if (fabs(ret + 1 - x) < 5e-16l) ++ret;
    return ret;
}

class linearsieve {
    vector<bool> isprime;
    vector<int> low;
    vector<ll> f;
    vector<int> primes;
    function<ll(ll)> fp;
    size_t _size;
public:
    linearsieve(int n, function<ll(ll)> fp) : fp(fp) {
        _size = n + 1;
        isprime.resize(n + 1, true);
        low.resize(n + 1, 0);
        isprime[0] = isprime[1] = false;
        f[1] = fp(1);
        for (int i = 2; i <= n; ++i) {
            if (isprime[i]) {
                primes.push_back(i);
                low[i] = i;
                f[i] = fp(i);
            }
            for (int p : primes) {
                if (i * p > n) break;
                isprime[i * p] = false;
                if (i % p == 0) {
                    low[i * p] = low[i] * p;
                    if (low[i] == i) f[i * p] = fp(i);
                    else f[i * p] = f[i / low[i]] * f[low[i] * p];
                    break;
                } else {
                    low[i * p] = p;
                    f[i * p] = f[i] * f[p];
                }
            }
        }
    }
    bool operator[](int n) const { return isprime[n]; }
    ll operator()(ll n) const {
        if (n < _size) return f[n];
        ll ans = 1;
        for (ll pi = 0, p, pn = 0, m = cfloor(sqrtl(n)); n > 1 && pi <= primes.size(); ++pi) {
            p = primes[pi];
            if (n % p == 0) {
                while (n % p == 0) n /= p, ++pn;
                ans *= fp(cfloor(pow(p, pn)));
                m = cfloor(sqrtl(n)), pn = 0;
            } else if (p > m) {
                ans *= fp(n);
                break;
            }
        }
        return ans;
    }
    auto begin() -> decltype(primes.begin()) { return primes.begin(); }
    auto end() -> decltype(primes.end()) { return primes.end(); }
};

需要一个数论函数并预先定义$f_p(1), f_p(p), f_p(p^n)$，这些的计算需要$O(1)$

例如phi的定义为:

$$\varphi(xy)=\varphi(x)\varphi(y)\ \text{when}\ gcd(x,y)=1 $$$$\varphi(n)=\varphi(p_1^{a_1})\varphi(p_2^{a_2})...\varphi(p_k^{a_k}) $$$$\varphi(p^{a})=p^a-p^{a-1}$$

可以得到:

class phisieve {
    vector<bool> isprime;
    vector<int> phis;
    vector<int> primes;
public:
    phisieve(int n) {
        isprime.resize(n + 1, true);
        phis.resize(n + 1, 0);
        isprime[0] = isprime[1] = false;
        phis[1] = 1;
        for (int i = 2; i <= n; ++i) {
            if (isprime[i]) {
                primes.push_back(i);
                phis[i] = i - 1;
            }
            for (int p : primes) {
                if (i * p > n) break;
                isprime[i * p] = false;
                if (i % p) phis[i * p] = phis[i] * phis[p];
                else {
                    phis[i * p] = phis[i] * p;
                    break;
                }
            }
        }
    }
    bool operator[](ll n) const { return isprime[n]; }
    ll operator()(ll n) const {
        if (n < phis.size()) return phis[n];
        ll ans = 1;
        for (ll pi = 0, p, pn = 0, m = cfloor(sqrtl(n)); n > 1 && pi <= primes.size(); ++pi) {
            p = primes[pi];
            if (n % p == 0) {
                while (n % p == 0) n /= p, ++pn;
                ans *= cfloor(pow(p, pn) - pow(p, pn - 1));
                m = cfloor(sqrtl(n)), pn = 0;
            } else if (p > m) {
                ans *= n - 1;
                break;
            }
        }
        return ans;
    }
    auto begin() -> decltype(primes.begin()) { return primes.begin(); }
    auto end() -> decltype(primes.end()) { return primes.end(); }
};

注意到，之所以要用low是因为计算fp时未必能从fp的值反过来推出它是不是一个$p^a$

如果是phi这种函数，计算$\varphi(p^{a})=p^a-p^{a-1}$，对于phi来说作为素数幂是可以被确定的。如果对于单幂不那么确定的函数，采取low来标定按照线性筛的顺序未被筛掉的最小素数幂子项。

一般线性筛完全具备区间素分解的功能。其中观察phi和common的sieve过程，通过记录low或利用phi本身的性质依线性筛的顺序素分解区间的每个数的数论函数。

如果要取素分解的具体每个元素的值，即完全区间筛，可以考虑让fp作为素分解，f是一个记录素分解的表，还包括表的乘积（两个因子同时算进来）即可，并且由于f记录过low，可以不需要low完成。

样题:

Same GCDs —— codeforces