多是一篇（抄來的）min25學習筆記

時間  2019-11-24
標籤 一篇 來的 min25 min 學習 筆記 简体版
原文   原文鏈接

多是一篇（抄來的）min25學習筆記

一個要求不少的積性函數

咱們考慮有一個積性函數，這個函數知足能夠快速計算質數處的值html

且質數能夠寫成一個多項式的形式……並且這個多項式若是強行套在合數上，知足積性，我也不知道有沒有除了\(x^{k}\)別的多項式惹node

假如\(F(x) = x^{k}\)吧c++

咱們想要計算這個東西數組

\(g(n,j)\)表示前\(n\)個數裏，質數的和，加上合數中最小質因子大於\(P_{j}\)的和函數

那麼，怎麼求呢oop

咱們考慮已經求好\(g(n,j - 1)\)這個數組學習

那麼若是\(P_{j}^{2} > n\)spa

那麼很妙，就是\(g(n,j) = g(n,j - 1)\)，由於沒有以\(P_{j}\)爲最小質因子的合數code

若是\(P_{j}^{2} \leq N\)htm

那麼考慮從\(g(n,j - 1)\)到\(g(n,j)\)減小了什麼

例如，我顯然減小了合數中最小質因子爲\(P_{j}\)的值，由於定義裏是大於，\(g(n,j)\)中的合數質因子至少是\(P_{j + 1}\)

那麼最小質因子爲\(P_{j}\)的值怎麼算呢

咱們先要提出一個\(P_{j}\)出來，而後在剩下的\(\lfloor \frac{n}{P_{j}} \rfloor\)中找出全部最小質因子大於等於\(P_{j}\)的個數

好像\(g(n,j - 1)\)就很不錯！然而其中有一些小於\(P_{j}\)的質數，考慮把它們扣掉，這個能夠預處理，由於\(P_{j}\)的範圍在\(\sqrt{N}\)之內

記\(sum[j] = \sum_{i = 1}^{j} F(P_{i})\)

