int solve(int n)
{
    int sqrtn = sqrt(n + 0.5), cnt, ans = 1;
    for (int i = 2; i <= sqrtn; ++i) if (n % i == 0)
        {
            for (n /= i, cnt = 1; n % i == 0; n /= i) ++cnt;
            if ((i & 3) == 1) ans *= cnt << 1 | 1;
        }
    if (n > 1 && (n & 3) == 1) ans += ans << 1;
    return (ans - 1) >> 1;
}

使用 Pollard's rho 算法可以加速到,代码就不贴了。。

参考 A046080 中的 Mathematica 代码:

cnt[n_] := If[n == 1, 0, With[{expList = Select[FactorInteger[n], Mod[#[[1]], 4] == 1 &][[All, 2]]}, ((Times @@ (2 * expList + 1)) - 1) / 2]];

证明待填。。

Comments

comments powered by Disqus