LOJ #6052. 「雅礼集训 2017 Day11」DIV

定义复数 $ a + bi$ 为整数 $k$ 的约数,当且仅当 $a$ 和 $b$ 为整数且存在整数 $c$ 和 $d$ 满足 $ (a + bi)(c + di) = k $,给定 $n $,求出 $1$ 到 $n$ 的所有满足 $a > 0 $ 的约数 $ a + bi $ 的 $a$ 的和。答案模 $ 1004535809 $ 输出。

首先容易发现,对于整数 $n$ 来说,如果存在$(a + bi)\times (c + di) = $n,那么一定有

$$\begin{equation}\begin{cases}{} ac - bd = n,\newline ad + bc = 0 \end{cases}\end{equation}$$

即 $\frac{a}{c} = \frac{b}{-d}$,令

$$\frac{a}{c} = \frac{b}{-d} = \frac{p}{q},(p, q) = 1, a = xp, c = xq, b = yp, -d = yq$$

得到 $xy(p^2 + q^2) = n,(p, q) = 1$

以上可用的条件基本用完了,我们考虑列出答案的求和式。

下面如果不作特殊说明,一切除法都是下取整的。

$$ans = \sum\limits_{x = 1}^n \sum\limits_{(p, q) = 1,(p^2 + q^2) | x} p f(\frac{x}{p^2 + q^2})$$

其中,$f(x)$代表$x$的约数和。

发现其中比较强的限制是$p, q$ 的,我们把它提到前面,再令$F(t) = \sum\limits_{i = 1}^t f(i)$,得:

$ans = \sum\limits_{(p, q) = 1, (p^2 + q^2) \leq n} pF(\frac{n}{p^2 + q^2})$

虽然化简到了这一步,但是容易发现仍然不好计算,我们考虑枚举$(p^2 + q^2)$,

同时令$g(x) = \sum\limits_{(p^2 + q^2) = x,(p, q) = 1} p$

$ans = \sum\limits_{x} g(x) F(\frac{n}{x})$

$\frac{n}{x}$的取值不超过$\sqrt {n}$个,如果能够快速的求出$G(x) = \sum\limits_{i = 1}^x g(x)$,就能够在可以接受的时间内得到答案。

因为$(p, q) = 1$ 的限制比较强,所以我们考虑去掉这个限制,令$H(x) = \sum\limits_{p^2 + q^2 \leq x} p$

观察一下可以发现$H(x) = \sum\limits_{d = 1}^x d\times G(\frac{x}{d^2})$

得到 $G(x) = H(x) - \sum\limits_{d = 2}^{\sqrt x} G(\frac{x}{d^2}) x$

利用杜教筛在 $O(n^{\frac{2}{3}} \log n)$ 得到所有的$G(x)$,同时对于$F(x)$,也利用杜教筛的思想优化到 $O(n^{\frac{2}{3}})$。

我们发现最后得到的答案仍然不对,因为我们上面计算的是 $c > 0$,根据对称性我们把答案乘二即可,对于$c = 0$ 的,直接加上$F(n)$即可。

于是我们在$O(n^{\frac{2}{3}} \log n )$ 的时间内得到了本题的答案。

#include <bits/stdc++.h>

using namespace std;

#define ll             long long
#define db            double
#define up(i,j,n)        for (int i = j; i <= n; i++)
#define down(i,j,n)    for (int i = j; i >= n; i--)
#define cadd(a,b)        a = add (a, b)
#define cpop(a,b)        a = pop (a, b)
#define cmul(a,b)        a = mul (a, b)
#define pr            pair<int, int>
#define fi            first
#define se            second
#define SZ(x)        (int)x.size()
#define bin(i)        (1 << (i))
#define Auto(i,node)    for (int i = LINK[node]; i; i = e[i].next)

template<typename T> inline bool cmax(T & x, T y){return y > x ? x = y, true : false;}
template<typename T> inline bool cmin(T & x, T y){return y < x ? x = y, true : false;}

const int mod = 1004535809;
const int inv2 = (mod + 1) / 2;
const int MAXN = 5e6 + 5;
const int L = 5e6;

inline int add(int a, int b){a += b; return a >= mod ? a - mod : a;}
inline int pop(int a, int b){a -= b; return a < 0 ? a + mod : a;}
inline int mul(int a, int b){return 1LL * a * b % mod;}

int gcd(int a, int b) { return !b ? a : gcd(b, a % b); }

ll n;

int dd[MAXN][2], ff[MAXN][2];
bool chk[MAXN][2];

int &D(ll x) {return x <= L ? dd[x][0] : dd[n / x][1];}
int &F(ll x) {return x <= L ? ff[x][0] : ff[n / x][1];}
bool &vis(ll x) {return x <= L ? chk[x][0] : chk[n / x][1];}

int S(ll a, ll b) { return mul(mul(add(a % mod, b % mod), (b - a + 1) % mod), inv2); }

int calc(ll x){
    if (x <= L || vis(x)) return F(x);
    vis(x) = 1; 
    int sum = 0;
    for (ll i = 1; i * i <= x; i++) cadd(sum, mul(i, (int)sqrt(x - i * i)));
    for (ll i = 2; i * i <= x; i++) cpop(sum, mul(calc(x / (i * i)), i));
    F(x) = sum;
    return F(x);
}

int cal(ll x){
    if (x <= L) return D(x);
    int sum = 0;
    for (ll i = 1, o; i <= x; i = o + 1) {
        o = x / (x / i);
        cadd(sum, mul((x / i) % mod, S(i, o)));
    }
    return sum;
}

int main(){
    scanf("%lld", &n);
    up (i, 1, L) for (int j = i; j <= L; j += i) cadd(D(j), i);
    up (i, 1, L) cadd(D(i), D(i - 1));
    for (int i = 1; i * i <= L; i++) for (int t = i * i, j = 1; t + j * j <= L; j++) if (gcd(i, j) == 1)
        cadd(F(t + j * j), i);
    up (i, 1, L) cadd(F(i), F(i - 1));
    int ans = 0;
    for (ll i = 1, o; i <= n; i = o + 1) {
        o = n / (n / i);
        cadd(ans, mul(pop(calc(o), calc(i - 1)), cal(n / i)));
    }
    cmul(ans, 2);
    cadd(ans, cal(n));
    printf("%d\n", ans);
    printf("%d\n", cal(n));
    return 0;
}