51nod1687 方差求和

题目

51nod1687 方差求和
题目链接

解题思路

简化公式
∑ i = 1 n ( a i − a ˉ ) 2 = ∑ i = 1 n ( a i 2 − 2 a i a ˉ + a ˉ 2 ) = ∑ i = 1 n a i 2 − 2 a ˉ ∑ i = 1 n a i + n a ˉ 2 = ∑ i = 1 n a i 2 − 2 n a ˉ 2 + n a ˉ 2 = ∑ i = 1 n a i 2 − n a ˉ 2 \sum_{i=1}^n(a_i-\bar a)^2\\=\sum_{i=1}^n(a_i^2-2a_i\bar a+\bar a^2)\\=\sum_{i=1}^na_i^2-2\bar a\sum_{i=1}^nai+n\bar a^2\\=\sum_{i=1}^na_i^2-2n\bar a^2+n\bar a^2\\=\sum_{i=1}^na_i^2-n\bar a^2 ∑i=1n​(ai​−aˉ)2=∑i=1n​(ai2​−2ai​aˉ+aˉ2)=∑i=1n​ai2​−2aˉ∑i=1n​ai+naˉ2=∑i=1n​ai2​−2naˉ2+naˉ2=∑i=1n​ai2​−naˉ2
考虑增量法
n a ˉ 2 = ( ∑ i = 1 n a i ) 2 n n\bar a^2=\frac{(\sum_{i=1}^na_i)^2}{n} naˉ2=n(∑i=1n​ai​)2​
( x + a ) 2 n + 1 − x 2 n = n x 2 + 2 n x a + n a 2 − ( n + 1 ) x 2 n ( n + 1 ) = 2 n a x + n a 2 − x 2 n ( n + 1 ) = 2 a x + a 2 n + 1 − x 2 n ( n + 1 ) \frac{(x + a)^2}{n+1}-\frac{x^2}{n}\\=\frac{nx^2+2nxa+na^2-(n+1)x^2}{n(n+1)}\\=\frac{2nax+na^2-x^2}{n(n+1)}\\=\frac{2ax+a^2}{n+1}-\frac{x^2}{n(n+1)} n+1(x+a)2​−nx2​=n(n+1)nx2+2nxa+na2−(n+1)x2​=n(n+1)2nax+na2−x2​=n+12ax+a2​−n(n+1)x2​
无法摆脱与长度的关系,不可行。
考虑枚举长度,同时计算长度相同的区间
我们知道区间长度和左端点L就可以计算右端点R了。
a n s = ∑ l e n = 1 n ∑ L = 1 n − l e n + 1 ( ∑ i = L R a i 2 − l e n a ˉ 2 ) ans=\sum_{len=1}^n\sum_{L=1}^{n-len+1}(\sum_{i=L}^Ra_i^2-len\bar a^2) ans=∑len=1n​∑L=1n−len+1​(∑i=LR​ai2​−lenaˉ2)
其中 ∑ l e n = 1 n ∑ L = 1 n − l e n + 1 ∑ i = L R a i 2 = ∑ i = 1 n a i 2 ∗ i ∗ ( n − i + 1 ) \sum_{len=1}^n\sum_{L=1}^{n-len+1}\sum_{i=L}^Ra_i^2=\sum_{i=1}^na_i^2*i*(n-i+1) ∑len=1n​∑L=1n−len+1​∑i=LR​ai2​=∑i=1n​ai2​∗i∗(n−i+1)可以简单求出。
∑ l e n = 1 n ∑ L = 1 n − l e n + 1 l e n a ˉ 2 \sum_{len=1}^n\sum_{L=1}^{n-len+1}len\bar a^2 ∑len=1n​∑L=1n−len+1​lenaˉ2
= ∑ l e n = 1 n ∑ L = 1 n − l e n + 1 ( ∑ i = L R a i ) 2 l e n =\sum_{len=1}^n\frac{\sum_{L=1}^{n-len+1}(\sum_{i=L}^Ra_i)^2}{len} =∑len=1n​len∑L=1n−len+1​(∑i=LR​ai​)2​
= ∑ l e n = 1 n ∑ L = 1 n − l e n + 1 ( ∑ i = L R a i 2 + 2 ∑ i = L R ∑ j = i + 1 R a i a j ) l e n =\sum_{len=1}^n\frac{\sum_{L=1}^{n-len+1}(\sum_{i=L}^Ra_i^2+2\sum_{i=L}^R\sum_{j=i+1}^Ra_ia_j)}{len} =∑len=1n​len∑L=1n−len+1​(∑i=LR​ai2​+2∑i=LR​∑j=i+1R​ai​aj​)​
其中 ∑ l e n = 1 n ∑ L = 1 n − l e n + 1 ∑ i = L R a i 2 l e n \sum_{len=1}^n\frac{\sum_{L=1}^{n-len+1}\sum_{i=L}^Ra_i^2}{len} ∑len=1n​len∑L=1n−len+1​∑i=LR​ai2​​也能够简单维护。
令 f ( l e n ) = ∑ L = 1 n − l e n + 1 ∑ i = L R ∑ j = i + 1 R a i a j f(len)=\sum_{L=1}^{n-len+1}\sum_{i=L}^R\sum_{j=i+1}^Ra_ia_j f(len)=∑L=1n−len+1​∑i=LR​∑j=i+1R​ai​aj​
我们发现 f ( l e n ) f(len) f(len)的二阶差分可以分为两部分,一部分可以用NTT维护,一部分可以简单维护。
整体复杂度 O ( T n l o g 2 n ) O(Tnlog_2n) O(Tnlog2​n)。

