部分引用自：http://blog.csdn.net/v5zsq/article/details/77255048

所以假設方程 x^2+x+1=0 在模p意義下的解為d，則答案就是滿足(ai/aj) mod p = d的數對(i,j)的數量(i<j)。

現在把問題轉化為解這個模意義下的二次方程。

x^2+x+1=0

配方：x^2+x+1/4+3/4=0

(x+1/2)^2+3/4=0

同乘4：(2x+1)^2+3=0

即(2x+1)^2=-3 (mod p)

換句話說，我們必須保證-3+p是p的一個二次剩余。

倘若-3+p不是p的一個二次剩余的話，無解。

如果是的話，我們就可以通過Cipolla算法得到2x+1的兩個可能值，也就把二次同余方程變成了兩個線性同余方程。

然後在一般情況下可以通過擴展歐幾裏得算法求得這兩個個線性同余方程的最小非負整數解，也就得到了這個二次同余方程的兩個解。

不過這題的情況有點特別，形式非常簡單，根據大神的總結，在我們得到2x+1的兩個值d₁、d₂之後，可以這樣得到x1、x2。

於是我們可以通過枚舉+map來得到點對的個數。

註意特殊情況：①倘若x1或者x2為零的話，意味著(ai/aj)mod p=0，這是不可能的，註意不要統計這種。

②p為2時無解。

③p為3時有解，但是由於p-3等於零，所以不能直接用上面的方法。不過經過簡單的推導，我們發現唯一合法的情況是ai=aj,且ai、aj均不為零，直接特判掉。

求解二次剩余的Cipolla算法：

http://blog.csdn.net/a_crazy_czy/article/details/51959546

http://blog.csdn.net/philipsweng/article/details/50000903

一句話：通過勒讓德符號來判斷n是不是模p的二次剩余，然後random一個a，使得(a^2-n)不是二次剩余，然後通過復數快速冪來求二次剩余的值。

#include<cstdio>
#include<map>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;
ll p,a[200005];
struct Complex{
	ll a,b;
	Complex(const ll &a,const ll &b){
		this->a=a;
		this->b=b;
	}
	Complex(){}
};
map<ll,int>ma;
int T,n;
ll Quick_Mul(ll x,ll p,ll mod){
	if(!p){
		return 0ll;
	}
	ll res=Quick_Mul(x,p>>1,mod);
	res=(res+res)%mod;
	if((p&1ll)==1ll){
		res=(x%mod+res)%mod;
	}
	return res;
}
ll Quick_Pow(ll x,ll p,ll mod){
	if(!p){
		return 1ll;
	}
	ll res=Quick_Pow(x,p>>1,mod);
	res=Quick_Mul(res,res,mod);
	if((p&1ll)==1ll){
		res=Quick_Mul(x%mod,res,mod);
	}
	return res;
}
ll aa,nn;
Complex Complex_Mul(const Complex &A,const Complex &B,const ll &p){
	return Complex((Quick_Mul(A.a,B.a,p)+Quick_Mul(Quick_Mul(A.b,B.b,p),(Quick_Mul(aa,aa,p)-nn+p)%p,p))%p,
	(Quick_Mul(A.a,B.b,p)+Quick_Mul(A.b,B.a,p))%p);
}
Complex Complex_Quick_Pow(Complex x,ll p,ll mod){
	if(!p){
		return Complex(1ll,0ll);
	}
	Complex res=Complex_Quick_Pow(x,p>>1,mod);
	res=Complex_Mul(res,res,mod);
	if((p&1ll)==1ll){
		res=Complex_Mul(x,res,mod);
	}
	return res;
}
ll ran1(){
	return ((rand()<<16)|rand());
}
ll ran(){
	return ((ran1()<<16)|ran1());
}
pair<ll,ll> work(ll n,ll p){
	aa=ran()%p;
	nn=n;
	while(Quick_Pow((Quick_Mul(aa,aa,p)-n+p)%p,(p-1ll)/2ll,p)!=p-1ll){
		aa=ran()%p;
	}
	ll res=Complex_Quick_Pow(Complex(aa,1ll),(p+1ll)/2ll,p).a;
	return make_pair(res,p-res);
}
ll ans;
ll f(ll x,ll p){
	return (x&1ll)==1ll ? (x-1ll)/2ll : (x-1ll+p)/2ll;
}
int main(){
//	freopen("1009.in","r",stdin);
//	freopen("hdu6128.out","w",stdout);
	srand(233);
	scanf("%d",&T);
	for(;T;--T){
		ans=0;
		ma.clear();
		scanf("%d%lld",&n,&p);
		for(int i=1;i<=n;++i){
			scanf("%lld",&a[i]);
		}
		if(p==2ll || Quick_Pow(p-3ll,(p-1ll)/2ll,p)==p-1ll){
			puts("0");
			continue;
		}
		if(p==3ll){
			for(int i=1;i<=n;++i){
				if(a[i]){
					ans+=ma[a[i]];
				}
				++ma[a[i]];
			}
			printf("%lld\n",ans);
			continue;
		}
		pair<ll,ll> d=work(p-3ll,p);
		d.first=f(d.first,p);
		d.second=f(d.second,p);
		for(int i=1;i<=n;++i){
			if(d.first){
				ans+=ma[Quick_Mul(a[i],d.first%p,p)];
			}
			if(d.second){
				ans+=ma[Quick_Mul(a[i],d.second%p,p)];
			}
			if(a[i]){
				++ma[a[i]];
			}
		}
		printf("%lld\n",ans);
	}
	return 0;
}

【數論】【二次剩余】【map】hdu6128 Inverse of sum