[ICPC 2018 Shenyang Onsite E][Gym101955E]The Kouga Ninja Scrolls

Link

Solution

曼哈顿转切比雪夫之后x,y独立了。
问题变成从区间中选取两个颜色不同的数字使得差值最大。
维护区间最大次大、最小次小值,满足最大次大、最小次小值分别是不同的颜色,最大差值一定来源于这四个数字。枚举合法情况判一判就好了。

Code

const int XN=1e5+11;
const long long INF=1e18;

struct SegTree {

    struct Data {
        std::pair<long long,int> mn[2],mx[2];

        Data() {
            mn[0]=mn[1]=std::pair<long long,int>(INF,-1);
            mx[0]=mx[1]=std::pair<long long,int>(-INF,-1);
        }

        long long Solve() {
            long long res=0;
            for(int i=0;i<2;++i)
                for(int j=0;j<2;++j)
                    if(mn[i].second!=mx[j].second)
                        Enlarge(res,mx[j].first-mn[i].first);
            return res;
        }

        friend Data Merge(Data const &a,Data const &b) {
            Data res;
            auto Update=[&](Data const &x) {
                for(int i=0;i<2;++i) {
                    if(x.mn[i].first<res.mn[0].first) {
                        if(x.mn[i].second!=res.mn[0].second)
                            res.mn[1]=res.mn[0];
                        res.mn[0]=x.mn[i];
                    } else if(x.mn[i].first<res.mn[1].first && x.mn[i].second!=res.mn[0].second)
                        res.mn[1]=x.mn[i];

                    if(x.mx[i].first>res.mx[0].first) {
                        if(x.mx[i].second!=res.mx[0].second)
                            res.mx[1]=res.mx[0];
                        res.mx[0]=x.mx[i];
                    } else if(x.mx[i].first>res.mx[1].first && x.mx[i].second!=res.mx[0].second)
                        res.mx[1]=x.mx[i];

                }
            };
            Update(a);Update(b);
            return res;
        }
    };

    struct Node {
        Node *son[2];
        Data v;
        Node(std::pair<long long,int> const &x) {
            v.mn[0]=v.mx[0]=x;
        }

        void *operator new(size_t flag) {
            static Node *pool=(Node*)malloc(XN*10*sizeof(Node)),*mem=0;
            return flag?mem++:mem=pool;
        }

        void Up() {
            v=Merge(son[0]->v,son[1]->v);
        }
    }*root;

    int L,R;

    SegTree(std::pair<long long,int> a[],int n):L(1),R(n) {
        std::function<void(Node*&,int,int)> Build=
            [&](Node *&pos,int L,int R) {
                pos=new Node(a[L]);
                if(L==R)
                    return;
                else {
                    int M=(L+R)/2;
                    Build(pos->son[0],L,M);
                    Build(pos->son[1],M+1,R);
                    pos->Up();
                }

            };
        Build(root,1,n);
    }

    long long Query(int l,int r) {
        std::function<Data(Node*,int,int)> Query=
            [&](Node *pos,int L,int R) {
                if(l<=L && R<=r)
                    return pos->v;
                else {
                    int M=(L+R)/2;
                    Data res;
                    if(l<=M)
                        res=Merge(res,Query(pos->son[0],L,M));
                    if(M+1<=r)
                        res=Merge(res,Query(pos->son[1],M+1,R));
                    return res;
                }
            };
        return Query(root,L,R).Solve();
    }

    template <class T>
    void Modify(int goal,T F) {
        std::function<void(Node*,int,int)> Modify=
            [&](Node *pos,int L,int R) {
                if(L==R) {
                    F(pos->v.mn[0]);
                    F(pos->v.mx[0]);
                } else {
                    int M=(L+R)/2;
                    goal<=M?Modify(pos->son[0],L,M):Modify(pos->son[1],M+1,R);
                    pos->Up();
                }
            };
        Modify(root,L,R);
    }
};

void Work() {
    SegTree::Node::operator new(0);

    auto Convert=[](int x,int y)->std::pair<int,int> { return std::pair<int,int>(x+y,x-y); };

    int n,m;fin>>n>>m;
    static std::pair<long long,int> cx[XN],cy[XN];
    for(int i=1;i<=n;++i) {
        int x,y,c;fin>>x>>y>>c;
        std::tie(cx[i].first,cy[i].first)=Convert(x,y);
        cx[i].second=cy[i].second=c;
    }

    SegTree segX(cx,n),segY(cy,n);

    for(int q=1;q<=m;++q) {
        int op;fin>>op;
        switch(op) {
            case 1 : {
                int id,x,y;fin>>id>>x>>y;
                std::tie(x,y)=Convert(x,y);
                segX.Modify(id,[&](std::pair<long long,int> &a) { a.first+=x; });
                segY.Modify(id,[&](std::pair<long long,int> &a) { a.first+=y; });
            } break;
            case 2 : {
                int id,c;fin>>id>>c;
                auto Change=[&](std::pair<long long,int> &a) { a.second=c; };
                segX.Modify(id,Change);
                segY.Modify(id,Change);
            } break;
            case 3 : {
                int l,r;fin>>l>>r;
                fout<<std::max(segX.Query(l,r),segY.Query(l,r))<<'\n';
            } break;
        }
    }
}