那麼咱們有
\[ g(n,j) =\left\{\begin{matrix} g(n,j - 1) & P_{j}^{2} > n\\ g(n,j - 1) - F(P_{j})(g(\lfloor \frac{n}{P_{j}} \rfloor,j - 1) - sum(j - 1)) & P_{j}^{2} \leq n \end{matrix}\right. \]

同時咱們發現，咱們只須要用到\(\lfloor \frac{n}{i} \rfloor\)的取值，一共\(2\sqrt{N}\)個位置，並且質數的話也只須要前\(\sqrt{N}\)個

初始化的話，對於\(g(n,0)\)，咱們把這個\(F(x)\)多項式套在每一個數上面，求一個\(\sum_{i = 2}^{n} x^{i}\)這個在k固定的時候有k + 1次多項式的公式能夠快速計算

通常不會太大吧，若是拉格朗日插值的話可能和\(\sqrt{n}\)相乘會有點涼

g的用處

求完\(g\)就有很神奇的事情發生了！若是你認爲\(F(x) = 1\)的話，咱們甚至能夠求出\(n\)之內質數的個數！在我不知道這些質數是什麼的狀況下！這簡直是我這個菜雞之前想都不敢想的事

不過，咱們的征途是星辰大海1到n的前綴和

那麼咱們相似的定義一個\(S(n,j)\)表示前\(n\)個數裏，全部數的最小的質因子大於等於\(P_{j}\)的方案數（注意這裏就不包括小於\(P_{j}\)的質數了）

咱們先搞定\(S(n,j)\)裏的全部質數的\(F(x)\)的和

方法是\(g(n,|P|) - sum[j - 1]\)

那麼咱們枚舉每一個數的最小質因子

直接列個式子就是
\[ S(n,j) = g(n,|P|) - sum[j - 1] + \sum_{k = j}^{P_{k}^{2} \leq n}\sum_{e = 1}^{P_{k}^{e + 1} \leq n} [F(P_{k}^{e}) S(\lfloor \frac{n}{P_{k}^{e}}\rfloor,j + 1) + F(P_{k}^{e + 1})] \]
意義就是咱們枚舉一個數的最小質因子是\(P_{k}^{e}\)，而後乘上\(\lfloor\frac{n}{P_{k}^{e}}\rfloor\)中最小質因子大於\(P_{j + 1}\)的

可是全部的\(F(P_{k}^{e}),e\geq 2\)都會漏掉，因此要重算一下

這樣的話複雜度就是……我抄的那我的沒寫……好像是\(O(\frac{n^{3/4}}{\log n})\)

好吧，這不重要……畢竟我……也不會證複雜度

這麼看起來彷彿挺簡單的好難

幾道例題

LOJ#6053. 簡單的函數

恍然以爲WC好像見過這道題……咋就變成板子了……板子都不會我也要廢了qwq

考慮質數的函數值怎麼表示，因爲只有2是偶數，除了2之外剩下的質數值都是\(f(p) = p - 1\)，而後\(f(2) = 3\)

這個時候咱們要大膽想象，咱們就認爲全部質數都是\(f(p) = p - 1\)，遇到2和1這個難受位置特判一下就行了

可是這樣不知足積性

不過一個多項式嘛，我把每一塊都拆開不就知足積性了嗎

設\(F(x) = x\)，統計一個\(g(n,|P|)\)

再來一個\(F(x) = 1\)，統計一個\(h(n,|P|)\)

能夠同時進行，方法都是同樣的

因而咱們的\(S(n,j)\)初始的時候統計質數要加上\(g(n,|P|) - h(n,|P|) - \sum_{i = 1}^{j - 1}(P_{i} - 1)\)

若是\(j = 1\)的話，咱們把2算成了1，那就再加上2

因而照着剛纔的討論，問題就完美解決了

最後加上1處的函數值是1

#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int,int>
#define mp make_pair
#define pb push_back
#define space putchar(' ')
#define enter putchar('\n')
#define eps 1e-10
#define ba 47
#define MAXN 200005
//#define ivorysi
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
template<class T>
    void read(T &res) {
    res = 0;T f = 1;char c = getchar();
    while(c < '0' || c > '9') {
        if(c == '-') f = -1;
        c = getchar();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 +c - '0';
        c = getchar();
    }
    res *= f;
}
template<class T>
    void out(T x) {
    if(x < 0) {x = -x;putchar('-');}
    if(x >= 10) {
        out(x / 10);
    }
    putchar('0' + x % 10);
}
const int MOD = 1000000007;
int inc(int a,int b) {
    return a + b >= MOD ? a + b - MOD : a + b;
}
int mul(int a,int b) {
    return 1LL * a * b % MOD;
}
int fpow(int x,int c) {
    int res = 1,t = x;
    while(c) {
        if(c & 1) res = mul(res,t);
        t = mul(t,t);
        c >>= 1;
    }
    return res;
}
void update(int &x,int y) {
    x = inc(x,y);
}
int64 N,pos[MAXN],sqr;
int prime[MAXN],tot,sum[MAXN];
bool nonprime[MAXN];
int cnt,g[MAXN],h[MAXN];
int id[2][MAXN];
int getsum(int64 x) {
    int u = x % MOD;
    int res = 1LL * u * (u + 1) % MOD;
    res = mul(res,(MOD + 1) / 2);
    return res;
}
void process(int64 S) {
    for(int i = 2 ; i <= S ; ++i) {
        if(!nonprime[i]) {
            prime[++tot] = i;
            sum[tot] = sum[tot - 1];
            update(sum[tot],i); 
        }
        for(int j = 1 ; j <= tot ; ++j) {
            if(prime[j] > S / i) break;
            nonprime[i * prime[j]] = 1;
            if(i % prime[j] == 0) break;
        }
    }
}
int S(int64 x,int y) {
    if(x <= 1 || prime[y] > x) return 0;
    int res = 0;
    if(y == 1) res += 2;
    int k = x <= sqr ? id[0][x] : id[1][N / x];
    update(res,inc(g[k],MOD - inc(inc(sum[y - 1],MOD - y + 1),h[k])));
    for(int j = y ; j <= tot ; ++j) {
        int64 t1 = prime[j],t2 = 1LL * prime[j] * prime[j];
        if(t2 > x) break;
        for(int e = 1 ; t2 <= x ; t1 = t2,t2 = t2 * prime[j],++e) {
            int t = inc(mul(prime[j] ^ e,S(x / t1,j + 1)),prime[j] ^ (e + 1));
            update(res,t);
        }
    }
    return res;
}
int main(){
    #ifdef ivorysi
    freopen("f1.in","r",stdin);
    #endif
    read(N);
    sqr = sqrt(N);
    process(sqr);
    for(int64 i = 1 ; i <= N ; ++i) {
        int64 r = N / (N / i);
        pos[++cnt] = N / i;
        h[cnt] = (N / i - 1) % MOD;
        g[cnt] = getsum(N / i);update(g[cnt],MOD - 1);
        if(N / i <= sqr) id[0][N / i] = cnt;
        else id[1][N / (N / i)] = cnt;
        i = r;
    }
    for(int j = 1 ; j <= tot ; ++j) {
        for(int i = 1 ; i <= cnt && 1LL * prime[j] * prime[j] <= pos[i] ; ++i) {
            int k = (pos[i] / prime[j]) <= sqr ? id[0][pos[i] / prime[j]] : id[1][N / (pos[i] / prime[j])];
            int t = inc(g[k],MOD - sum[j - 1]);
            t = mul(t,prime[j]);
            update(g[i],MOD - t);
            update(h[i],MOD - inc(h[k],MOD - (j - 1)));
        }
    }
    int ans = inc(S(N,1),1);
    out(ans);enter;
    return 0;
}

LOJ#572. 「LibreOJ Round #11」Misaka Network 與求和

當時我打的時候不會min25

一年過去了

我居然如今才學min25？！

這個，好像是相似min25的一種求和方法

首先咱們熟練的把gcd轉換成互質數對
\[ ans = \sum_{d = 1}^{n} f(d)^{k} \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor}\sum_{t|i,j}\mu(t) \\ \sum_{d = 1}^{n}f(d)^{k} \sum_{t = 1}^{\lfloor \frac{n}{d} \rfloor}\mu(t) \lfloor \frac{n}{dt} \rfloor^{2} \\ \sum_{T = 1}^{n} \lfloor \frac{n}{T} \rfloor^{2} \sum_{d | T}f(d)\mu(\frac{T}{d}) \]
這個時候咱們發現要求的就是\(f * \mu(T)\)

