初探FFT（快速傅里葉變換）

第一次接觸省選的知識點呢！zrf大佬在課堂上講的很是清楚，但因爲本蒟蒻實在太菜了，直接掉線了。今天趕忙惡補一下。

那麼這篇博客將分爲兩塊，第一塊是FFT的推導和實現，第二塊則是FFT在OI上的應用

由於博主是蒟蒻，不免有些寫錯的地方，還請各位大佬不吝指正。

目標是可以讓像博主這樣的蒟蒻都能學會FFT(都有耐心看完這篇博客)

 

1、FFT的推導與實現

一、多項式的表示

最多見的表示方式天然是係數表示

誒誒誒，別走啊，我說清楚點還不行嗎？

 其實就是咱們常見的表達方式

這種表達式的優點在於咱們能夠用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;
}

View Code

可是當咱們信心滿滿的提交上去時，咱們會發現

爲何會這樣呢？

 是算法錯了嗎？

咱們去數據比較友好的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); } }

原本還想再多寫幾道題的，可是立刻要省選了，實在沒有時間，那就到這吧……

話說都看到這裏了，那就點個推薦吧！