C++ version of http://www.zhihu.com/question/27714609

beal.cc

#include <algorithm>
#include <map>
#include <string>
#include <vector>
#include <ext/numeric>
#include <assert.h>
#include <stdint.h>

class UnsignedInt // copyable
{
 public:
  enum Radix
  {
    kDec,
    kHex,
  };

  UnsignedInt(uint32_t x = 0)
  {
    if (x > 0)
    {
      limbs_.push_back(x);
    }
  }

  UnsignedInt(const std::string& x, Radix r = kDec);

  std::string toHex() const;
  std::string toDec() const;
  bool isZero() const { return limbs_.empty(); }
  bool isNormal() const { return isZero() || limbs_.back() != 0; }
  bool lessThan(const UnsignedInt& x) const
  {
    const value_type& rhs = x.limbs_;
    return rhs.size() > limbs_.size() ||
           (rhs.size() == limbs_.size() &&
            std::lexicographical_compare(limbs_.rbegin(),
                                         limbs_.rend(),
                                         rhs.rbegin(),
                                         rhs.rend()));
  }
  bool operator<(const UnsignedInt& x) const
  {
    return lessThan(x);
  }

  UnsignedInt operator+(const UnsignedInt& x) const
  {
    UnsignedInt tmp(*this);
    tmp.add(x);
    return tmp;
  }

  void swap(UnsignedInt& rhs) { limbs_.swap(rhs.limbs_); }

  void assign(const uint32_t x)
  {
    if (x == 0)
    {
      limbs_.clear();
    }
    else
    {
      limbs_.resize(1);
      limbs_[0] = x;
    }
  }

  void add(const uint32_t x)
  {
    if (limbs_.empty())
    {
      limbs_.push_back(x);
      return;
    }

    uint64_t sum = limbs_[0] + static_cast<uint64_t>(x);
    limbs_[0] = sum & kMask_;
    uint64_t carry = sum > kMask_;

    for (size_t i = 1; i < limbs_.size() && carry; ++i)
    {
      sum = limbs_[i] + carry;
      limbs_[i] = sum & kMask_;
      carry = sum > kMask_;
    }

    if (carry)
    {
      limbs_.push_back(carry);
    }
  }

  void add(const UnsignedInt& x)
  {
    const value_type& rhs = x.limbs_;
    if (rhs.size() > limbs_.size())
    {
      limbs_.resize(rhs.size());
    }
    size_t len = std::min(limbs_.size(), rhs.size());
    uint64_t carry = 0;
    for (size_t i = 0; i < len; ++i)
    {
      uint64_t sum = limbs_[i] + static_cast<uint64_t>(rhs[i]) + carry;
      limbs_[i] = sum & kMask_;
      carry = sum > kMask_;
    }
    if (carry)
    {
      for (size_t i = len; i < limbs_.size() && carry; ++i)
      {
        uint64_t sum = limbs_[i] + carry;
        limbs_[i] = sum & kMask_;
        carry = sum > kMask_;
      }
    }
    if (carry)
    {
      limbs_.push_back(carry);
    }
  }

  // void sub(const UnsignedInt& x);

  void multiply(const uint32_t x)
  {
    const uint64_t multiplier = x;
    if (multiplier == 0)
    {
      limbs_.clear();
      return;
    }
    uint64_t carry = 0;
    for (size_t i = 0; i < limbs_.size(); ++i)
    {
      uint64_t prod = limbs_[i] * multiplier + carry;
      limbs_[i] = prod & kMask_;
      carry = prod / 0x10000 / 0x10000;
    }
    if (carry)
    {
      limbs_.push_back(carry);
    }
  }

  void multiply(const UnsignedInt& x)
  {
    const value_type& rhs = x.limbs_;
    if (limbs_.empty() || rhs.empty())
    {
      limbs_.clear();
      return;
    }
    value_type result(limbs_.size() + rhs.size());
    for (size_t i = 0; i < rhs.size(); ++i)
    {
      if (rhs[i] == 0)
        continue;
      uint64_t carry = 0;
      for (size_t j = 0; j < limbs_.size(); ++j)
      {
        uint64_t prod = uint64_t(rhs[i]) * limbs_[j]
                        + result[i+j]
                        + carry;
        result[i+j] = prod & kMask_;
        carry = prod / 0x10000 / 0x10000;
      }
      result[i + limbs_.size()] = carry;
    }
    if (result.back() == 0)
      result.pop_back();
    result.swap(limbs_);
  }