int main() {
    int T;fin>>T;

    for(int ks=1;ks<=T;++ks) {
        printf("Case #%d:\n",ks);
        Work();
    }
    return 0;
}

[ICPC 2018 Shenyang Onsite D][Gym 101955D]Diameter of a tree

Link

感觉很有意思,非常想A掉这道题,但是找了好久都没有任何的题解或者AC代码。无奈看到过掉这道题的一位KUT小哥哥CF在线,于是讨教了一发,朝鲜小哥哥马上就回给我了AC代码!心里有点小激动
KUT小哥哥的function对象用的实在6…感觉这样写特别紧凑易读!

Solution

任何一个边数为c边权和为w的路径,在x点的值就是cx+w
问题可以转化为一堆直线,求x点处最高的直线的纵坐标。
半平面交可以转化为凸包,因此维护一个全局凸包就行了。
点分治,对于分出来的子树的凸包,需要求出从不同子树选出两点加和产生新点集的凸包。
为了保证复杂度,KUT小哥哥用了类似线段树一样的分治思想。
这是这道题的精华部分,然而我现在还是不懂为什么KUT小哥哥AddHull这么做是对的。先A了过几天再看看吧。

Code

const int XN=2e5+11;

struct Point {
    long long x,y;

    long long operator ()(int a) {
        return x*a+y;
    }

    friend Point operator +(Point const &a,Point const &b) {
        return Point{a.x+b.x,a.y+b.y};
    }   

    friend Point operator -(Point const &a,Point const &b) {
        return Point{a.x-b.x,a.y-b.y};
    }   

    friend long long Outer(Point const &a,Point const &b) {
        return a.x*b.y-a.y*b.x;
    }

    friend bool operator <(Point const &a,Point const &b) {
        return a.x<b.x || (a.x==b.x && a.y<b.y);
    }
};

std::vector<Point> ConvexHull(std::vector<Point> a) {
    if(!std::is_sorted(a.begin(),a.end()))
        std::sort(a.begin(),a.end());
    std::vector<Point> hull;
    for(Point p : a) {
        while(hull.size()>1 && Outer(hull.back()-hull[hull.size()-2],p-hull.back())>=0)
            hull.pop_back();
        hull.push_back(p);
    }
    return hull;
}

std::vector<Point> MergeHull(std::vector<Point> const &a,std::vector<Point> const &b) {
    std::vector<Point> c(a.size()+b.size());
    std::merge(a.begin(),a.end(),b.begin(),b.end(),c.begin());
    return ConvexHull(c);
}

std::vector<Point> AddHull(std::vector<Point> const &v1,std::vector<Point> const &v2) {
    std::vector<Point> hull;
    for(size_t p1=0,p2=0;p1<v1.size() && p2<v2.size();) {
        hull.push_back(v1[p1]+v2[p2]);
        if(p1+1<v1.size() && p2+1<v2.size()) {
            if(Outer(v1[p1+1]-v1[p1],v2[p2+1]-v2[p2])<0)
                ++p1;
            else
                ++p2;
        } else if(p1+1<v1.size())
            ++p1;
        else
            ++p2;
    }
    return hull;
}

struct Edge {
    int to,v;
    Edge *pre;
    void *operator new(size_t flag) {
        static Edge *pool=(Edge*)malloc(XN*2*sizeof(Edge)),*mem;
        return flag?mem++:mem=pool;
    }
}*G[XN];

void AddEdge(int x,int y,int v) {
    G[x]=new Edge{y,v,G[x]};
    G[y]=new Edge{x,v,G[y]};
}

long long mxw[XN];

bool ud[XN];

void DC(int pos=1) {
    static int size[XN];
    std::function<int(int)> Centroid=[&](int pos) {
        std::function<void(int,int)> GetSize=[&](int pos,int fa) {
            size[pos]=1;
            for(Edge *e=G[pos];e;e=e->pre)
                if(fa!=e->to && !ud[e->to]) {
                    GetSize(e->to,pos);
                    size[pos]+=size[e->to];
                }
        };
        GetSize(pos,0);
        int tot=size[pos];
        bool found;
        do {
            found=0;
            for(Edge *e=G[pos];e;e=e->pre)
                if(!ud[e->to] && size[e->to]<size[pos] && size[e->to]*2>=tot) {
                    pos=e->to;
                    found=1;
                    break;
                }
        } while(found);
        return pos;
    };

    static std::vector<Point> v[XN];

    ud[pos=Centroid(pos)]=1;
    std::vector<std::vector<Point>> all;
    for(Edge *e=G[pos];e;e=e->pre)
        if(!ud[e->to]) {
            v[e->to].clear();
            std::function<void(int,int,std::vector<Point>&,int,long long)> DFS=
                [&](int pos,int fa,std::vector<Point> &v,int dep,long long length) {
                Enlarge(mxw[dep],length);
                v.push_back({dep,length});
                for(Edge *e=G[pos];e;e=e->pre)
                    if(e->to!=fa && !ud[e->to])
                        DFS(e->to,pos,v,dep+1,length+e->v);
            };
            DFS(e->to,pos,v[e->to],1,e->v);
            all.push_back(ConvexHull(v[e->to]));
        }

    std::function<void(int,int,int)> Solve=[&](int pos,int L,int R) {
        static std::vector<Point> T[XN*4];
        if(L==R)
            T[pos]=all[L];
        else {
            int M=(L+R)/2;
            Solve(pos<<1,L,M);Solve(pos<<1|1,M+1,R);
            T[pos]=MergeHull(T[pos<<1],T[pos<<1|1]);
            std::vector<Point> sum=AddHull(T[pos<<1],T[pos<<1|1]);
            for(auto p : sum)
                Enlarge(mxw[p.x],p.y);
        }
    };

    if(all.size())
        Solve(1,0,all.size()-1);

    for(Edge *e=G[pos];e;e=e->pre)
        if(!ud[e->to])
            DC(Centroid(e->to));
}

