SDOI2013 R1 Day2

目录

• 2018.3.25 Test
  • 总结
  • T1 BZOJ.3129.[SDOI2013]方程(扩展Lucas 容斥)
  • T2 洛谷.3305.[SDOI2013]费用流(最大流ISAP 二分)
  • T3 BZOJ.3131.[SDOI2013]淘金(数位DP 堆)
  • 考试代码
    • T1
    • T2
    • T3

---

2018.3.25 Test

时间：7:30~11:30 (最后半小时不做了)
期望得分：50+100+20=170
实际得分：40+44+20=104

总结

T1：善用容斥。
T2：要求输出小数当然有它的道理。。
T3：行列分开考虑。看好怎么取模。如要求最大的答案，然后对答案取模，那么在计算过程中不能取模，最后求答案时才可以！

T1 BZOJ.3129.[SDOI2013]方程(扩展Lucas 容斥)

题目链接

如果n个未知数是正整数且没有限制，那么方案数是C(m-1,n-1).(把m个数分割成n份)
不需要把x变为非负整数解，然后去求C(n-k+1,k).(当然无所谓吧。。按标程来)
(geq) 的限制很好做，减去即可。那么对于 (xileq Ai) 的我们容斥，去掉x的限制，然后把 (x>Ai) 的减掉。
利用扩展Lucas计算组合数.复杂度 (O(2^{n1}*Ex\_Lucas)).
复杂度有点高了。。扩展Lucas预处理阶乘！(预处理后飞快啊，以前从来没预处理。。)

//1620kb	1364ms
//预处理前:840kb	10552ms
#include <cmath>
#include <cstdio>
#include <cctype>
#include <algorithm>
#define gc() getchar()
typedef long long LL;
const int N=1e4+23;

int n,P,n1,n2,m,A[10],cnt,pi[666],pk[666];
LL Ans,fac[100005]/*这个开多大还是不清楚*/;

inline int read()
{
	int now=0;register char c=gc();
	for(;!isdigit(c);c=gc());
	for(;isdigit(c);now=now*10+c-'0',c=gc());
	return now;
}
LL FP(LL x,int k,LL p)
{
	LL t=1ll;
	for(; k; k>>=1,x=x*x%p)
		if(k&1) t=t*x%p;
	return t;
}
void Exgcd(int a,int b,int &x,int &y)
{
	if(!b) x=1,y=0;
	else Exgcd(b,a%b,y,x),y-=a/b*x;
}
inline LL inv(int a,int p)
{
	int x,y; Exgcd(a,p,x,y);
	return (LL)((x%p+p)%p);
}
void Init_P(int p)
{
	for(int i=2,lim=sqrt(p+1); i<=lim; ++i)
		if(!(p%i))
		{
			int _pk=1;
			while(!(p%i)) p/=i,_pk*=i;
			pi[++cnt]=i, pk[cnt]=_pk;
			if(p==1) break;
		}
	if(p!=1) pi[++cnt]=pk[cnt]=p;
}
LL Fact(int n,int pi,int pk)
{
	if(!n) return 1ll;//边界！ 
	LL res=1ll;
	if(n>=pk) res=FP(fac[pk],n/pk,pk);
	(res*=fac[n%pk])%=pk;
	return res*Fact(n/pi,pi,pk)%pk;
}
LL C(int n,int m,int pi,int pk,int mod)
{
	fac[0]=1;
	for(int i=1; i<=pk; ++i)
		if(i%pi) fac[i]=fac[i-1]*i%pk;
		else fac[i]=fac[i-1];
	LL a=Fact(n,pi,pk),b=Fact(m,pi,pk),c=Fact(n-m,pi,pk);int k=0;
	for(int i=n; i; i/=pi) k+=i/pi;
	for(int i=m; i; i/=pi) k-=i/pi;
	for(int i=n-m; i; i/=pi) k-=i/pi;
	LL res=a*inv(b,pk)%pk*inv(c,pk)%pk*FP(pi,k,pk)%pk;
	return 1ll*res*(mod/pk)%mod*inv(mod/pk,pk)%mod;
}
LL Ex_Lucas(int n,int m)
{
	if(n<m) return 0ll;
	LL res=0;
	for(int i=1; i<=cnt; ++i)
		(res+=C(n,m,pi[i],pk[i],P))%=P;
	return res;
}
void DFS(int p,int cnt,int sum)
{
	if(p>n1)
		if(cnt&1) (Ans-=Ex_Lucas(m-sum-1,n-1))%=P;//前边的数减去A[i]，即满足>A[i]。
		else (Ans+=Ex_Lucas(m-sum-1,n-1))%=P;//C(m-sum-1,n-1) not C(n-sum-1,m-1)!
	else DFS(p+1,cnt,sum),DFS(p+1,cnt+1,sum+A[p]);
}

