hud4048
核武的思路:http://hi.baidu.com/aekdycoin/blog/item/69e075f5e200bcf47709d769.html
参考完核武报告,没完全懂~翻了翻组合数学书,这实质上就是求元素可重复的圆排列--许多组合数学书里讲莫比乌斯反演时唯一的例题~
先求 r 种不同元素,放在圆的n个点上 , 有几种放法, 旋转后相同算一种.
做法是:设周期为为d,d|n,的方案数是Td(r) 则 每个这样的圆排列,对应着d个线排列
线排列总数是r^n,所以:
莫比乌斯反演后:
总的方案数M_n就是各种周期的和:
因 d|k|n , 设k'=k/d . n'=n/d
注意到后面一部分和等于n'的欧拉函数,故Mn可以写成:
问题基本解决,剩下就是核武报告里说差不多,把最小公约数是i的倍数的方法数求出,容斥原理搞搞.
最后的答案应该是:
其中r_n表示重量是n的倍数的石子种数
mu(n) 和 phi(n) 都可以在预处理时筛出来,
n的分解只需要搞一遍就够了~
遗憾的是,还是很慢,board的第一名时空复杂度都比我写的好~ym
SPOJ3871 GCDEX
这题的描述对我这种语言能力残缺的人来说,复述起来太难:
现在证明,T(n)是一个积性函数
证明:
设 m , n 满足 (m,n)=1:
则:
由于(m,n)=1,故dn与dm互素,所有dn,dm的组合正好取遍mn的约数
到了这里,只需考虑n等于素数的幂时T(n)的值:
根据T(n)的定义,显然有:
有了这个性质,就有可以在O(n)时间内用这个 牛X的筛法 把表打出来,然后求和还是O(n),查询什么的显然O(1)咯lol...
代码:
#include <stdlib.h> #include <stdio.h> #include <string.h> #include <math.h> #define FF(i,b) for(i=0;i<(b);i++) #define FOR(i,a,b) for(i=(a);i<=(b);i++) #define ZERO(x) do{memset(x,0,sizeof(x));}while(0) #define SWAP(a,b) do{int _tmp=a;a=b;b=_tmp;}while(0) #define INFI 2000000000 #define INFF 1e100 #define maxn 1000000 typedef long long LL; int p[maxn+1],ph[maxn+1],pk[maxn+1],T[maxn+1]; LL G[maxn+1]; void sevie() { memset(p,0,sizeof(p)); int i,j; p[0]=ph[0]=0;ph[1]=1;T[1]=1; for(i=2;i<=maxn;i++) { if(!p[i]){p[++p[0]]=i;ph[i]=i-1;pk[i]=i;T[i]=(i<<1)-1;} for(j=1;p[j]<=maxn/i;j++) { p[p[j]*i]=1; if(i%p[j]){ ph[i*p[j]]=ph[i]*p[j]; pk[i*p[j]]=p[j]; T[i*p[j]]=T[i]*T[p[j]]; } else{ ph[i*p[j]]=ph[i]*p[j]; pk[i*p[j]]=pk[i]*p[j]; if(i==pk[i])T[i*p[j]]=p[j]*T[i]+ph[i*p[j]]; else T[i*p[j]]=T[i/pk[i]]*T[pk[i]*p[j]]; break; } } } G[0]=0; for(i=1;i<=maxn;i++)G[i]=G[i-1]+T[i]-i; } int main ( int argc, char *argv[] ) { sevie(); int n; while(scanf("%d",&n)!=EOF) { if(!n)break; printf("%lld\n",G[n]); } return EXIT_SUCCESS; }