021.cpp

Go to the documentation of this file.

#include <iostream>
using namespace std;

/* Find smallest prime factor*/
unsigned
factor (unsigned *N, unsigned k0 = 3)   
{
  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;
}

/*calculate p^0+p^1+...+p^k*/
unsigned
term (unsigned p, unsigned 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 */
unsigned
d (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;
}

/*memorize impossible integers*/
class lookup_table
{
  bool x[10000];
public:
    lookup_table ()
  {
    bool *p = x;
    for (unsigned i = 0; i < 10000; i++, p++)
       *p = 0;
  }
  bool exception (unsigned k)
  {
    return x[k];
  }
  void push (unsigned k)
  {
    x[k] = 1;
  }
};

int
main ()
{

  unsigned Sum = 0;
  lookup_table table;

  for (unsigned a = 3; a < 10000; a++)
    {
      if (table.exception (a))
        continue;
      unsigned b = d (a);
      if (b == a || b >= 10000)
        continue;
      if (d (b) == a)
        {
          Sum += a + b;
          table.push (b);       //so b will never appear again
          //cout << a << ' '<<b << ' ';
        }
      else if (b > a && b < 10000)      //add an exception to reduce computation time
        table.push (b);
    }

  cout << Sum << endl;
  return 0;
}