Problem: F - Point Distance

Problem Statement

You work for invention center as a part time programmer. This center researches movement of protein molecules. It needs how molecules make clusters, so it will calculate distance of all pair molecules and will make a histogram. Molecules’ positions are given by a N x N grid map. Value C_{xy} of cell (x, y) in a grid map means that C_{xy} molecules are in the position (x, y).

You are given a grid map, please calculate histogram of all pair molecules.

Input

Input is formatted as follows.


N
C_{11} C_{12} ... C_{1N}
C_{21} C_{22} ... C_{2N}
...
C_{N1} C_{N2} ... C_{NN}


First line contains the grid map size N (1 \leq N \leq 1024). Each next N line contains N numbers. Each numbers C_{xy} (0 \leq C_{xy} \leq 9) means the number of molecule in position (x,y). There are at least 2 molecules.

Output

Output is formatted as follows.


D_{ave}
d_{1} c_{1}
d_{2} c_{2}
...
d_{m} c_{m}


Print D_{ave} which an average distance of all pair molecules on first line. Next, print histogram of all pair distance. Each line contains a distance d_i and the number of pair molecules c_i(0 \lt c_i). The distance d_i should be calculated by squared Euclidean distance and sorted as increasing order. If the number of different distance is more than 10,000, please show the first 10,000 lines. The answer may be printed with an arbitrary number of decimal digits, but may not contain an absolute or relative error greater than or equal to 10^{-8}.

Sample Input 1


2
1 0
0 1


Output for the Sample Input 1


1.4142135624
2 1


Sample Input 2


3
1 1 1
1 1 1
1 1 1


Output for the Sample Input 2


1.6349751825
1 12
2 8
4 6
5 8
8 2


Sample Input 3


5
0 1 2 3 4
5 6 7 8 9
1 2 3 2 1
0 0 0 0 0
0 0 0 0 1


Output for the Sample Input 3


1.8589994382
0 125
1 379
2 232
4 186
5 200
8 27
9 111
10 98
13 21
16 50
17 37
18 6
20 7
25 6






Example

Input

2
1 0
0 1


Output

1.4142135624
2 1

Solution (in CPP):

Solution:
#include<bits/stdc++.h>
typedef long long int ll;
typedef unsigned long long int ull;
#define BIG_NUM 2000000000
#define HUGE_NUM 99999999999999999
#define MOD 1000000007
#define EPS 0.000000001
using namespace std;




#define MAX 25
typedef complex<double> COMPLEX;
typedef vector<COMPLEX> V;
typedef vector<V> MATRIX;


ll ans[3000005];

struct Info{
	Info(int arg_left,int arg_right){
		left = arg_left;
		right = arg_right;
	}
	int left,right;
};

int POW[MAX];
int LR_rev[3000005];
vector<Info> info;
vector<COMPLEX> MULT[MAX];


void calc_swap_pair(int N,int num_pow){

	LR_rev[0] = 0;

	for(int i = 1; i < N; i++){

		if(i%2 == 1){

			LR_rev[i] += LR_rev[i-1]+POW[num_pow-1]; //1を足す代わりにPOW[num_pow-1]を足す

		}else{

			LR_rev[i] = LR_rev[i/2]/2; //半分のインデックスの値に、2を掛ける代わりに2で割る
		}

		if(i < LR_rev[i]){
			info.push_back(Info(i,LR_rev[i]));
		}
	}
}

void calc_COMPLEX(int N){

	for(int LEN = 2,num_pow = 0; LEN <= N; LEN *= 2,num_pow++){

		MULT[num_pow].resize(LEN/2);

		MULT[num_pow][0] = COMPLEX(cos(0),sin(0));

		COMPLEX mult = COMPLEX(cos((2*M_PI)/LEN),sin((2*M_PI)/LEN));

		for(int i = 1; i < LEN/2; i++){

			if(i%2 == 1){
				MULT[num_pow][i] = MULT[num_pow][i-1]*mult;
			}else{
				MULT[num_pow][i] = MULT[num_pow-1][i/2];
			}
		}
	}
}

//離散フーリエ変換
vector<COMPLEX> DFT(vector<COMPLEX> A) {

	int N = A.size();

	/*偶数項を左に、奇数項を右に仕分けるルールで要素数1まで仕分けした場合の、
	 * 最終結果を先に作っておく(ビットが左右反転したものがその位置に来る)*/
	for(int i = 0; i < info.size(); i++){

		swap(A[info[i].left],A[info[i].right]);
	}

	COMPLEX a,add;

	for(int LEN = 2,num_pow = 0; LEN <= N; LEN *= 2,num_pow++){

		for(int start = 0; start < N; start += LEN){
			for(int i = 0; i < LEN/2; i++){
				//マージソートの容量で、配列を更新していく
				a = A[start+i];
				add = MULT[num_pow][i]*A[start+LEN/2+i];

				A[start+i] = a+add;
				A[start+LEN/2+i] = a-add;
			}
		}
	}
	return A;
}