void Work() {
    int n,m;
    fin>>n>>m;
    Edge::operator new(0);
    for(int i=1;i<=n;++i) {
        G[i]=0;
        mxw[i]=-1;
        ud[i]=0;
    }
    for(int i=1;i<=n-1;++i) {
        int x,y,v;fin>>x>>y>>v;
        AddEdge(x,y,v);
    }
    DC();
    std::vector<Point> poly;
    for(int i=1;i<=n-1;++i)
        if(mxw[i]>=0)
            poly.push_back({i,mxw[i]});
    poly=ConvexHull(poly);

    for(int q=1;q<=m;++q) {
        int x;fin>>x;
        int L=0,R=poly.size()-1;
        while(R-L>2) {
            int M1=(L*2+R)/3,M2=(L+R*2)/3;
            if(poly[M1](x)>poly[M2](x))
                R=M2;
            else
                L=M1;
        }
        long long Ans=0;
        for(int i=L;i<=R;++i)
            Enlarge(Ans,poly[i](x));
        fout<<Ans<<'\n';
    }
}

int main() {
    int T;fin>>T;
    for(int ks=1;ks<=T;++ks) {
        printf("Case #%d:\n",ks);
        Work();
    }
    return 0;
}

[BZOJ 4176]Lucas的数论

Link

Solution

比较套路地画出来

\sum_{d}\mu(d)F([\frac nd])^2

其中

F(n)=\sum_{i=1}^n[\frac ni]

主要是求F函数部分的复杂度如何估计

\begin{aligned}
T(n)\approx&\sum_{d=1}^{\sqrt n}(\sqrt d+\sqrt{\frac nd})\\
\approx&\int_{1}^{\sqrt n}(\sqrt x+\sqrt {\frac nx}){\rm d}x\\
=&(2(nx)^{\frac 12}+\frac {2x^{\frac 32}}3)\bigg|_1^{\sqrt n}\\
=&O(n^{\frac 34})
\end{aligned}

Code

const int P=1000000007,N=1e6,XN=N+11;

int prime[XN],mu[XN],pcnt,smu[XN];
void Prep() {
    static bool notPrime[XN];
    mu[1]=1;
    for(int i=2;i<=N;++i) {
        if(!notPrime[i]) {
            prime[++pcnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0) {
                mu[i*prime[j]]=0;
                break;
            } else
                mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;++i)
        smu[i]=smu[i-1]+mu[i];
}

int Smu(int n) {
    if(n<=N)
        return smu[n];
    else {
        static std::unordered_map<int,int> map;
        if(map.count(n))
            return map[n];
        long long res=1;
        for(int l=2,r;l<=n;l=r+1) {
            r=n/(n/l);
            (res-=(long long)(r-l+1)*Smu(n/l))%=P;
        }
        return map[n]=res;
    }
}

int F(int n) {
    long long res=0;
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        (res+=(long long)(r-l+1)*(n/l))%=P;
    }
    return res;
}

int main() {
    Prep();
    int n;fin>>n;
    long long Ans=0;
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        long long x=F(n/l);
        (x*=x)%=P;
        (Ans+=(long long)(Smu(r)-Smu(l-1))*x)%=P;
    }
    Ans<0 && (Ans+=P);
    fout<<Ans<<'\n';
    return 0;
}

BZOJ 3512 Dzy loves math IV

Link

Solution

n很小,不妨枚举。
对于一个固定的n
\sum_{i=1}^m\varphi(ni)。不妨设n为无平方因子数,那么

\begin{aligned}
S(n,m)=&\sum_{i=1}^m\varphi(ni)\\
=&\sum_{i=1}^m\varphi(i)\varphi(\frac n{\gcd(i,n)})\gcd(i,n)\\
=&\sum_{i=1}^m\varphi(i)\varphi(\frac n{\gcd(i,n)})\sum_{d|\gcd(n,i)}\varphi(d)\\
=&\sum_{i=1}^m\varphi(i)\sum_{d|\gcd(n,i)}\varphi(\frac {\gcd(i,n)}d)\varphi(\frac n{\gcd(i,n)})\\
\end{aligned}

