未分类 - ShinriiTin's Blog

codechef Dynamic GCD

题意:给出一棵点权树，每次或者询问两点之间路径上所有点点权的gcd，或者把两点之间路径上所有点的点权同时加上一个数

考虑链上的做法，得到了链上做法之后在树上就只需要树链剖分就行了。

由更相减损术得，序列a₁,a₂,..,a_n的gcd等于序列a₁,a₂-a₁,a₃-a₂,...,a_n-a_n-1的gcd。

那么对于链上我们就维护一阶差分的区间gcd，以及原数列的增量就好。

前一个信息我们用线段树维护，后一个信息我们用bit维护。

在树上的推广:用线段树维护每条重链的区间gcd，用bit维护每个点的增量。

#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

#define g() getchar()
template<class Q>void Scan(Q&x){
	char c;int f=1;
	while(c=g(),c<48||c>57)if(c=='-')f=-1;
	for(x=0;c>47&&c<58;c=g())x=10*x+c-48;
	x*=f;
}

#define lch x->l,l,mid
#define rch x->r,mid+1,r
#define rep(i,a,b) for(int i=a;i<b;++i)

inline int gcd(int a,int b){
	for(;b;swap(a,b),b%=a);
	return a>0?a:-a; 
}

const int MAXN=50005;

typedef struct node*star;

struct node{
	int x;
	star l,r,p;
	inline void update(){
		x=gcd(l->x,r->x);
	}
}pool[MAXN<<1];
star tail=pool;

int a[MAXN],b[MAXN];

star pos[MAXN];

void build(star&x,int l,int r){
	x=tail++;
	if(l==r){
		x->x=b[l];
		pos[l]=x;
		return;
	}
	int mid=(l+r)>>1;
	build(lch),build(rch);
	x->l->p=x->r->p=x;
	x->update();
}

int query(star x,int l,int r,int ll,int rr){
	if(ll>rr)return 0;
	if(l==ll&&r==rr)return x->x;
	int mid=(l+r)>>1;
	if(rr<=mid)return query(lch,ll,rr);
	if(ll>mid) return query(rch,ll,rr);
	return gcd(query(lch,ll,mid),query(rch,mid+1,rr));
}

int n;

struct edge{
	int t,n;
	edge(int t=0,int n=0):t(t),n(n){}
}e[MAXN<<1];
int head[MAXN],tot;
inline void AddEdge(int s,int t){
	e[++tot]=edge(t,head[s]),head[s]=tot;
	e[++tot]=edge(s,head[t]),head[t]=tot;
}

int ch[MAXN],p[MAXN];

int dep[MAXN],siz[MAXN];

void dfs(int x){
	siz[x]=1;
	for(int i=head[x],y;i;i=e[i].n)
		if((y=e[i].t)!=p[x]){
			p[y]=x,dep[y]=dep[x]+1;
			dfs(y);
			siz[x]+=siz[y];
			if(siz[y]>siz[ch[x]])ch[x]=y;
		}
}

int dfn[MAXN],cnt; 

int bel[MAXN],top[MAXN],bel_cnt;

void dfs1(int x,int beln){
	dfn[x]=++cnt;
	bel[x]=beln;
	if(ch[x])dfs1(ch[x],beln);
	for(int i=head[x],y;i;i=e[i].n)
		if((y=e[i].t)!=p[x]){
			if(y==ch[x])continue;
			top[++bel_cnt]=y;
			dfs1(y,bel_cnt);
		}
	if(ch[x])b[dfn[ch[x]]]=a[ch[x]]-a[x];
}

star rt;

inline void Init(){
	Scan(n);
	rep(i,1,n){
		int x,y;
		Scan(x),Scan(y);
		AddEdge(++x,++y);
	}
	rep(i,1,n+1)Scan(a[i]);
	dfs(1);
	top[++bel_cnt]=1;
	dfs1(1,bel_cnt);
	build(rt,1,n); 
}

int Q;

int bit[MAXN];

