利用容斥原理求尤拉計劃565題的S(1E9,2017)

題目介紹

找出能被2017整除的小於或等於ubound平方根的質數的冪(1-9行)                            
加總能被2017整除的大於ubound平方根的質數的倍數，因為這些不同質數的倍數沒有公因數，因此它們之間不會重複                            
將ubound除以1-9行的值，得出倍數的個數，用等差數列公式求和                            
因為這些和包含重複計算，因此要扣除                            
第一列標黃部分是大於ubound平方根，但小於12101的最大倍數(82367)的能被2017整除的質數，因為它們在計算第5列標藍的資料時用過，需要扣除                            
只要質數冪乘以另一個質數冪的倍數小於ubound，都需要扣除                            
本題恰好只有1次方符合扣除條件                            
本題不存在3個質數冪連乘積小於ubound的情況                            

利用這個思路，不需要儲存已用過標記。利用篩法分步計算求可以求得S(1E11,2017)                            

難看的程式碼

int d[100000]= {0};
char b[100000]= {0};
int k=0;

long long in_exclude()
{
    //include
    int min_p=1e6;
    long long p_times[1000];
    long long n=0;
    long long x=0;
    for(int i=0; i<k; i++)
    {
        x=pow(d[i],b[i]-1);
        if(x<min_p)min_p=x;
        if(x>NN)break;
        n=NN/x;
        sum0+=x*n*(n+1)/2;

    }
    printf("sum0=%lld\n",sum0);
    //exclude
    int j=0;
    for(int i=0; i<=NN/min_p; i++)
        if(a[i]==1 && (i+1)%mmm==0)
        {
            p_times[j++]=i;
        }
    int j1=j;

    for(int i=0; i<k; i++)
    {
        x=pow(d[i],b[i]-1);
        if(x<=sqrt(NN))
            for(j=0; j<j1; j++)
            {
                if(NN/p_times[j]<x)break;
                if(x>p_times[j])continue;
                long long xx=x*p_times[j];
                long long nn=NN/xx;
                if(nn<1)break;
                sum0-=xx*nn*(nn+1)/2;
            }
    }

    printf("sum=%lld\n",sum0);
    return sum0;
}