而由于n为无平方因子数,\gcd(\frac {n}{\gcd(i,n)},\gcd(i,n))=1
因此上式可以化为

\begin{aligned}
&\sum_{i=1}^m\varphi(i)\sum_{d|\gcd(n,i)}\varphi(\frac {\gcd(i,n)}d)\varphi(\frac n{\gcd(i,n)})\\
=&\sum_{i=1}^m\varphi(i)\sum_{d|\gcd(n,i)}\varphi(\frac nd)\\
=&\sum_{d|n}\varphi(\frac nd)\sum_{i=1}^{[\frac md]}\varphi(id)\\
=&\sum_{d|n}\varphi(\frac nd)\sum_{i=1}^{[\frac md]}\varphi(id)\\
=&\sum_{d|n}\varphi(\frac nd)S(d,[\frac md])
\end{aligned}

边界情况S(1,m)就用杜教筛求\varphi的前缀和。
而S的第一维状态数为n,第二维状态数为O(\sqrt m),因为第二维会出现的状态都是m经过多次除法然后下取整的结果,只有O(\sqrt m)种。而每次计算枚举n的约数复杂度为O(\sqrt n),因此总复杂度就是O(n(\sqrt m+\sqrt n)+T),其中T为调用杜教筛求\varphi的耗时。

Code

const int P=1e9+7;

int Add(int x,int const &y) {
    return (x+=y)>=P?x-P:x;
}

int Minus(int x,int const &y) {
    return (x-=y)<0?x+P:x;
}

int Mul(long long x,int const &y) {
    return x*y%P;
}

const int N=1e6,XN=N+11;

int prime[XN],pcnt,phi[XN],dnum[XN],sphi[XN];
void Prep() {
    static bool notPrime[XN];
    phi[1]=dnum[1]=1;
    for(int i=2;i<=N;++i) {
        if(!notPrime[i]) {
            prime[++pcnt]=i;
            dnum[i]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0) {
                dnum[i*prime[j]]=dnum[i];
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            } else {
                dnum[i*prime[j]]=dnum[i]*prime[j];
                phi[i*prime[j]]=phi[i]*phi[prime[j]];
            }
        }
    }
    for(int i=1;i<=N;++i)
        sphi[i]=Add(sphi[i-1],phi[i]);
}

int Sphi(int n) {
    if(n<=N)
        return sphi[n];
    else {
        static std::unordered_map<int,int> map;
        if(map.count(n))
            return map[n];
        int &res=map[n]=(long long)n*(n+1)/2%P;
        for(int l=2,r;l<=n;l=r+1) {
            r=n/(n/l);
            res=Minus(res,Mul(r-l+1,Sphi(n/l)));
        }
        return res;
    }
}

int Phi(int n) {
    if(n<=N)
        return phi[n];
    else {
        static std::unordered_map<int,int> map;
        if(map.count(n))
            return map[n];
        int &res=map[n]=n;
        for(int d=2,m=sqrt(n)+0.5;d<=m;++d)
            if(n%d==0) {
                do n/=d;
                while(n%d==0);
                (res/=d)*=(d-1);
            }
        if(n!=1)
            (res/=n)*=(n-1);
        return res;
    }
}

int S(int n,int m) {
    if(n==1)
        return Sphi(m);
    else {
        static std::map<std::pair<int,int>,int> map;
        std::pair<int,int> pid(n,m);
        if(map.count(pid))
            return map[pid];
        int &res=map[pid]=0;
        for(int d=1;d*d<=n;++d)
            if(n%d==0) {
                res=Add(res,Mul(Phi(n/d),S(d,m/d)));
                if(d*d!=n)
                    res=Add(res,Mul(Phi(d),S(n/d,m/(n/d))));
            }
        return res;
    }
}

int main() {
    Prep();
    int n,m,Ans=0;fin>>n>>m;
    for(int i=1;i<=n;++i)
        Ans=Add(Ans,Mul(i/dnum[i],S(dnum[i],m)));
    fout<<Ans<<'\n';
    return 0;
}

总结一些积性函数的性质

  1. \sum_{d|n}\mu(d)=[n=1]
    x=\prod p_i^{k_i}f(x)=\sum_{d|x}\mu(d)
    f(x)=\prod f(p_i^{k_i})
    单独考虑f(p_i^{k_i})=f(p_i)=\mu(1)+\mu(p_i)=0
    因此f(x)=0
  2. \mu^2(n)=\sum_{d^2|n}\mu(d)
    \mu^2(n)=1\Leftrightarrow n是无平方因子数
    i是无平方因子数时,\sum_{d^2|n}\mu(d)=\mu(1)=1
    否则,\sum_{d^2|n}\mu(d)。现在限定d是无平方因子数。
    n的次数\ge 2的质因子集合为S
    那么,\sum_{d^2|n}\mu(d)=\sum_{P\subset S}(-1)^{|P|}=0
  3. n=\sum_{d|n}\varphi(d)
    不会证明
  4. \sum_{i=1}^n[(n,i)=1]i=\frac {n\varphi(n)+[n=1]}{2}
    如果in互质,则(n,i)=(n,n-i)=1,即与n互质的数字成对出现且和为n。因此,对于n\ge 2的情况,有\frac {\varphi(n)}2对,结论就显然了。
  5. D(i)i的约数个数,则D(ij)=\sum_{x|i}\sum_{y|j}[(x,y)=1]
    建立ij的约数d到有序数对(\gcd(i,d),\frac {d}{\gcd(i,d)})的映射。
    那么,一个与一个因数d构成映射的有序数对(a,b)满足a|i,b|j,\gcd(\frac ia,b)=1。而对于任意满足上述条件的数对,都可以反解出d来。
    因此,d与上述数对构成一一映射。
    \begin{aligned} D(ij)=&\sum_{x|i}\sum_{y|j}[(\frac ix,y)=1]\\ =&\sum_{x|i}\sum_{y|j}[(x,y)=1] \end{aligned}

