P3768 简单的数学题 (杜教筛+推式子）

题目描述

由于出题人懒得写背景了，题目还是简单一点好。

输入一个整数n和一个整数p，你需要求出(∑i=1n∑j=1nijgcd(i,j)) mod p(sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p(∑i=1n​∑j=1n​ijgcd(i,j)) mod p，其中gcd(a,b)表示a与b的最大公约数。

刚才题面打错了，已修改

输入输出格式

输入格式：

一行两个整数p、n。

输出格式：

一行一个整数(∑i=1n∑j=1nijgcd(i,j)) mod p(sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p(∑i=1n​∑j=1n​ijgcd(i,j)) mod p。

输入输出样例

输入样例#1： 复制

998244353 2000

输出样例#1： 复制

883968974

说明

对于20%的数据，n≤1000n leq 1000n≤1000。

对于30%的数据，n≤5000n leq 5000n≤5000。

对于60%的数据，n≤106n leq 10^6n≤106，时限1s。

对于另外20%的数据，n≤109n leq 10^9n≤109，时限3s。

对于最后20%的数据，n≤1010n leq 10^{10}n≤1010，时限6s。

对于100%的数据，5×108≤p≤1.1×1095 imes 10^8 leq p leq 1.1 imes 10^95×108≤p≤1.1×109且p为质数。

---

---

---

---

---

_{这题首先不解释： #define int usigned long long   qwq}

少mod一个这题就得wa，好毒瘤啊。

---

---

---

---

---

---

---

#include<bits/stdc++.h>
using namespace std;

#define int unsigned long long

int read() {
    char cc = getchar(); int cn = 0, flus = 1;
    while(cc < '0' || cc > '9') {  if( cc == '-' ) flus = -flus;  cc = getchar();  }
    while(cc >= '0' && cc <= '9')  cn = cn * 10 + cc - '0', cc = getchar();
    return cn * flus;
}

const int N = 5e6 + 5;
#define inf 9999999
#define g(x) (((( (x % mod) * ((x + 1) % mod) / 2 ) % mod) * ( (x % mod) * ((x + 1) % mod) / 2 % mod)) % mod)
#define g2(x) (( (x % mod) * ((x + 1) % mod) % mod * ((2 * x + 1) % mod) % mod * inv6 % mod) % mod)

int n, phi[N], pr[N], f[N], top, mod, inv6, inv2;
bool isp[N];

map<int, int> map1;
void init() {
    phi[1] = f[1] = 1;
    for ( int i = 2; i <= n; ++ i ) {
        if( !isp[i] )  pr[++top] = i, phi[i] = i - 1;
        for ( int j = 1; j <= top; ++ j ) {
            if( i * pr[j] > n ) break;
            isp[i * pr[j]] = 1;
            if( i % pr[j] == 0 ) {
                phi[i * pr[j]] = phi[i] * pr[j];
                break;
            }
            phi[i * pr[j]] = phi[i] * phi[pr[j]];
        }
        f[i] = f[i - 1] + ((1ll * phi[i] * i) % mod * i) % mod; f[i] %= mod;
    }
} 

int fpow( int x, int k ) {
    int base = x, ans = 1;
    while( k ) {
        if( k & 1 ) ans *= base, ans %= mod;
        base *= base, base %= mod;
        k >>= 1;
    }
    return ans;
}

int phi_sum( int x ) {
    if( x <= n ) return f[x];
    if( map1[x] ) return map1[x];
    int sumId = g(x % mod) % mod, ans = 0;
    int l, r;
    for ( l = 2; l <= x; l = r + 1 ) {
        r = ( x / ( x / l ) );
        ans += ((1ll * (( g2(r) - g2((l - 1)) + mod ) % mod ) * phi_sum( x / l )) % mod);
        ans %= mod;
    }
    return map1[x] = ( (sumId - ans + mod) % mod );
}

int solve1( int x ) {
    int l, r, ans = 0;
    for ( l = 1; l <= x; l = r + 1 ) {
        r = ( x / ( x / l ) );
        int ret = g( x / l ) % mod;
        ret = ret * ((phi_sum(r) - phi_sum(l - 1) + mod) % mod), ret %= mod; 
        ans += ret, ans %= mod;
    }
    return (ans + mod) % mod;
}

signed main()
{
    mod = read(); 
//  inv2 = fpow( 2, mod - 2 );
    inv6 = fpow( 6, mod - 2 ); 
    int x = read();
    n = 5000000;  init();
    printf("%lld", solve1(x));
    return 0;
}