杜教ふるい学習報告(訓練に伴って更新)


目次
杜教ふるいについての簡単な説明
トレーニング記録:
51 nod 1244モビウス関数の和
51 nod 1239オーラ関数の和
bzoj3944 sum
hdu5608 Function
 
 
杜教ふるいについての簡単な説明
skywalkertのブログを見て大体分かりました.
author: skywalkert  original article: http://blog.csdn.net/skywalkert/article/details/50500009  last update time : 2017-04-01
一般的な考え方は、ディリクレイボリュームによって接頭辞和をよりよく計算する関数を構築することができ、ボリュームに使用される別の関数も計算しやすい場合、計算プロセスを簡略化することができる.
関数を構築した後の前n^ \frac{2}{3}項の接頭辞と前処理を行い,残りの部分は前に構築した関数によって再帰される.f(n) = g(i) - \sum_{i = 2}^{n}f(\frac{n}{i})
この形式では,f(n)は求められる接頭辞和であり,g(n)は構造関数である.そして再帰する.
特に注意!!!サブリニアふるいの複雑さがO(n)より低いため、一般的なデータ範囲は1 e 8-1 e 11で、型取り!!!必ず型取りに注意!!!
前置技能:モビウス反転、常用積性関数、線形ふるい
トレーニング記録:
51 nod 1244モビウス関数の和
題意:モビウス関数の接頭辞和を求める.
構想:データ範囲は1 e 10であるため、明らかに線形のやり方を採用することができなくて、杜教ふるいを採用して、1つの関数を構築してf(n) = \sum_{i = 1}^{n}g(i) - \sum_{i = 2}^{n}f(\frac{n}{i})を満たして、モビウス関数の性質[n==1]=\sum_{d|n}^{n}\mu(d)に気づいて、1=\sum_{i = 1}^{n}[n==1]=\sum_{i=1}^{n}\sum_{d|n}^{n}\mu(d)=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor \frac{n}{i} \rfloor}\mu(d)があって、そこで、私達はg(n)=1を構築します.
その後式はf(n)=g(n)-\sum_{i=2}^{n}f(\frac{n}{i})=1-\sum_{i=2}^{n}f(\frac{n}{i})となり,4 e 6範囲のモビウス関数のプレフィックスと再帰を前処理すると終了し,計算した項をmapで記録して複雑さを低減する.
複雑度:O(n^\frac{2}{3})
#include 
using namespace std;
typedef long long ll;
const int N = 4641600;
const int inf = 0x3f3f3f3f;
const int mod = 1e6 + 3;
//const double pi = acos(-1.0);
const double eps = 1e-6;
const int inv = 5e8 + 4;
int gcd(int a,int b){return b ? gcd(b,a % b) : a;}
ll n,m;
int mu[N],prime[N],cnt;
bool vis[N];
inline ll read()
{
    ll x=0,f=1;char ch=getchar();
    while(ch'9'){if(ch=='-') f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0'; ch=getchar();}
    return x*f;
}
void Mobius()
{
    cnt = 0;
    mu[1] = 1;mu[0] = 0;
    for(int i = 2; i < N; ++i){
        if(!vis[i]){
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1; j <= cnt; ++j){
            if(i * prime[j] >= N) break;
            vis[i * prime[j]] = true;
            if(i % prime[j]== 0){
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
        mu[i] += mu[i - 1];
    }
}
mapmp;
int work(ll jq)
{
    if(jq < N) return mu[jq];
    if(mp[jq]) return mp[jq];
    int res = 1;ll last;
    for(ll i = 2;i <= jq;i = last + 1){
        last = jq / (jq / i);
        res -= (last - i + 1) * work(jq / i);
    }
    mp[jq] = res;
    return res;
}
int main()
{
    Mobius();//cout << pow((ll)10000000000,2.0/ 3);
    n = read(),m = read();
    printf("%d
",work(m) - work(n - 1)); return 0; }

51 nod 1239オーラ関数の和
題意:オラ関数の接頭辞と
考え方:オーラ関数の性質n=\sum_{n|d}^{n}\phi(d)=\sum_{d = 1}^{\lfloor \frac{n}{d} \rfloor}\phi(d)に注目するため、\sum_{i=1}^{n}i=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor \frac{n}{\frac{i}{d}} \rfloor }\phi(d)=\sum_{i=1}^{n}\phi(\lfloor \frac{n}{i} \rfloor )があるので、g(n)=\sum_{i=1}^{n}iを構成するので、f(n)=g(n)-\sum_{i=2}^{n}\phi(\lfloor \frac{n}{i} \rfloor )を解くことになる.
複雑度:O(n^\frac{2}{3})
#include 
using namespace std;
typedef long long ll;
const int N = 4641600;
const int inf = 0x3f3f3f3f;
const int mod = 1e9 + 7;
//const double pi = acos(-1.0);
const double eps = 1e-6;
const int inv = 5e8 + 4;
int gcd(int a,int b){return b ? gcd(b,a % b) : a;}
ll phi[N];
int prime[N],cnt;
bool isprime[N];
ll n;
inline ll add(ll x)
{
    return x >= mod ? x - mod : x;
}
inline ll sub(ll x)
{
    return x < 0 ? x + mod : x;
}
void get_phi()
{
   phi[1] = 1;
   for(int i = 2;i < N;++i){
       if(!isprime[i]){
             prime[++cnt] = i;
             phi[i] = i - 1;
        }
       for(int j = 1;j <= cnt && i * prime[j] < N;++j)
       {
          isprime[i * prime[j]] = 1;
          if(i % prime[j] == 0){
             phi[i * prime[j]] = phi[i] * prime[j];
             break;
          }
          else phi[i * prime[j]] = phi[i] * (prime[j] - 1);
       }
   }
   for(int i = 2;i < N;++i) phi[i] = add(phi[i] + phi[i - 1]);
}

mapmp;
ll work(ll jq)
{
    if(jq < N) return phi[jq];
    if(mp[jq]) return mp[jq];
    ll res = jq % mod * ((jq + 1) % mod) % mod * inv % mod,last;
    for(ll i = 2;i <= jq;i = last + 1){
        last = jq / (jq / i);
        res = sub(res - (last - i + 1) * work(jq / i) % mod);
    }
    mp[jq] = res;
    return res;
}
int main()
{
    get_phi();
    scanf("%lld",&n);
    printf("%lld
",work(n)); return 0; }

bzoj3944 sum
題意:オラ接頭辞とモビウス関数接頭辞と
アイデア:同じように、スタックのオーバーフローを避けるために構造体を開くことに注意します!!!
#pragma comment(linker, “/STACK:1024000000,1024000000”
#include 
using namespace std;
typedef long long ll;
const int N = 5e6 + 7;
const int inf = 0x3f3f3f3f;
const int mod = 1e9 + 7;
//const double pi = acos(-1.0);
const double eps = 1e-6;
const int inv = 5e8 + 4;
int gcd(int a,int b){return b ? gcd(b,a % b) : a;}
ll phi[N];
int prime[N],cnt,mu[N];
bool isprime[N];
ll n;
void get_phi()
{
   phi[1] = 1;mu[1] = 1;mu[0] = 0;
   for(int i = 2;i < N - 2;++i){
       if(!isprime[i]){
             prime[++cnt] = i;
             phi[i] = i - 1;
             mu[i] = -1;
        }
       for(int j = 1;j <= cnt && i * prime[j] < N - 2;++j){
          isprime[i * prime[j]] = 1;
          if(i % prime[j] == 0){
             phi[i * prime[j]] = phi[i] * prime[j];
             mu[i * prime[j]] = 0;
             break;
          }
          phi[i * prime[j]] = phi[i] * (prime[j] - 1);
          mu[i * prime[j]] = -mu[i];
       }
       mu[i] += mu[i - 1];
   }
   for(int i = 2;i < N - 2;++i) phi[i] = phi[i] + phi[i - 1];
}
struct node{
    ll p;
    int m;
};
mapmp,mp1;
node work(ll jq)
{
    if(jq < N - 3) return (node){phi[jq],mu[jq]};
    if(mp[jq]) return (node){mp[jq],mp1[jq]};
    node res = (node){(jq & 1 ? jq : jq / 2) * (jq & 1 ? (jq + 1) / 2 : jq + 1),1},tmp;
    ll last;
    for(ll i = 2;i <= jq;i = last + 1){
        last = jq / (jq / i),tmp = work(jq / i);
        res.p -= (last - i + 1) * tmp.p;
        res.m -= (last - i + 1) * tmp.m;
    }
    if(!mp[jq]) mp[jq] = res.p;
    if(!mp1[jq]) mp1[jq] = res.m;
    return res;
}
int main()
{
    get_phi();
    int t;
    scanf("%d",&t);
    while(t--){
        scanf("%lld",&n);
        if(!n){
            printf("0 0
"); continue; } node r = work(n); printf("%lld %d
",r.p,r.m); } return 0; }

hdu5608 Function
标题:既知f(x)がn^2-3n+2=\sum_{d|n}^{n}f(d),F(x)=\sum_{i=1}^{n}f(n)を満たすF(x)を求める.
構想:モビウスの再演を思い浮かべる.令f(n) = \sum_{d|n}^{n}g(d)*\mu(\frac{n}{d})、だからF(n)=\sum_{i=1}^{n}f(n)=\sum_{i=1}^{n}\sum_{d=1}^{\lfloor \frac{n}{d} \rfloor }g(i)*\mu(d)、前のモビウスが押した経験があって、私たちは直接G(n)=\sum_{i=1}^{n}g(x)*1を構築して、突然はっきりしているのではないでしょうか.この問題は解決しました.
//#pragma comment(linker, “/STACK:1024000000,1024000000”
#include 
using namespace std;
typedef long long ll;
const int N = 1e6 + 7;
const int inf = 0x3f3f3f3f;
const int mod = 1e9 + 7;
const int inv2 = 5e8 + 4;
//const double pi = acos(-1.0);
const double eps = 1e-6;
const int inv = 333333336;
int gcd(int a,int b){return b ? gcd(b,a % b) : a;}
int prime[N],cnt,mu[N];
ll g[N];
bool isprime[N];
ll n;
inline ll gao(ll x)
{
    return (x - 1) * (x - 2) % mod * x % mod * inv % mod;
}
void get_phi()
{
   mu[1] = 1;
   for(int i = 2;i < N;++i){
       if(!isprime[i]){
             prime[++cnt] = i;
             mu[i] = -1;
        }
       for(int j = 1;j <= cnt && i * prime[j] < N;++j){
          isprime[i * prime[j]] = 1;
          if(i % prime[j] == 0){
             mu[i * prime[j]] = 0;
             break;
          }
          mu[i * prime[j]] = -mu[i];
       }
       //mu[i] += mu[i - 1];
   }
   for(int i = 1;i < N;++i)
       for(int j = i;j < N;j += i)
           g[j] = ((g[j] + (ll)mu[j / i] * (i - 1) % mod * (i - 2) % mod) % mod + mod) % mod;
    for(int i = 2;i < N;++i) g[i] = (g[i] + g[i - 1]) % mod;
}
mapmp;
ll solve(ll pos)
{
    if(pos < N) return g[pos];
    if(mp[pos]) return mp[pos];
    ll res = gao(pos),last;
    for(ll i = 2;i <= pos;i = last + 1){
        last = pos / (pos / i);
        res = ((res - (last - i + 1) * solve(pos / i) % mod) % mod + mod) % mod;
    }
    mp[pos] = res;
    return res;
}
int main()
{
    get_phi();
    int t;
    scanf("%d",&t);
    while(t--){
        scanf("%lld",&n);
        printf("%lld
",solve(n)); } return 0; }