然而……

回憶起學過的杜教篩，咱們能夠捲上一個\(I(x) = 1\)

\(f * \mu * I = f * \varepsilon = f\)

而後套到杜教篩裏

能夠獲得
\[ S(n) = \sum_{i = 1}^{n}f(n) - \sum_{t = 2}^{n}S(\lfloor \frac{n}{t}\rfloor) \]
而後咱們要作的就是求\(f(n)\)的前綴和了

用一種相似min25的方法

設\(s(n,j)\)爲\([1,n]\)裏面，質數大於等於\(j - 1\)的和，\(f(n) = S(n,1)\)，默認\(prime[0] = 1\)

再枚舉第二大的質數的時候才統計貢獻

考慮min25的時候，咱們分了質數，和合數兩部分統計

若是是質數的部分，統計了\([1,\lfloor \frac{n}{prime[j - 1]} \rfloor]\)中，大於\(prime[j - 1]\)的一個質數，咱們用\(prime[j - 1]\)乘上這個質數，這個數第二大的質因子即爲\(prime[j - 1]\)

若是是合數的部分，從\(j\)到\(tot\)枚舉最小質數的時候，顯然前面加的是已經統計完次大質因子貢獻的合數，因此係數是1，惟一漏掉的是\(prime[j]^{e + 1}\)，加上便可

列一個式子就是
\[ s(n,j) = prime[j - 1]^{k} * (g(n,j) - j + 1)\\+\sum_{j = 1}^{prime[j]^{2} \leq n} \sum_{e = 1}^{prime[j]^{e + 1} \leq n}s(\lfloor \frac{n}{prime[j]^{e}}\rfloor,j + 1) + prime[j]^{k} \]
unordered_map比手寫哈希慢了兩倍不止……感受不行……

#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int, int>
#define mp make_pair
#define pb push_back
#define space putchar(' ')
#define enter putchar('\n')
#define eps 1e-10
#define ba 47
#define MAXN 200005
//#define ivorysi
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
template <class T>
void read(T &res) {
    res = 0;
    T f = 1;
    char c = getchar();
    while (c < '0' || c > '9') {
        if (c == '-')
            f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = getchar();
    }
    res *= f;
}
template <class T>
void out(T x) {
    if (x < 0) {
        x = -x;
        putchar('-');
    }
    if (x >= 10) {
        out(x / 10);
    }
    putchar('0' + x % 10);
}

u32 N, sqr, K;
u32 prime[MAXN], tot, h[MAXN], pos[MAXN];
u32 id[2][MAXN], cnt, pw[MAXN], rec;
bool nonprime[MAXN];
int c[10000005];
u32 getpos(u32 x) { return x <= sqr ? id[0][x] : id[1][N / x]; }
struct HASH {
    int sumE,head[974711 + 5];u32 mo = 974711;
    struct node {
    pair<u32,u32> x;u32 val;int next;
    }E[10000005];
    u32 hash_key(pair<u32,u32> t) {
    return t.fi << 16 | t.se;
    }
    void add(pair<u32,u32> x,u32 v) {
    u32 u = hash_key(x) % mo;
    E[++sumE].x = x;E[sumE].val = v;
    E[sumE].next = head[u];head[u] = sumE;
    }
    bool count(pair<u32,u32> x) {
    u32 u = hash_key(x) % mo;
    for(int i = head[u] ; i ; i = E[i].next) {
        if(E[i].x == x) return true;
    }
    return false;
    }
    u32 query(pair<u32,u32> x) {
    u32 u = hash_key(x) % mo;
    for(int i = head[u] ; i ; i = E[i].next) {
        if(E[i].x == x) return E[i].val;
    }
    return 0;
    }
}hsh;
u32 rf[MAXN];
bool vis[MAXN];
u32 fpow(u32 x, u32 c) {
    u32 res = 1, t = x;
    while (c) {
        if (c & 1)
            res = res * t;
        t = t * t;
        c >>= 1;
    }
    return res;
}