[ICPC 2018 Shenyang Online]Convex Hull

Link

Solution

(n+1)\sum_{i=1}^n\mu^2(i)i^2-\sum_{i=1}^n\mu^2(i)i^3

Min25筛过不去……
根据\mu^2(n)=\sum_{d^2|n}\mu(d)

\begin{aligned}&\sum_{i=1}^ni^2\mu^2(i)\\=&\sum_{i=1}^ni^2\sum_{d^2|i}\mu(d)
\\=&\sum_{d}\mu(d)d^4\sum_{i=1}^{[\frac n{d^2}]}i^2\end{aligned}
\begin{aligned}&\sum_{i=1}^ni^3\mu^2(i)\\=&\sum_{i=1}^ni^3\sum_{d^2|i}\mu(d)
\\=&\sum_{d}\mu(d)d^6\sum_{i=1}^{[\frac n{d^2}]}i^3\end{aligned}

Code

const int N=1e5,XN=N+11;

long long b20=(1ll<<20)-1,P;

long long Add(long long x,long long const &y) {
    return (x+=y)>=P?x-P:x;
}

void AddBy(long long &x,long long const &y) {
    (x+=y)>=P && (x-=P);
}

long long Minus(long long x,long long const &y) {
    return (x-=y)<0?x+P:x;
}

void MinusBy(long long &x,long long const &y) {
    (x-=y)<0 && (x+=P);
}

long long Mul(long long const &x,long long const &y,long long P=::P) {
    return ((((((x>>20)*(y>>20)<<20)+(x>>20)*(y&b20)+(y>>20)*(x&b20))%P)<<20)+(x&b20)*(y&b20))%P;
}

void MulBy(long long &x,long long const &y) {
    x=Mul(x,y);
}

int prime[XN],pcnt,mu[XN];
void Prep() {
    mu[1]=1;
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i]) {
            prime[++pcnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0) {
                mu[i*prime[j]]=0;
                break;
            } else
                mu[i*prime[j]]=-mu[i];
        }
    }
}

long long Sum2(long long x) {
    return Mul(Mul(x,x+1,P*6),x*2+1,P*6)/6;
}

long long Sum3(long long x) {
    x=Mul(x,x+1,P*2)/2;
    return Mul(x,x);
}

int main() {
    Prep();
    long long n;
    while(std::cin>>n>>P) {
        long long Ans,a=0,b=0;
        for(long long d=1;d*d<=n;++d)
            if(mu[d]) {
                long long sqd=d*d%P,tbd=sqd*d%P;
                AddBy(a,Mul(mu[d]==1?mu[d]:P-1,Mul(Mul(sqd,sqd),Sum2(n/(d*d)))));
                AddBy(b,Mul(mu[d]==1?mu[d]:P-1,Mul(Mul(tbd,tbd),Sum3(n/(d*d)))));
            }
        Ans=Minus(Mul(n+1,a),b);
        fout<<Ans<<'\n';
    }
    return 0;
}

Code(TLE)

const int N=1e5,XN=N+11;

long long b20=(1ll<<20)-1,P;

long long Add(long long x,long long const &y) {
    return (x+=y)>=P?x-P:x;
}

void AddBy(long long &x,long long const &y) {
    (x+=y)>=P && (x-=P);
}

long long Minus(long long x,long long const &y) {
    return (x-=y)<0?x+P:x;
}

void MinusBy(long long &x,long long const &y) {
    (x-=y)<0 && (x+=P);
}

long long Mul(long long const &x,long long const &y,long long P=::P) {
    return ((((((x>>20)*(y>>20)<<20)+(x>>20)*(y&b20)+(y>>20)*(x&b20))%P)<<20)+(x&b20)*(y&b20))%P;
}

void MulBy(long long &x,long long const &y) {
    x=Mul(x,y);
}

