Language: C++

#include<vector>
#include<cmath>
#include<map>
#include<cstdlib>
#include<iostream>
#include<sstream>
#include<fstream>
#include<string>
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<set>
#include<stack>
#include<bitset>
#include<functional>
#include<ctime>
#include<queue>
#include<deque>
#include<complex>
using namespace std;
#define pb push_back
#define pf push_front
typedef long long lint;
typedef complex<double> P;
#define mp make_pair
#define fi first
#define se second
typedef pair<int,int> pint;
#define All(s) s.begin(),s.end()
#define rAll(s) s.rbegin(),s.rend()
#define REP(i,a,b) for(int i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)
lint mo=1000000007;
lint ma[4100][4100];
int A[252521],B[252521];
lint extgcd(lint a, lint b, lint &x, lint &y) {
  lint g = a; x = 1; y = 0;
  if (b != 0) g = extgcd(b, a % b, y, x), y -= (a / b) * x;
  return g;
}
lint invMod(lint a, lint m) {
  lint x, y;
  if (extgcd(a, m, x, y) == 1) return (x + m) % m;return 0;
}

#define N 4096
const double pi = 4.0*atan(1.0);
const P I(0, 1);
vector<P> fft(int n, double theta, vector<P> a) {
  for (int m = n; m >= 2; m >>= 1) {
    int mh = m >> 1;
    for (int i = 0; i < mh; i++) {
      P w = exp(i*theta*I);
      for (int j = i; j < n; j += m) {
        int k = j + mh;
        P 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]);
  }
	return a;
}
vector<vector<P> > fft2(int n,double theta, vector<vector<P> > a){
	for (int m = n; m >= 2; m >>= 1) {
    int mh = m >> 1;
    for (int i = 0; i < mh; i++) {
      P w = exp(i*theta*I);
      for (int j = i; j < n; j += m) {
        int k = j + mh;
        vector<P> x;
      	rep(l,n) x.pb(a[j][l] - a[k][l]);
      	rep(l,n) a[j][l] += a[k][l];
      	rep(l,n) a[k][l] = w * x[l];
      }
    }
    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]);
  }
	return a;
}
int main()
{
	int n,a,b;lint out=0;
	memset(ma,0,sizeof(ma));
	scanf("%d",&n);
	rep(i,n){
		scanf("%d %d",&A[i],&B[i]);ma[A[i]][B[i]]++;
	}
	vector<P> vec(4100,P(0,0));
	vector<vector<P> > mat(4100,vec);
	rep(i,2050) rep(j,2050) mat[i][j].real()+=P(ma[i][j],0);
	//cout<<a<<endl;
	rep(i,4100) mat[i]=fft(N,2*pi/N,mat[i]);
	mat=fft2(N,2*pi/N,mat);
	rep(i,N) rep(j,N) mat[i][j]*=mat[i][j];
	mat=fft2(N,-2*pi/N,mat);
	rep(i,4100) mat[i]=fft(N,-2*pi/N,mat[i]);
	//cout<<a<<endl;
	
	vector<lint> zyo,rz;
	zyo.pb(1);
	rep(i,N) zyo.pb((zyo[i]*(i+1))%mo);
	rep(i,N) rz.pb(invMod(zyo[i],mo));
	rep(i,4098) rep(j,4098){
		lint x=(lint)(mat[i][j].real()+0.01);
		x/=N*N;
		ma[i][j]=x;
	}
	rep(i,n) ma[A[i]*2][B[i]*2]--;
	rep(i,4098) rep(j,4098){
		lint x=ma[i][j]/2;
		if(x>0){
			//cout<<i<<' '<<j<<' '<<x<<endl;
			out+=((x*zyo[i+j]%mo)*rz[i])%mo*rz[j];
			out%=mo;
		}
	}
	cout<<out<<endl;
}


Compiler output:
a.cc: In function 'int main()':
a.cc:102:49: error: no match for 'operator+=' (operand types are 'double' and 'P' {aka 'std::complex<double>'})
  102 |         rep(i,2050) rep(j,2050) mat[i][j].real()+=P(ma[i][j],0);
      |                                 ~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~