SRM 335 Div1Hard

機械的に問題を典型で処理するというのはかなり楽だと思うんですよね。

概要
n項ある整数列が2つある。(n<=50000)各項の数は0~49998.
持ち点は初め0である。
以下のことをn回繰り返す。

2つある数列からそれぞれ一回も選ばれていないものをランダムに一つ選ぶ。
1つめの数列から選ばれたものをA,2つめの数列から選ばれたものをBとする。
A=Bのとき、持ち点に(A-B)^2を足す。

最終状態での期待値を求めよ。


解法
まず、各項に対して独立で考えることはできるので、O(n^2)は自明。
さて、オーダーを落とそう。
このオーダーを落とすところにおいて、想定解法ではそれっぽい思考をしているが、機械的に典型に持ち込める。
上の「」の中にある、A-B=iとなるものの個数を数えることが出来ればよい。
ここで、以下のスライドを思い出す。

http://acm-icpc.aitea.net/index.php?plugin=attach&refer=2012%2FPractice%2F%E6%98%A5%E3%82%B3%E3%83%B3%E3%83%86%E3%82%B9%E3%83%88%2F%E8%AC%9B%E8%A9%95&openfile=Point%20Distance.pdf

上のものの完全に劣化版じゃん…

FFTを書く。以上。

注意:
doubleの桁落ち誤差に注意しましょう。long doubleか、整数の処理で済ませるところはlong longが無難です。(それでずっとはまっていました)

// I like wolves!!
#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <stack>
#include <cassert>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <queue>
#include <complex>
#include <string.h>
using namespace std;
const long double PI = 4.0*atan(1.0);
const long double EPS = 1e-6;
typedef complex<long double> Complex;
const Complex I(0, 1);

void fft(int n,long double theta, Complex a[]) {
  for (int m = n; m >= 2; m >>= 1) {
    int mh = m >> 1;
    for (int i = 0; i < mh; i++) {
      Complex w = exp(i*theta*I);
      for (int j = i; j < n; j += m) {
        int k = j + mh;
        Complex x = a[j] - a[k];
        a[j] += a[k];
        a[k] = w * x;
      }
    }
    theta *= 2;
  }
  int i = 0;
  for (int j = 1; j < n - 1; j++) {
    for (int k = n >> 1; k > (i ^= k); k >>= 1);
    if (j < i) swap(a[i], a[j]);
  }
}
class RandomFights {
	public:
	int a[131072];
	int b[131072];
	Complex ax[131072];
	Complex bx[131072];
	Complex cx[131072];
	double expectedNrOfPoints(vector <int> A, vector <int> B, int n) {
		int at=0;
		int sz=131072;
		for(int i=0;i<sz;i++){
			a[i]=b[i]=0;
			ax[i]=bx[i]=cx[i]=Complex(0.0,0.0);
		}
		int m=A.size();
		for(int i=0;i<n;i++){
			int val=A[at];
			a[val]++;
			int s=(at+1)%m;
			A[at]=((A[at]^A[s])+13)%49999;
			at=s;
		}
		at=0;
		m=B.size();
		for(int i=0;i<n;i++){
			int val=B[at];
			b[65535-val]++;
			int s=(at+1)%m;
			B[at]=((B[at]^B[s])+13)%49999;
			at=s;
		}
		//for(int i=0;i<65536;i++)if(a[i]!=b[65535-i])printf("%d: %d %d\n",i,a[i],b[65535-1]);
		for(int i=0;i<sz;i++){
			ax[i]=Complex((long double)a[i],0);
			bx[i]=Complex((long double)b[i],0);
		}
		/*for(int i=0;i<sz;i++){
			ax[i]/=10.0;
			bx[i]/=10.0;
		}*/
		fft(sz,PI*2/sz,ax);
		fft(sz,PI*2/sz,bx);
		/*for(int i=0;i<sz;i++){
			ax[i]*=10.0;
			bx[i]*=10.0;
		}*/
		for(int i=0;i<sz;i++){
			cx[i]=ax[i]*bx[i];
		}
		fft(sz,-PI*2/sz,cx);
		//for(int i=0;i<sz;i++)cx[i]/=sz;
		long long ret=0;  //こいつをdoubleにして数時間ハマってた
		long double gosa=0;
		for(int i=0;i<sz-1;i++){
			int v=65535-i;
			long double val=cx[i].real()/sz;
			long long q=(long long)(val+EPS);
			assert((abs(val-q)<0.2||val+EPS>0));
		//	gosa=max(gosa,abs(val-q));
		//	if(q)printf("%d: %lld\n",i,q);
		//	if(abs(val-q)>0.1)printf("%d: %f\n",i,val);
			if(v>0){
				ret-=q*v*v;
			}else{
				ret+=q*v*v;
			}
		}
		//printf("%.12Lf\n",gosa);
		return (double)ret/n;
	}
};

//Powered by [KawigiEdit] 2.0!