[Gym 102114F]Fireflies

Link

第一步的题意转化就并没想到(根本没往子集的方向去想)
Consider all the lattice points in a hypercube \ { (x_1,x_2,…,x_n)|0\le x_i < p_i \ } . Find a maximal subset such that there are no two points (x_1,x_2,…,x_n),(y_1,y_2,…,y_n) meeting the condition x_i≤y_i for all i. Report its size modulo 10^9+7.

Solution

对于任意集合S,如果它的一个子集集合F满足\forall X,Y\in F,X\not\subset Y,Y\not\subset X,则我们把这样的子集F称为Sperner族。
一个集合的极大Sperner族的元素个数为\binom{|S|}{\lfloor\frac {|S|}{2}\rfloor},是为Sperner定理。达到这个数目只需要取出所有\lfloor\frac {|S|}{2}\rfloor元子集。
Sperner定理可以直接用于解决问题在p_i=2时的情形,将某一位坐标的0,1转化为在集合中的选与不选。

Sperner 定理的一种推广可以给出本问题的答案:选择满足所有坐标之和等于 M = \lfloor \frac{1}{2} \sum_{i = 1}^{n}(p_i-1) \rfloor的位置即可达到最大的数量。

最终问题变成了,求\sum x_i=M,x_i\in[0,p_i)的解数,也就是\sum x_i=M+n,x_i\in[1,p_i]的解数。

由容斥原理,答案为\sum_{T\subset S}(-1)^{|T|}\binom{M-1-\sum_{i\subset T}p_i}{n-1}

使用meet-in-middle策略。
由于\binom{a-x}{b}可以写成一个b次的多项式,而n-1又很小,因此可以处理出左侧集合的关于右侧集合选出元素之和x对答案贡献的多项式的前缀和,枚举右侧集合的元素,然后利用左侧集合多项式的前缀和来处理答案。

\begin{aligned}
&\sum_{T\subset S}(-1)^{|T|}\binom{M-1-\sum_{i\subset T}p_i}{n-1}\\
=&\sum_{T_1\in S_1}\sum_{T_0\in S_0}(-1)^{|T_0|+|T_1|}\binom{M-1-\sum_{i\in T_0}p_i-\sum_{i\in T_1}p_i}{n-1}\\
=&\sum_{T_1\in S_1}(-1)^{|T_1|}\sum_{T_0\in S_0}(-1)^{|T_0|}F(\sum_{i\in T_1}p_i)\\
\end{aligned}

Code

#include <bits/stdc++.h>

const int XN=40,P=1e9+7;

int Pow(long long base,int v) {
    long long res;
    for(res=1;v;v>>=1,(base*=base)%=P)
        if(v&1)
            (res*=base)%=P;
    return res;
}

std::vector<int> Conv(std::vector<int> const &a,std::vector<int> const &b) {
    std::vector<int> c(a.size()+b.size()-1,0);
    for(size_t i=0;i<a.size();++i)
        for(size_t j=0;j<b.size();++j)
            (c[i+j]+=1ll*a[i]*b[j]%P)%=P;
    return c;
}

int main() {
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);

    freopen("input","r",stdin);

    int ifact[XN]={1};
    for(int i=1;i<XN;++i)
        ifact[i]=1ll*ifact[i-1]*Pow(i,P-2)%P;

    int T;std::cin>>T;

    while(T--) {
        int n;std::cin>>n;
        std::vector<int> p(n),s0,s1;
        long long M=0;
        for(int i=0;i<n;++i) {
            std::cin>>p[i];
            M+=p[i]-1;
            (i<(n+1)/2?s0:s1).push_back(p[i]);
        }
        (M/=2)+=n;
        std::vector<std::pair<long long,int>> v0;
        for(int S=0;S<(1<<s0.size());++S) {
            long long sum=0;
            for(size_t i=0;i<s0.size();++i)
                if(S>>i&1)
                    sum+=s0[i];
            v0.push_back({sum,__builtin_popcount(S)&1?-1:1});
        }
        std::sort(v0.begin(),v0.end());
        std::vector<std::vector<int>> f(v0.size(),std::vector<int>(n,0));
        for(size_t i=0;i<v0.size();++i) {
            std::vector<int> p{1};
            for(int j=0;j<n-1;++j)
                p=Conv(p,{(int)((M-1-v0[i].first-j)%P),-1});//XXX
            for(int &x : p)
                x=1ll*x*ifact[n-1]*v0[i].second%P;
            for(int j=0;j<n;++j)
                f[i][j]=((i?f[i-1][j]:0)+p[j])%P;
        }
        //需要枚举右集合的状态 一下子统计上所有左集合中和+和<=M-1的状态的
        //(-1)^J0 * (-1)^J1*F(X0+X1)
        int Ans=0;
        for(int S=0;S<(1<<s1.size());++S) {
            long long sum=0;
            for(size_t i=0;i<s1.size();++i)
                if(S>>i&1)
                    sum+=s1[i];
            if(sum<=M-1) {
                int p=std::upper_bound(v0.begin(),v0.end(),std::make_pair(M-1-sum,1))-v0.begin()-1,
                    c=__builtin_popcount(S)&1?-1:1;
                int res=0;
                sum%=P;
                for(int i=n-1;i>=0;--i)
                    res=(res*sum+f[p][i])%P;
                (Ans+=c*res)%=P;
            }
        }
        if(Ans<0)
            Ans+=P;
        std::cout<<Ans<<'\n';
    }
    return 0;
}

