zoukankan      html  css  js  c++  java
  • 【51nod】2026 Gcd and Lcm

    题解

    话说LOJ说我今天宜学数论= =看到小迪学了杜教筛去蹭了一波小迪做的题

    标解的杜教筛的函数不懂啊,怎么推的毫无思路= =
    所以写了个复杂度稍微高一点的??

    首先,我们发现f是个积性函数,那么我们就有……
    (prod_{i = 1}^{k}f(p_{i}^{a_{i}}))
    我们发现,对于每个质因子,gcd是取较小值,lcm取较大值
    (f(lcm(x,y)) * f(gcd(x,y)) = prod_{i = 1}^{k} f(p_{i}^{max(a_{i},b_{i}) + min(a_{i},b_{i})}))
    (max(a,b) + min(a,b) = a + b)
    那么就有
    (f(lcm(x,y)) * f(gcd(x,y)) = prod_{i = 1}^{k} f(p_{i}^{max(a_{i},b_{i}) + min(a_{i},b_{i})}) = f(x) * f(y))

    所以我们只要求出([sum_{i = 1}^{n} f(i)]^2)就是答案了!

    怎么求呢
    (S(n) = sum_{i = 1}^{n}sum_{d | i} mu(d)cdot d)
    (S(n) = sum_{d = 1}^{n}sum_{d | i} mu(d)cdot d)
    (S(n) = sum_{d = 1}^{n} mu(d)cdot d cdot lfloor frac{n}{d} floor)
    我们可以数论分块处理(lfloor frac{n}{d} floor)
    那么我们考虑计算(sum_{d = 1}^{n} mu(d)cdot d)
    我们发现这个函数卷上一个(Id(x))等于(e)

    (sum_{i = 1}^{n} [i = 1] = sum_{i = 1}^{n} sum_{d |i} mu(d) cdot d cdot frac{i}{d} = sum_{k = 1}^{n} k sum_{d}^{frac{n}{k}} mu(d) cdot d)
    所以最后就是
    (S(n) = 1 - sum_{i = 2}^{n} S(lfloor frac{n}{i} floor))

    代码

    #include <iostream>
    #include <cstdio>
    #include <vector>
    #include <algorithm>
    #include <cmath>
    #include <cstring>
    #include <map>
    //#define ivorysi
    #define pb push_back
    #define space putchar(' ')
    #define enter putchar('
    ')
    #define mp make_pair
    #define pb push_back
    #define fi first
    #define se second
    #define mo 974711
    #define MAXN 1000000
    #define RG register
    using namespace std;
    typedef long long int64;
    typedef double db;
    template<class T>
    void read(T &res) {
        res = 0;char c = getchar();T f = 1;
    	    while(c < '0' || c > '9') {
    			if(c == '-') f = -1;
    			c = getchar();
    	    }
    	    while(c >= '0' && c <= '9') {
    		res = res * 10 + c - '0';
    		c = getchar();
        }
        res *= f;
    }
    template<class T>
    void out(T x) {
        if(x < 0) {putchar('-');x = -x;}
        if(x >= 10) {
    		out(x / 10);
        }
        putchar('0' + x % 10);
    }
    int N;
    const int MOD = 1000000007;
    struct node {
    	int x,v,next;
    }E[2000006];
    int head[mo + 5],sumE;
    int prime[MAXN + 5],tot,S[MAXN + 5],mu[MAXN + 5];
    bool nonprime[MAXN + 5];
    int inc(int a,int b) {
    	a = a + b;
    	if(a >= MOD) a -= MOD;
    	return a;
    }
    void add(int u,int x,int v) {
    	E[++sumE].x = x;E[sumE].v = v;E[sumE].next = head[u];
    	head[u] = sumE;
    }
    void Insert(int x,int v) {
    	add(x % mo,x,v);
    }
    int Query(int x) {
    	int u = x % mo;
    	for(int i = head[u] ; i ; i = E[i].next) {
    		if(E[i].x == x) return E[i].v;
    	}
    	return -1;
    }
    int f(int x) {
    	if(x <= MAXN) return S[x];
    	int c = Query(x);
    	if(c != -1) return c;
    	int res = 0;
    	for(int i = 2 ; i <= x ; ++i) {
    		int r = x / (x / i);
    		res = inc(res,1LL * (r - i + 1) * (r + i) / 2 % MOD * f(x / i) % MOD);
    		i = r;
    	}
    	res = inc(1,MOD - res);
    	Insert(x,res);
    	return res;
    }
    void Solve() {
    	mu[1] = 1;S[1] = 1;
    	for(int i = 2 ; i <= MAXN ; ++i) {
    		if(!nonprime[i]) {
    			mu[i] = -1;
    			prime[++tot] = i;
    		}
    		for(int j = 1 ; j <= tot ; ++j) {
    			if(prime[j] > MAXN / i) break;
    			nonprime[i * prime[j]] = 1;
    			if(i % prime[j] == 0) break;
    			else mu[i * prime[j]] = -mu[i];
    		}
    		S[i] = (S[i - 1] + mu[i] * i + MOD) % MOD;
    	}
    	read(N);
    	int res = 0;
    	for(int i = 1 ; i <= N ; ++i) {
    		int r = N / (N / i);
    		res = inc(1LL * (f(r) + MOD - f(i - 1)) * (N / i) % MOD,res);
    		i = r;
    	}
    	out(1LL * res * res % MOD);enter;
    }
    int main() {
    #ifdef ivorysi
        freopen("f1.in","r",stdin);
    #endif
       	Solve();
        return 0;
    }
    
  • 相关阅读:
    Python进阶08 生成器
    Python进阶07 迭代
    Python进阶06 列表推导
    Python进阶05 函数property
    Python基础14 集合
    支付宝支付和微信消息推送
    Redis
    django之contenttype
    数据分析系列
    IPython
  • 原文地址:https://www.cnblogs.com/ivorysi/p/9157781.html
Copyright © 2011-2022 走看看