HDU6128 二次剩余/二次域求二次剩余解/LL快速乘法取模

LINK

题意：求满足模p下$frac{1}{a_i+a_j}equivfrac{1}{a_i}+frac{1}{a_j}$的对数，其中$n,p(1leq nleq10^5,2leq pleq10^{18})$

思路：推式子，两边同乘$(a_i + a_j)^3$,得$a_i^2+a_j^2 equiv {a_i·a_j} mod{p}$,进一步$a_i^2+a_j^2+a_i·a_jequiv {0} mod{p}$，然后？然后会点初中数竞,或者数感好会因式分解就能看出来，两边再乘个$a_i-a_j$就是$a_i^3-a_j^3$了，直接得到了$a_i$和$a_j$的关系,可喜可贺。然而显然我这种蠢得的人看不出来的，一般性的我们使用一元二次求根公式，视$a_j$为常数,得在模p下的解$a_i=a_j·frac{-1 pm sqrt{-3}}{2}$，那么得到了$a_i$和$a_j$的关系，但是有$sqrt{-3}$的存在，也就是要想办法求到$x^2 equiv{-3} mod{p}$的解，来替代原式中的根号，这个头次见到还是浙大校赛的ZOJ3774，当时还奇怪为什么有人能直接看出xx数就是$sqrt{5}$的二次剩余。

这里用到二次剩余的相关知识，前置知识是看二次剩余有关的教学ppt(欧拉判别，勒让德符号之类的)，求二次剩余的解有种二次域的求法（参见wiki上的Cipolla's algorithm），还有这位菊苣的博客，acdream菊苣的博客上也有模板和一些介绍......

因为半桶水，证明或算法在此不表。

其他要注意给出数据有可能为0，且要进行欧拉判别模p下是否存在二次剩余为-3的解，以及p=3的情况，还有就是因为LL下乘法溢出的问题，注意使用O（1）的$2^{64}$LL的取模乘法。

/** @Date    : 2017-08-17 20:07:47
  * @FileName: 1009.cpp
  * @Platform: Windows
  * @Author  : Lweleth (SoungEarlf@gmail.com)
  * @Link    : https://github.com/
  * @Version : $Id$
  */
#include <bits/stdc++.h>
#define LL long long
#define PII pair<int ,int>
#define MP(x, y) make_pair((x),(y))
#define fi first
#define se second
#define PB(x) push_back((x))
#define MMG(x) memset((x), -1,sizeof(x))
#define MMF(x) memset((x),0,sizeof(x))
#define MMI(x) memset((x), INF, sizeof(x))
using namespace std;

const int INF = 0x3f3f3f3f;
const int N = 1e5+20;
const double eps = 1e-8;

LL mod;

struct T 
{  
    LL p, d;  
};  
LL w;

//O1乘法取模黑科技
LL mul(LL x,LL y)
{
    return (x * y-(LL)(x /(long double)mod * y + 1e-3) * mod + mod) % mod;
}

//二次域乘法
T multi_er(T a, T b)  
{  
    T ans;  
    ans.p = (mul(a.p, b.p) + mul(mul(a.d,b.d), w)) % mod;  
    ans.d = (mul(a.p, b.d) + mul(a.d, b.p)) % mod;  
    return ans;  
}  



LL quick_mod(LL a, LL b)
{  
    LL ans = 1;  
    a %= mod;  
    while(b)  
    {  
        if(b & 1)  
        {  
            ans = mul(ans , a);  
            b--;  
        }  
        b >>= 1;  
        a = mul(a , a);  
    }  
    return ans;  
}  
  
//二次域上快速幂  
T power(T a, LL b)  
{  
    T ans;  
    ans.p = 1;  
    ans.d = 0;  
    while(b)  
    {  
        if(b & 1)  
        {  
            ans = multi_er(ans, a);  
            b--;  
        }  
        b >>= 1;  
        a = multi_er(a, a);  
    }  
    return ans;  
}  
  
//求勒让德符号  
LL Legendre(LL a, LL p)  
{  
    return quick_mod(a, (p-1)>>1);  
}  

LL QuadraticResidue()
{
    LL rem = (-3 % mod + mod) % mod;
    if(rem == 0)//特判mod==3
        return 0;
    if(mod == 2)//特判非奇素数
        return 1;
    if(Legendre(rem, mod) + 1 == mod)//欧拉判别条件 非剩余
        return -1;
    LL b;
    while(1)//找一个非剩余求二次域上的单位w=sqrt(b^2 - rem)
    {
        b = rand() % mod;
        w = (mul(b, b) - rem + mod) % mod;
        if(quick_mod(w, (mod - 1)/2) + 1 == mod)//cipolla
            break;
    }
    T tmp;
    tmp.p = b;
    tmp.d = 1;
    T ans = power(tmp, (mod + 1) / 2);
    return ans.p;
}


vector<LL>a;
int main()
{
    int T;
    cin >> T;
    while(T--)
    {
        a.clear();
        LL n, p;
        scanf("%lld%lld", &n, &mod);
        for(int i = 0; i < n; i++)
        {
            LL t;
            scanf("%lld", &t);
            if(t > 0)//注意有0的...
                a.PB(t);
        }
        LL cnt = a.size();
        sort(a.begin(), a.end());

        ///////////
        LL ans = 0;
        if(mod == 2)//特殊情况无解
            ans = cnt * (cnt - 1) / 2LL;
        else 
        {
            LL t = QuadraticResidue();
            if(t == -1)
            {
                printf("0
");
                continue;
            }
            LL inv = (mod + 1) >> 1;
            LL x = mul((mod + t - 1LL)%mod, inv);
            LL y = mul((mod - t - 1LL)%mod, inv);
            for(int i = 0; i < cnt; i++)
            {
                LL tmp = mul(x , a[i]);
                ans += upper_bound(a.begin(), a.begin() + i, tmp) 
                - lower_bound(a.begin(), a.begin() + i, tmp);
            }    
            if(x != y)//两个解
            {
                for(int i = 0; i < cnt; i++)
                {
                    LL tmp = mul(y, a[i]);
                    ans += upper_bound(a.begin(), a.begin() + i, tmp) 
                    - lower_bound(a.begin(), a.begin() + i, tmp);
                }
            }
        }
        printf("%lld
", ans);
    }
    return 0;
}
/*
99
5 7
1 2 3 4 5
6 7
1 2 3 4 5 6
*/