Streaming_4_noip_day2 距离统计 （本原解，勾股数，欧拉筛）

数据
题意：给出$n\ast m$的点阵，有T组数据，对于每一组数据给出一个长度l，有多少对点阵上的点能构成长度为l的线段。
$n,m\le 10^9,T\le 10^3,l\in [1,2max(n,m)]$
思路：只需要求解有多少对正整数a,b满足$a^2+b^2=l^2$。该方程的解至少是一个本原解的k倍$(k\ge 1)$。则只需枚举l的因子，然后利用本原解的性质判断该因子是否满足本原解。
n的因子个数为$O(lnn)$，知乎大佬证明
枚举因子的方法：先用欧拉筛筛出$1\sim \sqrt{n}$内的所有素数（$1\sim 5\cdot 10^4$中素数不超过6000个）然后就可以在O(素数个数)下做出l的标准分解式。然后利用dfs即可求出l的所有因子了。
再利用$z=r^2+s^2$的性质，在$O(\sqrt{n})$下完成r，s的枚举，（判断一个数是否为完全平方数，可以直接用std自带的sqrt函数计算，再带回检验）并且r，s满足$(r,s)=1且2\nmid (r+s)$
则总复杂度为$O(6000+lnn\cdot \sqrt{n})$

const int N = 1e5 + 10;
int prim[6000], tot;
bool vis[N];
void euler(int n) {//欧拉筛
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) prim[++tot] = i;
        for (int j = 1; j <= tot; j++) {
            int t = i * prim[j];
            if (t > n) break;
            vis[t] = 1;
            if (i % prim[j] == 0) break;
        }
    }
}
int n, m, T, l;
int p[N], alp[N];//p[]存储标准分解的素数,alp为该素数的次幂数
LL ans;
int gcd(int a, int b) {return (!b) ? a : gcd(b, a % b);}
void calc(int ing) {
    int z = l / ing;
    for (int i = 1; i * i < z; i++) {
        int a = i;
        int b = z - a * a;
        LL tmp = sqrt(b);
        if (tmp * tmp != b) continue;
        b = tmp;
        if (b < a) swap(a, b);
        if (a == 0 || b == 0 || a == b || (!((a + b) & 1)) || gcd(a, b) > 1) continue;//判断a,b是否满足上述两条性质，并且保证x,y为正整数
        LL x = (b * b - a * a) * ing;
        LL y = 2 * a * b * ing;
        if (x < n && y < m) ans += (LL)(n - x) * (m - y);
        if (x < m && y < n) ans += (LL)(n - y) * (m - x);
    }
}

void dfs(int now, int id) {//枚举因数
    if (id > p[0]) {
        calc(now);
        return;
    }
    int pri = p[id];
    dfs(now, id + 1);
    for (int i = 1; i <= alp[id]; i++)  {
        now *= pri;
        dfs(now, id + 1);
    }
}
int main() {
    freopen("dist.in", "r", stdin);
    freopen("dist.out", "w", stdout);
    euler(5e4);
    read(n), read(m), read(T);
    while (T--) {
        read(l);
        p[0] = 0; ans = 0;
        if (l < n) ans += (LL)(n - l) * m;//计算x=l,y=0
        if (l < m) ans += (LL)(m - l) * n;//计算x=0,y=l
        int x = l;
        for (int i = 1; i <= tot; i++) {//对l进行标准分解
            if (prim[i] * prim[i] > x) break;
            alp[i] = 0;
            if (x % prim[i] == 0) p[++p[0]] = prim[i];
            while (x % prim[i] == 0) {
                ++alp[p[0]];
                x /= prim[i];
            }
        }
        if (x != 1) p[++p[0]] = x, alp[p[0]] = 1;//注意最后x也可能为素数
        dfs(1, 1);
        printf("%lld ", ans);
    }
    return 0;
}