[LOJ2720][NOI2018]君の名は。

Link

Solution

首先求出TS的对应区间内每个前缀能匹配到的最长长度,用后缀自动机维护right集合比较好做。
然后将T反转求一下后缀数组,每个后缀(原串的前缀)给答案带来贡献就是n-sa[i]+1-std::max(height[i],match[sa[i]])

Code

#include <bits/stdc++.h>

template <class T>
bool Enlarge(T &a,T const &b) {
    return a<b?a=b,1:0;
}

const int XN=5e5+11;

struct SegTree {
    static int n;

    struct Node {
        Node *son[2];
        int sum;

        Node(int sum):sum(sum) {
            son[0]=son[1]=0;
        }

        void Up() {
            sum=0;
            for(int d=0;d<2;++d)
                if(son[d])
                    sum+=son[d]->sum;
        }
    }*root;

    SegTree(Node *root=0):root(root) {}

    void Add(int goal,int v) {
        std::function<void(Node*&,int,int)> Add=[&](Node *&pos,int L,int R) {
            pos=pos?new Node(*pos):new Node(0);
            if(L==R)
                pos->sum+=v;
            else {
                int M=(L+R)/2;
                goal<=M?Add(pos->son[0],L,M)
                       :Add(pos->son[1],M+1,R);     
                pos->Up();
            }
        };

        Add(root,1,n);
    }

    int Sum(int l,int r) {
        int res=0;

        std::function<void(Node*,int,int)> GetSum=[&](Node *pos,int L,int R) {
            if(!pos)
                return;
            if(l<=L && R<=r)
                res+=pos->sum;
            else {
                int M=(L+R)/2;
                if(l<=M)
                    GetSum(pos->son[0],L,M);
                if(M+1<=r)
                    GetSum(pos->son[1],M+1,R);
            }
        };

        GetSum(root,1,n);
        return res;
    }

    friend SegTree Merge(SegTree const &x,SegTree const &y) {
        std::function<Node*(Node*,Node*)> Merge=[&](Node *x,Node *y)->Node* {
            if(!x || !y)
                return x?x:y;
            else {
                Node *pos=new Node(0);
                pos->son[0]=Merge(x->son[0],y->son[0]);
                pos->son[1]=Merge(x->son[1],y->son[1]);
                pos->Up();
                return pos;
            }
        };

        return SegTree(Merge(x.root,y.root));
    }
};

int SegTree::n;

struct SuffixAutomata {
    struct Node {
        Node *son[26];
        Node *par;
        int deg,maxRight;
        SegTree T;

        Node(int maxRight=0):par(0),deg(0),maxRight(maxRight) {
            memset(son,0,sizeof(son));
        }

        void *operator new(size_t) {
            static Node *pool=(Node*)malloc(XN*2*sizeof(Node)),*mem=pool;
            nodes.push_back(mem);
            return mem++;
        }
    }*root,*last;

    static std::vector<Node*> nodes;

    std::vector<Node*> loc;

    SuffixAutomata():root(new Node),last(root) {}

    void Extend(int x) {
        Node *p=last,*nx=new Node(last->maxRight+1);
        for(;p && !p->son[x];p->son[x]=nx,p=p->par);
        if(p==0) {
            nx->par=root;
        } else {
            Node *ox=p->son[x];
            if(p->maxRight+1==ox->maxRight) {
                nx->par=ox;
            } else {
                Node *o=new Node(*ox);
                o->maxRight=p->maxRight+1;
                ox->par=nx->par=o;
                for(;p && p->son[x]==ox;p->son[x]=o,p=p->par);
            }
        }
        loc.push_back(last=nx);
    }

    void Solve() {
        SegTree::n=loc.size();
        for(size_t i=0;i<loc.size();++i)
            loc[i]->T.Add(i+1,1);
        for(auto p : nodes)
            if(p->par)
                p->par->deg++;
        std::queue<Node*> Q;
        for(auto p : nodes)
            if(!p->deg)
                Q.push(p);
        while(!Q.empty()) {
            Node *p=Q.front();Q.pop();
            if(p->par) {
                p->par->T=Merge(p->par->T,p->T);
                if(--p->par->deg==0)
                    Q.push(p->par);
            }
        }
    }

    std::vector<int> Match(int l,int r,std::string const &s) {
        Node *pos=root;
        std::vector<int> res;
        int len=0;
        for(char c : s) {
            int x=c-'a';
            while(pos!=root && !pos->son[x])
                len=(pos=pos->par)->maxRight;
            if(pos->son[x]) {
                ++len;
                pos=pos->son[x];
            }
            while(len && !pos->T.Sum(l+len-1,r))
                if(--len==pos->par->maxRight)
                    pos=pos->par;
            res.push_back(len);
        }
        return res;
    }
};

std::vector<SuffixAutomata::Node*> SuffixAutomata::nodes;

