【FFT-快速傅立葉變換】


HDU-1402  A * B Problem Plus

題意:給定兩個整數,整數的長度最多能達到50000位,輸出兩個整數的乘積。

分析:題意非常的明了,一個驚世駭俗的想法是使用兩個數組將整數保留起來,然后模擬我們平常手算時的乘法,不過這樣一來時間復雜度將是O(N^2),由於N過大,因此此題因尋求更加快速的解法。

  對於任何一個N位的整數都可以看作是An*10^(n-1) + An-1*10^(n-2) + ... + A2*10^2 + A1*10 + A0。如果把10看作是一個自變量,那么任何一個整數就可以視作為一個多項式,兩個整數相乘也便可以看作是兩個多項式相乘。對於一個多項式,我們平時所接觸到的多是其系數表示法,普通的相乘也就是建立在兩個整數均采用系數表示法的基礎上進行的。那么要使得計算多項式相乘的復雜度下降的另一種方式就是尋找一種新的表示多項式的方法......

  若一個多項式的最高階位為N-1,那么取N個點對(xi, yi)就能夠唯一確定這個多項式,可以想象成有N個系數需要N個方程去求解。那么在此就可以尋找點對來表示一個多項式,對於一個大的數,看作多項式之后,那么舍棄掉原來以10為自變量的取值,而選取其他值,再通過計算多項式An*xi^(n-1) + An-1*xi^(n-2) + ... + A2*xi^2 + A1*xi + A0來保存這個多項式的信息。需要選取N個xi形成N對(xi, yi)方可唯一確定原來各個項前的系數,通過選取1的N次單位復根即可,並且利用單位復根的性質,可以使得計算量下降。

  通過點值法表示多項式后,計算乘法也就是O(N)的時間了,由於兩個數相乘使得項數變多,因此需要在之前盡可能多取點。FFT算法能夠在O(NlogN)時間內將系數法轉化為點值法,相乘后再有點值法轉為系數法,該題就是使用的這個方法。

  順便說下一FFT過程中,計算葉子DTF時采用的二進制平攤反轉置換,其作用是為了避免算法的遞歸而實現自底向上的計算方式。回顧一下在計算原串DFT的時候,假設離散點數為0-7,那么有以下過程:

(0 1 2 3 4 5 7) = (0 2 4 6) + (1 3 5 7)                           1
(0 2 4 6) = (0 4) + (2 6)                                               2
(1 3 5 7) = (1 3) + (5 7)                                               2
(0 4) = (0) + (4)                                                          3
(2 6) = (2) + (6)                                                          3
(1 3) = (1) + (3)                                                          3
(5 7) = (5) + (7)                                                          3

分析這些分組的二進制位會發現,第1次分組是根據第0位是否為1來划分的,即奇偶性;第2次分組是根據第1位是否為1來划分的;第三次分組是根據第2位是否為1來划分的。這個特性與與一般的按照大小划分的數很類似(首先按照最高為是否為1划分,然后是次高位...),因此就可以通過一個算法來使得00...00 - 11...11這樣遞增的序列中的每一個數實現高位和低位的翻轉,二進制平攤反轉置換就是用於達到這個目的。

算法從1開始到N-2(FFT算法要求N必須是2的冪,保證每次折半之后不會出現奇數),因此0和N-1翻轉后還是本身,接着維護好一個下標 j ,這個數就是與遞增中的第 i 個數翻轉之后對應的數,初始化 j 的下標示 N/2,這個數要是從后往前來定義二進制數的話,就會是1,例如若N=8,那么4的二進制位為100,假定只有3位2進制位組成2進制數,從右往左看,其值為1。接下來就是要找 j 的下一個數了,這個數從右往左看應該是2才能夠滿足要求,於是從右往左尋找,遇到1變為0,知道遇到0就跳出,並且將該位賦值為1,這個偉大的過程的作用僅僅只是給在右往左定義的二進制數 j 加了一個1。

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;

struct C {
    double r, i;
    C() {}
    C(double _r, double _i) : r(_r), i(_i) {}
    inline C operator + (const C & a) const {
        return C(r + a.r, i + a.i);
    }
    inline C operator - (const C & a) const {
        return C(r - a.r, i - a.i);
    }
    inline C operator * (const C & a) const {
        return C(r*a.r - i*a.i, r*a.i + i*a.r);
    }
};

typedef long long LL;
const double pi = acos(-1.0);
const int N = 50005;
C a[N<<2], b[N<<2];
char num1[N], num2[N];
LL ret[N<<2];

