题目描述

本题保证数据随机生成

有一个 $n$ 个点的完全图，结点编号为 $1,2……n$，结点 $i,j$ 之间有权值为 $\mathrm{lcm}(i,j)+\mathrm{gcd}(i,j)$ 的无向边。

请你求出这个图的最大生成树的大小，对 $2^{32}$ 取模。

输入格式

本题有多组测试数据，第一行一个整数 $T$ 代表数据组数。

接下来 $T$ 行，每行一个整数 $n$ 代表点数。

输出格式

共 $T$ 行，每行一个整数代表一组数据的答案。

输入输出样例 #1

输入 #1

9
10
1000
100000
10000000
1000000000
100000000000
10000000000000
1000000000000000
100000000000000000

输出 #1

422
499008694
4172096327
3128649679
2692599804
194024000
2969759816
505684415
3052141644

说明/提示

有子任务限制

• 对于 $5\%$ 的数据，$T=1,n\leq 2000$。
• 对于 $15\%$ 的数据，$T=1,n\leq 10^6$。
• 对于 $30\%$ 的数据，$T=1,n\leq 10^8$。
• 对于 $50\%$ 的数据，$T=1,n\leq 10^{10}$。
• 对于 $75\%$ 的数据，$T=1,n\leq 10^{15}$。
• 对于 $100\%$ 的数据，$T\leq 10,n\leq 10^{18}$。

下发文件中提供了一份高速的 Pollard_Rho 分解质因数代码。