std::vector<int> SAIS(std::vector<int> &s,int sig) {
    enum Type { L,S };

    int n=s.size();
    std::vector<Type> type(n);
    std::vector<int> sa(n),cnt(sig),lp(sig),sp(sig);

    auto LMS=[&](int x) {
        return x>0 && type[x]==S && type[x-1]==L;
    };

    auto Induce=[&]() {
        for(int i=0;i<sig;++i)
            lp[i]=i?cnt[i-1]:0;
        for(int i=0;i<n;++i)
            if(sa[i]>0 && type[sa[i]-1]==L)
                sa[lp[s[sa[i]-1]]++]=sa[i]-1;
        for(int i=0;i<sig;++i)
            sp[i]=cnt[i]-1;
        for(int i=n-1;i>0;--i)
            if(sa[i]>0 && type[sa[i]-1]==S)
                sa[sp[s[sa[i]-1]]--]=sa[i]-1;
    };

    type[n-1]=S;
    for(int i=n-2;i>=0;--i)
        type[i]=s[i]<s[i+1] || (s[i]==s[i+1] && type[i+1]==S)?S:L;
    for(int x : s)
        cnt[x]++;
    for(int i=1;i<sig;++i)
        cnt[i]+=cnt[i-1];
    std::vector<int> lms;
    for(int i=0;i<n;++i)
        if(LMS(i))
            lms.push_back(i);
    for(int i=0;i<sig;++i)
        sp[i]=cnt[i]-1;
    std::fill(sa.begin(),sa.end(),-1);
    for(int i : lms)
        sa[sp[s[i]]--]=i;
    Induce();
    std::vector<int> id(n,-1);
    int c=0;
    for(int i=0,last=-1;i<n;++i)
        if(LMS(sa[i])) {
            auto Equal=[&](int x,int y)->bool {
                do if(s[x++]!=s[y++]) return 0;
                while(!LMS(x) && !LMS(y));
                return s[x]==s[y];
            };
            id[sa[i]]=last!=-1 && Equal(sa[i],last)?c-1:c++;
            last=sa[i];
        }

    std::vector<int> s1;
    for(int i=0;i<n;++i)
        if(id[i]!=-1)
            s1.push_back(id[i]);
    std::vector<int> sa1(c);
    if(c==(int)lms.size())
        for(int i=0;i<c;++i)
            sa1[s1[i]]=i;
    else
        sa1=SAIS(s1,c);
    for(int i=0;i<sig;++i)
        sp[i]=cnt[i]-1;
    std::fill(sa.begin(),sa.end(),-1);
    for(int i=lms.size()-1;i>=0;--i)
        sa[sp[s[lms[sa1[i]]]]--]=lms[sa1[i]];
    Induce();
    return sa;
}

long long Calc(std::string const &s,std::vector<int> &match) {
    int n=s.size();

    std::vector<int> a(s.size());
    for(int i=0;i<n;++i)
        a[i]=s[i]-'a'+1;
    std::reverse(a.begin(),a.end());
    match.push_back(-1);
    std::reverse(match.begin(),match.end());
    a.push_back(0);

    auto sa=SAIS(a,26+1);
    a.insert(a.begin(),-1);
    std::vector<int> rank(n+1),height(n+1);
    for(int i=1;i<=n;++i)
        rank[++sa[i]]=i;
    for(int i=1,len=0;i<=n;++i)
        if(rank[i]!=1) {
            int j=sa[rank[i]-1];
            while(a[i+len]==a[j+len]) ++len;
            height[rank[i]]=len;
            if(len) len--;
        }
    long long res=0;
    for(int i=1;i<=n;++i)
        res+=n-sa[i]+1-std::max(height[i],match[sa[i]]);

    return res;
}

int main() {
    freopen("name.in","r",stdin);
    freopen("name.out","w",stdout);
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);
    SuffixAutomata sam;
    std::string s;
    std::cin>>s;
    for(char c : s)
        sam.Extend(c-'a');
    sam.Solve();
    int q;std::cin>>q;
    while(q--) {
        std::string t;
        int l,r;
        std::cin>>t>>l>>r;
        auto v=sam.Match(l,r,t);
        long long Ans=Calc(t,v);
        std::cout<<Ans<<'\n';
    }
    return 0;
}

胡扯SAIS

对于每一个后缀\text{suffix}(i),当\text{suffix}(i) < \text{suffix}(i + 1)时,是S型后缀,否则为L型后缀。
显然,\text{suffix}(i)S型后缀的充要条件为:
S_i < S_{i+1}S_i=S_{i + 1}\text{suffix}(i + 1)S型后缀。由此可以倒着一遍把每个后缀的类型推出来。为了方便处理,在字符串末尾加一个极小的字符,并且令只包含这一个字符的后缀类型为S型。
对于满足\text{suffix}(i-1)L\text{suffix}(i)S的后缀,记为LMS(Leftmost-S)后缀,并且记相邻两个LMS后缀之间(左右端点都包含)的子串为LMS子串。
popoqqq举例

p o p o q q q #
L S L S L L L S
* * *

LMS子串为opooqqq##
LMS子串排好字典序并标号,得到一个新的串。在当前例子中,LMS子串排序之后#opooqqq#,依次标号为0,1,2并且替换到原串的对应位置,得到的就是1,2,0。在这个新串中,后缀的排序,就是每个后缀在原串的左端点对应的LMS后缀的顺序。比如,在当前例子中,1,2,0的后缀排序之后为01,2,02,0,对应回原串,代表这三个LMS后缀的顺序为#opoqqq#oqqq#
得到LMS子串压成的新串之后就可以递归了,现在就考虑在已知LMS后缀的顺序的情况下如何求出所有后缀的顺序,也就是整个算法的精髓——诱导排序。
对于首字母相同的所有后缀,在最终的排序序列中,一定是先出现所有L型后缀,然后出现所有S型后缀。因为,对于S_i=S_j,\text{suffix}(i)\in \text{S-Type},\text{suffix}(j)\in \text{L-Type},一定有,\text{suffix}(i)>\text{suffix}(j)
证明:如果S_i=S_{i+1}S_j=S_{j+1},那么可以比较\text{suffix}(i+1),\text{suffix}(j+1)。否则,一定有S_i < S_{i+1}S_j>S_{j+1},可以推出结论。

