数字表格
(T)次询问，每次给出(n,m(n,mle 10^6))，(f)为斐波那契数列，(f_0=0,f_1=1)，求：
[prod _{i=1}^nprod _{j=1}^mf[gcd(i,j)]
]
题解做法
这个题的思路很妙啊。
看到那个(gcd)，又注意到前面是个连乘，我们可以构造函数(g(x))，使得：
[f(x)=prod _{d|x}g(d)
]
那么：
[egin{aligned}
ans&=prod _{i=1}^nprod _{j=1}^mprod _{d|i,d|j}g(d) \
&=prod _{d=1}^ng(d)^{lfloor frac{n}{d}
floorlfloorfrac{m}{d}
floor}
end{aligned}
]
分段计算即可。
(g)的求法很简单：
[g(x)=frac{f(x)}{prod {_{d|x,d
e x}g(d)}}
]
预处理复杂度为(O(nsqrt n+nlogn))，询问总复杂度为(O(Tn(sqrt n+sqrt m)))。
想法
其实这个东西如果想不到，还有一个方向可以走：
[egin{aligned}
ans&=prod _{i=1}^nprod _{j=1}^mf[gcd(i,j)] \
&=prod _{d=1}^nf_d^{prod _{i=1}^{lfloorfrac{n}{d}
floor}prod _{j=1}^{lfloor frac{m}{d}
floor}[gcd(i,j)=1]} \
&=prod _{d=1}^nf_d^{prod _{i=1}^{lfloorfrac{n}{d}
floor}prod _{j=1}^{lfloor frac{m}{d}
floor}prod _{k|i,k|j}mu (k)} \
&=prod _{d=1}^nf_d^{prod _{k=1}^{lfloorfrac{n}{d}
floor}mu (k)lfloorfrac{n}{kd}
floorlfloorfrac{m}{kd}
floor} \
令t&=kd \
ans&=prod _{t=1}^nprod _{k|t}f_k^{mu (frac{t}{k})lfloorfrac{n}{t}
floorlfloorfrac{m}{t}
floor}
end{aligned}
]
同样可以分段计算。其实上面的(g)就是下面的$$prod _{k|t}f_k^{mu (frac{t}{k})}$$。
这就是连乘意义下的莫比乌斯反演：
[f(n)=prod _{d|n}g(d) Leftrightarrow g(n)=prod _{d|n}f(d)^{mu (frac{n}{d})}
]
两边取对数即得证。
代码
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long giant;
int read() {
int x=0,f=1;
char c=getchar();
for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
for (;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
const int maxn=1e6+1;
const int q=1e9+7;
int f[maxn],g[maxn],iv[maxn],ig[maxn];
int multi(int x,int y) {
return (giant)x*y%q;
}
int Plus(int x,int y) {
return ((giant)x+y)%q;
}
int mi(int x,int y) {
int ret=1;
for (;y;y>>=1,x=multi(x,x)) if (y&1) ret=multi(ret,x);
return ret;
}
int inv(int x) {
return mi(x,q-2);
}
int pre[maxn],ipre[maxn];
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
#endif
iv[1]=1;
for (int i=2;i<maxn;++i) iv[i]=q-multi(iv[q%i],q/i);
f[0]=0,f[1]=g[1]=ig[1]=1;
for (int i=2;i<maxn;++i) {
f[i]=Plus(f[i-1],f[i-2]);
int &in=ig[i]=1;
for (int j=1;j*j<=i;++j) if (i%j==0) {
in=multi(in,ig[j]);
if (j*j!=i) in=multi(in,ig[i/j]);
}
g[i]=multi(f[i],in);
in=inv(g[i]);
}
pre[0]=ipre[0]=1;
for (int i=1;i<maxn;++i) ipre[i]=inv(pre[i]=multi(pre[i-1],g[i]));
int T=read();
while (T--) {
int ret=1,n=read(),m=read();
if (n>m) swap(n,m);
for (int l=1,r,c;l<=n;l=r+1) {
r=min(n/(n/l),m/(m/l));
int up=(giant)(n/l)%(q-1)*(m/l)%(q-1);
ret=multi(ret,mi(multi(pre[r],ipre[l-1]),up));
}
printf("%d
",ret);
}
return 0;
}
树点涂色
给出一颗(n)个点1为根的树，初始时每个点都有个不一样的颜色。我们定义一条路径的权值为路径上不同颜色的个数。有三个操作：
- 将一个点到根的的路径上的所有点都染成一个新的颜色
- 询问一条路径((x,y))的权值
- 询问一个点的子树中所有点，他们到根路径的权值最大的那个
(n,mle 10^5)。
分析
这个问题也很巧妙！我们注意到每次操作都修改一个点到根路径上的所有点，每次的颜色都是新的，这就可以发现和link-cut tree的access操作很像。并且有一个性质，由于每次修改都是修改到根的路径，并且都变成新的颜色，所以颜色一定是竖直链状的。这样我们可以想到，如果对每个点维护一个(f)值，如果它的颜色与父亲一样，那么这个值为1，否则为0，那么一个点到根路径的权值就是根到这个点路径上的所有点的(f)值和，我们称它为(g)。
我们可以用线段树来维护这个东西，由于询问中有子树询问，所以就使用dfs序的线段树。每次修改用一个access操作，每断掉一条虚边，我们就把这个虚子树的所有点的(g)值加1，表示下面的点到根路径上多了一种颜色；每连上一条实边，我们就把这个子树的所有点(g)值减1。这样层层往上，我们就可以维护所有点的(g)信息了。
询问路径((x,y))的权值，设(x,y)的( ext{lca})为(l)，那么答案就是(g_x+g_y-2g_l+1)，想象一下这个情况就知道它减掉两倍(g_l)会少掉一种颜色，所以要补回来。
代码
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
int read() {
int x=0,f=1;
char c=getchar();
for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
for (;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
const int maxn=1e5+10;
const int maxj=17;
struct edge {
int v,nxt;
} e[maxn<<2];
int h[maxn],tot=0;
int f[maxn][maxj],dep[maxn],first[maxn],second[maxn],dft=0,n;
int inv[maxn];
void add(int u,int v) {
e[++tot]=(edge){v,h[u]};
h[u]=tot;
}
int lca(int x,int y) {
if (dep[x]<dep[y]) swap(x,y);
for (int j=maxj-1;j>=0;--j) if (dep[f[x][j]]>=dep[y]) x=f[x][j];
if (x==y) return x;
for (int j=maxj-1;j>=0;--j) if (f[x][j]!=f[y][j]) x=f[x][j],y=f[y][j];
return f[x][0];
}
struct SGT {
int t[maxn<<2],tag[maxn<<2];
void update(int x) {
t[x]=max(t[x<<1],t[x<<1|1]);
}
void build(int x,int l,int r) {
if (l==r) {
t[x]=dep[inv[l]];
tag[x]=0;
return;
}
int mid=l+r>>1;
build(x<<1,l,mid);
build(x<<1|1,mid+1,r);
update(x);
}
void doit(int x,int d) {
t[x]+=d;
tag[x]+=d;
}
void pushdown(int x) {
doit(x<<1,tag[x]);
doit(x<<1|1,tag[x]);
update(x);
tag[x]=0;
}
void add(int x,int L,int R,int l,int r,int d) {
if (L==l && R==r) {
doit(x,d);
return;
}
int mid=L+R>>1;
pushdown(x);
if (r<=mid) add(x<<1,L,mid,l,r,d); else
if (l>mid) add(x<<1|1,mid+1,R,l,r,d); else
add(x<<1,L,mid,l,mid,d),add(x<<1|1,mid+1,R,mid+1,r,d);
update(x);
}
int query(int x,int L,int R,int l,int r) {
if (L==l && R==r) return t[x];
pushdown(x);
int mid=L+R>>1;
if (r<=mid) return query(x<<1,L,mid,l,r);
if (l>mid) return query(x<<1|1,mid+1,R,l,r);
return max(query(x<<1,L,mid,l,mid),query(x<<1|1,mid+1,R,mid+1,r));
}
int one(int p) {
p=first[p];
return query(1,1,n,p,p);
}
} sgt;
struct node {
int ch[2],fa;
};
struct LCT {
node t[maxn];
void link(int x,int y) {
t[y].fa=x;
}
int left(int x) {
while (t[x].ch[0]) x=t[x].ch[0];
return x;
}
bool isroot(int x) {
return !x || t[t[x].fa].ch[rson(x)]!=x;
}
bool rson(int x) {
return t[t[x].fa].ch[1]==x;
}
void rotate(int x) {
int f=t[x].fa,d=rson(x),c=t[x].ch[d^1];
if (!isroot(f)) t[t[f].fa].ch[rson(f)]=x;
int inv[maxn];
if (c) t[c].fa=f;
t[x].fa=t[f].fa,t[f].fa=x,t[x].ch[d^1]=f,t[f].ch[d]=c;
}
void splay(int x) {
while (!isroot(x)) {
if (isroot(t[x].fa)) rotate(x); else {
if (rson(x)==rson(t[x].fa)) rotate(t[x].fa),rotate(x); else
rotate(x),rotate(x);
}
}
}
void access(int x) {
for (int last=0;x;x=t[last=x].fa) {
splay(x);
int p=t[x].ch[1]?left(t[x].ch[1]):0;
if (p) sgt.add(1,1,n,first[p],second[p],1);
t[x].ch[1]=last;
p=t[x].ch[1]?left(t[x].ch[1]):0;
if (p) sgt.add(1,1,n,first[p],second[p],-1);
}
}
} lct;
void dfs(int x,int fa) {
inv[first[x]=++dft]=x;
f[x][0]=fa;
if (x==fa) dep[x]=1; else dep[x]=dep[fa]+1,lct.link(fa,x);
for (int i=h[x],v=e[i].v;i;i=e[i].nxt,v=e[i].v) if (v!=fa) dfs(v,x);
second[x]=dft;
}
int line(int x,int y) {
int l=lca(x,y);
int fx=sgt.one(x),fy=sgt.one(y),fl=sgt.one(l);
return fx+fy-fl-fl+1;
}
int sub(int x) {
int ret=sgt.query(1,1,n,first[x],second[x]);
return ret;
}
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
freopen("my.out","w",stdout);
#endif
n=read();
int m=read();
for (int i=2;i<=n;++i) {
int x=read(),y=read();
add(x,y),add(y,x);
}
dfs(1,1);
for (int j=1;j<maxj;++j) for (int i=1;i<=n;++i) f[i][j]=f[f[i][j-1]][j-1];
sgt.build(1,1,n);
while (m--) {
int op=read();
if (op==1) {
int x=read();
lct.access(x);
} else if (op==2) {
int x=read(),y=read();
int ans=line(x,y);
printf("%d
",ans);
} else if (op==3) {
int x=read();
int ans=sub(x);
printf("%d
",ans);
}
}
return 0;
}
序列计数
有多少长度为(n)，其中的数都小于等于(m)，和模(p)为0，且至少含有一个质数的数列？
(nle 10^9,mle 2 imes 10^7,ple 100)。
分析
至少含有一个质数的数列，就是随便乱选的数列个数减去不含质数的数列个数。对于一些可选的数，由于选的次数不限，所以每次我们都可以从当前的选法加上一个数，即(f_{ij})表示选了(i)个数，当前模(p)为(j)的方案数，每次通过一个固定的线性转移转移过去。所以可以用矩阵乘法优化解决。我写的是生成函数的方法，本质和矩阵是一样的。
代码
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long giant;
const int q=20170408;
int multi(int x,int y) {
return ((giant)x*y)%q;
}
int Plus(int x,int y) {
return ((giant)x+y)%q;
}
int sub(int x,int y) {
return Plus(x,q-y);
}
const int maxn=250;
int a[maxn],b[maxn],n,m,p,c[maxn];
void mul(int a[],int b[]) {
memset(c,0,sizeof c);
for (int i=0;i<p;++i) for (int j=0;j<p;++j) c[(i+j)%p]=Plus(c[(i+j)%p],multi(a[i],b[j]));
for (int i=0;i<p;++i) a[i]=c[i];
}
int first() {
memset(b,0,sizeof b),b[0]=1;
memset(a,0,sizeof a);
for (int i=1;i<=m;++i) ++a[i%p];
int y=n;
for (;y;y>>=1,mul(a,a)) if (y&1)
mul(b,a);
return b[0];
}
const int maxm=2e7+1;
const int pss=6e6+1;
bool np[maxm];
int pm[pss],ps=0;
int second() {
np[1]=true;
for (int i=2;i<=m;++i) {
if (!np[i]) pm[++ps]=i;
for (int j=1;j<=ps && (giant)pm[j]*i<(giant)maxm;++j) {
int tmp=pm[j]*i;
np[tmp]=true;
if (i%pm[j]==0) break;
}
}
memset(b,0,sizeof b),b[0]=1;
memset(a,0,sizeof a);
for (int i=1;i<=m;++i) if (np[i]) ++a[i%p];
int y=n;
for (;y;y>>=1,mul(a,a)) if (y&1) mul(b,a);
return b[0];
}
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
#endif
scanf("%d%d%d",&n,&m,&p);
int f=first();
int s=second();
int ans=sub(f,s);
printf("%d
",ans);
}
新生舞会
一个二分图匹配，每一个匹配((i,j))有两个属性(a_{ij},b_{ij})，求一个完备匹配，使得：
[C=frac{sum a_i}{sum b_i}
]
最大。(nle 100)。
分析
二分答案(C)，我们要求的是是否有一种匹配使得(frac{sum a_i}{sum b_i})大于(C)，即是否存在一个(sum a_i-Cb_i)大于0的匹配。直接做最大匹配即可，这里用的是KM算法。
代码
#include<cstdio>
#include<cctype>
#include<cstring>
#include<cmath>
#include<queue>
#include<algorithm>
using namespace std;
int read() {
int x=0,f=1;
char c=getchar();
for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
for (;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
const double eps=1e-10;
const double up=1e8;
const double inf=1e100;
const int maxn=205;
int a[maxn][maxn],b[maxn][maxn],n,match[maxn],pre[maxn];
double f[maxn][maxn],ex[maxn],ey[maxn],slack[maxn];
bool vy[maxn];
void augment(int now) {
int y;
fill(slack,slack+n+1,inf);
memset(vy,0,sizeof vy);
match[y=0]=now;
do {
vy[y]=true;
int to,x=match[y];
double d=inf;
for (int i=1;i<=n;++i) if (!vy[i]) {
if (ex[x]+ey[i]-f[x][i]<slack[i]) slack[i]=ex[x]+ey[i]-f[x][i],pre[i]=y;
if (slack[i]<d) d=slack[i],to=i;
}
for (int i=0;i<=n;++i) if (vy[i]) ex[match[i]]-=d,ey[i]+=d; else slack[i]-=d;
y=to;
} while (match[y]);
for (;y;y=pre[y]) match[y]=match[pre[y]];
}
double km() {
for (int i=1;i<=n;++i)
augment(i);
double ret=0;
for (int i=1;i<=n;++i) ret+=ex[i]+ey[i];
return ret;
}
bool ok(double c) {
memset(f,0,sizeof f),memset(ey,0,sizeof ey),memset(match,0,sizeof match),memset(pre,0,sizeof pre);
for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) f[i][j]=(double)a[i][j]-c*(double)b[i][j];
for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) ex[i]=(j==1?f[i][j]:max(ex[i],f[i][j]));
double cost=km();
return cost>=0;
}
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
freopen("my.out","w",stdout);
#endif
n=read();
for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) a[i][j]=read();
for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) b[i][j]=read();
double l=0,r=1e4,mid,ans;
while (fabs(l-r)>eps) {
mid=(l+r)/2;
if (ok(mid)) ans=mid,l=mid; else r=mid;
}
printf("%.6lf
",ans);
return 0;
}
相关分析
(n,mle 10^5)。
分析
这是一个区间修改，区间加法，区间查询奇怪东西的问题，可以看看它要维护的是什么东西，从查询入手。
[egin{aligned}
a&=frac{sum _{i=L}^R(x_i-overline x)(y_i-overline y)}{sum _{i=L}^R(x_i-overline x)^2} \
&=frac{sum _{i=L}^Rx_iy_i-overline xsum _{i=L}^Ry_i-overline ysum _{i=L}^Rx_i+(R-L+1)overline xoverline y}{sum _{i=L}^R x_i^2+2overline xsum _{i=L}^Rx_i+(R-L+1)overline x^2}
end{aligned}
]
所以维护区间(x)，(y)，(xy)，(x^2)的和就好啦，记得修改的优先级比区间加要高。
代码
#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long giant; // remember to change long long
typedef long double lb;
giant read() {
giant x=0,f=1;
char c=getchar();
for (;!isdigit(c);c=getchar()) if (c=='-') f=-1;
for (;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*f;
}
const giant maxn=1e5+10;
lb basex[maxn],basey[maxn];
int n;
giant pf[maxn];
struct Data {
/*
* x is the sum of xi in this range
* y is the sum of yi in this range
* xy is the sum of xi*yi in the range
* fang is the sum of xi*xi in this range
*/
lb x,y,xy,fang;
Data (lb x=0,lb y=0,lb xy=0,lb fang=0):x(x),y(y),xy(xy),fang(fang) {}
void init(lb _x,lb _y) {
x=_x,y=_y;
xy=x*y;
fang=x*x;
}
};
Data operator + (Data a,Data b) {
return Data(a.x+b.x,a.y+b.y,a.xy+b.xy,a.fang+b.fang);
}
Data operator += (Data &a,Data b) {
a=a+b;
}
struct Tag {
lb s,t;
int op;
Tag (lb s=0,lb t=0,int op=0):s(s),t(t),op(op) {}
void clear() {
s=t=op=0;
}
};
Tag operator + (Tag a,Tag b) {
return Tag(a.s+b.s,a.t+b.t,max(a.op,b.op));
}
struct SGT {
Data dat[maxn<<3];
Tag tag[maxn<<3];
void build(int x,int l,int r) {
if (l==r) {
dat[x].init(basex[l],basey[l]);
tag[x].clear();
return;
}
giant mid=l+r>>1;
build(x<<1,l,mid);
build(x<<1|1,mid+1,r);
dat[x]=dat[x<<1]+dat[x<<1|1];
}
void doit(giant x,giant l,giant r,lb s,lb t,int op) {
if (op==1) {
Data &d=dat[x];
d.xy+=t*d.x+s*d.y+(r-l+1)*s*t;
d.fang+=2*s*d.x+(r-l+1)*s*s;
d.x+=(r-l+1)*s;
d.y+=(r-l+1)*t;
tag[x]=tag[x]+Tag(s,t,op);
} else if (op==2) {
Data &d=dat[x];
giant qj=(l+r)*(r-l+1)/2,f=pf[r]-pf[l-1];
d.xy=f+qj*(s+t)+(r-l+1)*s*t;
d.x=qj+(r-l+1)*s;
d.y=qj+(r-l+1)*t;
d.fang=f+s*qj*2+(r-l+1)*s*s;
tag[x]=Tag(s,t,op);
}
}
void pushdown(giant x,giant l,giant mid,giant r) {
if (tag[x].op==0) return;
if (tag[x].op==2) {
doit(x<<1,l,mid,tag[x].s,tag[x].t,tag[x].op);
doit(x<<1|1,mid+1,r,tag[x].s,tag[x].t,tag[x].op);
} else if (tag[x].op==1) {
doit(x<<1,l,mid,tag[x].s,tag[x].t,tag[x].op);
doit(x<<1|1,mid+1,r,tag[x].s,tag[x].t,tag[x].op);
}
tag[x].clear();
}
void modify(giant x,giant L,giant R,giant l,giant r,lb s,lb t,int op) {
giant mid=L+R>>1;
pushdown(x,L,mid,R);
if (L==l && R==r) {
doit(x,L,R,s,t,op);
return;
}
if (r<=mid) modify(x<<1,L,mid,l,r,s,t,op); else
if (l>mid) modify(x<<1|1,mid+1,R,l,r,s,t,op); else {
modify(x<<1,L,mid,l,mid,s,t,op);
modify(x<<1|1,mid+1,R,mid+1,r,s,t,op);
}
dat[x]=dat[x<<1]+dat[x<<1|1];
}
void modify(giant l,giant r,lb s,lb t,int op) {
modify(1,1,n,l,r,s,t,op);
}
Data query(giant x,giant L,giant R,giant l,giant r) {
if (L==l && R==r) return dat[x];
giant mid=L+R>>1;
pushdown(x,L,mid,R);
if (r<=mid) return query(x<<1,L,mid,l,r);
if (l>mid) return query(x<<1|1,mid+1,R,l,r);
Data r1=query(x<<1,L,mid,l,mid);
Data r2=query(x<<1|1,mid+1,R,mid+1,r);
Data ret=r1+r2;
return ret;
}
Data query(giant l,giant r) {
return query(1,1,n,l,r);
}
} sgt;
lb calc(int thel,int ther) {
giant len=ther-thel+1;
Data d=sgt.query(thel,ther);
lb up=d.xy,xp=(lb)d.x/len,yp=(lb)d.y/len;
up-=xp*d.y+yp*d.x;
up+=len*xp*yp;
lb down=d.fang;
down-=2*xp*d.x;
down+=len*xp*xp;
lb ret=up/down;
return ret;
}
void radd(int thel,int ther,lb s,lb t) {
sgt.modify(thel,ther,s,t,1);
}
void rmodi(int thel,int ther,lb s,lb t) {
sgt.modify(thel,ther,s,t,2);
}
int main() {
#ifndef ONLINE_JUDGE
freopen("test.in","r",stdin);
freopen("my.out","w",stdout);
#endif
n=read();
for (int i=1;i<=n;++i) pf[i]=pf[i-1]+(giant)i*i;
int m=read();
for (int i=1;i<=n;++i) scanf("%Lf",basex+i);
for (int i=1;i<=n;++i) scanf("%Lf",basey+i);
sgt.build(1,1,n);
while (m--) {
int op=read();
if (op==1) {
int l=read(),r=read();
lb ans=calc(l,r);
printf("%.10Lf
",ans); // change back to .10lf
} else if (op==2) {
int l=read(),r=read();
lb s,t;
scanf("%Lf%Lf",&s,&t);
radd(l,r,s,t);
} else if (op==3) {
int l=read(),r=read();
lb s,t;
scanf("%Lf%Lf",&s,&t);
rmodi(l,r,s,t);
}
}
return 0;
}