计算catlan数时无法找出错误

Not able to figure out the mistake while calculating catlan number

在这段代码中,我使用 catlan 数计算唯一二叉搜索树的数量。 答案一直是正确的,直到输入值 n=13 来自 n>=14 答案被证明少了一个。例如。对于 n=14,我的答案是 2674439,而实际答案是 2674440。 这是代码。

#include "stdafx.h"
#include "iostream"
using namespace std;
class Solution {
public:
    double arr[20000] = {};
    double fact(int n) {
        if (n == 1 || n == 0)
          return 1;
        if (arr[n] != 0)
          return arr[n];
        return arr[n] = ceil(n*fact(n - 1));
}

int numTrees(int n) {
    int res = fact(2 * n) / ((fact(n + 1))*fact(n));
    return res;
}
};


int main()
{    
  Solution s;
  cout << s.numTrees(14)<<endl;
  return 0;
}

这个函数有问题:

int numTrees(int n) {
   int res = fact(2 * n) / ((fact(n + 1))*fact(n));
   return res;
}

基本上,将 double 值转换为 int 会损失一些精度,因此无法得出正确的值。如果你把它改成double,问题就消失了。

更正函数:

double numTrees(int n) 
{
   double res = fact(2 * n) / ((fact(n + 1))*fact(n));
   return res;
}

你的一个中间值,28!需要 98 位精度。

双精度有 52-53 位精度。

令人惊讶的部分不是在 14 处出现错误,而是在 14 之前没有出现错误。这是因为 采取了一些努力 double 来减少累积误差,基本就走运了。

在这种情况下,我们用乘法做很多数学运算,紧接着 none 用加法做数学运算。与主要力量一起工作是一个很好的举动:

struct product {
  std::map<std::size_t, std::ptrdiff_t> powers;
  product& operator*=( product const& rhs ) {
    for (auto&& e:rhs.powers)
      powers[e.first] += e.second;
    return tidy(*this);
  }
  product& operator/=( product const& rhs ) {
    for (auto&& e:rhs.powers)
      powers[e.first] -= e.second;
    return tidy(*this);
  }
  friend product operator*( product lhs, product const& rhs ) {
    lhs *= rhs;
    return lhs;
  }
  friend product operator/( product lhs, product const& rhs ) {
    lhs /= rhs;
    return lhs;
  }
  // 1/x overload:
  friend product operator~( product p ) {
    for (auto& e:p.powers)
      e.second = -e.second;
    return p;
  }
  product() = default;
  product(product const&) = default;
  product(product &&) = default;
  product& operator=(product const&) = default;
  product& operator=(product &&) = default;

  product( std::size_t in ); // TODO

  bool is_integral() const {
    for (auto& e:powers)
      if (e.second < 0) return false;
    return true;
  }
  template<class Scalar=std::size_t>
  Scalar numerator() const {
    Scalar r = 1;
    for( auto& e: powers )
      for (std::ptrdiff_t i = 0; i < e.second; ++i)
        r *= e.first;
    return r;
  }
  template<class Scalar=std::size_t>
  Scalar denominator() const {
    Scalar r = 1;
    for( auto& e: powers )
      for (std::ptrdiff_t i = 0; i > e.second; --i)
        r *= e.first;
    return r;
  }
  friend product& tidy(product& p) {
    for (auto it = p.powers.begin(); it != p.powers.end();) {
      if (!it->second)
        it = p.powers.erase(it);
      else
        ++it;
    }
    return p;
  }
};

这是一个小阶乘引擎:

struct factorial_t {
  std::vector<product> values;
  factorial_t():values(2) {}
  product const& operator()( std::size_t in ) {
    if (values.size() > in) {
      return values[in];
    }
    values.push_back( (*this)(in-1)*product{in} );
    return values.back();
  }
};
factorial_t factorial;

这对于荒谬的值来说是非常精确的。

numTrees 则变为:

template<class Scalar=std::size_t>
Scalar numTrees(std::size_t n) {
  auto res = factorial(2 * n) / ((factorial(n + 1))*factorial(n));
  return res.numerator<Scalar>();
}

剩下要做的就是编写质因数 a std::size_t

struct factor_t {
    std::vector<std::size_t> primes;
    factor_t():primes{2,3,5,7,11,13,17} {}

    bool make_primes( std::size_t up_to ) {
        if ((primes.back()+2) > up_to)
            return false;
        bool added_prime = false;
        for (std::size_t x = primes.back()+2; x < up_to; x += 2) {
            bool found = false;
            for (auto p:primes)
            {
                if (p*p > x) break;
                if (x%p) continue;
                found = true;
                break;
            }
            if (found)
                primes.push_back(x);
            added_prime = added_prime || found;
        }
        return added_prime;
    }
    product operator()( std::size_t in ) {
        product r;
        for (auto&& prime:primes)
        {
            while (!(in%prime)) {
                r.powers[prime]++;
                in /= prime;
            }
        }
        // are there more primes to apply?
        if (make_primes(std::sqrt(in)))
        {
            r *= (*this)(in);
        }
        else if (in != 1)
        {
            // in is a prime
            r.powers[in]++;
        }
        return r;
    }
};
factor_t factor;

product::product( std::size_t in ):
  product(factor(in))
{}

bob is your uncle.

(有趣的是,我写的产品代码 "works" 不小心用了 0,因为它 factor0 误认为是质数,而 product 具有正 0 因子的 .numerator()0。如果你除法,你最终会在分母中得到 0,如果它们取消我假装它永远不会会发生。但是,一般来说,不要给产品 class 一个 0。)