对于每个首字母,维护两个桶,分别用于存储以此为首字母的LS型后缀。
首先,将LMS后缀按顺序放在S型桶内,在同一个桶内放的具体位置无所谓(?),只要保持有序就好了。
然后依次对于所有的字符,先遍历S型桶,再遍历L型桶,如果桶中的\text{suffix}(i)满足\text{suffix}(i-1)L型,就将\text{suffix}(i-1)放到对应的L型桶中。所有的L型后缀都会被按照顺序依次插入到L型的桶内。
L型后缀都有序了,那么按照类似的操作可以再将S型后缀都依次插入桶中。
最终就排好了顺序。

还有LMS子串的排序没有解决,其思路也是使用诱导排序。做法是在一开始将LMS后缀按照任意顺序插入到对应S型桶内,诱导完之后LMS后缀相对之间就是有序的了。正确性还没咋理解,先坑着…

const int XN=1e6+11;

std::vector<int> SAIS(std::vector<int> const &s,int sig) {
    int n=s.size();
    std::vector<bool> st(n);
    std::vector<int> sa(n),cnt(sig),lp(sig),sp(sig);

    auto Induce=[&]() {
        for(int i=0;i<sig;++i) {
            lp[i]=i?cnt[i-1]:0;
            sp[i]=cnt[i]-1;
        }
        for(int i=0;i<n;++i)
            if(sa[i]>0 && !st[sa[i]-1])
                sa[lp[s[sa[i]-1]]++]=sa[i]-1;
        for(int i=n-1;i>0;--i)
            if(sa[i]>0 && st[sa[i]-1])
                sa[sp[s[sa[i]-1]]--]=sa[i]-1;
    };

    st[n-1]=1;
    for(int i=n-2;i>=0;--i)
        st[i]=s[i]<s[i+1] || (s[i]==s[i+1] && st[i+1]);
    for(int i=0;i<n;cnt[s[i++]]++);
    for(int i=1;i<sig;cnt[i]+=cnt[i-1],++i);
    for(int i=0;i<sig;sp[i]=cnt[i]-1,++i);
    std::fill(sa.begin(),sa.end(),0);
    std::vector<int> lms,id(n,-1);
    for(int i=1;i<n;++i)
        if(!st[i-1] && st[i]) {
            id[i]=lms.size();
            lms.push_back(sa[sp[s[i]]--]=i);
        }
    Induce();
    int c=0;
    static int t[XN];
    for(int i=0,j,k=-1;i<n;++i)
        if((j=id[sa[i]])!=-1) {
            t[j]=k!=-1 && std::equal(s.begin()+lms[j],s.begin()+lms[j+1],s.begin()+lms[k])?c-1:c++;
            k=j;
        }

    std::vector<int> sa1(c);
    if(c==(int)lms.size())
        for(int i=0;i<c;sa1[t[i]]=i,++i);
    else 
        sa1=SAIS(std::vector<int>(t,t+lms.size()),c);
    for(int i=0;i<sig;sp[i]=cnt[i]-1,++i);
    std::fill(sa.begin(),sa.end(),0);
    for(int i=lms.size()-1;i>=0;--i)
        sa[sp[s[lms[sa1[i]]]]--]=lms[sa1[i]];
    Induce();
    return sa;
}

参考

用后缀数组构造后缀树

后缀树就是把所有后缀插入到trie之后,把只有一个分支的单链合并到一起,使得状态数减小至O(n)
下图为bananas的后缀树

定义边的长度为边上字符串的长度,点的深度为点到根的边权和。
容易发现,叶节点的DFS序就是按照叶节点对应的后缀字典序排列的,且相邻叶节点的LCA的深度为相邻两个后缀LCP的长度。
那么构造的时候就将后缀从大到小加入,同时使得这棵树满足上述性质即可。代码一看就懂了。

int n;
std::vector<int> sa,rank,height;

struct Edge {
    int to,v;
};

std::vector<Edge> G[XN];
int root,pcnt,last,dis[XN],fa[XN],size[XN];

int NewNode() {
    G[++pcnt].clear();
    size[pcnt]=fa[pcnt]=dis[pcnt]=0;
    return pcnt;
}

void AddEdge(int x,int y,int v) {
    G[x].push_back({y,v});
    dis[y]=dis[fa[y]=x]+v;
}

void Build() {
    pcnt=0;
    root=NewNode();
    AddEdge(root,last=NewNode(),n-sa[1]+1);
    size[last]=1;
    for(int i=2;i<=n;++i) {
        int p=last,np=NewNode();
        while(dis[p]>height[i])
            p=fa[p];
        if(dis[p]!=height[i]) {
            int d=height[i]-dis[p],o=NewNode();
            auto e=G[p].back();
            G[p].pop_back();
            AddEdge(p,o,d);
            AddEdge(o,e.to,e.v-d);
            p=o;
        }
        AddEdge(p,np,n-sa[i]+1-height[i]);
        size[last=np]=1;
    } 
}

