整数的计算方法？
首先，我们定义整数开平方为非负整数映射至非负整数的函数：可利用乘法线性搜寻或二分搜寻，得到最大而平方不超过 的根。通过完全平方数（square number）数列，我们还可以在线性搜寻中只用加法，因为两个完全平方数的差为奇数数列：uint32_t isqrt0(uint32_t n) { uint32_t delta = 3; for (uint32_t square = 1; square <= n; delta += 2) square += delta; return delta / 2 - 1; } 因为问题是关于大整数的，我们要把大整数的位数（）也考虑在内。线性搜寻需要次迭代，每次迭代的加法需时间，合共。而二分搜寻最坏情况需要次迭代，每次的乘法需时间，合共。而一些数值方法（如牛顿迭代）只适合计算近似值，而且当中也涉及除法。我们换一个思路，参考 Integer Square Roots 这篇文章，开平方可以用类似长除法的方式计算，在二进制中只需要用比较和减法，32位无号整数的C实现如下：uint32_t isqrt1(uint32_t n) { uint32_t remainder = 0, root = 0, divisor; for (size_t i = 0; i < 16; i++) { root <<= 1; remainder <<= 2; remainder |= n >> 30; n <<= 2; // Extract 2 MSB from n divisor = (root << 1) + 1; if (divisor <= remainder) { remainder -= divisor; ++root; } } return root; }这个方法的迭代次数是 次（整数有多少位） ，每次迭代的加法、减法、移位、比较都是，合共时间，时间复杂度比线性和二分搜寻都要低。由于 divisor 和 root 的关系是固定的，如果空间是考虑因素（考虑到大整数或硬件实现），可以改为这种形式，省下 divisor 的存储：uint32_t isqrt2(uint32_t n) { uint32_t remainder = 0, root = 0; for (size_t i = 0; i < 16; i++) { root <<= 1; ++root; remainder <<= 2; remainder |= n >> 30; n <<= 2; // Extract 2 MSB from n if (root <= remainder) { remainder -= root; ++root; } else --root; } return root >>= 1; }接下来，我们把这算法实现写成 C++11 泛形形式，接受任何无号整数类型：template T isqrt(const T& n) { T remainder{}, root{}; auto bitCount = isqrt_traits::bitCount(n); for (size_t i = bitCount; i > 0; ) { i -= 2; root <<= 1; ++root; remainder <<= 2; remainder |= isqrt_traits::extractTwoBitsAt(n, i); if (root <= remainder) { remainder -= root; ++root; } else --root; } return root >>= 1; } T 需要支持 <<=、>>=、<=、前置++、前置--、|= uint8_t，还需要提供一个 isqrt_traits 去抽像两个额外操作，对于内建的无符号整数类型，它的通用 isqrt_traits 是这样的：template struct isqrt_traits { static_assert(std::is_unsigned::value, "generic isqrt only on unsigned types"); // Number of bits in multiples of two static size_t bitCount(const T& n) { T a(n); size_t count = 0; while (a > 0) { a >>= 2; count += 2; } return count; } // Extract the i+1, i bits static uint8_t extractTwoBitsAt(const T& n, size_t i) { return static_cast((n >> i) & 3); } }; 在 isqrt2 的每个迭代中，我们是通过移位来取得 的两个位，而在 isqrt 中，我们用 extractTwoBitsAt(n, i) 去取得第 i+1 和 第 i 位。这种改动是因为大整数中可直接取得某个位，而不需另外复制一个大整数来做移位操作。这里的 bitCount() 其实可简单返回 sizeof(T) * 8，但这里加上简单优化，循环找出最高的非零两位。然后，我们只需要设计一个支持上述操作的大整数类型，以 std::vector 储存，U 一般可设置为 uint32_t 或 uint64_t，并加入十六进制流输出：template class biguint { public: biguint() : v{0} {} biguint(std::initializer_list init) : v(init) {} biguint& operator<<=(size_t shift) { assert(shift <= unitBitCount); U inBits = 0; for (auto& x : v) { U outBits = x >> (unitBitCount - shift); x = (x << shift) | inBits; inBits = outBits; } if (inBits) v.push_back(inBits); return *this; } biguint& operator>>=(size_t shift) { assert(shift <= unitBitCount); U inBits = 0; for (auto itr = v.rbegin(); itr != v.rend(); ++itr) { U outBits = *itr << (unitBitCount - shift); *itr = (*itr >> shift) | inBits; inBits = outBits; } if (v.back() == 0) v.pop_back(); return *this; } biguint& operator|=(uint8_t rhs) { v[0] |= rhs; return *this; } biguint& operator-=(const biguint& rhs) { assert(rhs <= *this); U inBorrow = 0; for (size_t i = 0; i < v.size(); i++) { U r = i < rhs.v.size() ? rhs.v[i] : 0; U previous = v[i]; v[i] -= r + inBorrow; inBorrow = v[i] > previous ? 1 : 0; } assert(inBorrow == 0); while (v.size() > 1 && v.back() == 0) v.pop_back(); return *this; } biguint& operator++() { for (auto& x : v) if (++x != 0) return *this; v.push_back(1); return *this; } biguint& operator--() { assert(!(v.size() == 1 && v[0] == 0)); // non-zero for (auto& x : v) if (x-- != 0) return *this; return *this; } bool operator<=(const biguint& rhs) const { if (v.size() == rhs.v.size()) { for (auto i = v.size(); i-- > 0; ) if (v[i] < rhs.v[i]) return true; else if (v[i] > rhs.v[i]) return false; return true; } else return v.size() < rhs.v.size(); } friend std::ostream& operator<<(std::ostream& os, const biguint& x) { auto f(os.flags()); os << "0x" << std::hex; for (auto itr = x.v.rbegin(); itr != x.v.rend(); ++itr) os << *itr; os.flags(f); return os; } friend struct isqrt_traits; private: static const size_t unitBitCount = sizeof(U) * 8; std::vector v; }; 并为 biguint 提供一个 isqrt_traits：template struct isqrt_traits> { static size_t bitCount(const biguint& n) { return biguint::unitBitCount * (n.v.size() - 1) + isqrt_traits::bitCount(n.v.back()); } static uint8_t extractTwoBitsAt(const biguint& n, size_t i) { return static_cast((n.v[i / biguint::unitBitCount] >> (i % biguint::unitBitCount)) & 3); } }; 我简单测试了一下 45765 和 50! 的开平方：int main() { // floor(sqrt(45765)) = 213 std::cout << isqrt1(45765) << std::endl; std::cout << isqrt2(45765) << std::endl; std::cout << isqrt(45765) << std::endl; // 50! = 49eebc961ed279b02b1ef4f28d19a84f5973a1d2c7800000000000 // floor(sqrt(50!)) = 899310e94a8b185249821ebce70 std::cout << isqrt(biguint{0x00000000, 0xd2c78000, 0x4f5973a1, 0xf28d19a8, 0xb02b1ef4, 0x961ed279, 0x49eebc}) << std::endl; } 输出$ g++ -std=c++11 -o isqrt isqrt.cpp && ./isqrt 213 213 213 0x899310e94a8b185249821ebce7050! 开平方的结果和 https://www.wolframalpha.com/input/?i=floor(sqrt(50!))+in+hex 吻合（知乎插入 URL 有 bug）。原整代码在 Big integer square root ?· GitHub注意：未经完整测试。---更新1：按@算海无涯 的提示，时间复杂度的次序应为---更新2：isqrt0() 之前有錯，謝 @LOOP 反馈
