HDU 6128 : Inverse of sum

题目大意：有$n$ 个小于质数$p$ 的非负整数$a_i$ ，你想知道有多少对$i,j(1 \leqslant i<j\leqslant n)$ ，使得模$p$ 意义下$\frac{1}{a_i+a_j} \equiv \frac{1}{a_i}+\frac{1}{a_j}$ ，即这两数的和的逆元等于它们逆元的和，注意零元没有逆元。

二次剩余

$$
(a+b)^{-1} \equiv a^{-1}+b^{-1}(mod , p)\\
1 \equiv (a+b)(a^{-1}+b^{-1})(mod , p)\\
a^{-1}b+ab^{-1}+1 \equiv 0(mod , p)\\
a^2+b^2+ab\equiv 0 (mod , p)\\
(ab^{-1})^2+(ab^{-1})+1\equiv 0 (mod , p)\\
ab^{-1}\equiv \frac{-1\pm\sqrt{-3}}{2}(mod , p)
$$

故原式等价于模$p$ 意义下，两数之比为$\frac{-1\pm\sqrt{-3}}{2}$ ，根号可以用二次剩余解决。

下面分类考虑几种情况：

• $p=2$ ，此时可以直接求得$-3$ 的二次剩余是$1$
• $(\frac{-3}{p})=1$ ，套用上述公式遍历一遍数组可求得
• $(\frac{-3}{p})=-1$ ，即$-3$ 在模$p$ 下不存在二次剩余，故不存在满足条件的对
• $(\frac{-3}{p})=0$ ，即$p|-3$ ，此时$\sqrt{-3} \equiv 0 (mod , p)$

需要注意的是，此题会出现long long相乘溢出的情况，可以用快速乘处理。

代码如下：

#include <iostream>
#include <cstdlib>
#include <ctime>
#include <map>
using namespace std;
typedef long long ll;
const int N=100005;
ll T,n,p,w,a[N];
map<ll,ll>mp;
ll Min(ll a,ll b){
    return a<b?a:b;
}
ll Mul(ll a,ll b){
    ll r=0;
    while(b){
        if(b&1)r=(r+a)%p;
        a=(a+a)%p;
        b>>=1;
    }
    return r;
}
struct Complex{
    ll x,y;
    Complex friend operator*(Complex x,Complex y){
        Complex z;
        z.x=(Mul(x.x,y.x)+Mul(Mul(x.y,y.y),w))%p;
        z.y=(Mul(x.y,y.x)+Mul(x.x,y.y))%p;
        return z;
    }
};
ll powmod(ll a,ll n){
    ll r=1;
    while(n){
        if(n&1)r=Mul(r,a);
        a=Mul(a,a);
        n>>=1;
    }
    return r;
}
ll ComPow(Complex a,ll n){
    Complex r;
    r.x=1;r.y=0;
    while(n){
        if(n&1)r=r*a;
        a=a*a;
        n>>=1;
    }
    return Min(r.x%p,(p-r.x%p)%p);
}
ll Cipolla(ll a,ll p){
    ll x=rand()%p;
    while(powmod((Mul(x,x)-a+p)%p,(p-1)/2)!=p-1)
        x=rand()%p;
    w=(Mul(x,x)-a+p)%p;
    Complex t;
    t.x=x;t.y=1;
    return ComPow(t,(p+1)/2);
}
ll solve(ll rd,ll p){
    mp.clear();
    ll ans=0;
    ll x=Mul(((-1+rd)%p+p)%p,powmod(2,p-2));
    ll y=Mul(((-1-rd)%p+p)%p,powmod(2,p-2));
    if(x==y){
        for(int i=n-1;i>=0;--i)if(a[i]%p){
            ll b=Mul(a[i],x);
            if(mp.count(b))ans+=mp[b];
            mp[a[i]%p]++;
        }
    }else{
        for(int i=n-1;i>=0;--i)if(a[i]%p){
            ll bx=Mul(a[i],x),by=Mul(a[i],y);
            if(mp.count(bx))ans+=mp[bx];
            if(mp.count(by))ans+=mp[by];
            mp[a[i]%p]++;
        }
    }
    return ans;
}
int main(void){
    cin>>T;
    while(T--){
        cin>>n>>p;
        for(int i=0;i<n;++i)
            cin>>a[i];
        ll ans;
        if(p==2){
            ans=solve(1,p);
        }else if(powmod(p-3,(p-1)/2)==1){
            ll rd=Cipolla(p-3,p);
            ans=solve(rd,p);
        }else if(powmod(p-3,(p-1)/2)==0){
            ans=solve(0,p);
        }else ans=0;
        cout<<ans<<"\n";
    }
}

立方公式

除了用二次剩余还有另一种更优雅的解法。
$$
(a+b)^{-1} \equiv a^{-1}+b^{-1}(mod , p)\\
1 \equiv (a+b)(a^{-1}+b^{-1})(mod , p)\\
a^{-1}b+ab^{-1}+1 \equiv 0(mod , p)\\
a^2+b^2+ab\equiv 0 (mod , p)\\
$$
当$a-b \not \equiv 0 (mod , p)$ 时，可以在上述等式两边同乘上$a-b$ ，得到
$$
(a-b)(a^2+ab+b^2)\equiv 0 (mod ,p)\\
a^3-b^3\equiv0(mod , p)
$$
如此只需记录当$a \not\equiv b(mod ,p)$ 时，所有满足$a^3-b^3\equiv 0 (mod ,p)$ 的对数即可。

查看代码 。