//離散逆フーリエ変換
vector<COMPLEX> inverseDFT(vector<COMPLEX> A){

	int N = A.size();

	vector<COMPLEX> TMP = DFT(A);

	//係数を入れ替える&サイズNで各要素を割る
	vector<COMPLEX> RET(N);

	for(int i = 0; i < N; i++){

		RET[i] = TMP[(N-i)%N];
		RET[i] = RET[i]/COMPLEX(N,0);
	}

	return RET;
}

//AとBのサイズは、大きい方に揃えておく
MATRIX convolution2D(MATRIX A,MATRIX B){

	int SIZE = A.size();

	for(int i = 0; i < MAX; i++){

		if(SIZE == POW[i]){
			calc_swap_pair(SIZE,i);
					break;
		}
	}
	calc_COMPLEX(SIZE);

	MATRIX RET(SIZE,V(SIZE));

	//行についてのフーリエ変換
	for(int row = 0; row < SIZE; row++){

		A[row] = DFT(A[row]);
		B[row] = DFT(B[row]);
	}

	//列についてのフーリエ変換を行単位で処理するため、転置
	for(int col = 1; col < SIZE; col++){
		for(int row = 0; row < col; row++){
			swap(A[row][col],A[col][row]);
			swap(B[row][col],B[col][row]);
		}
	}

	//列についてのフーリエ変換
	for(int row = 0; row < SIZE; row++){

		A[row] = DFT(A[row]);
		B[row] = DFT(B[row]);
	}

	//元に戻すため、再転置
	for(int col = 1; col < SIZE; col++){
		for(int row = 0; row < col; row++){
			swap(A[row][col],A[col][row]);
			swap(B[row][col],B[col][row]);
		}
	}

	//ここまでの結果を計算
	for(int row = 0; row < SIZE; row++){
		for(int col = 0; col < SIZE; col++){

			RET[row][col] = A[row][col]*B[row][col];
		}
	}

	//各行をフーリエ逆変換
	for(int row = 0; row < SIZE; row++){

		RET[row] = inverseDFT(RET[row]);
	}

	//転置
	for(int col = 1; col < SIZE; col++){
		for(int row = 0; row < col; row++){
			swap(RET[row][col],RET[col][row]);
		}
	}

	for(int row = 0; row < SIZE; row++){

		RET[row] = inverseDFT(RET[row]);
	}

	//元に戻すため、再転置
	for(int col = 1; col < SIZE; col++){
		for(int row = 0; row < col; row++){
			swap(RET[row][col],RET[col][row]);
		}
	}

	return RET;
}

int main(){

	POW[0] = 1;
	for(int i = 1; i < MAX; i++){

		POW[i] = POW[i-1]*2;
	}

	ll N;
	scanf("%lld",&N);

	ll SIZE = 1;
	while(SIZE < 2*N)SIZE *= 2;

	MATRIX A(SIZE,V(SIZE));
	MATRIX B(SIZE,V(SIZE)); //上下左右反転

	ll tmp;
	ll first_sum = 0;

	for(int row = 0; row < N; row++){
		for(int col = 0; col < N; col++){

			scanf("%lld",&tmp);
			first_sum += tmp;

			A[row][col] += tmp;
			B[(N-1)-row][(N-1)-col] += tmp;
		}
	}

	MATRIX ret = convolution2D(A,B);

	ll maximum = (N-1)*(N-1)*2;

	for(ll i = 0; i <= maximum; i++){

		ans[i] = 0;
	}

	double sum = 0;
	ll all_num = 0;

	for(ll row = 0; row <= 2*N-2; row++){
		for(ll col = 0; col <= 2*N-2; col++){

			ll dist = ((N-1)-row)*((N-1)-row)+((N-1)-col)*((N-1)-col);
			ll num = (ll)round(ret[row][col].real());

			all_num += num;

			sum += num*sqrt((double)dist);

			ans[dist] += num;
		}
		printf("\n");
	}

	all_num -= first_sum;

	printf("%.10lf\n",sum/all_num);

	int count = 0;

	ans[0] -= first_sum;

	for(ll i = 0; i <= maximum; i++){

		if(ans[i] == 0)continue;

		printf("%lld %lld\n",i,ans[i]/2);

		count++;
		if(count == 10000)break;
	}

	return 0;
}