023.cpp

Go to the documentation of this file.

#include <iostream>

unsigned
factor (unsigned *N, unsigned k0 = 3)   //destructive, find smallest prime factor only
{
        if (!(*N & 1)) {        //multiple of 2
                *N >>= 1;                       //faster than *N/=2
                return 2;
        }

        for (unsigned k = k0; k <= *N; k++)
                if (*N % k == 0) {
                        *N /= k;
                        return k;
                }
        return 1;
}

unsigned
term (unsigned p, unsigned k)   //calculate p^0+p^1+...+p^k
{
        unsigned t = 1;
        for (unsigned i = 1; i <= k; i++) {
                t = t * p + 1;
        }
        return t;
}

//sum of all divisors of n except n
//copied from 021.cpp
unsigned
sum_of_divisors (unsigned N)
{
        unsigned answer = 1;

        unsigned previous_prime = 2;    //modified from 012.cpp
        unsigned power = 0;

        unsigned r = N;
        unsigned div = 2;
        do {
                div = factor (&r, div);
                if (div == previous_prime)
                        power++;
                else {          //determined the power of the previous prime factor
                        answer *= term (previous_prime, power);

                        previous_prime = div;
                        power = 1;
                }
                //cout << div<<','<<answer<<' ';
        } while (div > 1);

        answer -= N;                    //exclude N itself from the sum
        return answer;
}

const unsigned max = 28123;

class abundant
{
        bool array[max+1];
public:
        abundant() {
                for (unsigned i=0;i<12;i++) {
                        array[i]=0;
                }
                for (unsigned i=12;i<=max;i++) {
                        if (sum_of_divisors(i)>i)
                                array[i]=1;
                        else
                                array[i]=0;
                }
        }
        unsigned search_back(unsigned x) {
                unsigned b;
                for (b=x-1;!array[b] && b>=12;b--);
                return b;
        }
        bool exist(unsigned x) {
                return array[x];
        }
};

int main()
{
        abundant abundant_nos;

        unsigned long sum = 0;

        //any numbers smaller than 24 cannot be expressed as a sum of two abundant numbers
        for (unsigned i=1;i<=23;i++)
                sum+=i;

        for (unsigned i=24;i<=max;i++) {
                bool flag = 0;
                //search until b<12
                for (unsigned b=abundant_nos.search_back(i);
                        b>=12;b=abundant_nos.search_back(b)) {
                        //find a = i-b
                        unsigned a=i-b;

                        //if i is the sum of a,b
                        if (abundant_nos.exist(a)) {
                                //std::cout << i << '=' << a<<'+'<<b<<std::endl;
                                flag=1;
                                break;
                        }
                        if (flag)break;
                }
                //i cannot be written as the sum of two abundant numbers.
                if (!flag) {
                        //std::cout << i << ' ';
                        sum+=i;
                }
        }

        std::cout << sum << std::endl;

        return 0;
}