【POI08 PER】模意義下的的除法

模數這種東西到底還是好，將廣大OIER從高精度的深淵中解放出來

但是，模數也不總是那麼方便的

比方說，我們的操作中不可避免地出現了除法

大部分時候題目會告訴我們模數的數是質數，或者至少除數和模數的數互質，這樣我們就能用歐拉定理來解決這個問題了

然而，不幸的事終還是發生了

這道題所要求的模數的數就是任意的~~

題目本身不難，從後往前按位統計+樹狀數組即可，但由於可重排列公式中不可避免地出現了除法，怎麼辦？

盾哥從這個題的標程中摳出了下面這個不太優美的方法

先將要模數的數分解質因數，然後每次要乘或要除一個數的時候，就將那個數的質因子中跟模數的數的相同的提出來，這部分因子“暴力”（就是有除法就直接約分，要用的時候再快速冪算值），另外的部分就和要模數的數互質了，可以直接求其乘法逆元

由於一個數的質因數個數是log層級的，所以通過一大堆預先處理可以讓模數過程均攤O(logN)

可是，演算法常數還是非常大，怎麼辦呢？

我想了一個比較有效常數最佳化，就是算一個數的乘法逆元可以“記憶化”，加入這個最佳化程式能快上將近一倍

感覺演算法不太優美，求更漂亮的實現~~

 

代碼：

program syj;<br />const<br /> maxn=300005;<br />type<br /> arr=array[0..43]of longint;<br />var<br /> n,max,i,j,k,phm,mo,dt:longint; s,ans:int64;<br /> a,b,c,sh,ne,u:array[0..maxn]of longint;<br /> ph:array[0..maxn]of longint;<br /> d:array[0..43]of longint;<br /> z:array[0..43,0..maxn]of longint;<br /> f,g:arr;<br />procedure update(i:longint);<br />begin<br /> while i<=max do begin<br /> inc(b[i]);inc(i,i and-i);<br /> end;<br />end;<br />procedure ask(i:longint;var s:longint);<br />begin<br /> s:=0;<br /> while i>0 do begin<br /> inc(s,b[i]);dec(i,i and-i);<br /> end;<br />end;<br />function calc(i:longint):longint;<br /> procedure mul(n:longint);<br /> begin<br /> if n=1 then ph[i]:=i else begin<br /> mul(n>>1);<br /> ph[i]:=int64(ph[i])*ph[i] mod mo;<br /> if odd(n) then ph[i]:=int64(ph[i])*i mod mo;<br /> end;<br /> end;<br />begin<br /> if ph[i]>0 then exit(ph[i]);<br /> mul(phm);<br /> calc:=ph[i];<br />end;<br />function work(i,j:longint;var f:arr):int64;<br />var k:longint;<br />begin<br /> if i=0 then exit(0);<br /> work:=1;<br /> while i>1 do begin<br /> k:=sh[i];<br /> if u[k]>0 then inc(f[u[k]],j) else<br /> if j>0 then work:=work*k mod mo<br /> else work:=work*calc(k) mod mo;<br /> i:=ne[i];<br /> end;<br />end;<br />begin<br /> assign(input,'per.in');reset(input);<br /> assign(output,'per.out');rewrite(output);<br /> readln(n,mo);<br /> k:=mo;phm:=1;<br /> for i:=2 to trunc(sqrt(k)) do<br /> if k mod i=0 then begin<br /> inc(dt); d[dt]:=i; u[i]:=dt; j:=1;<br /> while k mod i=0 do begin<br /> k:=k div i;j:=j*i;<br /> end;<br /> phm:=phm*(j div i*(i-1));<br /> end;<br /> if k>1 then phm:=phm*(k-1);<br /> dec(phm,1);<br /> if (k>0)and(k<=n) then begin<br /> inc(dt);d[dt]:=k;u[k]:=dt;<br /> end;<br /> for i:=1 to dt do begin<br /> z[i,0]:=1;<br /> for j:=1 to n do z[i,j]:=int64(z[i,j-1])*d[i] mod mo;<br /> end;<br /> for i:=1 to n do begin<br /> read(a[i]);<br /> if a[i]>max then max:=a[i];<br /> end;<br /> if n>max then max:=n;<br /> for i:=2 to max do<br /> for j:=1 to max div i do<br /> if sh[i*j]=0 then begin<br /> sh[i*j]:=i; ne[i*j]:=j;<br /> end;<br /> s:=1; ans:=1;<br /> for i:=n downto 1 do begin<br /> ask(a[i]-1,k);<br /> fillchar(g,sizeof(g),0);<br /> j:=work(k,1,g);<br /> inc(c[a[i]]);<br /> s:=s*work(c[a[i]],-1,f) mod mo;<br /> j:=s*j mod mo;<br /> if j>0 then<br /> for k:=1 to dt do<br /> j:=int64(j)*z[k,g[k]+f[k]] mod mo;<br /> inc(ans,j);<br /> update(a[i]);<br /> s:=s*work(n-i+1,1,f) mod mo;<br /> end;<br /> writeln(ans mod mo);<br /> close(input);close(output);<br />end.<br /> 

速度比盾哥快，哈哈~