  // returns remainder
  uint32_t devide(const uint32_t x)
  {
    uint64_t carry = 0;
    for (size_t i = limbs_.size(); i > 0; --i)
    {
      uint64_t prod = carry << 32;
      prod += limbs_[i-1];

      limbs_[i-1] = prod / x;
      carry = prod % x;
    }

    if (!limbs_.empty() && limbs_.back() == 0)
      limbs_.pop_back();
    assert(isNormal());
    return carry;
  }

  // pow(this, n)
  void power(uint32_t n)
  {
    if (n == 0)
    {
      assign(1);
    }
    else if (n > 1)
    {
      // FIXME: if n is a power of 2, we can do this in-place.
      UnsignedInt result(1);
      while (n)
      {
        if (n & 1)
        {
          result.multiply(*this);
        }
        multiply(*this);
        n /= 2;
      }
      swap(result);
    }
  }

  typedef std::vector<uint32_t> value_type;

  const value_type& getValue() const
  {
    return limbs_;
  }

  void setValue(int n, uint32_t x)
  {
    limbs_.assign(n, x);
  }

 private:
  void parseDec(const std::string& str);
  void parseHex(const std::string& str);

  uint32_t highest() const
  {
    return limbs_.empty() ? 0 : limbs_.back();
  }

  value_type limbs_;

  const static uint32_t kMask_ = 0xFFFFFFFF;
};

std::string UnsignedInt::toDec() const
{
  const uint32_t segment = 1000000000;
  std::string result;
  if (limbs_.empty())
  {
    result.push_back('0');
    return result;
  }

  result.reserve(9*limbs_.size() + 7*limbs_.size()/11 + 1);
  // log10(2**32) = 32*log10(2) = 9.633 digits per word
  UnsignedInt copy = *this;
  while (copy.limbs_.size() > 1)
  {
    uint32_t x = copy.devide(segment);
    for (int i = 0; i < 9; ++i)
    {
      int lsd = x  % 10;
      x /= 10;
      result.push_back(lsd + '0');
    }
  }

  uint32_t x = copy.limbs_[0];
  do
  {
    int lsd = x % 10;
    x /= 10;
    result.push_back(lsd + '0');
  } while (x != 0);

  std::reverse(result.begin(), result.end());

  return result;
}

/////////////////////////////////////////////////////////////////////////////

const int g_min_base = 2;
const int g_min_power = 3;
int g_max_base = 2500;
int g_max_power = 200;
std::vector<std::vector<UnsignedInt>> g_pow;
std::map<UnsignedInt, int> g_table;

void initial_data()
{
  g_pow.resize(g_max_base+1);
  for (int z = g_min_base; z <= g_max_base; ++z)
  {
    g_pow[z].resize(g_max_power+1);
    UnsignedInt zr(z);
    zr.power(g_min_power);
    for (int r = g_min_power; r <= g_max_power; ++r)
    {
      g_pow[z][r] = zr;
      g_table[zr] = r;
      zr.multiply(z);
    }
  }
  printf("initialized %zd table elements\n", g_table.size());
  /*
  for (int x = 0; x < g_pow.size(); ++x)
  {
    const auto& powx = g_pow[x];
    for (int r = 0; r < powx.size(); ++r)
    {
      printf("%d ** %d = %s\n", x, r, powx[r].toDec().c_str());
    }
  }
  */
}

bool ispower(int a, int b)
{
  while (a > 1)
  {
    if (a % b == 0)
      a /= b;
    else
      return false;
  }
  return true;
}

int gcd(int x, int y)
{
    while (x > 0)
    {
      int old = x;
      x = y % x;
      y = old;
    }
    return y;
}

void report()
{
  abort();
}

void beal0()
{
  for (int x = g_min_base; x <= g_max_base; ++x)
  {
    bool flag2 = false, flag3 = false;
    if (ispower(x, 2))
      flag2 = true;
    else if (ispower(x, 3))
      flag3 = true;
    const auto& powx = g_pow[x];
    if (x % 10 == 0)
      printf("x = %d\n", x);
    for (int y = g_min_base; y <= g_max_base; ++y)
    {
      if (y > x)
        break;
      if (gcd(x, y) > 1)
        continue;
      if (flag2) {
        if (ispower(y, 3))
          continue;
      } else if (flag3 && ispower(y, 2))
        continue;
      const auto& powy = g_pow[y];
      for (int m = g_min_power; m <= g_max_power; ++m)
      {
        if (x == g_max_base && m == g_max_power)
          continue;
        const auto& xm = powx[m];
        for (int n = g_min_power; n <= g_max_power; ++n)
        {
          UnsignedInt sum = xm + powy[n];
          auto it = g_table.find(sum);
          if (it != g_table.end())
            report();
        }
      }
    }
  }
}

int main(int argc, char* argv[])
{
  initial_data();
  beal0();
}