namespace pollardrho
{
    using i64 = long long;
    using u64 = unsigned long long;
    using u32 = unsigned int;
    using u128 = __uint128_t;
    namespace
    {
        template <class U0, class U1>
        struct Montgomery
        {
            constexpr static unsigned B0 = sizeof(U0) * 8U;
            U0 n, nr, rs, np;
            constexpr Montgomery(const U0 &Mod) { SetMod(Mod); }
            constexpr U0 GetMod() const noexcept { return n; }
            constexpr void SetMod(const U0 &Mod)
            {
                // assert(Mod >= 2), assert(Mod % 2 == 1);
                // assert((Mod >> (B0 - 2)) == 0);
                n = nr = Mod, rs = -static_cast<U1>(n) % n;
                for (u32 i = 0; i < __lg(B0); i++)
                {
                    nr *= 2 - n * nr;
                }
                np = Reduce(static_cast<U0>(1), rs);
            }
            constexpr U0 Reduce(const U0 &x) const noexcept
            {
                const U0 q = x * nr;
                const U0 m = (static_cast<U1>(q) * n) >> B0;
                return n - m;
            }
            constexpr U0 Reduce(const U0 &x, const U0 &y) const noexcept
            {
                const U1 t = static_cast<U1>(x) * y;
                const U0 c = t, d = t >> B0;
                const U0 q = c * nr;
                const U0 m = (static_cast<U1>(q) * n) >> B0;
                return d + n - m;
            }
            constexpr U0 Reduce(const U0 &x, const U0 &y, const U0 &z) const noexcept
            {
                const U1 t = static_cast<U1>(x) * y;
                const U0 c = t, d = t >> B0;
                const U0 q = c * nr;
                const U0 m = (static_cast<U1>(q) * n) >> B0;
                return z + d + n - m;
            }
            constexpr U0 val(const U0 &x) const noexcept
            {
                const u64 t = Reduce(x);
                return (t == n) ? static_cast<U0>(0) : t;
            }
            constexpr U0 zero() const noexcept { return static_cast<U0>(0); }
            constexpr U0 one() const noexcept { return np; }
            constexpr U0 raw(const U0 &x) const noexcept { return Reduce(x, rs); }
            template <class U, typename std::enable_if<std::is_unsigned<U>::value, int>::type = 0>
            constexpr U0 trans(const U &x) const noexcept
            {
                if (__builtin_expect(x < n, 1))
                    return raw((U0)x);
                return Reduce((U0)(x % n), rs);
            }
            template <class S, typename std::enable_if<std::is_signed<S>::value, int>::type = 0>
            constexpr U0 trans(S x) const noexcept
            {
                if (__builtin_expect(0 <= x && x < static_cast<S>(n), 1))
                    return raw((U0)x);
                x %= static_cast<S>(n);
                if (x < 0)
                    x += static_cast<S>(n);
                return Reduce((U0)x, rs);
            }
            constexpr U0 neg(const U0 &x) const noexcept { return (x != 0) ? (2 * n - x) : x; }
            constexpr U0 inc(const U0 &x) const noexcept { return add(x, np); }
            constexpr U0 dec(const U0 &x) const noexcept { return sub(x, np); }
            constexpr U0 add(const U0 &x, const U0 &y) const noexcept { return (x + y >= 2 * n) ? (x + y - 2 * n) : (x + y); }
            constexpr U0 sub(const U0 &x, const U0 &y) const noexcept { return (x < y) ? (x - y + 2 * n) : (x - y); }
            constexpr U0 mul(const U0 &x, const U0 &y) const noexcept { return Reduce(x, y); }
            constexpr U0 mul_add(const U0 &x, const U0 &y, const U0 &z) const noexcept { return Reduce(x, y, z); }
            constexpr bool same(const U0 &x, const U0 &y) const noexcept
            {
                const U0 dif = x - y;
                return (dif == 0) || (dif == n) || (dif == -n);
            }
        };
        constexpr bool Is_Prime(u64 x) noexcept
        {
            if (x <= 1)
            {
                return false;
            }
            if (x % 2 == 0)
            {
                return x == 2;
            }
            constexpr array<u64, 10> Base{2, 7, 61, 2, 325, 9375, 28178, 450775, 9780504, 1795265022};
            const u32 s = __builtin_ctzll(x - 1);
            const u64 d = (x - 1) >> s;
            const int q = 63 ^ __builtin_clzll(d);
            const Montgomery<u64, u128> Mod(x);
            const u32 l = (x >> 32) ? 3 : 0;
            const u32 r = (x >> 32) ? 10 : 3;
            for (u32 _ = l; _ < r; _++)
            {
                u64 base = Base[_];
                if (base % x == 0)
                {
                    continue;
                }
                base = Mod.trans(base);
                u64 a = base;
                for (int i = q - 1; ~i; i--)
                {
                    a = Mod.mul(a, a);
                    if ((d >> i) & 1)
                    {
                        a = Mod.mul(a, base);
                    }
                }
                if (Mod.same(a, Mod.one()))
                {
                    continue;
                }
                for (u32 t = 1; t < s && !Mod.same(a, x - Mod.one()); ++t)
                {
                    a = Mod.mul(a, a);
                }
                if (!Mod.same(a, x - Mod.one()))
                {
                    return false;
                }
            }
            return true;
        }
        vector<u64> Factorize(u64 n)
        {
            vector<u64> ans;
            if (n % 2 == 0)
            {
                u32 z = __builtin_ctzll(n);
                ans.push_back(2ULL), n >>= z;
            }
            auto upd = [&](const u64 &x)
            {
                for (auto p : ans)
                {
                    if (x == p)
                        return;
                }
                ans.push_back(x);
            };
            auto Pollard_Rho = [&](const u64 &n)
            {
                if (n % 2 == 0)
                    return 2ULL;
                const Montgomery<u64, u128> Mod(n);
                const u64 C1 = 1, C2 = 2, M = 512;
                u64 Z1 = 1, Z2 = 2, ans = 0;
                auto find = [&]()
                {
                    u64 z1 = Z1, z2 = Z2;
                    for (u64 k = M;; k *= 2)
                    {
                        const u64 x1 = z1 + n, x2 = z2 + n;
                        for (u64 j = 0; j < k; j += M)
                        {
                            const u64 y1 = z1, y2 = z2;
                            u64 q1 = 1, q2 = 2;
                            z1 = Mod.mul_add(z1, z1, C1), z2 = Mod.mul_add(z2, z2, C2);
                            for (u64 i = 0; i < M; ++i)
                            {
                                u64 t1 = x1 - z1, t2 = x2 - z2;
                                z1 = Mod.mul_add(z1, z1, C1), z2 = Mod.mul_add(z2, z2, C2);
                                q1 = Mod.mul(q1, t1), q2 = Mod.mul(q2, t2);
                            }
                            q1 = Mod.mul(q1, x1 - z1), q2 = Mod.mul(q2, x2 - z2);
                            const u64 q3 = Mod.mul(q1, q2), g3 = __gcd(n, q3);
                            if (g3 == 1)
                            {
                                continue;
                            }
                            if (g3 != n)
                            {
                                ans = g3;
                                return;
                            }
                            const u64 g1 = __gcd(n, q1);
                            const u64 g2 = __gcd(n, q2);
                            const u64 C = g1 != 1 ? C1 : C2;
                            const u64 x = g1 != 1 ? x1 : x2;
                            u64 z = g1 != 1 ? y1 : y2;
                            u64 g = g1 != 1 ? g1 : g2;
                            if (g == n)
                            {
                                do
                                {
                                    z = Mod.mul_add(z, z, C);
                                    g = __gcd(n, x - z);
                                } while (g == 1);
                            }
                            if (g != n)
                            {
                                ans = g;
                                return;
                            }
                            Z1 += 2, Z2 += 2;
                            return;
                        }
                    }
                };
                do
                {
                    find();
                } while (!ans);
                return ans;
            };
            auto DFS = [&](auto &&self, const u64 &n) -> void
            {
                if (Is_Prime(n))
                    return upd(n);
                u64 d = Pollard_Rho(n);
                self(self, d), self(self, n / d);
            };
            if (n > 1)
                DFS(DFS, n);
            sort(ans.begin(), ans.end());
            return ans;
        }
    }
}