HDU 5307 He is Flying【FFT】

题意：
有n段路，每段路长度为$s_i$，你从节点i到节点j，可以获得一个开心值$j-i+1$，然后问你，主人公走过了所有总长度为$s$的段，问你有多少开心值。
思路：
首先想到的就是类似于组合数学那样，于是就想到了母函数：

$$(\sum ix^{\sum s_i})(\sum x^{-\sum s_{i-1}})-(\sum x^{\sum s_i})(\sum (i-1)x^{-\sum s_{i-1}})$$
预处理出前缀和，用FFT做。
当然要预处理出$s=0$的结果，用$O(n)$扫描一遍就行了。
官方的题解说，FFT可能会被卡常数，但是用了long double之后就可以顺利过了
另外，负数的处理方法有几种，一种是加上一个偏移量，另一种是将整个复数数组翻转过来。

//      whn6325689
//      Mr.Phoebe
//      http://blog.youkuaiyun.com/u013007900
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <cstring>
#include <climits>
#include <complex>
#include <fstream>
#include <cassert>
#include <cstdio>
#include <bitset>
#include <vector>
#include <deque>
#include <queue>
#include <stack>
#include <ctime>
#include <set>
#include <map>
#include <cmath>
#include <functional>
#include <numeric>
#pragma comment(linker, "/STACK:1024000000,1024000000")

using namespace std;

#define eps 1e-9
#define PI acos(-1.0)
#define INF 0x3f3f3f3f
#define LLINF 1LL<<62
#define speed std::ios::sync_with_stdio(false);

typedef long long ll;
typedef long double ld;
typedef pair<ll, ll> pll;
typedef complex<ld> point;
typedef pair<int, int> pii;
typedef pair<pii, int> piii;
typedef vector<int> vi;

#define CLR(x,y) memset(x,y,sizeof(x))
#define CPY(x,y) memcpy(x,y,sizeof(x))
#define clr(a,x,size) memset(a,x,sizeof(a[0])*(size))
#define cpy(a,x,size) memcpy(a,x,sizeof(a[0])*(size))

#define mp(x,y) make_pair(x,y)
#define pb(x) push_back(x)
#define lowbit(x) (x&(-x))

#define MID(x,y) (x+((y-x)>>1))
#define ls (idx<<1)
#define rs (idx<<1|1)
#define lson ls,l,mid
#define rson rs,mid+1,r
#define root 1,1,n

template<class T>
inline bool read(T &n)
{
    T x = 0, tmp = 1;
    char c = getchar();
    while((c < '0' || c > '9') && c != '-' && c != EOF) c = getchar();
    if(c == EOF) return false;
    if(c == '-') c = getchar(), tmp = -1;
    while(c >= '0' && c <= '9') x *= 10, x += (c - '0'),c = getchar();
    n = x*tmp;
    return true;
}
template <class T>
inline void write(T n)
{
    if(n < 0)
    {
        putchar('-');
        n = -n;
    }
    int len = 0,data[20];
    while(n)
    {
        data[len++] = n%10;
        n /= 10;
    }
    if(!len) data[len++] = 0;
    while(len--) putchar(data[len]+48);
}
//-----------------------------------

const int MAXN=50010;
const ld pi=acosl(-1.0);

struct Complex
{
    ld r,i;
    Complex(){}
    Complex(ld r ,ld i):r(r),i(i) {}
    Complex operator + (const Complex& t) const
    {
        return Complex(r+t.r,i+t.i) ;
    }
    Complex operator - (const Complex& t) const
    {
        return Complex(r-t.r,i-t.i);
    }
    Complex operator * (const Complex& t) const
    {
        return Complex(r*t.r-i*t.i,r*t.i+i*t.r);
    }
}las[MAXN*5],pre[MAXN*5],c[MAXN*5];

void FFT(Complex y[],int n,int rev)//rev=-1表示逆变换
{
    for(int i=1,j,k,t; i<n; i++)  //进行蝶型变换
    {
        for(j=0,k=n>>1,t=i; k; k>>=1,t>>=1) j=j<<1|t&1;
        if(i<j ) swap(y[i],y[j]);
    }
    for ( int s = 2 , ds = 1 ; s <= n ; ds = s , s <<= 1 )
    {
        Complex wn=Complex(cosl(rev*2*pi/s),sinl(rev*2*pi/s)),w=Complex(1,0),t;
        for(int k=0; k<ds; k++,w=w*wn)
        {
            for(int i=k; i<n; i+=s)
            {
                y[i+ds]=y[i]-(t=w*y[i+ds]);
                y[i]=y[i]+t;
            }
        }
    }
    if(rev==-1) for(int i=0; i<n; i++) y[i].r/=n;
}

int n;
ll val[MAXN*5],ans[MAXN*5];
ll pu[MAXN*5];