void Process() {
    prime[0] = 1;
    pw[0] = 1;
    for (u32 i = 2; i <= sqr; ++i) {
        if (!nonprime[i]) {
            prime[++tot] = i;
            pw[tot] = fpow(i, K);
        }
        for (int j = 1; j <= tot; ++j) {
            if (prime[j] > sqr / i)
                break;
            nonprime[i * prime[j]] = 1;
            if (i % prime[j] == 0)
                break;
        }
    }
    for (u32 i = 1; i <= N; ++i) {
        u32 r = N / (N / i);
        pos[++cnt] = N / i;
        h[cnt] = N / i - 1;
        if (N / i <= sqr)
            id[0][N / i] = cnt;
        else
            id[1][N / (N / i)] = cnt;
        i = r;
    }
    for (int j = 1; j <= tot; ++j) {
        for (int i = 1; i <= cnt && 1LL * prime[j] * prime[j] <= pos[i]; ++i) {
            int k = getpos(pos[i] / prime[j]);
            h[i] = h[i] - (h[k] - (j - 1));
        }
    }
}
u32 SF(u32 x, u32 y) {
    if (x <= 1 || prime[y] > x)
        return 0;
    if (hsh.count(mp(x, y)))
        return hsh.query(mp(x, y));
    ++rec;
    u32 res = pw[y - 1] * (h[getpos(x)] - (y - 1));
    for (u32 j = y; j <= tot && 1LL * prime[j] * prime[j] <= x; ++j) {
        int64 t1 = prime[j], t2 = 1LL * prime[j] * prime[j];
        for (int e = 1; t2 <= x; t1 = t2, t2 *= prime[j], ++e) {
            res += SF(x / t1, j + 1) + pw[j];
        }
    }
    // out(x);space;out(y);enter;
    // printf("%.3lf\n",(db)clock() / CLOCKS_PER_SEC);
    hsh.add(mp(x,y),res);
    return res;
}
u32 FF(u32 x) {
    if (vis[getpos(x)])
        return rf[getpos(x)];
    u32 res = SF(x, 1);
    for (u32 i = 2; i <= x; ++i) {
        u32 r = x / (x / i);
        res -= (r - i + 1) * FF(x / i);
        i = r;
    }
    vis[getpos(x)] = 1;
    rf[getpos(x)] = res;
    return res;
}
int main() {
#ifdef ivorysi
    freopen("f1.in", "r", stdin);
#endif
    read(N);
    read(K);
    sqr = sqrt(N);
    Process();
    u32 ans = 0;
    u32 pre = 0;
    for (u32 T = 1; T <= N; ++T) {
        u32 r = N / (N / T);
        u32 nw = FF(r);
        ans += (N / T) * (N / T) * (nw - pre);
        pre = nw;
        T = r;
    }
    out(ans);
    enter;
    // out(rec);enter;
    // printf("%.3lf\n",(db)clock() / CLOCKS_PER_SEC);
    return 0;
}

【51nod】1847 奇怪的數學題

這個怎麼跟上面那個題好像啊！

因此直接給出杜教篩的部分
\[ S(n) = \sum_{i = 1}^{n} sgcd(n)^{k} - \sum_{i = 1}^{n}S(\lfloor \frac{n}{i} \rfloor) \]
em……徹底同樣好吧

發現\(sgcd(n) = \frac{n}{min_{p}(n)}\)，就是n除以最小質數

這個和g函數就已經很相似了……

考慮到\(g(n,j)\)求的是\([1,n]\)中的質數和最小質因子大於\(p_{j}\)的\(x^{k}\)和

而後只要枚舉一個質數\(p_{i} * g(\lfloor \frac{n}{p_{i}}\rfloor,i - 1)\)就搞定了，最後再統計一遍純質數的答案

不須要後面的那部分了

至於天然數冪和的那部分能夠用乘方轉斯特林數推一個公式出來

（還傻乎乎求了一遍後面的我真是zz）

#pragma GCC optimize("Ofast,unroll-loops,no-stack-protector,fast-math")
#include <bits/stdc++.h>
#define fi first
#define se second
#define pii pair<int,int>
#define mp make_pair
#define pb push_back
#define space putchar(' ')
#define enter putchar('\n')
#define eps 1e-10
#define ba 47
#define MAXN 100005
//#define ivorysi
using namespace std;
typedef long long int64;
typedef unsigned int u32;
typedef double db;
template<class T>
void read(T &res) {
    res = 0;T f = 1;char c = getchar();
    while(c < '0' || c > '9') {
    if(c == '-') f = -1;
    c = getchar();
    }
    while(c >= '0' && c <= '9') {
    res = res * 10 +c - '0';
    c = getchar();
    }
    res *= f;
}
template<class T>
void out(T x) {
    if(x < 0) {x = -x;putchar('-');}
    if(x >= 10) {
    out(x / 10);
    }
    putchar('0' + x % 10);
}
u32 N,sqr,K;
u32 h[MAXN],g[MAXN],f[MAXN],pos[MAXN],id[2][MAXN],sum[MAXN],cnt;
u32 prime[MAXN],tot,S[55][55],rec[55];
bool nonprime[MAXN];
u32 getpos(u32 x) {
    return x <= sqr ? id[0][x] : id[1][N / x];
}
u32 getsum(u32 x) {
    u32 res = 0;
    for(u32 j = 0 ; j <= K ; ++j) {
    if(x < j) break;
    u32 t = S[K][j];
    bool f = 0;
    for(u32 h = 0 ; h < j + 1 ; ++h) {
        u32 p = (x + 1 - h);
        if(!f && p % (j + 1) == 0) {f = 1;p /= j + 1;}
        t = t * p;
    }
    res += t;
    }
    return res;
}
u32 fpow(u32 x,u32 c) {
    u32 res = 1,t = x;
    while(c) {
    if(c & 1) res = res * t;
    t = t * t;
    c >>= 1;
    }
    return res;
}
void Process() {
    for(u32 i = 2 ; i <= sqr ; ++i) {
    if(!nonprime[i]) {
        prime[++tot] = i;
        sum[tot] = sum[tot - 1] + fpow(i,K);
    }
    for(int j = 1 ; j <= tot ; ++j) {
        if(prime[j] > sqr / i) break;
        nonprime[prime[j] * i] = 1;
        if(i % prime[j] == 0) break;
    }
    }
    S[0][0] = 1;
    for(int i = 1 ; i <= K ; ++i) {
    for(int j = 1 ; j <= i ; ++j) {
        S[i][j] = S[i - 1][j - 1] + j * S[i - 1][j];
    }
    }
    for(u32 i = 1 ; i <= N ; ++i) {
    u32 r = N / (N / i);
    pos[++cnt] = N / i;
    if(N / i <= sqr) id[0][N / i] = cnt;
    else id[1][N / (N / i)] = cnt;
    h[cnt] = N / i - 1;
    g[cnt] = getsum(N / i) - 1;
    i = r;
    }
    for(u32 j = 1 ; j <= tot ; ++j) {
    for(u32 i = 1 ; i <= cnt && 1LL * prime[j] * prime[j] <= pos[i] ; ++i) {
        u32 k = getpos(pos[i] / prime[j]);
        g[i] -= (sum[j] - sum[j - 1]) * (g[k] - sum[j - 1]);
        h[i] -= h[k] - (j - 1);
        f[i] += g[k] - sum[j - 1];
    }
    }
    for(int i = 1 ; i <= cnt ; ++i) f[i] += h[i];
}
bool vis[MAXN];
u32 rs[MAXN];
u32 UT(u32 n) {
    if(vis[getpos(n)]) return rs[getpos(n)];
    u32 res = f[getpos(n)];
    for(u32 i = 2 ; i <= n ; ++i) {
    u32 r = n / (n / i);
    res -= (r - i + 1) * UT(n / i);
    i = r;
    }
    vis[getpos(n)] = 1;rs[getpos(n)] = res;
    return res;
}
void Solve() {
    read(N);read(K);sqr = sqrt(N);
    Process();
    u32 ans = 0,pre = 0,nw;
    for(u32 i = 1 ; i <= N ; ++i) {
    u32 r = N / (N / i);
    nw = UT(r);
    ans += (N / i) * (N / i) * (nw - pre);
    pre = nw;
    i = r;
    }
    out(ans);enter;
}
int main(){
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    Solve();
    return 0;
}
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程