BZOJ3309 DZY Loves Math

參考自這位大佬的推導思路以及這位大佬的實現思路

題目簡述

定義\(f(n)\)為\(n\)所含質因子的最大冪指數，如\(f(18)=f(2*3^2)=2\)
輸入T組數據，求\(\sum_{i=1}^{n} \sum_{j=1}^{m} f(gcd(i,j))\)的值
\(T<=10000,1<=a,b<=10^7\)

公式推導(略)

\[\sum_{i=1}^{n} \sum_{j=1}^{m} f(gcd(i,j))\]
\[=\sum_{T=1}^{min(n,m)} \lfloor \frac nT \rfloor \lfloor \frac mT \rfloor \sum_{d|T} f(d) \mu(\frac Td)\]

• 令\(g(T)=\sum_{d|T}f(d) \mu(\frac Td)\)
  則現在的目標是在線性時間內求出\(g(x)\)的值
  設\(T=p_1^{a_1}p_2^{a_2}…p_k^{a_k}\)，\(d=p_1^{b_1}p_2^{b_2}…p_k^{b_k}\)，\(c_i=a_i-b_i\)
  則\(\frac Td=p_1^{c_1}p_2^{c_2}…p_k^{c_k}\)
  要使\(\mu(\frac Td) \neq0\)，則\(c_1,c_2,…,c_k \leq1\)

分類討論\(g(T)\)的值

1. \(a_1,a_2,…,a_k\)並不全都相等
   此時設\(a_x\)最大，對於所有包含\(p_x^{a_x}\)

   的\(d\)均有\(f(d)=a_x\)
   對於這些\(d\)，\(p_i(i \neq x)\)的指數的選擇均有兩個，為\(a_i-1,a_i\)，則一共能組合出\(2^{k-1}\)個不同的\(d\)
   設其中\(d_1,d_2\)的質因子中除\(p_y\)一項之外的指數均相等，則\(\mu(d_1)=-\mu(d_2)\)
   而這樣\(\mu\)值相反的\(d\)共有\(2^{k-2}\)組，各自相加後的和就為0即
   \[[\sum_{d|T,p_x^{a_x} |d} f(d) \mu(\frac Td) ]=0\]
   之後對於剩下取\(p_x^{a_x-1}\)的\(d\) 執行同樣的分析，發現所有以\(a_q\)
   為最大值的的\(d\)的和均為0，所以此時\(g(T)=0\)

2. 令\(A=a_1=a_2=…=a_k\)即\(a_1,a_2,…,a_k\)均相等
   此時\(f(d)=\begin {cases} A(b_1,b_2,…,b_k中至少有一個值為A)\\ A-1(b_1=b_2=…=b_k=A-1) \end {cases}\)
   假設\(b_1=b_2=…=b_k=A-1\)的時\(d\) 值\(f(d)\)也為A
   則\(g(T)\)類似於上一種情況的推導，其的和也為0
   當然，最後\(g(T)\)要減去某特殊\(d\)多加的\(1*\mu(\frac Td)\)
   所以此情況下\(g(T)=-(-1)^{k}=(-1)^{k+1}\)

而線性求g(T),則需要記錄每個點的最小質因子及個數，除去該最小質因子後的值
這裏就感性理解一下吧OvO

#include <iostream>
#include <cstdio>
#include <cstring>
#define N 10000010
using namespace std;
int n,prim[N],pcnt,at[N],la[N],cnt[N];//at[i]記錄i的最小質因子，la[i]=記錄i除以所有at[i]後的值，cnt[i]記錄at[i]的數量 
long long sum[N],g[N],f[N];
void eular()
{
    g[1]=0;
    for(int i=2;i<=N-10;i++)
    {
        if(!at[i])
        {
            prim[++pcnt]=i;
            at[i]=pcnt;
            g[i]=la[i]=cnt[i]=1;
        }
        for(int j=1;j<=at[i];j++)
        {
            if((long long)i*prim[j]>N-10)break;
            at[i*prim[j]]=j;
            if(j==at[i])//prim[j]為i的最小質因子，即也為prim[j]*i的最小質因子 
            {
                la[i*prim[j]]=la[i];
                cnt[i*prim[j]]=cnt[i]+1;
                if(la[i]==1)g[i*prim[j]]=1;
                else g[i*prim[j]]=(cnt[i]+1==cnt[la[i]]?-g[la[i]]:0);
            }
            else//prim[j]為新出現的質因子，個數為1 
            {
                la[i*prim[j]]=i;
                cnt[i*prim[j]]=1;
                g[i*prim[j]]=(cnt[i]==1?-g[i]:0);
            }
        }
        sum[i]=sum[i-1]+g[i];
    }
}
long long mb(int n,int m)
{
    long long ans=0;
    if(n>m)swap(n,m);
    for(int l=1,r;l<=n;l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        ans+=(sum[r]-sum[l-1])*(n/l)*(m/l);
    }
    return ans;
}
int main(){
    eular();
    int a,b,T;
    cin>>T;
    while(T--)
    {
        scanf("%d%d",&a,&b);
        printf("%lld\n",mb(a,b));   
    }
    return 0;
} 

BZOJ3309 DZY Loves Math