Project Eular 632/633

先讨论一下632

我们先求一遍,然后容斥

先求一遍的时候,我们对于每个数x

如果x的因子里面带平方,那么就忽略掉,否则就在ans[x的因子个数]上加(10¹⁶/x²)的答案

之后容斥,我们就可以求出这个答案

#include<set>
#include<map>
#include<list>
#include<queue>
#include<stack>
#include<string>
#include<time.h>
#include<math.h>
#include<memory>
#include<vector>
#include<bitset>
#include<fstream>
#include<stdio.h>
#include<utility>
#include<sstream>
#include<string.h>
#include<iostream>
#include<stdlib.h>
#include<algorithm>
using namespace std;
int ans[100000005];
int last_p[100000005];
const long long n=100000000ll*100000000;
long long final_ans[15];
const int modo=1000000007;
int c(int n,int m)
{
    long long i;
    long long p=1;
    for (i=1;i<=n;i++)
    {
        p*=i;
    }
    for (i=1;i<=m;i++)
    {
        p/=i;
    }
    for (i=1;i<=n-m;i++)
    {
        p/=i;
    }
    return p;
}
const int m=100000000;
int main()
{
    memset(ans,0,sizeof(ans));
    int i,j;
    ans[1]=0;
    for (i=2;i<=m;i++)
    {
        if (ans[i]==-1) continue;
        if (i<=10000)
        {
            for (j=i*i;j<=m;j+=i)
            {
                ans[j]=-1;
                last_p[j]=i;
            }
        }
    }
    for (i=2;i<=m;i++)
    {
        if (ans[i]==0)
        {
            ans[i]=1;
        }
        else
        {
            if ((i/last_p[i])%last_p[i]==0)
            {
                ans[i]=-10000;
            }
            else
            {
                ans[i]=ans[i/last_p[i]]+1;
            }
        }
    }
    for (i=1;i<=m;i++)
    {
        if (ans[i]<0) continue;
        //cout<<ans[i]<<" ";
        final_ans[ans[i]]+=(n/i/i);
    }
    //cout<<endl;
    for (i=14;i>=0;i--)
    {
        int j;
        for (j=14;j>i;j--)
        {
            final_ans[i]-=final_ans[j]*c(j,i);
        }
    }
    long long p=1;
    for (i=0;i<=14;i++)
    {
        cout<<final_ans[i]<<endl;
        if (final_ans[i]!=0)
        {
            p=(long long)p*(final_ans[i]%modo)%modo;
        }
    }
    cout<<endl<<endl<<p<<endl;
    system("pause");
    return 0;
}

然后633是一个近似的过程

由于要我们求7的,我们就求7+的

显然我们已经不满足于10¹⁶这么小的范围,那只好扩大

把程序改成搜索,开始爆搜

加个剪枝还是挺快的

质数可以少选一些,后面的不要了

跑出最终结果需要一段时间,还要耐心等待...

反正最后跑出来结果以后我们还可以再稍微加一加试一试对不对,我们精度4位就够,题目也只要求5位

#include<set>
#include<map>
#include<list>
#include<queue>
#include<stack>
#include<string>
#include<time.h>
#include<math.h>
#include<memory>
#include<vector>
#include<bitset>
#include<fstream>
#include<stdio.h>
#include<utility>
#include<sstream>
#include<string.h>
#include<iostream>
#include<stdlib.h>
#include<algorithm>
using namespace std;
const long long m=1000000000000ll;
int c(int n,int m)
{
    long long i;
    long long p=1;
    for (i=1;i<=n;i++)
    {
        p*=i;
    }
    for (i=1;i<=m;i++)
    {
        p/=i;
    }
    for (i=1;i<=n-m;i++)
    {
        p/=i;
    }
    return p;
}
bool prime[4000005];
int p[500005];
int cnt=0;
double ans[15];
double calc(long long x,int y,int c=0)
{
    if (x>m)
    {
        return 0;
    }
    if (y==0)
    {
        return (long long)((double)m*m/x/x);
    }
    int i;
    double sum=0;
    for (i=c;i<cnt;i++)
    {
        double k=calc(x*p[i],y-1,i+1);
        if (k==0)
        {
            return sum;
        }
        sum+=k;
    }
    return sum;
}
int main()
{
    memset(prime,true,sizeof(prime));
    int i,j;
    ans[1]=0;
    for (i=2;i<=10000;i++)
    {
        if (prime[i])
        {
            for (j=i*i;j<=4000000;j+=i)
            {
                prime[j]=false;
            }
        }
    }
    for (i=2;i<=4000000;i++)
    {
        if (prime[i])
        {
            p[cnt++]=i;
        }
    }
    ans[10]=max(0.0,calc(1,10));
    ans[9]=max(0.0,calc(1,9));
    ans[8]=calc(1,8);
    ans[7]=calc(1,7);
    ans[9]-=c(10,9)*ans[10];
    ans[8]-=c(10,8)*ans[10];
    ans[7]-=c(10,7)*ans[10];
    ans[8]-=c(9,8)*ans[9];
    ans[7]-=c(9,7)*ans[9];
    ans[7]-=c(8,7)*ans[8];
    cout<<ans[7]<<endl;
    system("pause");
    return 0;
}