---
title: Tsr and Variance
tags:
  - Fast Fourier Transform
id: "3604"
categories:
  - - OI
    - Number Theory
date: 2020-01-31 12:00:00
---

## 题目描述

Tsr is a cute boy with handsome moustache.

You are given a sequence with length $n$. Tsr wants you to calculate the sum of variance of each successive subsequence. Note: The variance in this problem should't be divided by length.

Recall $\overline a_{l, r}={1 \over r-l+1} \sum_{i=l}^r a_i$. Then you are supposed to calculate $\sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline a_{i, j})^2$.

## 题意概述

给定一个长度为$n$的序列，求$\sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline a_{i, j})^2$。有$T$组数据。

数据范围：$1 \le T \le 20, \ 1 \le n \le 10000, \ 1 \le a_i \le 10$。

## 算法分析

$$\begin{align} \sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k-\overline a_{i, j})^2 &= \sum_{i=1}^n \sum_{j=i}^n \sum_{k=i}^j (a_k^2-2a_k \overline a_{i, j}+\overline a_{i, j}^2) \\\\ &= \sum_{k=1}^n \sum_{i=1}^k \sum_{j=k}^n a_k^2-2\sum_{i=1}^n \sum_{j=i}^n \overline a_{i, j} \sum_{k=i}^j a_k+ \sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline a_{i, j}^2 \\\\ &= \sum_{k=1}^n k(n-k+1)a_k^2-\sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline a_{i, j}^2 \end{align}$$

前一项可以在$O(n)$时间内求出，只需考虑如何求后一项。令$s_i=\sum_{j=1}^i a_j$。

$$\begin{align} \sum_{i=1}^n \sum_{j=i}^n (j-i+1) \overline a_{i, j}^2  &= \sum_{j=1}^n \sum_{i=1}^j {(s_j-s_{i-1})^2 \over j-i+1} \\\\ &= \sum_{j=1}^n \sum_{i=1}^j {s_j^2-2s_js_{i-1}+s_{i-1}^2 \over j-i+1} \\\\ &= \sum_{j=1}^n s_j^2 \sum_{i=1}^j {1 \over i}-2\sum_{j=1}^n s_j \sum_{i=1}^j {s_{i-1} \over j-(i-1)}+\sum_{j=1}^n \sum_{i=1}^j {s_{i-1}^2 \over j-(i-1)} \end{align}$$

第一项也可以在$O(n)$时间内求出，而后面两项都是卷积的形式，可以用 FFT 在$O(n\log n)$时间内求出。

## 代码

```cpp
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>

int const N = 10005, T = 400005;
double const PI = acos(-1);

int a[N], rev[T];

struct Point {
    double x, y;
    Point(double _x = 0, double _y = 0) : x(_x), y(_y) {}
    friend Point const operator + (Point const &a, Point const &b) {
        return Point(a.x + b.x, a.y + b.y);
    }
    friend Point const operator - (Point const &a, Point const &b) {
        return Point(a.x - b.x, a.y - b.y);
    }
    friend Point const operator * (Point const &a, Point const &b) {
        return Point(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
    }
    friend Point const operator / (Point const &a, double const &n) {
        return Point(a.x / n, a.y / n);
    }
} wn[T], A[T], B[T], C[T];

int init(int n) {
    int m = n, l = 0;
    for (n = 1; n <= m; n <<= 1, ++l) ;
    for (int i = 1; i < n; ++i) {
        rev[i] = rev[i >> 1] >> 1 | (i & 1) << l - 1;
    }
    for (int i = 0; i < n >> 1; ++i) {
        wn[i] = Point(cos(2 * PI / n * i), sin(2 * PI / n * i));
    }
    return n;
}

void fft(Point *a, int n, int inv = 0) {
    for (int i = 0; i < n; ++i) {
        if (i < rev[i]) {
            std::swap(a[i], a[rev[i]]);
        }
    }
    for (int i = 1; i < n; i <<= 1) {
        for (int j = 0; j < n; j += i << 1) {
            for (int k = 0; k < i; ++k) {
                Point x = a[j + k], y = wn[n / (i << 1) * k] * a[j + k + i];
                a[j + k] = x + y;
                a[j + k + i] = x - y;
            }
        }
    }
    if (inv) {
        std::reverse(a + 1, a + n);
        for (int i = 0; i < n; ++i) {
            a[i] = a[i] / n;
        }
    }
}

int main() {
    int T;
    scanf("%d", &T);
    for (; T--;) {
        int n;
        scanf("%d", &n);
        double ans = 0;
        for (int i = 1; i <= n; ++i) {
            scanf("%d", &a[i]);
            ans += 1. * a[i] * a[i] * i * (n - i + 1);
            a[i] += a[i - 1];
        }
        double sum = 0;
        for (int i = 1; i <= n; ++i) {
            sum += 1. / i;
            ans -= sum * a[i] * a[i];
        }
        int len = init(n << 1);
        for (int i = 0; i < n; ++i) {
            A[i] = Point(a[i]);
            B[i] = Point(1. * a[i] * a[i]);
            C[i] = Point(1. / (i + 1));
        }
        for (int i = n; i < len; ++i) {
            A[i] = B[i] = C[i] = Point(0);
        }
        fft(A, len);
        fft(B, len);
        fft(C, len);
        for (int i = 0; i < len; ++i) {
            A[i] = A[i] * C[i];
            B[i] = B[i] * C[i];
        }
        fft(A, len, 1);
        fft(B, len, 1);
        for (int i = 0; i < n; ++i) {
            ans += 2 * a[i + 1] * A[i].x - B[i].x;
        }
        printf("%.10f\n", ans);
    }
    return 0;
}
```