int prime[XN*2],pcnt;
void Prep() {
    static bool notPrime[XN*2];
    for(int i=2;i<=N*2;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N*2;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
}

namespace Min25 {
    std::function<long long(int,int)> F;

    long long n;
    int lim,psz;

    struct Identifier {
        int id[2][XN],cnt;
        int &operator [](long long x) {
            int &res=x<=lim?id[0][x]:id[1][n/x];
            if(res==0)
                res=++cnt;
            return res;
        }
    }id;

    long long g[XN*2],fps[XN];

    long long H(long long n,int m) {
        if(n<=1 || m>psz)
            return 0;
        long long res=Minus(g[id[n]],fps[m-1]);
        for(int i=m;i<=psz && (long long)prime[i]*prime[i]<=n;++i) {
            long long pt=prime[i],pt1=pt*prime[i];
            for(int t=1;pt1<=n;++t,pt=pt1,pt1*=prime[i])
                AddBy(res,Add(Mul(F(prime[i],t),H(n/pt,i+1)),F(prime[i],t+1)));
        }
        return res;
    }

    long long Solve(long long n,std::function<long long(int,int)> F,std::function<long long(long long)> gInit) {
        static long long kp[XN*2];
        int kpc=0;
        lim=sqrt(n)+0.5,psz=std::upper_bound(prime+1,prime+1+pcnt,lim)-prime;
        for(int i=id.cnt=0;i<=lim;++i)
            id.id[0][i]=id.id[1][i]=0;
        Min25::F=F;
        Min25::n=n;
        for(long long l=1,r;l<=n;l=r+1) {
            r=n/(n/l);
            g[id[kp[++kpc]=n/l]]=gInit(n/l);
        }
        for(int i=1;i<=psz;++i)
            fps[i]=Add(fps[i-1],F(prime[i],1));
        for(int j=1;j<=psz;++j)
            for(int i=1;i<=kpc && (long long)prime[j]*prime[j]<=kp[i];++i)
                MinusBy(g[id[kp[i]]],Mul(F(prime[j],1),Minus(g[id[kp[i]/prime[j]]],fps[j-1])));//XXX
        return H(n,1);
    }
}

long long Sum2(long long x) {
    return Mul(Mul(x,x+1,P*6),x*2+1,P*6)/6;
}

long long Sum3(long long x) {
    x=Mul(x,x+1,P*2)/2;
    return Mul(x,x);
}

int main() {
    Prep();
    long long n;
    while(std::cin>>n>>P) {
        //(n+1)*mu^2(i)*i^2-mu^2(i)*i^3
        long long a=Add(1,Min25::Solve(n,[](int p,int t)->long long { return t==1?Mul(p,p):0; },
                                         [](long long x)->long long { return Minus(Sum2(x),1);})),
                  b=Add(1,Min25::Solve(n,[](int p,int t)->long long { return t==1?Mul(p,Mul(p,p)):0; },
                                         [](long long x)->long long { return Minus(Sum3(x),1);})),
                  Ans=Minus(Mul(n+1,a),b);
        fout<<Ans<<'\n';
    }
    return 0;
}

[CF1067E]Random Forest Rank

Link

给定一个森林,每条边有\frac 12的概率被删去,随机删掉一些边,求剩下的图的临界矩阵的秩的期望值。

Solution

满秩的矩阵的行列式不为0。结合行列式的定义|A|=\sum_{\pi_i}(-1)^{\delta(\pi_i)}\prod a_{\pi_{ij}}。对于一个行列式不为零的邻接矩阵,其必要条件是\exist \pi_i,\text{s.t.}\prod a_{\pi_{ij}}\not=0
\pi_i循环分解,循环内相邻的点都有边连接。然而,对于一个不存在环的森林,这样的循环大小只能为2——即一条边连接的两个点。因此,求出这个森林的最大匹配,设为m,那么就能选出一个大小为2m的点集,其邻接矩阵满足\exist \pi_i,\text{s.t.}\prod a_{\pi_{ij}}\not=0的条件。
另外,选出的这个点集,完美匹配是唯一的,因此,只会有一个\pi_i满足上式子。
由此,问题转化为求期望最大匹配。
然后就不会了 先坑着吧

Code

[loj 6053]简单的函数

Link

Solution

p xor 1并不是积性函数
但是除了p=2的时候,p^1=p-1
p和1都是积性函数,就可以筛出来然后加上2把p=2的情况修正好。

Code

const int P=1e9+7,inv2=500000004;

int Add(int x,int const &y) {
    return (x+=y)>=P?x-P:x;
}

int Minus(int x,int const &y) {
    return (x-=y)<0?x+P:x;
}

int Mul(long long x,int const &y) {
    return x*y%P;
}


const int N=2e5,XN=N+11;

int prime[XN],pcnt,psum[XN],fsum[XN];
void Prep() {
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
    for(int i=1;i<=pcnt;++i)
        psum[i]=Add(psum[i-1],prime[i]);
}

int lim,psz;long long n;
struct Identifier {
    int id[2][XN],cnt;
    int &operator [](long long x) {
        int &res=x<=lim?id[0][x]:id[1][n/x];
        if(res==0)
            res=++cnt;
        return res;
    }
}id;


int F(int p,int t) {
    return p^t;
}

int g[XN*2];

int H(long long n,int m) {
    if(n<=1 || m>psz)
        return 0;
    int res=Minus(g[id[n]],fsum[m-1]);
    for(int i=m;i<=psz && (long long)prime[i]*prime[i]<=n;++i) {
        long long pt=prime[i],pt1=pt*prime[i];
        for(int t=1;pt1<=n;++t,pt=pt1,pt1*=prime[i])
            res=Add(res,Add(Mul(F(prime[i],t),H(n/pt,i+1)),F(prime[i],t+1)));
    }
    return res;
}

int main() {
    Prep();
    static long long kp[XN*2];
    int kpc=0;
    fin>>n;
    lim=sqrt(n)+0.5,psz=std::upper_bound(prime+1,prime+1+pcnt,lim)-prime;
    static int g1[XN*2],g2[XN*2];
    for(long long l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        kp[++kpc]=n/l;
        g1[id[n/l]]=Minus((n/l%P)*(n/l%P+1)%P*inv2%P,1);
        g2[id[n/l]]=Minus(n/l%P,1);
    }
    for(int j=1;j<=psz;++j)
        for(int i=1;i<=kpc && (long long)prime[j]*prime[j]<=kp[i];++i) {
            g1[id[kp[i]]]=Minus(g1[id[kp[i]]],Mul(prime[j],Minus(g1[id[kp[i]/prime[j]]],psum[j-1])));
            g2[id[kp[i]]]=Minus(g2[id[kp[i]]],Minus(g2[id[kp[i]/prime[j]]],j-1));
        }
    for(int i=1;i<=psz;++i)
        fsum[i]=Add(fsum[i-1],prime[i]^1);
    for(int i=1;i<=kpc;++i)
        g[i]=Add(Minus(g1[i],g2[i]),kp[i]>=2?2:0);
    fout<<Add(H(n,1),1)<<'\n';
    return 0;
}

[ICPC 2018 Xuzhou Online]Easy Math

Link

Min25

Solution

\sum_{i=1}^n\mu(ix)=\mu(x)\sum_{i=1}^n\mu(i)[\gcd(i,x)=1]

考虑计算

\begin{aligned}
&\sum_{i=1}^nf(i)\\
&f(i)=\mu(i)[\gcd(i,n)=1]
\end{aligned}

首先,

f(i)=\begin{cases}
1,&i=1\\
-1,&i=p{\ \rm and\ }p\not|x\\
0,&i=p^k{\ \rm or\ }p|x
\end{cases}

然后,分\gcd(i,x)的不同情况讨论,容易知道f(i)是积性函数。
因此可以使用min25筛法了。

Tips

还是不必拘泥于g的形式……

Code

const int N=5e4,XN=N+11;

int prime[XN],pcnt;
void Prep() {
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
}

int n,lim,psz;long long x;
struct Identifier {
    int id[2][XN],cnt;
    int &operator [](int x) {
        int &res=x<=lim?id[0][x]:id[1][n/x];
        if(res==0)
            res=++cnt;
        return res;
    }
}id;

long long dv[100];int dt;
int Mu(long long x) {
    int res=1;long long cx=x;
    for(long long d=2;d*d<=cx && x!=1;++d)
        if(x%d==0) {
            x/=d;
            dv[++dt]=d;
            if(x%d==0)
                return 0;
            else
                res*=-1;
        }
    if(x!=1) {
        res*=-1;
        dv[++dt]=x;
    }
    return res;
}

int F(int p,int t) {
    return t==1 && x%p?-1:0;
}

int g[XN*2],fsum[XN];

long long H(int n,int m) {
    if(n<=1 || m>psz)
        return 0;
    long long res=g[id[n]]-fsum[m-1];
    for(int i=m;i<=psz && prime[i]*prime[i]<=n;++i) {
        long long pt=prime[i],pt1=pt*prime[i];
        for(int t=1;pt1<=n;++t,pt=pt1,pt1*=prime[i])
            res+=F(prime[i],t)*H(n/pt,i+1)+F(prime[i],t+1);
    }
    return res;
}

int main() {
    Prep();
    static int kp[XN*2],kpc;
    fin>>n>>x;
    int mun=Mu(x);
    if(mun==0)
        fout<<0<<'\n';
    else {
        lim=sqrt(n)+0.5,psz=std::upper_bound(prime+1,prime+1+pcnt,lim)-prime;
        for(int l=1,r;l<=n;l=r+1) {
            r=n/(n/l);
            g[id[kp[++kpc]=n/l]]=n/l-1;
        }
        for(int j=1;j<=psz;++j)
            for(int i=1;i<=kpc && prime[j]*prime[j]<=kp[i];++i)
                g[id[kp[i]]]-=g[id[kp[i]/prime[j]]]-(j-1);//(1-(j-Count(j)));
        for(int i=1;i<=psz;++i)
            fsum[i]=fsum[i-1]+F(prime[i],1);
        for(int i=1;i<=kpc;++i)
            g[i]=std::upper_bound(dv+1,dv+1+dt,kp[i])-dv-1-g[i];
        fout<<mun*(H(n,1)+1)<<'\n';
    }
    return 0;
}

递归求解+杜教筛

网赛的时候写出了这个式子,但是感觉不大能过就一直在自闭……

Solution

\begin{aligned}
&\sum_{i=1}^n\mu(im)\\=&\mu(m)\sum_{i=1}^n\mu(i)[(i,m)=1]\\=&\mu(m)\sum_{i=1}^n\mu(i)\sum_{d|i,d|m}\mu(d)\\=&\mu(m)\sum_{d|m}\mu(d)\sum_{i=1}^{[\frac nd]}\mu(id)\end{aligned}

S(n,m)=\sum_{d|m}\mu(d)S([\frac nd],d)

Tips

  1. 搞清楚哪个是要算的函数……一开始多乘了一个\mu愣是没看出来。
  2. 这种递归的式子,与其对着自闭,不如大胆写,说不定就过了呢。

Code

const int N=2e6,XN=N+11;

int mert[XN],prime[XN],pcnt;
void Prep() {
    static bool notPrime[XN];
    mert[1]=1;
    for(int i=2;i<=N;++i) {
        if(!notPrime[i]) {
            prime[++pcnt]=i;
            mert[i]=-1;
        }
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0) {
                mert[i*prime[j]]=0;
                break;
            } else
                mert[i*prime[j]]=-mert[i];
        }
    }
    for(int i=2;i<=N;++i)
        mert[i]+=mert[i-1];
}

