[BOJ 10916] Xtreme gcd sum

이 문제는 [BOJ 14862] 최대공약수 기댓값의 확장 버전이라고 할 수 있다. 최대공약수 기댓값 문제를 풀면서 klimmek55님의 코드를 보고 감명받아 이 문제를 푼 것이다. 굉장히 멋진 아이디어이니 참고하면 좋을 것 같다.

 

문제에서 주어진 코드는 다음과 같이 바꿔쓸 수 있다.

\[\sum_{x_{1}=a_{1}}^{b_{1}}\sum_{x_{2}=a_{2}}^{b_{2}}...\sum_{x_{n}=a_{n}}^{b_{n}} \gcd(x_{1},x_{2},...,x_{n})\]

이런 경우에 Mobius inversion을 이용해서 다음과 같이 바꿔쓸 수 있다. (Mobius inversion 포스팅을 참고하라.)

\[\sum_{k=1}^{\min(b_{1},b_{2},...,b{n})} \left(\left[\frac{b_{1}}{k}\right]-\left[\frac{a_{1}-1}{k}\right]\right) \left(\left[\frac{b_{2}}{k}\right]-\left[\frac{a_{2}-1}{k}\right]\right)...\left(\left[\frac{b_{n}}{k}\right]-\left[\frac{a_{n}-1}{k}\right]\right)\]

함수 $f$를 이용해서 다음과 같이 쓰자. 단, $f(0)=f(1)$

\[\sum_{k=1}^{\min(b_{1},b_{2},...,b{n})} f(k)\]

이때 아래의 식은 $f(k)$의 인수라고 부른다.

\[\left(\left[\frac{b_{i}}{k}\right]-\left[\frac{a_{i}-1}{k}\right]\right)\]

이때 $f(k)$의 인수의 집합을 $s_{k}$라고 하자.

 

먼저 $dx[i],dy[i]$를 정의하자. ($i$는 자연수)

\[dx[i]=\prod_{x \in s_{i},\ x \ne 0}x \prod_{x \in s_{i-1},\ x \ne 0} \frac{1}{x}\]

($f(i)$의 인수들 중 0이 아닌 인수들의 곱과 $f(i-1)$의 인수들 중 0이 아닌 인수들의 역수의 곱을 곱한 것)

\[dy[i]=\sum_{x \in s_{i}} [x =0]-\sum_{x \in s_{i-1}} [x =0]\]

($f(i)$의 인수들 중 0인 것들의 개수$-f(i-1)$의 인수들 중 0인 것들의 개수)

 

우리가 $dx, dy$를 모두 구했다고 하자. $tot=f(1)$라고 하자. 이때 다음이 성립한다.

\[tot \cdot dx[1] = f(1)\]

\[tot \cdot dx[1] \cdot dx[2] = f(2)\]

\[tot \cdot dx[1] \cdot dx[2] \cdot dx[3] = f(3)\]

위와 같은 방식으로 $f(i)$를 전부 계산할 수 있다.

하지만 $f(i)$에 0인 인수가 껴있다면 $f(i)$가 0이기 때문에 처리할 필요가 없다. 반면 $dx$는 항상 0보다 크므로 $f(i)$가 0인 것을 알려줄 수 없다는 것이다. 그래서 필요한 것이 $dy$이다.

만약 $dx[i]$의 분모에 0인 인수가 $n$개 포함되어 있다면 반드시 $dx[i+1]$의 분자에 0인 인수가 $n$개 포함되어 있다. $dy[i]=n$이고 $dy[i+1]=-n$인 것이다. $tot$에 $dx[i+1]$까지 곱했다면 0인 인수가 포함되어 있는 $f(i)$는 사라지고 $f(i+1)$만 남게 된다. 그리고 이때의 $dy$의 누적 합은 0이다. 단순하게 위의 경우를 생각했지만 실제로는 연쇄적으로 $0$인 인수가 생기고 사라지는 현상이 발생한다. 결론적으로, $dy$의 누적 합이 $0$이 되었다는 것은 $0$인 인수가 모두 소거되었음을 의미하므로 누적 합이 0이 될때마다 답을 갱신해주면 된다. 

 