inline void Add(int x,int w){
	for(;x<=n;x+=x&-x)bit[x]+=w;
}
inline void Add(int l,int r,int w){
	if(l>r)return;
	Add(l,w),Add(r+1,-w);
}

inline int Sum(int x){
	int res=0;
	for(;x;x-=x&-x)res+=bit[x];
	return res;
}

inline void Query(int u,int v){
	int ans=0;
	int bu,bv;
	while((bu=top[bel[u]])!=(bv=top[bel[v]])){
		if(dep[bu]>dep[bv]){
			swap(bu,bv),swap(u,v);
		}
		int x=query(rt,1,n,dfn[bv],dfn[v]);
		x=gcd(x,Sum(dfn[bv])+a[bv]);
		ans=gcd(ans,x);
		v=p[bv];
	}
	if(dep[u]>dep[v])swap(u,v);
	ans=gcd(ans,query(rt,1,n,dfn[u]+1,dfn[v]));
	int x=Sum(dfn[u])+a[u];
	printf("%d\n",gcd(ans,x));
}

inline void modify(int x,int c){
	star p=pos[dfn[x]];
	if(!p)return;
	p->x+=c;
	for(p=p->p;p;p=p->p)p->update();
}

inline void Modify(int u,int v,int w){
	int bu,bv;
	while((bu=top[bel[u]])!=(bv=top[bel[v]])){
		if(dep[bu]>dep[bv]){
			swap(bu,bv),swap(u,v);
		}
		Add(dfn[bv],dfn[v],w);
		modify(ch[v],-w);
		v=p[bv];
	}
	if(dep[u]>dep[v])swap(u,v);
	Add(dfn[u],dfn[v],w);
	if(u!=top[bel[u]])modify(u,w);
	modify(ch[v],-w);
}

inline void Solve(){
	Scan(Q);
	rep(i,0,Q){
		char c; 
		while(c=g(),c<'A'||c>'Z');
		int u,v,w;
		Scan(u),Scan(v);
		++u,++v;
		if(c=='F')Query(u,v);
		else Scan(w),Modify(u,v,w);
	}
}

int main(){
	Init(); Solve(); return 0;
}

codeforces round265 div1 D World of Darkraft - 2

详见zhj的15年集训队论文

Portal

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

#define rep(i,a,b) for(int i=a;i<b;++i)
typedef long double ld; 

const int MAXN=1000; 

int n,k,cur;

ld dp[2][MAXN];

int main(){
	scanf("%d%d",&n,&k);
	for(int i=n;i;--i){
		cur^=1;
		int B=701;
		rep(j,1,B){
			dp[cur][j]=(dp[cur^1][j]+(1.+j)/2)*j;
			dp[cur][j]+=(dp[cur^1][j+1]+j);
			dp[cur][j]/=j+1;
			dp[cur][j]+=dp[cur^1][j]*(k-1);
			dp[cur][j]/=k;
		}
	}
	ld ans=dp[cur][1]*k;
	printf("%.20f\n",(double)ans);
	return 0;
}

library

题意:给出一棵点权树，多次询问从一个点u走到另一个点v的路径中，第一个点权不小于w的点是什么

用lct就可以水过去了。。。

~~Portal~~ 我只有数据。。。

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

#define g() getchar()
template<class Q>void Scan(Q&x){
	char c;int f=1;
	while(c=g(),c<48||c>57)if(c=='-')f=-1;
	for(x=0;c>47&&c<58;c=g())x=10*x+c-48;
	x*=f;
}

#define Max(x,y) if(x<y)x=y
#define rep(i,a,b) for(int i=a;i<b;++i)

const int MAXN=2e5+5;

typedef struct node*star;

struct node{
	int x,mx;
	bool rev;
	star ch[2],p;
	inline int d(){
		rep(i,0,2)
			if(p->ch[i]==this)return i;
		return -1;
	}
	inline void set_ch(star x,int d){
		if(~d)ch[d]=x; x->p=this;
	}
	inline void update();
	inline void down(){
		if(!rev)return;
		swap(ch[0],ch[1]);
		ch[0]->rev^=1,ch[1]->rev^=1;
		rev=0;
	}
}pool[MAXN];
star tail=pool;
inline star null(){
	star x=tail++;
	x->ch[0]=x->ch[1]=x->p=x;
	return x;
}
star nu=null();