long long Mert(int n) {
    if(n<=N)
        return mert[n];
    else {
        static std::unordered_map<int,long long> map;
        if(map.count(n))
            return map[n];
        long long &res=map[n]=1;
        for(int l=2,r;l<=n;l=r+1) {
            r=n/(n/l);
            res-=Mert(n/l)*(r-l+1);
        }
        return res;
    }
}

int Mu(long long x) {
    if(x<=N)
        return mert[x]-mert[x-1];
    else {
        static std::unordered_map<long long,int> map;
        if(map.count(x))
            return map[x];
        int &mu=map[x]=1;
        long long m=x;
        for(int i=1;(long long)prime[i]*prime[i]<=m && x!=1;++i)
            if(x%prime[i]==0 && (mu=-mu,(x/=prime[i])%prime[i]==0))
                return mu=0;
        if(x!=1)
            mu=-mu;
        return mu;
    }
}

long long S(int n,long long m) {
    if(n==0)
        return 0;
    else if(m==1)
        return Mert(n);
    else {
        static std::map<std::pair<int,long long>,long long> map;
        if(map.count(std::pair<long long,int>(n,m)))
            return map[std::pair<long long,int>(n,m)];
        long long &res=map[std::pair<long long,int>(n,m)];
        if(Mu(m)==0)
            return res=0;
        for(long long d=1;d*d<=m;++d)
            if(m%d==0) {
                res+=Mu(d)*S(n/d,d);
                if(m!=d*d)
                    res+=Mu(m/d)*S(n/(m/d),m/d);
            }
        return res*=Mu(m);
    }
}