$dx$를 계산하는 법을 알아보자. 분모와 분자를 따로 계산한다.

$dx[1]*=\frac{b_{1}-a_{1}}{b_{1}-a_{1}}$로 시작한다. $dx[2]$의 분모는 $dx[1]$의 분자가 되고 $dx[2]$의 분자는 $dx[1]$의 분모의 각 항을 2로 나눈 값이 된다. 즉, $dx[2]*=\frac{b_{1}/2-a_{1}/2}{b_{1}-a_{1}}$. 이것을 $dx[b_{1}]$까지 반복하면 된다. 지금껏 $k=1$인 경우를 했으므로 $k=n$이 될 때까지 반복한다. 이해가 쉽도록 분수를 사용하였으나 계산할 때에는 모듈러 역원을 이용해 모두 정수로 계산해야 한다. 또한 $b_{1}/2$는 $b_{1}$을 2로 나눈 몫을 의미한다. $dy$는 $dx$를 계산하는 과정에서 얻을 수 있다.

 

모듈러 역원과 오일러 파이 함수는 미리 전처리를 1부터 $10^{6}$까지 해두는 것이 좋다.

$m$이 $b_{i}$의 최솟값이라고 한다면 $dx$계산에 $O(n\sqrt{m})$, 답을 계산하는 것이 $O(m)$이므로 전체 시간복잡도는 $O(n\sqrt{m}+m)$이다.

 

아래는 이 과정을 구현한 코드이다.

# 코드 보기

#include<stdio.h>
#define MAX 1000000
#define MAXN 10001
typedef long long ll;
const ll mod=1e9+7;
int a[MAXN],b[MAXN],pa[MAXN],pb[MAXN];
ll dx[MAX+1],dy[MAX+1],inv[MAX+1];
char sp[MAX+1];
ll phi[MAX+1];
ll fpow(int n,int k){
    ll s=1,h=n;
    for(;k;k>>=1){
        if(k&1) s=s*h%mod;
        h=h*h%mod;
    }
    return s;
}
int min(int a,int b){
    return a>b?b:a;
}
int main(){
    int i,j;
    for(i=1;i<=MAX;i++) phi[i]=1;
    for(i=2;i<=MAX;i+=2){
        if(sp[i]) continue;
        for(j=i;j<=MAX;j+=i){
            sp[j]=1;
            phi[j]=phi[j]*(i-1)%mod;
            ll h=j/i;
            while(h%i==0){
                h/=i;
                phi[j]=phi[j]*i%mod;
            }
        }
        if(i==2) i--;
    }
    for(i=1;i<=MAX;i++){
        inv[i]=fpow(i,mod-2);
        dx[i]=1;
    }
 
    int n,mini=1e9;
    scanf("%d",&n);
    for(i=0;i<n;i++){
        scanf("%d %d",&a[i],&b[i]);
        a[i]--;
        pa[i]=a[i];
        pb[i]=b[i];
        mini=min(mini,b[i]);
    }
    
    for(i=0;i<n;i++){
        for(j=1;j<=b[i];){
            int jump=1e9;
            if(pa[i]!=pb[i]) dx[j]=dx[j]*inv[pb[i]-pa[i]]%mod;
            else dy[j]--;
            
            if(a[i]/j!=b[i]/j) dx[j]=dx[j]*(b[i]/j-a[i]/j)%mod;
            else dy[j]++;
            
            if(a[i]>=j) jump=min(a[i]/(a[i]/j),jump);
            jump=min(b[i]/(b[i]/j),jump);
            
            pa[i]=a[i]/j;
            pb[i]=b[i]/j;
            
            j=jump+1;
        }
    }
    
    ll tot=1,ans=0;
    int cnt=0;
    for(i=0;i<n;i++) tot=tot*(b[i]-a[i])%mod;
    for(i=1;i<=mini;i++){
        tot=tot*dx[i]%mod;
        cnt+=dy[i];
        if(cnt==0){
            ans=(ans+tot*phi[i])%mod;
        }
    }
    printf("%lld",ans);
}

# 닫기