bzoj 2313 分组 2012-01-15
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了bzoj 2313 分组 2012-01-15相关的知识,希望对你有一定的参考价值。
http://www.zybbs.org/JudgeOnline/problem.php?id=2313
_______________________________________________
对于n要求方案数:
对于n的约数s,即将n个数分为s组,每组有n/s=p个数共有:
ans (p)=(C(p,n)*C(p,n-p)*C(p,n-2p)*...*C(p,p))/(s!)
=(n!*(n-p)!*...*p!)/((p!)^s*(n-p)!*(n-2p)!..1!*s!)
=(n!)/((p!)^s*s!)
种方案
所以总共的方案数就有 total=sigma( ans(si))
题目要求的是 m^total mod (10^9-401)
设prime=10^9-401
prime-1 分解质因数为 di[i]={2,13,5281,7283}
首先需要知道:①费尔马小定理: a^(p-1)=1(mod p)
② (k-1)! =-1 (mod k)
所以 m^(prime-1)=1 (1 mod prime);
m^ans=m^(ans mod (prime-1)) (mod prime);
所以现在需要求 ans mod (prime-1)
等价于求满足 ans =xi (mod di[i]) (1<=i<=4) 最后将每一个yi 用中国剩余定理合并即求出 x=ans( mod (prime-1))。
最后即可算出每一个x的和即为 total mod (prime-1)
对于 ans=(n!)/((p!)^s*s!) mod di
怎么求 n! mod di 呢?
因为 n!=1*2..(di-1)*di*(di+1)*..*(q*di)*(q*di+1)...*n (mod di)
=((di-1)! mod di)^q*(n-q*di)!*q!*(di)^q ( mod di)
=(-1)^q*(n-q*di)!*q!*di^q ( mod di)
因为(n-q*di)不是特别大(n-q*di)! mod di的值可以预处理出来
所以可以用递归求一个数的阶乘模di的值 求出来的值可以化简为 fac*di^pow
将分数ans 的分子存成 di^pow1*fac1
分母存成 di^pow2*fac2
易证pow1>=pow2 所以当 pow1>pow2时 ans= (fac1/fac2) (mod di);
当 pow1=pow2时 ans= 0 (mod di)
!!!程序合并ans =xi (mod di[i]) 用的是求线性模方程组的方法,具体百度
_______________________________________________
1 Program Stone; 2 const di:array[1..4]of longint=(2,13,5281,7283); 3 prime=999999599; 4 type data=record pow,fac:int64;end; 5 var i,j,t:longint; 6 n,m,ans:int64; 7 aj:array[1..4,0..8000]of int64; 8 s_up,s_low:array[1..4]of data; 9 ci:array[1..4]of int64; 10 function pow(c,x,p:int64):int64; //求 c^x mod p 的值 11 var i,k:int64; 12 begin 13 pow:=1; 14 i:=x; 15 k:=c; 16 while i>0 do 17 begin 18 if i mod 2=1 then pow:=k*pow mod p; 19 i:=i shr 1; 20 k:=k*k mod p; 21 end; 22 end; 23 24 function factorial(p,a:int64):data; //求 p! mod di[a] 的值,存储成 fac*di[a]^pow 25 var i,j:int64; 26 s:data; 27 begin 28 i:=p div di[a];j:=p mod di[a]; 29 s.fac:=1;s.pow:=0; 30 if i>0 then s:=factorial(i,a); 31 if i mod 2=1 then s.fac:=-s.fac; 32 s.fac:=(s.fac*aj[a,j]) mod di[a]; 33 s.pow:=(s.pow+i); 34 factorial:=s; 35 end; 36 function exgcd(a,b:int64;var x,y:int64):int64; //扩展欧几里德 37 38 var t:int64; 39 begin 40 if b=0 then begin x:=1;y:=0;exit(a);end; 41 exgcd:=exgcd(b,a mod b,x,y); 42 t:=x; 43 x:=y; 44 y:=t-(a div b)*y; 45 end; 46 function work(a,b,n:int64):int64; //求解 ax=n (mod b)的x 值 47 var y,g:int64; 48 begin 49 g:=exgcd(a,b,work,y); 50 work:=((work*(n div g) mod (b div g)+(b div g))mod (b div g)); //求出的解应该为最小整数解 51 end; 52 function lcm(a,b:int64):int64; //求最小公倍数 53 var i,j,k,gcd:int64; 54 begin 55 i:=a;j:=b; 56 while j<>0 do 57 begin 58 k:=i;i:=j;j:=k mod j; 59 end; 60 lcm:=a*b div i; 61 end; 62 Procedure doans(i:int64); //对于将n分成每组有i 个的方案数 63 var j:longint; 64 ai,yi,xi,k:int64; 65 s:data; 66 begin 67 for j:=1 to 4 do //s_up表示分子,s_low表示分母 68 begin 69 s_low[j]:=factorial(n div i,j); 70 s:=factorial(i,j); 71 s_low[j].pow:=s_low[j].pow+s.pow*(n div i); 72 s_low[j].fac:=s_low[j].fac*pow(s.fac,n div i,di[j]); 73 if s_up[j].pow<>s_low[j].pow then ci[j]:=0 74 else ci[j]:=work(s_low[j].fac,di[j],s_up[j].fac); 75 end; 76 ai:=di[1];yi:=ci[1]; //合并ans =xi (mod di[i]) 77 for j:=2 to 4 do 78 begin 79 xi:=work(ai,di[j],(ci[j]-yi)); 80 yi:=ai*xi+yi; 81 ai:=lcm(ai,di[j]); 82 end; 83 ans:=(ans+yi)mod (prime-1); 84 end; 85 Procedure main; 86 var i,j:longint; 87 s:data; 88 begin 89 if n>1 then ans:=2 else ans:=1; 90 for i:=1 to 4 do 91 s_up[i]:=factorial(n,i); 92 for i:=2 to round(sqrt(n)) do //枚举n的约数 93 if n mod i=0 then 94 begin 95 doans(i); 96 if n div i<>i then doans(n div i); 97 end; 98 end; 99 Begin 100 assign(input,‘input.in‘);reset(input); 101 readln(t); 102 for i:=1 to 4 do 103 begin 104 aj[i,0]:=1; 105 for j:=1 to di[i] do 106 aj[i,j]:=aj[i,j-1]*j mod di[i]; //预处理 i! mod di[i]的值 (i<=di) 107 end; 108 for i:=1 to t do 109 begin 110 readln(n,m); 111 main; 112 writeln(pow(m,ans,prime)); 113 end; 114 end. 115 116
以上是关于bzoj 2313 分组 2012-01-15的主要内容,如果未能解决你的问题,请参考以下文章