代码

#include <cstdio>
#include <algorithm>
#include <cstring>
using namespace std;

typedef long long ll;

void read(int &x) {
	x = 0; char c = getchar();
	while (c < '0' || c > '9') c = getchar();
	while (c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
}

const int N = 3e5 + 100;

namespace NTT {
#define pw(n) (1<<n)
	const int N = 262144, P = 998244353, g = 3;//或P=1004535809
	int n, m, bit, bitnum = 0, a[N + 5], b[N + 5], rev[N + 5];
	void getrev(int l) {
		for (int i = 0; i < pw(l); i++) {
			rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
		}
	}
	int fastpow(int a, int b) {
		int ans = 1;
		for (; b; b >>= 1, a = 1LL * a*a%P) {
			if (b & 1)ans = 1LL * ans*a%P;
		}
		return ans;
	}
	void NTT(int *s, int op) {
		for (int i = 0; i < bit; i++)if (i < rev[i])swap(s[i], s[rev[i]]);
		for (int i = 1; i < bit; i <<= 1) {
			int w = fastpow(g, (P - 1) / (i << 1));
			for (int p = i << 1, j = 0; j < bit; j += p) {
				int wk = 1;
				for (int k = j; k < i + j; k++, wk = 1LL * wk*w%P) {
					int x = s[k], y = 1LL * s[k + i] * wk%P;
					s[k] = (x + y) % P;
					s[k + i] = (x - y + P) % P;
				}
			}
		}
		if (op == -1) {
			reverse(s + 1, s + bit);
			int inv = fastpow(bit, P - 2);
			for (int i = 0; i < bit; i++)s[i] = 1LL * s[i] * inv%P;
		}
	}
	int solve(int *aa, int nn, int *bb, int mm, int *c) {
		n = nn; m = mm;
		bit = bitnum = 0;
		for (int i = 0; i <= n; i++) a[i] = aa[i];
		for (int i = 0; i <= m; i++) b[i] = bb[i];
		m += n;
		for (bit = 1; bit <= m; bit <<= 1)bitnum++;
		getrev(bitnum);
		NTT(a, 1);
		NTT(b, 1);
		for (int i = 0; i < bit; i++) a[i] = 1LL * a[i] * b[i] % P;
		NTT(a, -1);
		for (int i = 0; i < bit; i++) c[i] = a[i];
		for (int i = 0; i < bit; i++) a[i] = b[i] = 0;
		return bit;
	}
}

int n;
int sa[N], sum[N], pre[N], suf[N], aa[N], bb[N];

int main() {
	//freopen("0.txt", "r", stdin);
	int T; read(T);
	while (T--) {
		read(n);
		for (int i = 1; i <= n; i++) read(sa[i]);
		double ans = 0;
		for (int i = 1; i <= n; i++) {
			ans += 1LL * i * (n - i + 1) * sa[i] * sa[i];
			sum[i] = sum[i - 1] + sa[i] * sa[i];
			pre[i] = pre[i - 1] + sa[i];
		}
		for (int i = n; i >= 1; i--) suf[i] = suf[i + 1] + sa[i];
		ll r = sum[n];
		for (int len = 1; len <= n; len++) {
			ans -= r / double(len);
			int L = len + 1, R = n - len;
			if (L <= r) r += sum[R] - sum[L - 1];
			else r -= sum[L] - sum[R - 1];
		}
		for (int i = 1; i <= n; i++) {
			aa[i] = sa[i];
			bb[i] = sa[n - i + 1];
		}
		NTT::solve(aa, n + 1, bb, n + 1, aa);
		ll r1 = 0, r2 = 0;
		for (int len = 2; len <= n; len++) {
			r1 += aa[n - len + 2] - sa[n - len + 2] * suf[n - len + 3] - sa[len - 1] * pre[len - 2];
			r2 += r1;
			ans -= 2 * r2 / double(len);
		}
		printf("%.2lf\n", ans);
		for (int i = 0; i <= n + 1; i++) sum[i] = pre[i] = suf[i] = 0;
	}
	return 0;
}
上一篇:python中的运算符


下一篇:能源DEA--弱处置与强处置、自然处置与管理处置(disposability)