参考https://blog.csdn.net/u011542204/article/details/51231962

[Gym 102114J]Just So You Know

Link

字符集出这么大的范围卡后缀自动机有点令人不适。

Solution

题意为给定一个字符串,找出所有不同的子串并统计其出现次数,并以出现次数为频率将所有的子串做哈夫曼编码的编码长度。
第一步,要求出,出现次数为i的子串的个数c_i。如果使用后缀自动机,可以对于每个状态求出Right集合R的大小,那么一个状态p就会对i=|R|c_i带来p->maxRight-p->par->maxRight的贡献。
但是由于\Sigma特别大,这样会MLE。
tls给的做法是用后缀数组构建出后缀树然后在后缀树上做统计。考虑后缀树的一条长度为l的边,那么从根到这条边上的l个字符可以产生l个子串。假设这条边下方有k个后缀的终点,那么就会有k个后缀满足这l个子串是他们的前缀,也就是说,这l个子串都出现了k次。因此,就会对c_k带来l的贡献。
接下来,就是求出哈夫曼编码的权值。
哈夫曼编码,在频率升序存储的前提下,有线性的做法。开两个队列,一开始将所有字符的频率扔到Q_1,然后循环到\rm{size}(Q_1)+\rm{size}(Q_2)=1为止。每次循环,从两个队列的队首取两次最小值,然后将和插到Q_2末尾。显然,取出的两个值的和是不下降的,Q_2序列一定是有序的,算法正确性显然。
在本问题中,对于i\le n的部分,直接扫一遍c_i数组,模拟堆的操作就好了。对于合并产生的值v,如果v\le n,就还是将频率累加到c_v待以后处理,否则,就将v插入到Q_1中。
最终,Q_1中的元素是有序的。另外,因为Q_1中每个元素的权值都>n,且所有元素的权值和\le \frac {n(n+1)} 2,所以Q_1中的元素个数为O(n)
接下来套用如上算法就可以O(n)解决求最优编码部分。
后缀数组转后缀树的复杂度为O(n),采用O(n)的后缀数组实现,则总复杂度为O(n)

Code

#include <bits/stdc++.h>

const int XN=2e6+11;

//(v,c) v->O(n)
long long times[XN];

std::vector<int> SAIS(std::vector<int> &s,int sig) {
    enum Type { L,S };

    int n=s.size();
    std::vector<Type> type(n);
    std::vector<int> sa(n),cnt(sig),lp(sig),sp(sig);

    auto LMS=[&](int x) {
        return x>0 && type[x]==S && type[x-1]==L;
    };

    auto Induce=[&]() {
        for(int i=0;i<sig;++i)
            lp[i]=i?cnt[i-1]:0;
        for(int i=0;i<n;++i)
            if(sa[i]>0 && type[sa[i]-1]==L)
                sa[lp[s[sa[i]-1]]++]=sa[i]-1;
        for(int i=0;i<sig;++i)
            sp[i]=cnt[i]-1;
        for(int i=n-1;i>0;--i)
            if(sa[i]>0 && type[sa[i]-1]==S)
                sa[sp[s[sa[i]-1]]--]=sa[i]-1;
    };

    type[n-1]=S;
    for(int i=n-2;i>=0;--i)
        type[i]=s[i]<s[i+1] || (s[i]==s[i+1] && type[i+1]==S)?S:L;
    for(int x : s)
        cnt[x]++;
    for(int i=1;i<sig;++i)
        cnt[i]+=cnt[i-1];
    std::vector<int> lms;
    for(int i=0;i<n;++i)
        if(LMS(i))
            lms.push_back(i);
    for(int i=0;i<sig;++i)
        sp[i]=cnt[i]-1;
    std::fill(sa.begin(),sa.end(),-1);
    for(int i : lms)
        sa[sp[s[i]]--]=i;
    Induce();
    std::vector<int> id(n,-1);
    int c=0;
    for(int i=0,last=-1;i<n;++i)
        if(LMS(sa[i])) {
            auto Equal=[&](int x,int y)->bool {
                do if(s[x++]!=s[y++]) return 0;
                while(!LMS(x) && !LMS(y));
                return s[x]==s[y];
            };
            id[sa[i]]=last!=-1 && Equal(sa[i],last)?c-1:c++;
            last=sa[i];
        }

    std::vector<int> s1;
    for(int i=0;i<n;++i)
        if(id[i]!=-1)
            s1.push_back(id[i]);
    std::vector<int> sa1(c);
    if(c==(int)lms.size())
        for(int i=0;i<c;++i)
            sa1[s1[i]]=i;
    else
        sa1=SAIS(s1,c);
    for(int i=0;i<sig;++i)
        sp[i]=cnt[i]-1;
    std::fill(sa.begin(),sa.end(),-1);
    for(int i=lms.size()-1;i>=0;--i)
        sa[sp[s[lms[sa1[i]]]]--]=lms[sa1[i]];
    Induce();
    return sa;
}

int n;
std::vector<int> sa,rank,height;

struct Edge {
    int to,v;
};

std::vector<Edge> G[XN];
int root,pcnt,last,dis[XN],fa[XN],size[XN];

int NewNode() {
    G[++pcnt].clear();
    size[pcnt]=fa[pcnt]=dis[pcnt]=0;
    return pcnt;
}

void AddEdge(int x,int y,int v) {
    G[x].push_back({y,v});
    dis[y]=dis[fa[y]=x]+v;
}

void Build() {
    pcnt=0;
    root=NewNode();
    AddEdge(root,last=NewNode(),n-sa[1]+1);
    size[last]=1;
    for(int i=2;i<=n;++i) {
        int p=last,np=NewNode();
        while(dis[p]>height[i])
            p=fa[p];
        if(dis[p]!=height[i]) {
            int d=height[i]-dis[p],o=NewNode();
            auto e=G[p].back();
            G[p].pop_back();
            AddEdge(p,o,d);
            AddEdge(o,e.to,e.v-d);
            p=o;
        }
        AddEdge(p,np,n-sa[i]+1-height[i]);
        size[last=np]=1;
    } 
}

void DFS(int pos) {
    for(Edge &e : G[pos]) {
        DFS(e.to);
        size[pos]+=size[e.to];
        times[size[e.to]]+=e.v;
    }
}

void Work() {
    std::cin>>n;
    std::vector<int> a(n+1);
    for(int i=1;i<=n;++i) {
        std::cin>>a[i-1];
        a[i-1]++;
        times[i]=0;
    }
    a[n]=0;
    sa=SAIS(a,100+1);
    a.insert(a.begin(),0);
    rank.resize(n+1);
    height.resize(n+1);
    for(int i=1;i<=n;++i)
        rank[++sa[i]]=i;
    for(int i=1,len=0;i<=n;++i)
        if(rank[i]!=1) {
            int j=sa[rank[i]-1];
            while(a[i+len]==a[j+len]) ++len;
            height[rank[i]]=len;
            if(len) len--;
        }

    Build();
    DFS(root);

    long long sum=0,cnt=(long long)n*(n+1)/2;
    std::deque<long long> Q1,Q2;

    for(int i=1;i<=n;)
        if(times[i]) {
            if(times[i]==1) {
                int j=i+1;
                while(j<=n && !times[j])
                    ++j;
                if(j<=n) {
                    times[j]--;
                    if(i+j<=n)
                        times[i+j]++;   
                    else
                        Q1.push_back(i+j);
                    sum+=i+j;
                } else
                    Q1.push_front(i);
                times[i++]--;
            } else {
                sum+=i*2ll*(times[i]/2);
                if(i*2<=n)
                    times[i*2]+=times[i]/2;
                else
                    for(int t=times[i]/2;t--;Q1.push_back({i*2}));
                times[i]%=2;
            }
        } else
            i++;
    while(Q1.size()+Q2.size()>1) {
        long long a[2];
        for(int i=0;i<2;++i) {
            if(Q1.size() && (Q2.empty() || Q1.front()<Q2.front())) {//NAIVE
                a[i]=Q1.front();
                Q1.pop_front();
            } else {
                a[i]=Q2.front();
                Q2.pop_front();
            }
        }
        sum+=a[0]+a[1];
        Q2.push_back(a[0]+a[1]);
    }

    long long d=std::__gcd(sum,cnt);
    long long p=sum/d,q=cnt/d;
    if(q==1)
        std::cout<<p<<'\n';
    else
        std::cout<<p<<'/'<<q<<'\n';
}

int main() {
    freopen("input","r",stdin);
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);std::cout.tie(0);
    int T;std::cin>>T;
    while(T--)
        Work();
    return 0;
}

[ICPC 2019 Nanchang Invitational Online G]tsy’s number

这几周很忙…好久没有静心做点题了。。
北理校赛的期望DP题做不对
南昌网赛的积性函数题不会做
各种专题都有待加强
这几天先再看看数论吧。

Link

Solution

首先

\begin{aligned}
&\frac{\varphi(i)\varphi(j^2)\varphi(k^3)}{\varphi(i)\varphi(j)\varphi(k)}\\
=&\frac {j^2\prod_{p|j^2}(1-\frac 1p)k^3\prod_{p|k^3}(1-\frac 1p)}
{j\prod_{p|j}(1-\frac 1p)k\prod_{p|k}(1-\frac 1p)}\\
=&jk^2
\end{aligned}

式子化为

\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^njk^2\varphi((i,j,k))\\
=&\sum_{d}\varphi(d)\sum_{i=1}^{[\frac nd]}\sum_{j=1}^{[\frac nd]}\sum_{k=1}^{[\frac nd]}
(jd)(kd)^2[(i,j,k)=1]\\
=&\sum_{d}d^3\varphi(d)\sum_{i=1}^{[\frac nd]}\sum_{j=1}^{[\frac nd]}\sum_{k=1}^{[\frac nd]}
jk^2[(i,j,k)=1]\\
=&\sum_{d}d^3\varphi(d)\sum_{i=1}^{[\frac nd]}\sum_{j=1}^{[\frac nd]}\sum_{k=1}^{[\frac nd]}
jk^2\sum_{t|i,t|j,t|k}\mu(t)\\
=&\sum_{d}d^3\varphi(d)\sum_{t}t^3\mu(t)\sum_{i=1}^{[\frac n{dt}]}\sum_{j=1}^{[\frac n{dt}]}\sum_{k=1}^{[\frac n{dt}]}jk^2\\
=&\sum_{T}T^3\sum_{d|T}\varphi(d)\mu(\frac Td)\sum_{i=1}^{[\frac nT]}\sum_{j=1}^{[\frac nT]}\sum_{k=1}^{[\frac nT]}jk^2\\
=&\sum_{T}T^3\sum_{d|T}\varphi(d)\mu(\frac Td)[\frac nT](\sum_{j=1}^{[\frac nT]}j)(\sum_{k=1}^{[\frac nT]}k^2)
\end{aligned}

考虑筛出T^3(\varphi*\mu)(T)的前缀和。
都是积性函数,随便筛的。

\begin{aligned}
&(p^t)^3(\mu * \varphi)(p^t)\\
=&p^{3t}\sum_{d|p^t}\mu(d)\varphi(\frac {p^t}d)\\
=&p^{3t}(\mu(1)\varphi(p^t)+\mu(p)\varphi(p^{t-1}))\\
=&p^{3t}(\varphi(p^t)-\varphi(p^{t-1}))\\
\end{aligned}

Code

线性筛也写的不熟,写一次忘一次。

#include <bits/stdc++.h>

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

unsigned int Sum(unsigned int x) {
    return (long long)x*(x+1)/2;
}

unsigned int Sum2(unsigned int x) {
    return (unsigned long long)x*(x+1)%(6ull<<30)*(2*x+1)%(6ull<<30)/6;
}

unsigned int Pow3(unsigned int x) {
    return x*x*x;
}

unsigned int g[XN];
void Prep() {
    static bool notPrime[XN];
    static struct { int p,pt; }mp[XN];
    static int prime[XN],pcnt;

    g[1]=1;

    for(int i=2;i<=N;++i) {
        if(!notPrime[i]) {
            prime[++pcnt]=i;
            mp[i]={i,i};
            g[i]=Pow3(i)*(i-2);
        }
        for(int j=1,x;j<=pcnt && prime[j]<=mp[i].p && (x=i*prime[j])<=N;++j) {
            notPrime[x]=1;
            if(mp[i].p==prime[j])
                mp[x]={prime[j],mp[i].pt*prime[j]};
            else
                mp[x]={prime[j],prime[j]};
            auto Phi=[&](int pt)->int {
                return pt==1?1:(mp[x].p-1)*(pt/mp[x].p);
            };
            g[x]=g[x/mp[x].pt]*Pow3(mp[x].pt)*(Phi(mp[x].pt)-Phi(mp[x].pt/mp[x].p));
        }
    }
    for(int i=1;i<=N;++i)
        g[i]+=g[i-1];
}

void Work() {
    int n;std::cin>>n;
    unsigned int Ans=0;
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        Ans+=(n/l)*Sum(n/l)*Sum2(n/l)*(g[r]-g[l-1]);
    }
    std::cout<<(Ans%(1<<30))<<'\n';
}

int main() {
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);
    //freopen("input","r",stdin);
    Prep();
    /*
    for(int i=1;i<=100;++i)
        std::cerr<<Sum(i)-Sum(i-1)<<' '<<Sum2(i)-Sum2(i-1)<<'\n';
    */
    int T;std::cin>>T;
    while(T--)
        Work();
    return 0;
}

[ICPC 2017 ECNA B]Craters

Link

Solution

求出所有可能在凸包上的切线,然后半平面交。。。
二维量的eps比较微妙

Code

const int XN=1e5+11;
const double eps=1e-4,inf=1e100,pi=acos(-1);

int sgn(double const &x) {
    return (x>-eps)-(x<eps);
}

struct Point {
    double x,y;

    double Length() const {
        return sqrt(x*x+y*y);
    }

    double Length2() const {
        return x*x+y*y;
    }

    Point Unit() const {
        double len=sqrt(x*x+y*y);
        return {x/len,y/len};
    }

    Point TurnRight() const {
        return {y,-x};
    }

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

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

    friend Point operator *(Point const &p,double const &k) {
        return {p.x*k,p.y*k};
    }

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

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

    friend double Dist(Point const &a,Point const &b) {
        return sqrt(Inner(a-b,a-b));
    }

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

    friend bool operator ==(Point const &a,Point const &b) {
        return a.x==b.x && a.y==b.y;
    }
}p[XN];

struct Circle {
    Point o;
    double r;

    bool Include(Circle const &a) {
        return sgn(a.r+Dist(o,a.o)-r)<=0;
    }
}o[XN];

struct Line {
    Point p,v;

    friend bool operator <(Line const &a,Line const &b) {
        return a.p<b.p || (a.p==b.p && a.v<b.v);
    }
};

double Dist(Point const &a,Line const &l) {
    return fabs(Outer(a-l.p,l.v))/l.v.Length();
}

bool OnLeft(Point const &p,Line const &ln) {
    return sgn(Outer(ln.v,p-ln.p))>0;
}

bool Paral(Line const &ln1,Line const &ln2) {
    return sgn(Outer(ln1.v,ln2.v))==0;
}

Point Cross(Line const &ln1,Line const &ln2) {
    if(Paral(ln1,ln2)) {
        return ln1.p+ln1.v*inf;
    } else {
        double t=Outer(ln2.v,ln1.p-ln2.p)/Outer(ln1.v,ln2.v);
        return ln1.p+ln1.v*t;
    }
}

Point CutPoint(Point const &p,Circle const &a) {
    if(sgn(Dist(p,a.o)-a.r)==0)
        return p;
    else {
        double x=Dist(p,a.o),cosa=a.r/x,sina=sqrt(1-cosa*cosa);
        Point v=(p-a.o).Unit(),
              t={v.x*cosa-v.y*sina,v.y*cosa+v.x*sina};
        return a.o+t*a.r;
    }
}

