zoukankan      html  css  js  c++  java
  • [Luogu3768]简单的数学题

    Portal

    [求 sum_{i = 1}^{n}sum_{j = 1}^{n}(i,j)ij ]


    好像直接利用(varphi)很好做诶:

    [sum_{i = 1}^{n}sum_{j = 1}^{n}(i,j)ij\ = sum_{i = 1}^{n} isum_{j = 1}^{n}j sum_{d | (i,j)} varphi(d)\ = sum_{i = 1}^{n} isum_{j = 1}^{n}j sum_{d | i,d | j} varphi(d)\ = sum_{d = 1}^{n} d ^ 2varphi(d) sum_{i = 1}^{lfloor frac{n}{d} floor}isum_{j = 1}^{lfloor frac{n}{d} floor}j\ ]

    令$Sum(n) = sum_{i = 1}^{n} i $

    那么有:

    [原式 = sum_{d = 1}^{n}d^2varphi(d) Sum^2(lfloorfrac{n}{d} floor) ]

    后面一部分可以整除分块。考虑前面一部分:

    构造:(g(n) = n ^ 2)

    [(f * g)(n) = sum_{d | n} d^2varphi(d) * frac{n ^ 2}{d ^ 2} \ = n ^ 2 sum_{d | n} varphi(d) = n ^ 3 ]

    有:(sum{n ^ 3} = (sum{n})^2)

    然后杜教筛筛出前面的和,然后整除分块即可。


    这里有莫比乌斯反演推式子的办法,可以参考一下

    #include <bits/stdc++.h>
    #include <bits/extc++.h>
    using namespace std;
    #define rep(i, a, b) for(LL i = (a), i##_end_ = (b); i <= i##_end_; ++i)
    #define drep(i, a, b) for(LL i = (a), i##_end_ = (b); i >= i##_end_; --i)
    #define clar(a, b) memset((a), (b), sizeof(a))
    #define debug(...) fprintf(stderr, __VA_ARGS__)
    typedef long long LL;
    typedef long double LD;
    LL read() {
        char ch = getchar();
        LL x = 0, flag = 1;
        for (;!isdigit(ch); ch = getchar()) if (ch == '-') flag *= -1;
        for (;isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
        return x * flag;
    }
    
    __gnu_pbds :: gp_hash_table <LL, LL> Sum;
    const int Maxn = 10000009;
    LL isnprime[Maxn], prime[Maxn], tot, phi[Maxn];
    LL prefix[Maxn], Mod, n;
    
    LL fpm(LL a, LL tims) {
    	LL r = 1;
    	for (; tims; tims >>= 1, a = a * a % Mod) 
    		if (tims & 1) r = r * a % Mod;
    	return r;
    }
    
    LL inv(LL a) { return fpm(a, Mod - 2); }
    
    void init() {
    	Mod = read();
    	phi[1] = 1;
    	rep (i, 2, Maxn - 1) {
    		if (!isnprime[i]) 
    			prime[++tot] = i, phi[i] = i - 1;
    
    		for (int j = 1, k; j <= tot && (k = i * prime[j]) < Maxn; ++j) {
    			isnprime[k] = 1;
    			if (i % prime[j] == 0) {
    				phi[k] = phi[i] * prime[j];
    				break;
    			} else phi[k] = phi[prime[j]] * phi[i];
    		}
    	}
    
    	rep (i, 1, Maxn - 1) prefix[i] = (prefix[i - 1] + (phi[i] * i % Mod * i % Mod)) % Mod;
    }
    
    int inv2, inv6;
    inline LL singleSum(LL val) { return val % Mod * (val % Mod + 1) % Mod * inv2 % Mod; }
    inline LL doubleSum(LL val) { return (val % Mod * (val % Mod + 1) % Mod * (val % Mod * 2 + 1) % Mod) * inv6 % Mod; }
    inline LL tripleSum(LL val) { return singleSum(val) * singleSum(val) % Mod; }
    
    LL calcSum(LL val) {
    	if (val < Maxn) return prefix[val];
    	if (Sum[val]) return Sum[val];
    
    	LL ans = tripleSum(val);
    	for (LL l = 2, r; l <= val; l = r + 1) {
    		r = val / (val / l);
    		LL res = ((doubleSum(r) - doubleSum(l - 1)) % Mod + Mod) % Mod;
    		(ans += (Mod - res * calcSum(val / l) % Mod) % Mod) %= Mod;
    	}
    
    	return Sum[val] = (ans + Mod) % Mod;
    }
    
    void solve() {
    	inv2 = inv(2), inv6 = inv(6);
    	n = read();
    
    	LL ans = 0; 
    	for (LL l = 1, r; l <= n; l = r + 1) {
    		r = n / (n / l);
    		(ans += tripleSum(n / l) * ((calcSum(r) - calcSum(l - 1) + Mod) % Mod) % Mod) %= Mod;
    	}
    
    	cout << ans << endl;
    }
    
    int main() {
    	freopen("LG3768.in", "r", stdin);
    	freopen("LG3768.out", "w", stdout);
    
    	init();
    	solve();
    
    #ifdef Qrsikno
        debug("
    Running time: %.3lf(s)
    ", clock() * 1.0 / CLOCKS_PER_SEC);
    #endif
        return 0;
    }
    
  • 相关阅读:
    matplotlib 进阶之origin and extent in imshow
    Momentum and NAG
    matplotlib 进阶之Tight Layout guide
    matplotlib 进阶之Constrained Layout Guide
    matplotlib 进阶之Customizing Figure Layouts Using GridSpec and Other Functions
    matplotlb 进阶之Styling with cycler
    matplotlib 进阶之Legend guide
    Django Admin Cookbook-10如何启用对计算字段的过滤
    Django Admin Cookbook-9如何启用对计算字段的排序
    Django Admin Cookbook-8如何在Django admin中优化查询
  • 原文地址:https://www.cnblogs.com/qrsikno/p/10111165.html
Copyright © 2011-2022 走看看