inline void node::update(){
	mx=x;
	rep(i,0,2)
		if(ch[i]!=nu)Max(mx,ch[i]->mx);
}

inline void rot(star x){
	star y=x->p;
	int d=x->d();
	y->p->set_ch(x,y->d());
	y->set_ch(x->ch[d^1],d),y->update();
	x->set_ch(y,d^1);
}

void pd(star x){
	if(~x->d())pd(x->p); x->down();
}

inline void splay(star x){
	pd(x);
	for(int d1,d2;~(d1=x->d());rot(x))
		if(~(d2=x->p->d()))rot(d1^d2?x:x->p);
	x->update();
}

inline void Access(star x){
	star y=x,t=nu;
	for(;x!=nu;t=x,x=x->p){
		splay(x);
		x->ch[1]=t,x->update();
	}
	splay(y);
}

inline void Evert(star x){
	Access(x),x->rev^=1;
}

inline int Query(star x,star y,int w){
	Evert(x),Access(y);
	if(y->mx<w)return -1;
	while(y!=nu){
		y->down();
		if(y->ch[0]!=nu&&y->ch[0]->mx>=w)
			y=y->ch[0];
		else{
			if(y->x>=w)break;
			y=y->ch[1];	
		}
	}
	return y-pool;
}

int n,q;

struct edge{
	int t,n;
	edge(int t=0,int n=0):t(t),n(n){}
}e[MAXN<<1];
int head[MAXN],tot;
inline void AddEdge(int s,int t){
	e[++tot]=edge(t,head[s]),head[s]=tot;
	e[++tot]=edge(s,head[t]),head[t]=tot;
}

star p[MAXN];

void dfs(int x,int par=0){
	for(int i=head[x],y;i;i=e[i].n)
		if((y=e[i].t)!=par){
			p[y]->p=p[x];
			dfs(y,x);
		}
}

inline void Init(){
	Scan(n),Scan(q);
	rep(i,1,n+1){
		p[i]=tail++;
		*p[i]=*nu;
		Scan(p[i]->x);
		p[i]->mx=p[i]->x;
	}
	rep(i,1,n){
		int x,y;
		Scan(x),Scan(y);
		AddEdge(x,y);
	}
	dfs(1);
}

inline void Solve(){
	rep(i,0,q){
		int u,v,w;
		Scan(u),Scan(v),Scan(w);
		int res=Query(p[u],p[v],w);
		if(~res)printf("%d\n",res);
		else puts("sad");
	}
}

int main(){
	freopen("db.txt","r",stdin);
	freopen("st.out","w",stdout);
	Init(); Solve(); return 0;
}
/*
	4 2
	1 2 1 2
	1 2
	2 3
	3 4
	1 4 2
	3 4 2
*/

hdu4625 JZPTREE

题意:给出一棵树，边权均视为1，对于树上每一个点u，求∑_v dis(u,v)^k

我特意去买了本混凝土数学。。。

《具体数学》:xⁿ=∑_k S₂(n,k)*A(x,k) (S₂表示第二类斯特林数，我不会用Latex别D我。。。)

这个公式的正确性显然，把有x个元素的集合划分为k个集合，允许空集的方案数为xⁿ

再做个简单的变换可以得到:xⁿ=∑_k S₂(n,k)*k!*C(x,k)

则(x+1)ⁿ=∑_k S₂(n,k)*k!*C(x+1,k)

组合数有性质:C(n,k)=C(n-1,k)+C(n-1,k-1)

则(x+1)ⁿ=∑_k S₂(n,k)*k!*(C(x,k)+C(x,k-1))=xⁿ+∑_kS₂(n,k)*k!*C(n,k-1)

实际上，我们容易发现 S₂(n,k)*k! 这个算子一直没有发生改变，所以我们可以预处理出S₂(n,k)*k!的值

令dp[x][k]=∑_{y(y属于以x为根的子树)} C(dis(x,y),k)

则dp[x][k]=∑_{y(y是与x的儿子)} dp[y][k]+dp[y][k-1]

先dfs 一次求出以1为根的所有dp值，再dfs一次，使根为当前节点x并计算 ∑_k S₂(n,k)*k!*dp[x][k]

复杂度:预处理 O(k²)+dp O(nk)

Portal

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

#define g() getchar()
template<class Q>void Scan(Q&x){
	char c;int f=1;
	while(c=g(),c<48||c>57)if(c=='-')f=-1;
	for(x=0;c>47&&c<58;c=g())x=10*x+c-48;
	x*=f;
}

#define rep(i,a,b) for(int i=a;i<b;++i)

const int MAXN=50005;
const int MAXK=505;
const int Mod=10007;

struct edge{
	int t,n;
	edge(int t=0,int n=0):
		t(t),n(n){}
}e[MAXN<<1];
int head[MAXN],tot;
inline void AddEdge(int s,int t){
	e[++tot]=edge(t,head[s]),head[s]=tot;
	e[++tot]=edge(s,head[t]),head[t]=tot;
}

int n,k;

inline void Init(){
	Scan(n),Scan(k);
	rep(i,1,n+1)head[i]=0;
	tot=0;
	rep(i,1,n){
		int u,v;
		Scan(u),Scan(v);
		AddEdge(u,v); 
	}
}

int T;

int s[MAXK][MAXK],fac[MAXK];

inline void init(){
	s[0][0]=fac[0]=1;
	rep(i,1,MAXK){
		fac[i]=fac[i-1]*i%Mod;
		rep(j,1,MAXK){
			s[i][j]=(j*s[i-1][j]+s[i-1][j-1])%Mod;
		}
	}
	rep(i,1,MAXK)
		rep(j,1,MAXK){
			s[i][j]=s[i][j]*fac[j]%Mod;
		}
}

int dp[MAXN][MAXK];

void prev_dfs(int x,int p){
	dp[x][0]=1;
	rep(i,1,k+1)dp[x][i]=0;
	for(int i=head[x],y;i;i=e[i].n)
		if((y=e[i].t)!=p){
			prev_dfs(y,x);
			rep(j,0,k+1){
				dp[x][j]+=dp[y][j];
				if(j)dp[x][j]+=dp[y][j-1];
				dp[x][j]%=Mod;
			}
		}
}

int ans[MAXN];

void dfs(int x,int p){
	ans[x]=0;
	rep(i,0,k+1)
		ans[x]=(ans[x]+s[k][i]*dp[x][i])%Mod;
	for(int i=head[x],y;i;i=e[i].n)
		if((y=e[i].t)!=p){
			rep(j,0,k+1){
				dp[x][j]+=Mod-dp[y][j];
				if(j)dp[x][j]+=Mod-dp[y][j-1];
				dp[x][j]%=Mod;
			}
			rep(j,0,k+1){
				dp[y][j]+=dp[x][j];
				if(j)dp[y][j]+=dp[x][j-1];
				dp[y][j]%=Mod;
			}
			dfs(y,x);
			rep(j,0,k+1){
				dp[y][j]+=Mod-dp[x][j];
				if(j)dp[y][j]+=Mod-dp[x][j-1];
				dp[y][j]%=Mod;
			}
			rep(j,0,k+1){
				dp[x][j]+=dp[y][j];
				if(j)dp[x][j]+=dp[y][j-1];
				dp[x][j]%=Mod;
			}
		}
}

inline void Solve(){
	prev_dfs(1,0);
	dfs(1,0);
	rep(i,1,n+1)
		printf("%d\n",ans[i]);
}

int main(){
	for(init(),Scan(T);T--;Init(),Solve());
	return 0;
}

hdu4969 Just a Joke

题意:

有个点G在半径为R的圆O上做速率为v₁的匀速圆周运动，点B从圆心出发，向点G做速率为v₂的运动，已知O,B,G三点始终共线，求B和G相遇时走过的路程

由三点始终共线得B,G两点的角速度相同，w=v₁/R

把v₂分解为与半径平行的v₂y和与半径垂直的v₂x，设离圆心的距离为r，则v₂x=w*r=v₁*r/R

则v₂y=sqrt(v₂²-(v₁*r/R)²)

由v₂y=dr/dt得:dt=dr/v₂y

两边同时积分得:t(r)=R*arcsin(v₁*r/(v₂*R))/v₁+C

则B点运动时间T为R*arcsin(v₁/v₂)/v₁

知道了运动时间就好算路程了。。。s=v₂*T

Portal

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

int T;

double v1,v2,R,D;

int main(){
	for(scanf("%d",&T);T--;){
		scanf("%lf%lf%lf%lf",&v1,&v2,&R,&D);
		double len=R*v2*asin(v1/v2)/v1;
		if(len-D>1e-9)puts("Why give up treatment");
		else puts("Wake up to code");
	}
	return 0;
}

bzoj3695 滑行

有个东西叫光路最速。。。

所以我们二分出射角度，由折射定律算出入射位置验证。

Portal

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#define eps 1e-6
using namespace std;

const int MAXN=105;
typedef double db;
const db pi=acos(-1.);

int n;

db x,ans,h[MAXN],v[MAXN];

inline db check(db t){
	db x=0,y;ans=0;
	for(int i=1;i<=n;++i){
		y=h[i]*tan(t);
		x+=y;
		ans+=sqrt(h[i]*h[i]+y*y)/v[i];
		if(i<n)t=asin(v[i+1]*sin(t)/v[i]);
	}
	return x;
}

int main(){
	scanf("%d%lf",&n,&x);
	for(int i=n;i;--i)scanf("%lf",h+i);
	for(int i=n;i;--i)scanf("%lf",v+i);
	db l=0,r=pi/2.;
	for(int i=0;i<50;++i){
		db mid=(l+r)/2.;
		if(check(mid)>x)r=mid;
		else l=mid;
	}
	printf("%.3f\n",ans); 
	return 0;
}

codeforces 372D Choosing Subtree is Fun

题意:

选出一棵树的的一个子连通块，要求这个子连通块的大小不超过k

并且定义这个子连通块的价值为最长的一个区间[l，r]，满足区间中的点都在这个子连通块内，的长度。

枚举最大的编号，看最小的编号可以是多少。

具体做法，维护一个区间的点的虚树，每次加入一个新的点时，如果大小超过k，就删除编号最小的点。

用set维护dfs序来维护这个虚树，虚树大小=边长和+1

Portal

#include <set>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

#define g() getchar()
template<class Q>void Scan(Q&x){
	char c;int f=1;
	while(c=g(),c<48||c>57)if(c=='-')f=-1;
	for(x=0;c>47&&c<58;c=g())x=10*x+c-48;
	x*=f;
}

#define fi first
#define se second
#define rep(i,a,b) for(int i=a;i<b;++i)
typedef long long ll;
typedef double db;

const int MAXN=1e5+5;

struct edge{
	int t,n;
	edge(int t=0,int n=0):t(t),n(n){}
}e[MAXN<<1];
int head[MAXN],tot;
inline void AddEdge(int s,int t){
	e[++tot]=edge(t,head[s]),head[s]=tot;
	e[++tot]=edge(s,head[t]),head[t]=tot;
}

int n,k;

int dep[MAXN];

int dfn[MAXN],cnt;

int p[MAXN][20];

void dfs(int x){
	dfn[x]=++cnt;
	for(int i=0;p[x][i];++i)
		p[x][i+1]=p[p[x][i]][i];
	for(int i=head[x],y;i;i=e[i].n)
		if((y=e[i].t)!=*p[x]){
			*p[y]=x;
			dep[y]=dep[x]+1;
			dfs(y);
		}
}

