n以内与n互素的平方和（算数基本定理+容斥原理）

题目：多组输入，每组一个正整数n（n<=10^6）,求n以内与n互素的平方和

思路：先将n素数分解一下，求出所有素数因子，然后求出所有与n不互质的数的平方和（容斥原理），再用1-n（利用公式（n*（n+1）*（2n+1）/6）的平方和减去这些与n不互质的平方和。例如n为12，它的素数因子为2,3，所以所有与n不互质的数为2,3,4,6,8，9,10。其中素因子2 的贡献度为2^2+4^2+6^2+8^2+10^2=2^2（1^2+2^2+3^2+4^2+5^2）=2^2*(5*(5+1)*(2*5+1))，素因子3的贡献度为3^2+6^2+9^2=3^2*(1^2+2^3+3^2)=3^2(3*(3+1)*(2*3+1))，由于6（2*3）被重复算了一次再减去n以内6的倍数的平方和。（奇加偶减）容斥原理用代码实现一般有两种方式，一是用二进制0 1标记法，二是用dfs（深度优先遍历）。两种方法都会在下面代码中给出。

代码一：

 

/// 二进制0 1标记法
#include<algorithm>
#include<iostream>
#include<cstring>
#include<string>
using namespace std;
#define LL long long
const int maxn = 1e6+5;
int prim[maxn/10],cnt=0,m;

void getPrime(){
    int n=m;
    cnt=0;
    for(int i=2;i*i<=n;i++){
        if(n%i==0) {
            while(n%i==0) n/=i;
            prim[cnt++]=i;
        }
    }
    if(n>1) prim[cnt++]=n;
}

LL cal(int x){
    return (LL)x*(x+1)*(2*x+1)/6;///1-x的平方和公式 n较大所以要用long long
}

LL solve(){
    LL ans = 0;
    for(int i=1;i<(1<<cnt);i++){
        int f=0,s=1;
        for(int j=0;j<cnt;j++){
            if(i& (1<<j)){///<< 二进制左移运算符 左移n位就扩大2^n倍
                s*=prim[j];
                f++;
            }
        }
        ///奇加偶减
        if(f&1) ans += (LL)s*s*cal(m/s);
        else ans -= (LL)s*s*cal(m/s);
    }
    return ans;
}

int main(){
    while(cin>>m){
        getPrime();
        if(m==1) cout<<0<<endl;///注意1要特殊处理
        else cout<<cal(m)-solve()<<endl;
    }
}

代码二：

 

 

/// dfs
#include<algorithm>
#include<iostream>
#include<cstring>
#include<string>
using namespace std;
#define LL long long
const int maxn = 1e6+5;
int prim[maxn/10],cnt=0,m;
LL ans=0;

void getPrime(){///素数分解
    int n=m;
    cnt=0;
    for(int i=2;i*i<=n;i++){
        if(n%i==0) {
            while(n%i==0) n/=i;
            prim[++cnt]=i;
        }
    }
    if(n>1) prim[++cnt]=n;
}

LL cal(int x){
    return (LL)x*(x+1)*(2*x+1)/6;
}

void dfs(LL k,LL val,LL id){
    if(val>=m) return ;
    LL y,u;
    if(k%2==1){
        u=(m-1)/val;
        y=u*(u+1)*(2*u+1)/6;
        y=y*val*val;
        ans=ans-y;
    }
    else{
        u=(m-1)/val;///除去m的影响
        y=u*(u+1)*(2*u+1)/6;
        y=y*val*val;
        ans=ans+y;
    }
    for(int i=id+1;i<=cnt;i++)
        dfs(k+1,val*prim[i],i);
}

int main(){
    while(cin>>m){
        getPrime();
        ans=cal(m-1);///这里计算的是m-1以内的立方和 除去了m本身
        if(m==1) cout<<0<<endl;
        else {
            for(int i=1;i<=cnt;i++)
                dfs((LL)i,(LL)prim[i],(LL)i);
            cout<<ans<<endl;
        }
    }
}

这里再补充说明一下dfs是如何进行容斥的，以素因子2,3,5,7为例子结合代码二，如下图。