int main()
{
//    freopen("data.txt","r",stdin);
//    freopen("out.txt","w",stdout);
    for(int i=1;i<=MAXN*3;i++)
        pu[i]=pu[i-1]+1LL*i*(i+1)/2;
    int T,temp;
    read(T);
    while(T--)
    {
        read(n);temp=0;
        int L=1;CLR(ans,0);
        for(int i=1;i<=n;i++)
        {
            read(val[i]);
            val[i]+=val[i-1];
        }

        for (int q=1,h=0;q<=n;q=h+1)
            if(val[q]!=val[q-1])    h=q;
            else
            {
                for(h=q;h<n && val[h+1]==val[q];h++);
                    ans[0]+=pu[h-q+1];
            }
        while(L<=val[n]) L<<=1;
        L<<=1;
        CLR(las,0);CLR(pre,0);CLR(c,0);
        for(int i=1;i<=n;i++)   las[val[i]].r+=i;
        for(int i=1;i<=n;i++)   pre[-val[i-1]+val[n]].r+=1.0;
        FFT(las,L,1);FFT(pre,L,1);
        for(int i=0;i<L;i++)
            c[i]=las[i]*pre[i];
        FFT(c,L,-1);
        for(int i=1;i<L;i++)
            ans[i]+=ll(c[i].r+0.1);
        CLR(las,0);CLR(pre,0);CLR(c,0);
        for(int i=1;i<=n;i++)   las[val[i]].r+=1.0;
        for(int i=1;i<=n;i++)   pre[-val[i-1]+val[n]].r+=i-1;
        FFT(las,L,1);FFT(pre,L,1);
        for(int i=0;i<L;i++)
            c[i]=las[i]*pre[i];
        FFT(c,L,-1);
        for(int i=1;i<L;i++)
            ans[i]-=ll(c[i].r+0.1);
        printf("%lld\n",ans[0]);
        for(int i=1;i<val[n]+1;i++)
            printf("%lld\n",ans[i+val[n]]);
    }
    return 0;
}

下面是官方的题解
用的是两个$2^{31}$以内的模数分别FFT，最后用中国剩余定理合并。

#include <cstdio>
#include <string.h>
using namespace std;

#define P 50000000001507329LL
#define G 3
int T;
int n, s[110000];
long long ans[140000], a[140000], b[140000], c[140000], x[140000], w[140000];
long long pu[110000];
int nn;

long long Mul(long long x, long long y) {
    return (x*y-(long long)(x/(long double)P*y+1e-3)*P+P)%P;
}

long long Pow(long long x, long long y) {
    long long i, ans = 1;
    for (i = 1; i <= y; i *= 2, x = Mul(x, x)) if (y & i)  ans = Mul(ans, x);
    return ans;
}

void DFT(long long *a, int n) {
    int m, i, d, p, q;
    for (m = 1; (1 << m) <= n; m++){
        for (i = 0; i < (n >> m); i++)
        for (q = 0, d = p = i; d < n; q += (n >> m), d += (n >> (m - 1)), p += (n >> m)){
            x[p] = (Mul(a[d + (n >> m)], w[q]) + a[d]) % P;
            x[p + n / 2] = (Mul(a[d + (n >> m)], w[q + n / 2]) + a[d]) % P;
        }
        for (i = 0; i < n; i++) a[i] = x[i];
    }
}

void DFT1(long long *a, int n){
    int m, i, d, p, q;
    for (m = 1; (1 << m) <= n; m++) {
        for (i = 0; i < (n >> m); i++)
        for (q = 0, d = p = i; d < n; q += (n >> m), d += (n >> (m - 1)), p += (n >> m)){
            x[p] = (Mul(a[d + (n >> m)], w[n - q]) + a[d]) % P;
            x[p + n / 2] = (Mul(a[d + (n >> m)], w[n / 2 - q]) + a[d]) % P;
        }
        for (i = 0; i < n; i++) a[i] = x[i];
    }
}

void doit() {
    long long S = Pow(n, P - 2);
    DFT(a, n);
    DFT(b, n);
    for (int i = 0; i < n; i++)
        c[i] = Mul(a[i], b[i]);
    DFT1(c, n);
    for (int i = 0; i < n; i++)
        c[i] = Mul(c[i], S);
    for (int i = 0; i < n; i++)
        ans[i] = (ans[i] + c[i]) % P;
    // for (int i = 0; i <= s[n]; i++)
    //  for (int j = 0; j <= s[n]; j++)
    //      ans[i + j] += 1LL * a[i] * b[j];
}

int main() {
    // freopen("a.in", "r", stdin);
    // freopen("std.out", "w", stdout);
    scanf("%d", &T);
    n = 131072;
    for (int i = 0; i <= n; i++)
        w[i] = Pow(G, (P - 1) / n * i);
    while (T--) {
        scanf("%d", &nn);
        for (int i = 1; i <= nn; i++)
            scanf("%d", &s[i]), s[i] += s[i - 1];
        memset(ans, 0, sizeof ans);
        memset(a, 0, sizeof a);
        memset(b, 0, sizeof b);

        for (int i = 1; i <= nn; i++)
            a[s[i]] += i;
        for (int i = 1; i <= nn; i++)
            b[s[nn] - s[i - 1]]++;
        doit();
        memset(a, 0, sizeof a);
        memset(b, 0, sizeof b);
        for (int i = 0; i <= 50000; i++)
            a[i] = b[i] = 0;
        for (int i = 1; i <= nn; i++)
            a[s[i]] = (a[s[i]] + P - 1) % P;
        for (int i = 1; i <= nn; i++)
            b[s[nn] - s[i - 1]] += i - 1;
        doit();
        long long ans0 = 0;
        int q, h;
        for (int i = 1; i <= nn; i++)
            pu[i] = pu[i - 1] + 1LL * i * (i + 1) / 2;
        for (q = 1; q <= nn; q = h + 1)
            if (s[q] != s[q - 1]) h = q;
            else {
                for (h = q; h < nn && s[h + 1] == s[q]; h++);
                ans0 += pu[h - q + 1];
            }
        printf("%lld\n", ans0);
        for (int i = 1; i <= s[nn]; i++)
            printf("%lld\n", ans[i + s[nn]]);
    }
}