typedef pair<int,int>pr;
set<pr>s;
typedef set<pr>::iterator It;

int ans;

inline int swim(int x,int h){
	for(int i=0;h;++i,h>>=1)
		if(h&1)x=p[x][i];
	return x;
}

inline int lca(int x,int y){
	if(dep[x]>dep[y])swap(x,y);
	y=swim(y,dep[y]-dep[x]); 
	for(int i=19;~i;--i)
		if(p[x][i]!=p[y][i]){
			x=p[x][i];
			y=p[y][i];
		}
	return x==y?x:*p[x];
}

inline int dis(int x,int y){
	return dep[x]+dep[y]-dep[lca(x,y)]*2;
}

inline void insert(int x){
	It it=s.insert(pr(dfn[x],x)).fi;
	if(s.size()==1)return;
	It prev,next;
	It head=s.begin(),tail=--s.end();
	if(it==head)prev=tail;
	else prev=--it,++it;
	if(it==tail)next=head;
	else next=++it,--it;
	ans-=dis(prev->se,next->se);
	ans+=dis(prev->se,x);
	ans+=dis(x,next->se); 
}

inline void remove(int x){
	It it=s.find(pr(dfn[x],x));
	if(s.size()==1){
		s.erase(it);
		return;
	}
	It prev,next;
	It head=s.begin(),tail=--s.end();
	if(it==head)prev=tail;
	else prev=--it,++it;
	if(it==tail)next=head;
	else next=++it,--it;
	ans+=dis(prev->se,next->se);
	ans-=dis(prev->se,x);
	ans-=dis(x,next->se);
	s.erase(it); 
}

inline void Init(){
	Scan(n),Scan(k);
	rep(i,1,n){
		int x,y;
		Scan(x),Scan(y);
		AddEdge(x,y);
	}
	dfs(1);
}

int ans1;

inline void Solve(){
	int l=1;
	rep(r,1,n+1){
		insert(r);
		while((ans>>1)>=k)remove(l++);
		ans1=max(ans1,r-l+1);
	}
	printf("%d\n",ans1);
}

int main(){
	Init(); Solve(); return 0;
}

bzoj3570 Dzy Loves Physics I

看到题面就傻了。。。

完全弹性碰撞:m₁v₁+m₂v₂=m₁v₁'+m₂v₂'

由动能定理得:0.5*m₁*v₁²+0.5*m₂*v₂=0.5*m₁*v₁'²+0.5*m₂*v₂'²

当m₁=m₂时，v₁'=v₂,v₂'=v₁

又有条件a*v=C,即dv/dt *v=C

两边同时积分得: 0.5*v²=C*t

解得v=sqrt(2*C*t)

所以。。。求原来k小的速度，再算出现在的速度即可

Portal

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

const int MAXN=1e5+5;

int bit[MAXN],sz;

inline void Add(int x){
	for(;x<=sz;x+=x&-x)++bit[x];
}

inline int kth(int k){
	int ans=0;
	for(int i=19;~i;--i){
		ans+=1<<i;
		if(ans>=sz||bit[ans]>=k)ans-=1<<i;
		else k-=bit[ans];
	}
	return ans+1;
}

int n,q,C,a[MAXN],h[MAXN<<1];

int opt[MAXN],t[MAXN],x[MAXN];

inline int Hash(int x){
	return lower_bound(h+1,h+1+sz,x)-h;
}

int main(){
	scanf("%d%d",&n,&C);
	for(int i=1;i<=n;++i){
		scanf("%d%*d%*d",a+i);
		h[++sz]=a[i];
	}
	scanf("%d",&q);
	for(int i=1;i<=q;++i){
		scanf("%d",opt+i);
		if(opt[i])scanf("%d%d",t+i,x+i);
		else scanf("%d%*d%*d",x+i),h[++sz]=x[i];
	}
	sort(h+1,h+1+sz);
	sz=unique(h+1,h+1+sz)-(h+1);
	for(int i=1;i<=n;++i)Add(Hash(a[i]));
	for(int i=1;i<=q;++i){
		if(opt[i]){
			int v=h[kth(x[i])];
			printf("%.3f\n",sqrt(2.*C*t[i]+1.*v*v));
		}
		else Add(Hash(x[i]));
	}
	return 0;
}