void brc(C *y, int L) {
    int i, j, k;
    for (i=1,j=L>>1; i<L-1; ++i) { // 二進制平攤反轉置換 O(NlogN)
        if (i < j) swap(y[i], y[j]);
        k = L>>1;
        while (j >= k) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
}

void FFT(C *y, int L, int dir) {
    brc(y, L);
    for (int h = 2; h <= L; h <<= 1) { // 枚舉所需計算的點數 
        C wn(cos(dir*2*pi/h), sin(dir*2*pi/h)); // h次單位復根 
        for (int j = 0; j < L; j += h) { // 原序列被分成了L/h段h長序列 
            C w(1, 0); // 旋轉因子 
            for (int k = j; k < j+h/2; ++k) { // 因為折半定理,只需要計算枚舉一半的長度即可 
                C u = y[k];
                C t = w*y[k+h/2];
                y[k] = u + t;
                y[k+h/2] = u - t;
                w = w * wn; // 更新旋轉因子 
            }
        }
    }
    if (dir == 1) {
        for (int i = 0; i < L; ++i) {
            y[i] = y[i] * C(1.0/L, 0);
        }
    }
}

int main() {
    while (scanf("%s %s", num1, num2) != EOF) {
        memset(ret, 0, sizeof (ret));
        int len1 = strlen(num1), len2 = strlen(num2);
        int ML = len1+len2-1, L = 1;
        while (L < ML) L <<= 1;
        for (int i = len1-1, j = 0; i >= 0; --i, ++j) {
            a[j] = C(num1[i]-'0', 0);
        }
        for (int i = len2-1, j = 0; i >= 0; --i, ++j) {
            b[j] = C(num2[i]-'0', 0);
        }
        for (int i = len1; i < L; ++i) a[i] = C(0, 0);
        for (int i = len2; i < L; ++i) b[i] = C(0, 0);
        FFT(a, L, -1), FFT(b, L, -1);
        for (int i = 0; i < L; ++i) {
            a[i] = a[i] * b[i];
        }
        FFT(a, L, 1);
        for (int i = 0; i < L; ++i) {
            ret[i] = (LL)floor(a[i].r + 0.5);
        }
        for (int i = 0; i < L; ++i) {
            ret[i+1] += ret[i] / 10;
            ret[i] %= 10;
        }
        int p = L;
        while (!ret[p] && p) --p;
        while (p >= 0) printf("%d", (int)ret[p--]);
        puts("");
    }
    return 0;
} 
View Code

 

HDU-4609 3-idiots

題意:有N條線段,問從這N條線段中選出三條能過組成三角形的概率為多大?

分析:三條邊能夠組成三角形則滿足等式x+y < z。首先將所有的邊排序,然后將每條邊的長度構成多項式的指數,同一長度的邊的數量為系數。然后將這個多項式自己和自己作一個乘法,這里需要使用FFT來實現,去掉自己與自己的組合已經相互的組合情況,就能夠得到兩兩之間組合形成邊長和值為某一個值得方案數。使用這個方案數除以總方案數即可。

#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

const int N = 300005;

typedef long long LL;
struct C {
    double r, i;
    C() {}
    C(double _r, double _i) : r(_r), i(_i) {}
    inline C operator + (const C &a) const {
        return C(r + a.r, i + a.i);
    }
    inline C operator - (const C &a) const {
        return C(r - a.r, i - a.i);
    }
    inline C operator * (const C &a) const {
        return C(r*a.r-i*a.i, r*a.i+i*a.r);
    }
}a[N], b[N];

const double pi = acos(-1.0);
int n, num[N], cnt[N];
LL res[N], sum[N];

void brc(C *y, int l) {
    int i, j, k;
    for (i=1,j=l>>1; i<l-1; ++i) {
        if (i < j) swap(y[i], y[j]);
        k = l>>1;
        while (j >= k) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
}

void FFT(C *y, int l, int on) {
    int h, i, j, k;
    C u, t;
    brc(y, l); // 得到一個自底向上的序列 
    for (h = 2; h <= l; h <<= 1) { // 控制一個O(logn)的外層復雜度 
        C wn(cos(on*2*pi/h), sin(on*2*pi/h));
        for (j = 0; j < l; j+=h) { // 兩個for循環共組成O(n)的復雜度 
            C w(1, 0);
            for (k = j; k <j+h/2; ++k) {
                u = y[k];
                t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if (on == 1) {
        for (i = 0; i < l; ++i) {
            y[i]= y[i] * C(1.0/l, 0.0);
        }
    }
}

int main() {
    int T;
    scanf("%d", &T);
    while (T--) {
        scanf("%d", &n);
        int Max = 0;
        memset(cnt, 0, sizeof (cnt));
        for (int i = 0; i < n; ++i) {
            scanf("%d", &num[i]);
            Max = max(Max, num[i]);
            ++cnt[num[i]];
        }
        int L = 1;
        ++Max;
        while (L < (Max<<1)) L <<= 1;
        for (int i = 0; i < Max; ++i) {
            a[i] = C(cnt[i], 0);
        }
        for (int i = Max; i < L; ++i) {
            a[i] = C(0, 0);
        }
        FFT(a, L, -1);
        for (int i = 0; i < L; ++i) {
            a[i] = a[i] * a[i];
        }
        FFT(a, L, 1);
        for (int i = 0; i < L; ++i) {
            res[i] = (LL)floor(a[i].r + 0.5);
        }
        for (int i = 0; i < Max; ++i) {
            res[i<<1] -= cnt[i];
        }
        for (int i = 0; i < L; ++i) {
            res[i] >>= 1;
        }
        for (int i = 1; i < L; ++i) {
            sum[i] = sum[i-1] + res[i];
        }
        double ret = 0, den = 1.0*n*(n-1)*(n-2)/6.0;
        for (int i = 0; i < n; ++i) {
            ret += sum[num[i]] / den;
        }
        printf("%.7f\n", 1-ret);
    }
    return 0;
}
View Code

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM