第一次接觸省選的知識點呢！zrf大佬在課堂上講的非常清楚，但由於本蒟蒻實在太菜了，直接掉線了。今天趕緊惡補一下。
那么這篇博客將分為兩塊，第一塊是FFT的推導和實現，第二塊則是FFT在OI上的應用
因為博主是蒟蒻，難免有些寫錯的地方，還請各位大佬不吝指正。
目標是能夠讓像博主這樣的蒟蒻都能學會FFT(都有耐心看完這篇博客)
一、FFT的推導與實現
1、多項式的表示
最常見的表示方式自然是系數表示
誒誒誒，別走啊，我說清楚點還不行嗎？
其實就是我們常見的表達方式
這種表達式的優勢在於我們可以用O(n)的復雜度快速求值
這稱之為霍納法則，很明顯這是對的，因為實質上只是合並同類項。
但是遇到乘法系數表示的復雜度就不盡人意了
先等等，多項式乘法怎么乘呢？
就和一般的乘法一樣列豎式計算啊！
好的，我們繼續說系數表示下多項式乘法的復雜度。
比如說
令C(x)=A(x)*B(x)，則
這個公式看着非常復雜，但結合上面多項式乘法的原理和定義應該可以看得懂，最重要的是，兩個sigma有助於清晰直觀地看出復雜度是O(n^2)的。
如果有兩個長度為100000的多項式，讓我們求他們的乘積，那么用這種方法，我們只好說一句“告辭”了。
但如果一定要你求呢？雖然好像強人所難，但這仍然是可做的。
“當手里的每一張牌都是壞牌，想要贏一把的唯一辦法就是打破游戲規則。”--保羅·奧斯特《幻影書》《作文素材》
是的，我們發現系數表示已經阻礙住了我們的步伐，所以我們不得不舍棄它，使用一種新的表示方式--點值表示
什么是點值表示？
來來來，先一起迷糊一下：
一個次數界為n的多項式A(x)的點值表示就是一個由n個點值所組成的集合
{(x0,y0),(x1,y1),...,(xn-1,yn-1)}
對k=0,1,2...n-1，xk各不相同，yk=A(xk)
什么，你居然看懂了，好吧，我真是太菜了。
意思就是對於多項式函數，我們將n個不同的值代入該函數（n是指該多項式的最高次），將代入值與所得值組成一對，總共n對，來表示該函數。
比如說
我們代入0,1,2（當然，你要代什么198,328,648或者233,666,19260817之類的也沒人攔你）
則可以得出該函數的點值表示
{(0,1),(1,2),(2,5)}
好的，那么點值表示有什么優點呢？
假設
A(x)={(x0,y0),(x1,y1),...,(xn-1,yn-1)}
B(x)={(x0,g0),(x1,g1),...,(xn-1,gn-1)}
C(x)=A(x)*B(x)
噫，我記得A(x0)不就等於y0,B(x0)不就等於g0嗎？
因此
C(x0)=A(x0)*B(x0)=y0*g0
所以C的點值表達式是什么呢？
C(x)={(x0,y0*g0),(x1,y1*g1),...,(xn-1,yn-1*gn-1)}
復雜度竟然是驚人的O(n)!
可優越了！
但是不得不潑一盆冷水，沒有題目會選擇給出點值表達的多項式，然后讓你也還給他點值表達的多項式，而且就算給了，兩個多項式的x也不會都是分別相等的。
所以你不得不把讀進來的系數表達式轉換成點值，再在乘完以后把點值轉換回系數。
這樣問題就來了，點值表達式真的只對應一個系數表達式嗎？
首先我們先來介紹一個高科技的術語：插值
其實在這里插值什么的並不高深，可以理解為將一個多項式從點值表示轉換成系數表示。
然后就要來證明插值多項式的唯一性（點值表達式只對應一個系數表達式）了
用系數的方式表示y0如下
y0=a0+a1x1+a2x1^2+...+an-1x^n-1
然后你可以用矩陣乘法來表示y0
接着我們可以把y1,y2,y3.....全部代進去，放到一塊
就會成這樣。
左邊這個矩陣在各種矩陣中是特殊的存在，被稱之為范德蒙德矩陣
他的行列式值det(A)為
因為之前保證了xk各不相同，所以只要j不等於k，xk-xj一定不等於0
即det(A)不等於0
同時我們知道，如果有一個矩陣他沒有逆矩陣，他的det一定為0
所以這個矩陣一定有逆矩陣。
有逆矩陣代表着由唯一的矩陣使
即有唯一的a
所以在對k=0,1,2...n-1，xk各不相同的情況下，插值多項式具有唯一性。
證畢了，但是結束了嗎？
真遺憾並非如此，我們知道用求x所對應的點值復雜度是O(n),然后要算n個值，復雜度是O(n^2)，所以還是gg的！
那你還說用點值表示是優越的？！
先別激動，我們可以通過巧妙選取x的方式來優化這個算值的過程，使得復雜度成為O(nlogn)。
且聽我慢慢道來。
首先我們需要引入一個概念
單位復數根
n次單位復數根是指一個復數ω，滿足
那么我們可以做以下推斷：
所以n次單位根的解只有n個。
我們將n次復數根的第m個解定義為
其中m同樣表示為
的m次方，根據上面的定義，很顯然這是成立的。
所以
然后根據復平面上的乘法幾何意義滿足角度為輔角相加，長度為長度相乘的原則我們可以得出：
整個圓的度數為單位根與x軸正半軸夾角的n倍，因為單位長度所以長度沒有變，角度增加如前所述
然后自然
以此類推，每個單位根之間的夾角均相同，整個圓正好能排n個單位根。
即n個單位復數根均勻的分布在單位圓上。
然后再給出一個知名的公式（別嚇着了！這公式只是給你裝逼用的！）
其中u是弧度制的
那么
這玩意不求理解，只是為了求值用的~當然其實在寫程序的時候你會發現壓根用不着e
好的，那么推了這么多，為的其實是引出幾個引理（引理名稱以及定義均來自於《算法導論》）。
消去引理：
對於任何整數n≥0，k≥0，以及d＞0滿足
看着很高大上，但其實很簡單
請問下圖兩條邊有什么區別？
肯定沒區別的啦！
就跟1/4=2/8一個道理啊。
從而可以推出
折半引理：
如果n>0，那么2n個2n次單位復數根的平方的集合就是n個n次單位復數根的集合。
這可以從感官上來證明，
那么平方完以后就相當於輔角擴大到原來的兩倍。
所以左邊的就膨脹成了右邊的。
是的，如我們所見，獲得每個單位根正好兩次。
也就是（這個非常重要！）
然后
求和引理：
對於任意整數n≥1和不能被n整除的非負整數k，有（這個也很重要！）
這是怎么證的呢？
sigma看着有點難受，我們把它拆開來。
這不是一個等比序列嗎？
來來來，為您送上等比數列的求和公式
那么代入上式得：
那么如果k能被n整除呢？
DFT
好的，上面的引理都整出來了，我們開始找回我們很早很早之前提到的多項式的系數表示法
假設我們已經知道了所有的aj（廢話，不知道還做什么題啊！別忘了出題人給你的是系數表示法，然后我們要求點值表示法（滑稽））
然后我們根據之前的乘法，我們需要2n個點值，以保證獲得唯一插值的C
將我們已知的值
代入A(x)得
然后向量yk={y0,y1,y2...}就是向量a{a0,a1,a2...}的離散傅里葉變換（DFT），恭喜你get到了一個新的可以裝逼的專業術語了！
我們記
由此我們終於能夠引出FFT了！
FFT
假設我們有函數
我們可以把它拆成兩半
然后我們合並一下公因數
那么原式可表示為
這似乎還不夠明顯，讓我們再設兩個公式
然后不是非常明顯了嗎？
這有什么用呢？
還記得折半引理嗎？
然后我們只需要知道等號右邊那個玩意的值，我們就可以搞定左邊的了！
然后這就是整個FFT最巧妙的一步！
原式被化成了這個樣子：
用n/2的數據可以求出n個值，這就相當於縮小了一半的大小！
所以復雜度就變成了O(nlogn)!
當然，n必須是2的冪次，如果數量不夠的話，必須往高項補零！
這樣就可以寫出代碼了：
先別得意的太早！我們只是找到了將系數表示法通過快速求值變成點值表示法的方式，然后，應該怎么把點值表示法快速插值成系數表示法呢？
IDFT&IFFT
我們曾經用范德蒙德矩陣證明過插值唯一性，在那里我們有得到過這個玩意：
只需要給已知點值乘上一個范德蒙德矩陣的逆矩陣即可
那么現在我們的范德蒙德矩陣是什么樣的呢？
我們把這個矩陣單獨拎出來
這個矩陣我們略微找下規律
j行k列的值V(j,k)為：
那這個范德蒙德矩陣的逆矩陣是什么？
哦，對了，逆矩陣的定義
然后I是單位矩陣，長這樣：
只有i行i列的值為1，其他都為0
定理30.7
對j,k=0,1,...,n-1,此逆矩陣的(j,k)處元素為
這個定理其實規定了逆矩陣長這樣。
怎么證？
我們先把求和引理拖下來
對於任意整數n≥1和不能被n整除的非負整數k，有
那么范德蒙德矩陣與我們剛剛所設的矩陣的乘積的(j1,j2)處是多少呢？
很明顯，因為j1,j2=0,1,...,n-1
所以
顯然
很明顯范德蒙德矩陣與我們剛剛所設的矩陣的乘積就是一個單位矩陣
我們將逆矩陣代入公式，就得到了IDFT
這玩意和DFT很像不是嗎？
然后你可別問我
怎么求……
所以發現了嗎？
DFT是順時針旋轉，那么IDFT就是逆時針旋轉。
那么只需要在FFT上改動很小一點
把初始的Y坐標從正翻成負的，那么就搞定了！
別忘了我們的1/n！IFFT在最后還要將結果都除以n！
那么IFFT的代碼也就出來了。
既然這么像，不如整合一下，用一個kd函數表示DFT還是IDFT，如果kd為1，則做DFT，為-1做IDFT。
由此我們得到了遞歸版的FFT
那么多項式乘法的nlogn做法就寫出來了！
代碼如下：（遞歸版）
#include<cmath> #include<vector> #include<cstdio> #include<cstring> #include<iostream> #include<complex> #include<complex.h> #include<algorithm> #define hi puts("hi!"); using namespace std; double Pi=acos(-1.0); struct comp { double r,i; comp(){r=0,i=0;} comp(double a,double b):r(a),i(b){}; void print() { printf("%.12lf %.12lf\n",r,i); } }; inline comp operator +(const comp a,const comp b) { return comp(a.r+b.r,a.i+b.i); } inline comp operator -(const comp a,const comp b) { return comp(a.r-b.r,a.i-b.i); } inline comp operator *(const comp a,const comp b) { return comp(a.r*b.r-a.i*b.i,a.i*b.r+b.i*a.r); } int e,m; vector<comp>x,y; vector<comp> DFT(vector<comp>a,int kd) //disastrous fatal TLE { int n=a.size(); if(n==1) { return a; //Èç¹ûֻʣÏÂÒ»¸öÊý£¬ÄÇô¾ÍÖ±½Ó·µ»ØaÊý×é } vector<comp> a0,a1,y0,y1,ans; for(int i=0;i<n;i++) { if(i&1) { a1.push_back(a[i]); //ÆæÊýÏî·Öµ½a1Êý×é } else { a0.push_back(a[i]); //żÊýÏî·Öµ½a0Êý×é } } y0=DFT(a0,kd); y1=DFT(a1,kd); //·ÖÖÎ comp wn=comp((cos(2.0*Pi/(double) n)),(kd*sin(2.0*Pi/(double) n))); //Çó³öÒ»¸ön´Îµ¥Î»¸´Êý¸ùµÄÖµ comp w=comp(1.0,0.0); //¶¨ÒåµÚ0¸öµ¥Î»¸´Êý¸ùµÄÖµ ans.resize(a.size()); for(int i=0;i<n/2;i++,w=w*wn) //»ñÈ¡ÏÂÒ»¸ön´Îµ¥Î»¸´Êý¸ùµÄÖµ { ans[i]=y0[i]+y1[i]*w; //¶ÔÓ¦ÉÏÎÄһʽ ans[i+(n>>1)]=y0[i]-y1[i]*w; //¶ÔÓ¦ÉÏÎĶþʽ } return ans; } int main() { scanf("%d%d",&e,&m); for(int i=0;i<=e;i++) { int tmp; scanf("%d",&tmp); x.push_back(comp((double)tmp,0.0)); } for(int i=0;i<=m;i++) { int tmp; scanf("%d",&tmp); y.push_back(comp((double)tmp,0.0)); } int limit=1; while(limit<=e+m) { limit<<=1; } for(int i=e;i<limit-1;i++) { x.push_back(comp(0.0,0.0)); } for(int i=m;i<limit-1;i++) { y.push_back(comp(0.0,0.0)); } x=DFT(x,1); y=DFT(y,1); for(int i=0;i<limit;i++) { x[i]=x[i]*y[i]; } x=DFT(x,-1); for(int i=0;i<=e+m;i++) { printf("%d ",(int)(x[i].r/limit + 0.5)); } return 0; }
但是當我們信心滿滿的提交上去時，我們會發現
為什么會這樣呢？
是算法錯了嗎？
我們去數據比較友好的loj上測試一下
這說明算法沒錯。
但是耗時巨久，內存開銷巨大
想來也是，我們的FFT是遞歸版的，傳遞的是vector函數，不RE都奇怪啊！
那有沒有方法優化呢？
有的，首先是公用子表達式上的優化。
因為復數乘法比較慢，
這里y1[i]*w出現了兩遍，可以這么寫
可惜治標不治本，它的函數傳遞問題還是不能解決。
不妨把它從遞歸版改成迭代版吧！
FFT的迭代優化
要把它從遞歸改成迭代的，我們需要理解它的遞歸結構。
如圖為8個數FFT的分治過程。
發現了嗎？
沒發現我們再看一下程序顯示的結果
對於數據
1 2
1 2
1 2 1
第一組1 2來說
一開始我們的數組是這樣的（0,1,2,3）
但是我們實際遞歸的順序是（0,2,1,3）
所以遞歸的時候我們是按照上列排列合並的，
對於數列ans，a0是a0[0]，a4是a1[0]，並不是按照我們所理解的a0,a1
所以迭代時我們也需要按照上述格式現將數組排好。
那么怎么把原數組排成這樣子呢？
找一找規律
你需要萬能的二進制！
一般的數字1,2,3,4,5,6,7,8的二進制如下，如果反轉一下呢？
wow！這不就是我們需要排的順序嗎？
現在的問題變成了如何快速求出i的反轉數
這其實是個DP問題。
x的兩倍相當於x<<1那么反轉之后r[x]相當於r[x>>1]>>1
這還是好理解的，然后如果是奇數，它的末尾為一，反轉之后則開頭為一
於是DP公式就推出來了
r[x]=r[x>>1]>>1|(x&1)<<(l-1)
其中l是指當前情況二進制的最高位數
代碼如此寫：
之后的操作就和遞歸是一樣的了，一點一點的並上去。
那么最后的優化也做完了，我們可以開始寫FFT的迭代版了。
再往上一提交，舒爽！
可見常數已經大大減小了，但是由於數據范圍為1000000
內存占用還是很高的
但這也是沒辦法的事情不是嗎？（容♂許之心泛濫中）
以上，本篇博客的第一部分：FFT的推導與實現終於結束了，不記公式和代碼總共3600+，也挺不容易的……
然后之后就是狗尾續貂的第二部分：FFT在OI上的應用了，因為這篇博客是我邊學邊寫的，所以應用並不高深，主要是幾道例題，然后就全靠練習了……
二：FFT在OI上的應用
首先當然是我們最熟悉的多項式乘法了！
那么就先上模板吧！
多項式乘法（洛谷P3803/loj108/uoj34）
題目描述
給定一個n次多項式F(x)，和一個m次多項式G(x)。
請求出F(x)和G(x)的卷積。
輸入輸出格式
輸入格式：
第一行2個正整數n,m。
接下來一行n+1個數字，從低到高表示F(x)的系數。
接下來一行m+1個數字，從低到高表示G(x))的系數。
輸出格式：
一行n+m+1個數字，從低到高表示F(x)∗G(x)的系數。
輸入輸出樣例
輸入樣例#1：
1 2
1 2
1 2 1
輸出樣例#1：
1 4 5 2
原理剛才已經基本上顯示過了，這就是一個fft的板子
代碼如下：
#include<cmath> #include<vector> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define Pi acos(-1) using namespace std; int r[4000010]; int limit,n,m,l; struct comp { double r,i; comp(){} comp(double a,double b):r(a),i(b){} void print() { printf("%.12lf %.12lf\n",r,i); } }a[4000010],b[4000010]; inline comp operator +(const comp a,const comp b) { return comp(a.r+b.r,a.i+b.i); } inline comp operator -(const comp a,const comp b) { return comp(a.r-b.r,a.i-b.i); } inline comp operator *(const comp a,const comp b) { return comp(a.r*b.r-a.i*b.i,a.i*b.r+a.r*b.i); } void FFT(comp *x,int pd) { for(int i=0;i<limit;i++) { if(i<r[i]) { swap(x[i],x[r[i]]); } } for(int mid=1;mid<limit;mid<<=1) { comp wn=comp(cos(Pi/mid),pd*sin(Pi/mid)); for(int l=mid<<1,j=0;j<limit;j+=l) { comp w=comp(1.0,0.0); for(int k=0;k<mid;k++,w=w*wn) { comp u=x[k+j]; comp v=w*x[k+j+mid]; x[k+j]=u+v; x[k+j+mid]=u-v; } } } } int main() { scanf("%d%d",&n,&m); limit=1; while(limit<=n+m) { limit<<=1; l++; } for(int i=0;i<limit;i++) { r[i]=r[i>>1]>>1|((i&1)<<(l-1)); } for(int i=0;i<=n;i++) { scanf("%lf",&a[i].r); } FFT(a,1); for(int i=0;i<=m;i++) { scanf("%lf",&b[i].r); } FFT(b,1); for(int i=0;i<limit;i++) { a[i]=a[i]*b[i]; } FFT(a,-1); for(int i=0;i<=n+m;i++) { printf("%d ",(int)(a[i].r/limit+0.5)); } return 0; }
然后我們再來一道模板題過渡一下
A*Bproblem升級版(洛谷1919/RQNOJ314/BZOJ2179)
給出兩個n位10進制整數x和y，你需要計算x*y。
輸入輸出格式
輸入格式：
第一行一個正整數n。 第二行描述一個位數為n的正整數x。 第三行描述一個位數為n的正整數y。
輸出格式：
輸出一行，即x*y的結果。（注意判斷前導0）
輸入輸出樣例
輸入樣例#1：
1
3
4
輸出樣例#1：
12
說明
數據范圍：
n<=60000
這道題其實還是很巧妙的，把x按位數化成多項式，y也如法炮制，那么問題就變成了求多項式x與多項式y的乘積
為什么這樣可以呢？
設x=23，那么轉換成多項式就是{2,3}
y=24，轉換成多項式就是{2,4}
兩者的成績為
{4,14,12}
這可能看不出什么，但是進下位之后呢？
{5,5,2}
那啥，23*24不就等於552嗎？其實這本來就是乘法的豎式計算啊！
所以就可以用FFT錘爆了。
對於處理前導零，我選擇的是從后往前做x和y的多項式
如{0,0,2,3}
這樣最后的解也就自動對齊了！
代碼如下：
#include<cmath> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define Pi acos(-1) using namespace std; struct comp { double r,i; comp() {} comp(double a,double b):r(a),i(b) {} } a[250010],b[250010]; inline comp operator +(const comp a,const comp b) { return comp(a.r+b.r,a.i+b.i); } inline comp operator -(const comp a,const comp b) { return comp(a.r-b.r,a.i-b.i); } inline comp operator *(const comp a,const comp b) { return comp(a.r*b.r-a.i*b.i,a.i*b.r+a.r*b.i); } int l,limit,r[250010]; void FFT(comp *x,int pd) { for(int i=0; i<limit; i++) { if(i<r[i]) { swap(x[i],x[r[i]]); } } for(int mid=1; mid<limit; mid<<=1) { comp wn=comp(cos(Pi/mid),pd*sin(Pi/mid)); for(int l=mid<<1,j=0; j<limit; j+=l) { comp w=comp(1.0,0.0); for(int k=0; k<mid; k++,w=w*wn) { comp u=x[k+j]; comp v=w*x[k+j+mid]; x[k+j]=u+v; x[k+j+mid]=u-v; } } } } int main() { char s1[100000],s2[100000]; int ans[250010],n; scanf("%d\n",&n); limit=1; while(limit<=2*n) { limit<<=1; l++; } for(int i=0; i<limit; i++) { r[i]=r[i>>1]>>1|((i&1)<<(l-1)); } gets(s1); for(int i=0; i<n; i++) { a[limit-i-1].r=s1[i]-'0'; } FFT(a,1); gets(s2); for(int i=0; i<n; i++) { b[limit-i-1].r=s2[i]-'0'; } FFT(b,1); for(int i=0; i<limit; i++) { a[i]=a[i]*b[i]; } FFT(a,-1); for(int i=0; i<=n*2; i++) { ans[i]=(int)(a[limit-i].r/limit+0.5); } int x=0; for(int i=n*2; i>=0; i--) { ans[i]+=x; x=ans[i]/10; ans[i]%=10; } int len; for(int i=0; i<=n*2; i++) { if(ans[i]) { len=i; break; } } for(int i=len; i<=n*2; i++) { printf("%d",ans[i]); } }
來道應用題：
HDU4609 3-idiots
Problem Description
King OMeGa catched three men who had been streaking in the street. Looking as idiots though, the three men insisted that it was a kind of performance art, and begged the king to free them. Out of hatred to the real idiots, the king wanted to check if they were lying. The three men were sent to the king's forest, and each of them was asked to pick a branch one after another. If the three branches they bring back can form a triangle, their math ability would save them. Otherwise, they would be sent into jail.
However, the three men were exactly idiots, and what they would do is only to pick the branches randomly. Certainly, they couldn't pick the same branch - but the one with the same length as another is available. Given the lengths of all branches in the forest, determine the probability that they would be saved.
Input
An integer T(T≤100) will exist in the first line of input, indicating the number of test cases.
Each test case begins with the number of branches N(3≤N≤105).
The following line contains N integers a_i (1≤a_i≤105), which denotes the length of each branch, respectively.
Output
Output the probability that their branches can form a triangle, in accuracy of 7 decimal places.
Sample Input
2
4
1 3 3 4
4
2 3 3 4
Sample Output
0.5000000
1.0000000
題意大致為給你n（n<=100000）根棒子，每根棒子的長度為ai（ai<=100000）
求任意取出三根能組成三角形的概率。
這道題當然是要先用O(n^3)的方法TLE一遍了！
好吧，廢話不多說，我們開始懟正解（正經臉）
如果我們拿出一根長度為a[i]的棒子，假設另外兩根都比它要短，那么此時我們的方案數該怎么算呢？
首先兩根棒子的長度和要比a[i]長，每根的單一長度都比a[i]短。
如果我們將a數組排個序，那么a[i]之前的就都比它短了。
然后我們要取長度和大於它的
這里就可以用FFT優化了
如果把樣例
4
1 3 3 4
按照長度為i的數量轉化成新的數組，為：
{0,1,0,2,1}
然后組合數就是數組的卷積
什么是卷積？就是多項式的乘積啊！
比如說上面數組的卷積就是
{0,0,1,0,4,2,4,4,1}
其中的第i項就是指任取兩根組成的長度個數
但其實我們把兩根一樣的也加起來了，所以要去掉
{0,0,0,0,4,2,2,4,0}
然后a+b和b+a我們都算了
所以答案應該再除二
{0,0,0,0,2,1,1,2,0}
接着我們會發現長度和比a[i]大的不止都是長度比它小的
然后我們要減去這些情況
第一種是取了自己的
共有n-1種
第二種是取了一個比自己大的，一個比自己小的
共有（n-i-1）*i種
第三種是取了兩個都大於自己的
共有（n-i-1）*（n-i-1）/2種
然后就是O（n）的遍歷跑答案了！
這道題還是很好的，完整代碼如下
#include<cmath> #include<cstdio> #include<cstring> #include<iostream> #include<algorithm> #define Pi acos(-1) using namespace std; int n,a[400010],r[400010],l,maxx,limit; long long sum[400010],cnt[400010]; struct comp { double r,i; comp() {} comp(double a,double b):r(a),i(b) {} } num[400010]; inline comp operator +(const comp a,const comp b) { return comp(a.r+b.r,a.i+b.i); } inline comp operator -(const comp a,const comp b) { return comp(a.r-b.r,a.i-b.i); } inline comp operator *(const comp a,const comp b) { return comp(a.r*b.r-a.i*b.i,a.i*b.r+a.r*b.i); } void FFT(comp *a,int kd) { for(int i=0; i<limit; i++) { if(i<r[i]) { swap(a[i],a[r[i]]); } } for(int mid=1; mid<limit; mid<<=1) { comp wn=comp(cos(Pi/mid),kd*sin(Pi/mid)); for(int l=mid<<1,j=0; j<limit; j+=l) { comp w=comp(1.0,0.0); for(int k=0; k<mid; k++,w=w*wn) { comp u=a[k+j]; comp v=w*a[k+j+mid]; a[k+j]=u+v; a[k+j+mid]=u-v; } } } } int main() { int t; scanf("%d",&t); while(t--) { scanf("%d",&n); limit=1; l=0,maxx=0; for(int i=0; i<n; i++) { scanf("%d",&a[i]); maxx=max(maxx,a[i]); } sort(a,a+n); while(limit<=maxx*2) { limit<<=1; l++; } for(int i=0; i<limit; i++) { r[i]=r[i>>1]>>1|(i&1)<<(l-1); num[i].r=num[i].i=0.0; } for(int i=0; i<n; i++) { num[a[i]].r+=1.0; } FFT(num,1); for(int i=0; i<limit; i++) { num[i]=num[i]*num[i]; } FFT(num,-1); for(int i=0; i<limit; i++) { cnt[i]=(long long)(num[i].r/limit+0.5); } for(int i=0; i<n; i++) { cnt[a[i]<<1]--; } for(int i=0; i<limit; i++) { cnt[i]>>=1; } sum[0]=0; for(int i=1; i<limit; i++) { sum[i]=sum[i-1]+cnt[i]; } long long ans=0; for(int i=0; i<n; i++) { ans+=sum[limit-1]-sum[a[i]]; ans-=(n-1); ans-=(long long)(n-1-i)*i; ans-=(long long)(n-i-1)*(n-i-2)/2; } long long tot=(long long)n*(n-1)*(n-2)/6; printf("%.7lf\n",(double) ans/tot); } }
本來還想再多寫幾道題的，但是馬上要省選了，實在沒有時間，那就到這吧……
話說都看到這里了，那就點個推薦吧！