int main() {
    Prep();
    long long m;int n;
    fin>>n>>m;
    fout<<S(n,m)<<'\n';
    return 0;
}

SPOJ DIVCNTK

Link

Solution

使用min_25筛
由于求h的过程都只会用到[\frac ni]集合中的点,那么求g的过程中把[\frac ni]构成的集合求出来,对于这个集合的所有点分层递推(递推过程中需要用到的g也属于[\frac ni]集合),就可以求出g(i,|P|)
然后递归求h,并不需要记忆化。边界见代码。

Code

typedef unsigned long long qword;

const int N=2e5,XN=N+11;

int prime[XN],pcnt;
void Prep() {
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
}

long long n,k;
int lim,psz;

struct Identifier {
    int id[2][XN*2],cnt;
    int &operator [](long long x) {
        int &res=x<=lim?id[0][x]:id[1][n/x];
        if(res==0)
            res=++cnt;
        return res;
    }
}id;

qword g[XN*2];

qword H(long long n,int m) {
    if(n<=1 || m>psz)
        return 0;
    qword res=(k+1)*(g[id[n]]-m+1);
    for(int i=m;i<=psz && (long long)prime[i]*prime[i]<=n;++i) {
        long long pt=prime[i];
        for(int t=1;pt*prime[i]<=n;++t,pt*=prime[i])
            res+=H(n/pt,i+1)*(k*t+1)+(k*(t+1)+1);//id[n/pt]
    }
    return res;
}

void Work() {
    lim=sqrt(n)+0.5;
    id.cnt=0;
    for(int i=0;i<=lim;++i)
        id.id[0][i]=id.id[1][i]=0;
    psz=std::lower_bound(prime+1,prime+1+pcnt,lim)-prime;
    static long long kp[XN*2];int kpc=0;
    for(long long l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        g[id[kp[++kpc]=n/l]]=n/l-1;
    }
    for(int j=1;j<=psz;++j)
        for(int i=1;i<=kpc && (long long)prime[j]*prime[j]<=kp[i];++i)
            g[id[kp[i]]]-=g[id[kp[i]/prime[j]]]-(j-1);
    std::cout<<H(n,1)+1<<'\n';
}

int main() {
    Prep();
    int T;fin>>T;
    while(T--) {
        fin>>n>>k;
        Work();
    }
    return 0;
}