int main()
{
	int T=read(); P=read();
	Init_P(P);
	while(T--)
	{
		n=read(),n1=read(),n2=read(),m=read();
		LL tot=0;
		for(int i=1; i<=n1; ++i) A[i]=read();
		for(int i=1; i<=n2; ++i) m-=read()-1;//去下界 
		if(m<=1) {puts("0"); continue;}
		if(m==n) {puts("1"); continue;}
		Ans=0ll, DFS(1,0,0), printf("%lld
",(Ans+P)%P);//取模！
	}
	return 0;
}

T2 洛谷.3305.[SDOI2013]费用流(最大流ISAP 二分)

题目链接(BZOJ无SPJ)

题意即是求满足最大流条件下 使得流量最大的边尽量小。花费显然是最大的某条边上的流量*P。
但P是实数范围内！这说明流量也是实数，二分答案跑最大流检验。限制的话，直接重设流量。
贪心的话最大流过程中流量一直是整数，所以就GG了。。

#include <cmath>
#include <cstdio>
#include <cctype>
#include <algorithm>
#define gc() getchar()
const double eps=1e-6;
const int N=218,M=2018;

int n,m,P,Max,src,des,Enum,H[N],cur[N],fr[M],to[M],nxt[M],cost[M],num[N],lev[N],pre[N],que[N];
double cap[M],Cap[M];

inline int read()
{
	int now=0;register char c=gc();
	for(;!isdigit(c);c=gc());
	for(;isdigit(c);now=now*10+c-'0',c=gc());
	return now;
}
inline void AddEdge(int u,int v,int w)
{
	Max=std::max(Max,w);
	to[++Enum]=v, fr[Enum]=u, nxt[Enum]=H[u], H[u]=Enum, Cap[Enum]=cap[Enum]=w;
	to[++Enum]=u, fr[Enum]=v, nxt[Enum]=H[v], H[v]=Enum, Cap[Enum]=cap[Enum]=0;	
}
bool BFS()
{
	for(int i=src; i<des; ++i) lev[i]=des+1;
	lev[des]=0, que[0]=des; int h=0,t=1;
	while(h<t)
	{
		int x=que[h++];
		for(int i=H[x]; i; i=nxt[i])
			if(lev[to[i]]==des+1 && cap[i^1])
				lev[to[i]]=lev[x]+1, que[t++]=to[i];
	}
	return lev[src]<=des;
}
double Augment()
{
	double mn=1e8;
	for(int i=des; i!=src; i=fr[pre[i]])
		mn=std::min(mn,cap[pre[i]]);
	for(int i=des; i!=src; i=fr[pre[i]])
		cap[pre[i]]-=mn, cap[pre[i]^1]+=mn;
	return /*fabs(mn)<eps?0:*/mn;
}
double ISAP()
{
	if(!BFS()) return 0;
	for(int i=src; i<=des; ++i) ++num[lev[i]],cur[i]=H[i];
	int x=src; double res=0;
	while(lev[src]<=des)
	{
		if(x==des) x=src,res+=Augment();
		bool can=0;
		for(int i=cur[x]; i; i=nxt[i])
			if(lev[to[i]]==lev[x]-1 && cap[i])
			{
				can=1, cur[x]=i, pre[x=to[i]]=i;
				break;
			}
		if(!can)
		{
			int mn=des;
			for(int i=H[x]; i; i=nxt[i])
				if(cap[i]) mn=std::min(mn,lev[to[i]]);
			if(!--num[lev[x]]) break;
			++num[lev[x]=mn+1];
			cur[x]=H[x];
			if(x!=src) x=fr[pre[x]];
		}
	}
	return res;
}
void Init(double lim){
	for(int i=2; i<=Enum; ++i) cap[i]=std::min(Cap[i],lim);
}
double Solve(double mxflow)
{
	double l=1.0,r=Max,m;
	while(r-l>1e-4)
		if(Init(m=(l+r)*0.5),fabs(mxflow-ISAP())<eps) r=m;
		else l=m;
	return l;
}

int main()
{
	n=read(),m=read(),P=read(),Enum=src=1,des=n;
	for(int u,v,w,i=1; i<=m; ++i)
		u=read(),v=read(),w=read(),AddEdge(u,v,w);//f&=(w==1);
	double res1=ISAP();
	printf("%.0lf
%.4lf",res1,1.0*Solve(res1)*P);

	return 0;
}

T3 BZOJ.3131.[SDOI2013]淘金(数位DP 堆)

题目链接

40pts:
行列是独立的！即同在第i行的元素一定都会加到f(i)行，它们对所有列位置的贡献同其它行一样是固定的，只是加的行不同而已。列同理。而且行列的贡献是相同的。
所以求出(一维下)有多少个数能转移到某个位置i。找出两两相乘(二维)最大的K个。-> 这40pts我也不会写啊QAQ
其实是用堆 然后每个节点维护两个指针，拿一个就放两个？并不，这样在x==y时会有重。(其实也有办法 大概要判重之类？)
如果x是初始的s_x，则不断扩展 push(x,y+1)；如果x!=y，则push(x+1,y)(维护一个有序数对，感觉好6啊。。)
而且在x!=y时加两遍，因为行列是可以在反过来的！
(也可以先push进n个去再扩展吧。。)

100pts:
(n)很大无法直接算出某个位置的(val[])。
在(n)的范围内(f(x))的种类只有不是很多种(大概不到(20000)?)，且因子只会有(2.3.5.7)，可以通过数位DP算出(y=2^i*3^j*5^k*7^l)时有多少个x的数位乘积=(y(xleq n)).
(dp[len][0/1][i][j][k][l])表示当前是(n)的第(len)位，前面是否达到了上界，以及(y)的因子个数分别为(i,j,k,l)时的(x)个数。
最后统计所有(y)的(val[])。
注意dp的细节: 可以有前导(0)，除了最高位都要计算前导(0)的贡献，即(dp[now][0][0][0][0][0])++(不卡上界)；最高位初始化时是(dp[now][1][0][0][0][0]=1)(上界)。
（并没有自己写代码。。借鉴(抄)下xxy博客了）
另外这个取模是在取最多的情况下算出后取模！(val[])运算过程中不能取模。(是这样吧)

//12124kb	740ms
#include <queue>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
#define gc() getchar()
#define mod (1000000007)
typedef long long LL;
const int N=1e5+7;

int K,cnt,f[11][5]/*多一位数影响的数值*/,A[20];
//LL n,val[N],dp[2][2][40][26][18][15];
LL n,val[N],dp[2][2][41][27][19][16];
struct Node
{
	int x,y; LL v;
	Node(int _x,int _y): x(_x),y(_y),v(val[x]*val[y]) {}
	bool operator <(const Node &a)const{
		return v<a.v;//val升序 
	}
};
std::priority_queue<Node> q;

inline bool cmp(const int &a,const int &b){
	return a>b;
}
inline void Pre()
{//0:2 1:3 2:5 3:7
	f[2][0]=1, f[4][0]=2, f[6][0]=1, f[8][0]=3,
	f[3][1]=1, f[6][1]=1, f[9][1]=2,
	f[5][2]=1, f[7][3]=1;
}
inline LL FP(LL x,int k)
{
	LL t=1;
	for(; k; k>>=1,x=x*x)
		if(k&1) t=t*x;
	return t;
}
void DP()
{//2:40 3:26 5:18 7:15
	Pre();
	LL t=n; int len=0;
	while(t) A[len++]=t%10,t/=10;//我一定是sb了。。下标从0开始还用A[0]做长度 
	std::reverse(A,A+len);
	int now=0,nxt=1;
	for(int p=0; p<len; ++p,std::swap(now,nxt))//神tm还拿l做变量。。
	{
		
		if(!p) dp[now][1][0][0][0][0]=1;
		else ++dp[now][0][0][0][0][0];
		memset(dp[nxt],0,sizeof dp[nxt]);
		for(int s=0; s<=1; ++s)
			for(int i=0; i<=40; ++i)
				for(int j=0; j<=26; ++j)
					for(int k=0; k<=18; ++k)
						for(int l=0; l<=15; ++l)
							if(dp[now][s][i][j][k][l])//怎么这么。。
								for(int w=1,lim=s?A[p]:9; w<=lim; ++w)//不能取模  //另外不计算0！
									dp[nxt][s&&w==lim][i+f[w][0]][j+f[w][1]][k+f[w][2]][l+f[w][3]]+=dp[now][s][i][j][k][l];
	}
	LL tmp2=1,sum;
	for(int i=0; i<=40; ++i)
	{
		LL tmp3=tmp2;
		for(int j=0; j<=26; ++j)
		{
			LL tmp5=tmp3;
			for(int k=0; k<=18; ++k)
			{
				LL tmp7=tmp5;
				for(int l=0; l<=15; ++l)
				{
					if(sum=(dp[now][0][i][j][k][l]+dp[now][1][i][j][k][l])%mod) val[++cnt]=sum;
					if((tmp7*=7)>n) break;
				}
				if((tmp5*=5)>n) break;
			}
			if((tmp3*=3)>n) break;
		}
		if((tmp2<<=1)>n) break;
	}
}

int main()
{
	scanf("%lld%d",&n,&K);
	DP();
	std::sort(val+1,val+1+cnt,cmp);
	LL res=0;
	q.push(Node(1,1));
	while(!q.empty())
	{
		Node rt=q.top(); q.pop();
		(res+=rt.v)%=mod;
		if(!--K) break;
		if(rt.x!=rt.y){
			(res+=rt.v)%=mod;
			if(!--K) break;
			q.push(Node(rt.x+1,rt.y));
		}
		if(rt.x==1) q.push(Node(rt.x,rt.y+1));
	}
	printf("%lld",res);

	return 0;
}

考试代码

T1

恶心的Lucas+Ex_Lucas。

#include <cmath>
#include <cstdio>
#include <cctype>
#include <algorithm>
#define gc() getchar()
typedef long long LL;
const int N=2e4+23;

int n,n1,n2,m,A[10],B[10],fac[N],f[10][N],sum[N],cnt,pi[666],pk[666];

inline int read()
{
	int now=0;register char c=gc();
	for(;!isdigit(c);c=gc());
	for(;isdigit(c);now=now*10+c-'0',c=gc());
	return now;
}
LL FP(LL x,int k,LL p)
{
	LL t=1ll;
	for(; k; k>>=1,x=x*x%p)
		if(k&1) t=t*x%p;
	return t;
}
void Exgcd(int a,int b,int &x,int &y)
{
	if(!b) x=1,y=0;
	else Exgcd(b,a%b,y,x),y-=a/b*x;
}
inline LL inv(int a,int p)
{
	int x,y; Exgcd(a,p,x,y);
	return (LL)((x%p+p)%p);
}
LL Calc(int n,int m,int p)
{
	if(n<m) return 0ll;
	return 1ll*fac[n]*inv(fac[m],p)%p*inv(fac[n-m],p)%p;
}
LL Lucas(int n,int m,int p)//p=10007
{
//	printf("Lucas:C(%d,%d) mod %d = ",n,m,p);
	if(n<m) return 0ll;
	LL res=1;
	for(; m&&res; n/=p,m/=p)
		(res*=Calc(n%p,m%p,p))%=p;
//	printf("%d
",res);
	return res;
}
LL Fact(int n,int pi,int pk)
{
	if(!n) return 1ll;//边界！ 
	LL res=1ll;
	if(n>=pk)
	{
		for(int i=2; i<pk; ++i)
			if(i%pi) (res*=i)%=pk;
		res=FP(res,res/pk,pk);//这怎么写成res/pk了。。
	}
	for(int i=2; i<=n%pk; ++i)
		if(i%pi) (res*=i)%=pk;
	return res*Fact(n/pi,pi,pk)%pk;
}
LL C(int n,int m,int pi,int pk,int mod)
{
	LL a=Fact(n,pi,pk),b=Fact(m,pi,pk),c=Fact(n-m,pi,pk);int k=0;
	for(int i=n; i; i/=pi) k+=i/pi;
	for(int i=m; i; i/=pi) k-=i/pi;
	for(int i=n-m; i; i/=pi) k-=i/pi;
	LL res=a*inv(b,pk)%pk*inv(c,pk)%pk*FP(pi,k,pk)%pk;
	return res*(mod/pk)%mod*inv(mod/pk,pk)%mod;
}
LL Ex_Lucas(int n,int m,int p)
{
//	printf("Ex_Lucas:C(%d,%d) mod %d = ",n,m,p);
	if(n<m) return 0ll;
	if(n==m) return 1ll;
	int res=0;
	for(int i=1; i<=cnt; ++i)
		(res+=C(n,m,pi[i],pk[i],p))%=p;
//	printf("%d
",res);
	return (LL)res;
}
void Init_P(int p)
{
	for(int i=2,m=sqrt(p+1); i<=m; ++i)
		if(!(p%i))
		{
			int pk_=1;
			while(!(p%i)) p/=i,pk_*=i;
			pi[++cnt]=i, pk[cnt]=pk_;
			if(p==1) break;
		}
}

int main()
{
	freopen("equation.in","r",stdin);
	freopen("equation.out","w",stdout);

	int T=read(),p=read();
	fac[0]=1;
	if(p==10007)
		for(int x,y,i=1; i<p; ++i) fac[i]=fac[i-1]*i%p;
	else Init_P(p);
	while(T--)
	{
		n=read(),n1=read(),n2=read(),m=read()-n+n2;
		LL tot=0;
		for(int i=1; i<=n1; ++i) A[i]=read()-1,tot+=A[i];
		for(int i=1; i<=n2; ++i) B[i]=read(),m-=B[i];
//		printf("1: n:%d m:%d
",n,m);
		if(m<0) {puts("0"); continue;}
		if(!m) {puts("1"); continue;}
		if(m==1)
		{
			int res=n;
			for(int i=1; i<=n1; ++i) if(!A[i]) --res;
			printf("%d
",res);
			continue;
		}
		A[n1+1]=0;
		for(int j=0; j<=A[1]; ++j) f[1][j]=1;
		sum[0]=1;
		for(int j=1; j<=A[1]+A[2]; ++j) sum[j]=sum[j-1]+f[1][j];
		for(int tmp=A[1],tmp2=A[1]+A[2],i=2; i<=n1; ++i)//O(8m)
		{
			tmp+=A[i], f[i][0]=1;
			for(int j=1; j<=tmp; ++j)
			{
				if(j<std::min(j,A[i])) f[i][j]=sum[j];
				else f[i][j]=sum[j]-sum[j-std::min(j,A[i])-1], f[i][j]=(f[i][j]%p+p)%p;
//				printf("f[%d][%d]:%d %d bef:%d %d
",i,j,f[i][j],p,sum[j],sum[j-std::min(j-1,A[i])-1]);
			}
			tmp2+=A[i+1], sum[0]=1;
			for(int j=1; j<=tmp2; ++j)
			{
				sum[j]=sum[j-1]+f[i][j];
				if(sum[j]>=p) sum[j]-=p;
			}
		}
		LL res=0; f[n1][0]=1;
//		printf("2: n:%d m:%d
",n,m);
		tot=std::min(tot,(LL)m);
//		for(int i=0; i<=tot; ++i) printf("f:%d:%d
",i,f[n1][i]);
		if(p==10007) for(int i=0; i<=tot; ++i) (res+=1ll*f[n1][i]*Lucas(n+m-i-n1-1,m-i,p))%=p;
		else for(int i=0; i<=tot; ++i) (res+=1ll*f[n1][i]*Ex_Lucas(n+m-i-n1-1,m-i,p))%=p;
		printf("%I64d
",res);
	}
	return 0;
}

T2

答案都比标程正好大了一倍，气啊。

#include <cstdio>
#include <cctype>
#include <algorithm>
#define gc() getchar()
const int N=218,M=2018,INF=0x3f3f3f3f;

int n,m,P,Max,src,des,Enum,H[N],cur[N],fr[M],to[M],cap[M],nxt[M],cost[M],num[N],lev[N],pre[N],que[N];

inline int read()
{
	int now=0;register char c=gc();
	for(;!isdigit(c);c=gc());
	for(;isdigit(c);now=now*10+c-'0',c=gc());
	return now;
}
inline void AddEdge(int u,int v,int w)
{
	to[++Enum]=v, fr[Enum]=u, nxt[Enum]=H[u], H[u]=Enum, cap[Enum]=w;
	to[++Enum]=u, fr[Enum]=v, nxt[Enum]=H[v], H[v]=Enum, cap[Enum]=0;	
}
bool BFS()
{
	for(int i=src; i<des; ++i) lev[i]=des+1;
	lev[des]=0, que[0]=des; int h=0,t=1;
	while(h<t)
	{
		int x=que[h++];
		for(int i=H[x]; i; i=nxt[i])
			if(lev[to[i]]==des+1 && cap[i^1])
				lev[to[i]]=lev[x]+1, que[t++]=to[i];
	}
	return lev[src]<=des;
}
int Augment()
{
	int mn=INF;
	for(int i=des; i!=src; i=fr[pre[i]])
		mn=std::min(mn,cap[pre[i]]);
	for(int i=des; i!=src; i=fr[pre[i]])
		cap[pre[i]]-=mn, cap[pre[i]^1]+=mn;
//	printf("Augment:%d
",mn);
//	Max=std::max(Max,mn);
	return mn;
}
int ISAP()
{
	if(!BFS()) return 0;
	for(int i=src; i<=des; ++i) ++num[lev[i]],cur[i]=H[i];
	int x=src,res=0;
	while(lev[src]<=des)
	{
		if(x==des) x=src,res+=Augment();
		int way=0,mn=0;
		for(int i=cur[x]; i; i=nxt[i])
			if(lev[to[i]]==lev[x]-1 && cap[i])
				if(!mn) way=mn=i;
				else if(cap[way]>cap[i]) way=i;
		if(way) cur[x]=mn, pre[x=to[way]]=way;
		else
		{
			int mn=des;
			for(int i=H[x]; i; i=nxt[i])
				if(cap[i]) mn=std::min(mn,lev[to[i]]);
			if(!--num[lev[x]]) break;
			++num[lev[x]=mn+1];
			cur[x]=H[x];
			if(x!=src) x=fr[pre[x]];
		}
	}
	return res;
}

int main()
{
	freopen("costflow.in","r",stdin);
	freopen("costflow.out","w",stdout);

	n=read(),m=read(),P=read(),Enum=src=1,des=n;
	for(int u,v,w,i=1; i<=m; ++i)
		u=read(),v=read(),w=read(),AddEdge(u,v,w);//f&=(w==1);
	int res1=ISAP();
//	if(f) {printf("%d
%d.0000",res1,P); return 0;}
	for(int i=3; i<=Enum; i+=2) Max=std::max(Max,cap[i]);//应是从最大流后全局的反向弧找Max 
	printf("%d
%d.0000",res1,Max*P);

	return 0;
}

T3

#include <queue>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <algorithm>
#define gc() getchar()
#define mod (1000000007)
typedef long long LL;
const int N=1e6+7;

int K,f[N+3],sum[1005][1005];//size!
LL n;
std::priority_queue<int> q;

void Violence()
{
	for(int i=1; i<=n; ++i)
		for(int j=1; j<=n; ++j)
			++sum[f[i]][f[j]];
	for(int i=1; i<=n; ++i)
		for(int j=1; j<=n; ++j)
			if(sum[i][j]) q.push(sum[i][j]);
//	for(int i=1; i<=n; ++i,putchar('
'))
//		for(int j=1; j<=n; ++j)
//			printf("%d ",sum[i][j]);
	LL res=0;
	while(K-- && !q.empty()) res+=q.top(),q.pop();
	printf("%I64d",res%mod);
}

int main()
{
	freopen("gold.in","r",stdin);
	freopen("gold.out","w",stdout);

	scanf("%I64d%d",&n,&K);
	f[0]=1;
	for(int i=1; i<N; ++i)
		if(i%10) f[i]=f[i/10]*(i%10);
		else f[i]=0;
	f[0]=0;
//	for(n=100; n<=300; ++n,putchar('
'))
//		memset(sum,0,sizeof sum),printf("%I64d:	",n),Violence();
	if(n<=1000) {Violence(); return 0;}
	printf("%d",rand()%mod);

	return 0;
}