【BNU】39676 Point Distance【FFT求矩阵中点对的欧几里德距离】

传送门：【BNU】39676 Point Distance

题目分析：

纠结了好久，终于搞定了……

首先我们定义：
Dx,y=∑n−1i=0∑n−1j=0Ci,jCi+x,j+y

Dx,y表示两个点对之间行距为x，列距为y的对数。

然后我们将点映射到一维，得到表达式：
Dxn+y=∑n−1i=0∑n−1j=0Cin+jC(i+x)n+j+y
⇒   Dxn+y=∑n−1i=0∑n−1j=0Cin+jCin+j+(xn+y)
⇒   DA=∑n−1i=0∑n−1j=0CBCB+A

转化成多项式表达则为：

DA=∑n−1i=0∑n−1j=0CBx−B⋅CB+AxB+A

这样我们就可以用FFT优化了。

然而我们直接这么做得到的一个DA可能有多个值。原因是，出现了跳行的数据，避免的方法就是每一行元素个数翻倍，使得前一半或者后一半的元素为空，剩下的保持不变，这样FFT过程中的非法值就不会有了。（具体自行体会哈~）

上面的方法处理完以后，每个DA就对应一个确定的pair对：（0，0），（x，y）。然后这个（x，y）可能是非法的，需要注意一下。

一天就只做了这么一道水题，没救了，还有好多都不会做T T……

my  code:

#include <stdio.h>#include <string.h>#include <math.h>#include <algorithm>using namespace std ;typedef long long LL ;#define clr( a , x ) memset ( a , x , sizeof a )#define cpy( a , x ) memcpy ( a , x , sizeof a )const int MAXN = 1030 ;const int MAXM = 1030 * 1030 * 4 ;const double pi = acos ( -1.0 ) ;struct Complex {    double r , i ;    Complex () {}    Complex ( double r , double 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 ) ;    }} ;int n , m ;Complex x1[MAXM] , x2[MAXM] ;LL num[MAXM] ;LL ans[MAXM] ;void FFT ( Complex y[] , int n , int rev ) {    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 ( cos ( rev * 2 * pi / s ) , sin ( rev * 2 * pi / s ) ) ;        for ( int k = 0 ; k < n ; k += s ) {            Complex w ( 1 , 0 ) , t ;            for ( int i = k ; i < k + ds ; ++ i , w = w * wn ) {                y[i + ds] = y[i] - ( t = w * y[i + ds] ) ;                y[i] = y[i] + t ;            }        }    }    if ( rev < 0 ) {        for ( int i = 0 ; i < n ; ++ i ) {            y[i].r /= n ;        }    }}int dist ( int x , int y ) {    return x * x + y * y ;}void solve () {    LL tot = 0 ;    double ave = 0 ;    int x , n1 = 1 ;    clr ( num , 0 ) ;    clr ( ans , 0 ) ;    m = 2 * n * n ;    while ( n1 < 2 * m ) n1 <<= 1 ;    for ( int i = 0 ; i < n ; ++ i ) {        for ( int j = 0 ; j < n ; ++ j ) {            scanf ( "%d" , &x ) ;            num[i * 2 * n + j] = x ;            tot += x ;            ans[0] += x * ( x - 1 ) / 2 ;        }    }    /*for ( int i = 0 ; i < m ; ++ i ) {        printf ( "%d " , num[i] ) ;    }    printf ( "\n" ) ;*/    for ( int i = 0 ; i < m ; ++ i ) {        x1[i] = Complex ( num[i] , 0 ) ;        x2[i] = Complex ( num[m - i - 1] , 0 ) ;    }    for ( int i = m ; i < n1 ; ++ i ) {        x1[i] = x2[i] = Complex ( 0 , 0 ) ;    }    FFT ( x1 , n1 , 1 ) ;    FFT ( x2 , n1 , 1 ) ;    for ( int i = 0 ; i < n1 ; ++ i ) {        x1[i] = x1[i] * x2[i] ;    }    FFT ( x1 , n1 , -1 ) ;    /*for ( int i = 0 ; i < n1 ; ++ i ) {        printf ( "%d " , ( int ) ( x1[i].r + 0.5 ) ) ;    }    printf ( "\n" ) ;*/    for ( int i = m ; i < n1 ; ++ i ) {        num[i] = ( int ) ( x1[i].r + 0.5 ) ;        int x1 = 0 , y1 = 0 ;        int x2 = ( i - m + 1 ) / ( 2 * n ) ;        int y2 = ( i - m + 1 ) % ( 2 * n ) ;        if ( y2 >= n ) {            int tmp = n * 2 - y2 ;            y1 += tmp ;            x2 ++ ;            y2 = 0 ;        }        int dis = dist ( x1 - x2 , y1 - y2 ) ;        ans[dis] += num[i] ;        ave += sqrt ( 1.0 * dis ) * num[i] ;    }    tot = tot * ( tot - 1 ) / 2 ;    printf ( "%.10f\n" , ave / tot ) ;    int cnt = 0 ;    for ( int i = 0 ; i < MAXM ; ++ i ) if ( ans[i] ) {        printf ( "%d %lld\n" , i , ans[i] ) ;        ++ cnt ;        if ( cnt >= 10000 ) break ;    }}int main () {    while ( ~scanf ( "%d" , &n ) ) solve () ;    return 0 ;}