Point Mirror(Point const &p,Line const &l) {
    Point foot=l.p+l.v.Unit()*(Inner(p-l.p,l.v)/l.v.Length());
    return foot+(foot-p);
}

Line Cut(Circle const &a,Circle const &b) {
    /*a--b   b|
      ----   a|
      */
    if(sgn(a.r-b.r)==0) {
        return {a.o+(b.o-a.o).Unit().TurnRight()*a.r,b.o-a.o};
    } else {
        double d=Dist(a.o,b.o),
               x=a.r*d/(b.r-a.r);
        Point c=a.o+(a.o-b.o).Unit()*x,
              p=CutPoint(c,a);
        if(a.r>b.r) {
            p=Mirror(p,{a.o,b.o-a.o});
            return {p,c-p};
        } else
            return {c,p-c};
    }
}

std::vector<Line> Intersect(std::vector<Line> &ln) {
    std::sort(ln.begin(),ln.end(),[](Line const &a,Line const &b)->bool {
                auto Quad=[](Point const &a)->int {
                    return sgn(a.x)>0 && sgn(a.y)>=0?1:
                            (sgn(a.x)<=0 && sgn(a.y)>0?2:
                            (sgn(a.x)<0 && sgn(a.y)<=0?3:4));
                };
                int qa=Quad(a.v),qb=Quad(b.v);
                return qa==qb?sgn(Outer(b.v,a.v))<0:qa<qb;
            });
    static Point Qp[XN];static Line Qln[XN];
    int n=ln.size();
    int head,tail;Qln[head=tail=1]=ln[0];
    for(int i=1;i<n;++i) {
        while(tail-head>=1 && !OnLeft(Qp[tail-1],ln[i]))
            tail--;
        while(tail-head>=1 && !OnLeft(Qp[head],ln[i]))
            head++;
        Qln[++tail]=ln[i];
        if(Paral(Qln[tail-1],Qln[tail])){
            int d1=OnLeft(Qln[tail].p,Qln[tail-1]),d2=OnLeft(Qln[tail-1].p,Qln[tail]);
            if(!d1 || !d2) {
                tail--;
                if(d1)
                    Qln[tail]=ln[i];
            }
        }
        if(tail-head>=1)
            Qp[tail-1]=Cross(Qln[tail-1],Qln[tail]);
    }
    while(tail-head>=1 && !OnLeft(Qp[tail-1],Qln[head]))
        tail--;
    if(tail-head>=1)
        return std::vector<Line>(Qln+head,Qln+tail+1);
    else
        return {};
}

double ArcLen(double const &len,double const &r) {
    return r*asin(len/2/r)*2;
}

int main() {
    freopen("input","r",stdin);
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);
    srand(0x1f1f1f1f);
    int n;std::cin>>n;
    for(int i=1;i<=n;++i) {
        std::cin>>o[i].o.x>>o[i].o.y>>o[i].r;
        o[i].r+=10;
    }
    static bool ban[XN];
    for(int i=1;i<=n;++i)
        if(!ban[i]) {
            for(int j=1;j<=n;++j)
                if(i!=j && o[i].Include(o[j]))
                    ban[j]=1;
        }
    int cc=0;
    for(int i=1;i<=n;++i)
        if(!ban[i])
            o[++cc]=o[i];
    n=cc;
    double Ans=0;
    if(n==1)
        Ans=2*pi*o[1].r;
    else {
        std::map<Line,std::pair<int,int>> id;
        std::vector<Line> v;
        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j) 
                if(i!=j) {
                    Line cut=Cut(o[i],o[j]);
                    /*
                    for(auto x : {o[i],o[j]})
                        std::cerr<<Dist(x.o,cut)<<' '<<x.r<<'\n';
                    */
                    bool tak=1;
                    for(int k=1;k<=n;++k)
                        tak&=OnLeft(o[k].o,cut) && sgn(Dist(o[k].o,cut)-o[k].r)>=0;
                    if(tak)
                        v.push_back(cut);
                }
        v=Intersect(v);
        int s=v.size();
        std::vector<std::array<Point,2>> cp(s);
        std::vector<int> ids(s);
        for(int i=0;i<s;++i) {
            double mx=-inf,mn=inf;
            for(int j=1;j<=n;++j) {
                Point pts[2];
                pts[1]=Mirror(pts[0]=CutPoint(v[i].p,o[j]),{v[i].p,o[j].o-v[i].p});
                for(int d=0;d<2;++d)
                    if(sgn(Outer(pts[d]-v[i].p,v[i].v))==0) {
                        double t=Inner(pts[d]-v[i].p,v[i].v)/v[i].v.Length2();
                        if(Enlarge(mx,t)) {
                            cp[(i+1)%s][0]=pts[d];
                            ids[(i+1)%s]=j;
                        }
                        if(Reduce(mn,t)) {
                            cp[i][1]=pts[d];
                        }
                        break;
                    }
            }
        }
        for(int i=0;i<s;++i)
            Ans+=Dist(cp[i][1],cp[(i+1)%s][0])+ArcLen(Dist(cp[i][0],cp[i][1]),o[ids[i]].r);
    }
    std::cout<<std::fixed<<std::setprecision(10)<<Ans<<'\n';
    return 0;
}