bzoj2048 书堆

有n本厚度不计的完全相同的书，长度均为m。

问在能放稳的情况下，最上面一本书伸出桌面最长的长度是多少?

显然重心刚好在桌子边缘最优。

令f_n表示n本书时最长伸出桌面的长度，则有fn=∑_i(i>0&&i<=n) m/(2*i)

n较小的时候直接暴力计算

n较大的时候用1+1/2+1/3+...+1/n=ln(n+1)+r,r=0.5772156649

Portal

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

long long n,m;

int main(){
	scanf("%lld%lld",&n,&m);
	double ans=0;
	if(n>10000){
		ans=(log(1.+n)+0.57721566490153286060651209)*m/2.;
	}
	else{
		for(int i=1;i<=n;++i)ans+=1.*m/i/2.;
	}
	printf("%d\n",int(ans-1e-10));
	return 0;
}

bzoj3157 国王奇遇记 && bzoj 3516 国王奇遇记加强版

求∑_i(i>0&&i<=n)i^m*mⁱmod 1000000007

令f_k=∑_i(i>0&&i<=n)i^k*mⁱ

两边同乘(m-1),整理得:(m-1)f_k=n^k*mⁿ⁺¹+∑_i(i>0&&i<=n)mⁱ*[(i-1)^k-i^k]

用二项式定理展开(i-1)^k,整理得:(m-1)f_k=n^k*mⁿ⁺¹+∑_i(i>0&&i<=n)m_i*∑_j(j>=0&&j<k)C(k,j)*(-1)^k-j*i^j

交换一下两个循环的位置,整理得(m-1)f_k=n^k*mⁿ⁺¹+∑_j(j>=0&&j<k)C(k,j)*(-1)^k-j*∑_i(i>0&&i<=n)i^j*mⁱ

好像发现了什么熟悉的东西。。。

(m-1)f_k=n^k*mⁿ⁺¹+∑_j(j>=0&&j<k)C(k,j)*(-1)^k-j*f_j

这样就能O(m²)计算了

这两题做法一模一样。。。我就贴个3157的代码吧

Portal

#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;

const int MAXN=205;
const int Mod=1e9+7;
typedef long long ll;

int n,m;

inline ll qpow(ll x,int k){
	ll ans=1;
	for(;k;x=x*x%Mod,k>>=1)
		if(k&1)ans=ans*x%Mod;
	return ans;
}

int C[MAXN][MAXN];

inline void Init(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=m;++i){
		C[i][0]=C[i][i]=1;
		for(int j=1;j<i;++j){
			C[i][j]=(C[i-1][j]+C[i-1][j-1])%Mod;
		}
	}
}

int ans,s[MAXN];

inline void Solve(){
	int inv=qpow(m-1,Mod-2);
	int nn=1,mm=qpow(m,n+1);
	s[0]=(qpow(m,n+1)-m+Mod)*inv%Mod;
	for(int i=1;i<=m;++i){
		nn=(ll)nn*n%Mod;
		s[i]=(ll)nn*mm%Mod;
		for(int j=0;j<i;++j){
			int res=(ll)C[i][i-j]*s[j]%Mod;
			if((i-j)&1)s[i]=(s[i]+Mod-res)%Mod;
			else s[i]=(s[i]+res)%Mod;
		}
		s[i]=(ll)s[i]*inv%Mod;
	}
	printf("%d\n",s[m]);
}

int main(){
	Init();
	if(m==1){
		ans=((ll)n*(n+1)>>1)%Mod;
		printf("%d\n",ans);
	}
	else Solve();
	return 0;
}

其实还有个再加强版。。。太神，不会做