【集训队试题 礼物】数论、递归求值

题目见H8OJ 2142

http://www.zybbs.org/JudgeOnline/problem.php?id=2142

综合性比较强的求值题（这个名字好像不怎么样。。）

my solution：

 

本题有几个关键点，主要都是基于mod的快速求值

a.       列出目标式： 

b.       将p分解质因数p =（ai为质数）（共t个质因子）

c.       加速如下形式的式子的求值

S = K! mod (ai^bi) （a为质数）

即对p的t个质因子分别求值

目标式中在分子的阶乘式直接乘，在分母的则取其乘法逆元

d.       用中国剩余定理或直接累加将t个取模后的值合并得到答案

大概的步骤就是这样，但还有两个非常关键的问题还没有解决

A． 如何对分母的阶乘式取乘法逆元

由于式子中不是所有的项都与ai^bi互质，所以不能直接用欧拉定理，那我们就将所有分子分母中含ai的质因子提取出来，先约分再最后一次快速幂求出值来，提取完ai剩下的部分都与ai^bi互质，就能直接用欧拉定理求解

B． 如何加速下式

S = K! mod (ai^bi) （a为质数）

的求值

令mo = ai^bi

可以发现

1*2*…*(mo-1) mod mo

= (1 + mo)*(2 + mo)*…*(mo-1 + mo) mod mo

= (1 + mo*2)*(2 + mo*2)*…*(mo-1 + mo*2) mod mo

…

=(1 + mo*i)*(2 + mo*i)*…*(mo-1 + mo*i) mod mo

那么

1*2*…*K mod mo

=[1*2*…*(mo-1)]^(K div mo)

*mo*(2*mo)*(3*mo)*…*[(K div mo)*mo]

*(K div mo*mo+1)*(K div mo*mo+2)*…*K

其实就是将所以乘积mod mo相同的(1+i*mo)*(2+i*mo)*…*(mo-1 + i*mo)全部放到一起，i*mo放到一起，最后剩下的一部分放到一起

其中第一部分可以预处理+快速幂解决，第三部分可以朴素地算，第二部分提出一个mo后变成1*2*…*(K div mo)，可以递归地计算

由于每次递归都将问题规模缩为1/ai，所以递归次数为log级，每层的复杂度为   O(mo + logK)，总的复杂度就还是很低的，可以秒杀所有数据

至此，问题其实还并没有完全解决，我们还遗留了一些小问题没有解决。因为实际上由于我们要提出所有含ai的质因数，在上述快速求值的过程中，我们不能直接提出mo(实际上就是ai^bi)，而只能提出ai，虽然大体的递归计算过程差不多，但还是引发了一些细节问题，有如下几个

A． 我们求出的1*2*…*(mo-1) mod mo中mod ai = 0的项不能算进来，因为我们在对上面第二部分式子提出ai并递归计算的时候，实际上我们或直接或间接地将含有ai的项全都算进来了

B． 我们在上面的式子的第三部分朴素的计算部分中，mod ai = 0的项也不能算进来，理由同上

现在，问题才算基本完美解决，这里顺便提一下，最后那个合并答案的部分其实直接朴素地加是没有问题的，复杂度是有保证的，证明思路大概就是每次合并最多加当前这个ai^bi次，而即便是这样，最多也就加m*ai^bi次，这个次数其实是很少的，完全不会影响算法的效率

代码：

program syj; var lcm,k,x,t,n,m,p,s,ans,mo,r:int64; l,c,i,j,pt:longint; w:array[0..6]of int64; ph,a,q:array[0..20]of int64; function mul(a,n:int64):int64; begin if n=0 then mul:=1 else begin mul:=mul(a,n>>1); mul:=mul*mul mod mo; if odd(n) then mul:=mul*a mod mo; end; end; procedure work(n,z:int64); begin r:=1; if n<p then for j:=1 to n do r:=r*j mod mo else begin r:=mul(s,n div mo); for j:=n div mo*mo+1 to n do if j mod p<>0 then r:=r*j mod mo; inc(c,n div p*z); end; if z=1 then t:=t*r mod mo else t:=t*mul(r,ph[l]-1) mod mo; if n>=p then work(n div p,z); end; begin assign(input,'gift.in');reset(input); assign(output,'gift.out');rewrite(output); readln(p,n,m); for i:=1 to m do begin read(w[i]); x:=x+w[i]; if x>n then break; end; if x>n then begin writeln('Impossible'); close(input);close(output); halt; end; for i:=2 to trunc(sqrt(p)) do if p mod i=0 then begin inc(pt); a[pt]:=1; q[pt]:=i; while p mod i=0 do begin a[pt]:=a[pt]*i; p:=p div i; end; ph[pt]:=a[pt] div i*(i-1); end; if p>1 then begin inc(pt); a[pt]:=p; q[pt]:=p; ph[pt]:=p-1; end; ans:=0; lcm:=1; for l:=1 to pt do begin p:=q[l]; mo:=a[l]; t:=1; c:=0; s:=1; for i:=1 to mo-1 do if i mod p<>0 then begin j:=i; while j mod p=0 do j:=j div p; s:=s*j mod mo; end; work(n,1); work(n-x,-1); for i:=1 to m do work(w[i],-1); t:=t*mul(p,c)mod mo; while ans mod mo<>t do ans:=ans+lcm; lcm:=lcm*mo; end; writeln(ans); close(input);close(output); end.  

我一开始直接以ai分的段，居然也ac了。。