Projuect Euler 23

Non-abundant sums
Problem 23
A perfect number is a number for which the sum of its proper divisors is exactly equal to the number. For example, the sum of the proper divisors of 28 would be 1 + 2 + 4 + 7 + 14 = 28, which means that 28 is a perfect number.

A number n is called deficient if the sum of its proper divisors is less than n and it is called abundant if this sum exceeds n.

As 12 is the smallest abundant number, 1 + 2 + 3 + 4 + 6 = 16, the smallest number that can be written as the sum of two abundant numbers is 24. By mathematical analysis, it can be shown that all integers greater than 28123 can be written as the sum of two abundant numbers. However, this upper limit cannot be reduced any further by analysis even though it is known that the greatest number that cannot be expressed as the sum of two abundant numbers is less than this limit.

Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers.

Answer:
4179871
Completed on Fri, 30 Jun 2017, 06:40

新学约数和定理：
http://baike.baidu.com/link?url=_TZlbgFVQKEM8lxO3nM1cac_ewzsS0GqOGzrzNF2IaOaZcokgu__8q13zXU3-e8arXxalwntNek8ktPw6-450KIys3nF1n1pRPOrx-HLqBBF_nMn5UViMAEkyvmSNimZJbnGioJejvGOY-FrGvDbXK

#include<bits/stdc++.h>
using namespace std;
#define maxn 50005
int is_prime[maxn];  //素数标记数组，开始全为0，之后变为每个数最小素因子幂积
vector<int>prime;    //里面放的是素数
int abt[maxn];       //每个数的约数和
vector<int>num;  
void get_abount()   //线性筛模板改写为求素数和模板
{
    memset(is_prime,0,sizeof(is_prime));
    for(int i=2;i<maxn;i++)
    {
        if(!is_prime[i]) {      //如果是素数
            is_prime[i]=i;      //最小素数幂积就是本身
            abt[i]=i+1;         //因数和就是1+本身
            prime.push_back(i); 
        }
        for(int j=0;i*prime[j]<maxn;j++)
        {
           // cout<<j<<endl;
            if(i%prime[j]==0) {  //如果不互质，prime[j]一定是i的最小质因数，那么可以根据abt[i]推出abt[i*prime[j]]，公式太复杂，见百度
                is_prime[i*prime[j]]=is_prime[i]*prime[j];
                 abt[i*prime[j]]=abt[i]*(is_prime[i]*prime[j]*prime[j]-1)/(is_prime[i]*prime[j]-1);
                 if(abt[i*prime[j]]>2*i*prime[j])
                  num.push_back(i*prime[j]);
                break;
            }
            else {          //如果互质直接相乘
                is_prime[i*prime[j]]=prime[j];
                abt[i*prime[j]]=abt[i]*abt[prime[j]];
                if(abt[i*prime[j]]>2*i*prime[j]) {
                    num.push_back(i*prime[j]);
                }
            }
        }
    }
}
bool ans[maxn];
int main()
{
    get_abount();
    sort(num.begin(),num.end());
    memset(ans,0,sizeof(ans));
    for(int i=0; ;i++)
    {
        if(num[i]*2>28123) break;
        for(int j=i; ;j++)
        {
            if(num[i]+num[j]>28123)
                break;
            ans[num[i]+num[j]]=1;
        }
    }
    int sum=0;
    for(int i=1;i<=28123;i++)
    {
        if(!ans[i]) sum+=i;
    }
    cout<<sum<<endl;
}