整数的计算方法

首先，我们定义整数开平方为非负整数映射至非负整数的函数：可利用乘法线性搜寻或二分搜寻，得到最大而平方不超过 的根。通过完全平方数（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 <typename T> T isqrt(const T& n) { T remainder{}， root{} auto bitCount = isqrt_traits<T>::bitCount(n) for (size_t i = bitCount i > 0 ) { i -= 2 root <<= 1 ++root remainder <<= 2 remainder |= isqrt_traits<T>::extractTwoBitsAt(n， i) if (root <= remainder) { remainder -= root ++root } else --root } return root >>= 1 } T 需要支持 <<=、>>=、<=、前置++、前置--、|= uint8_t，还需要提供一个 isqrt_traits<T> 去抽像两个额外操作，对于内建的无符号整数类型，它的通用 isqrt_traits 是这样的：template <typename T> struct isqrt_traits { static_assert(std::is_unsigned<T>::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<uint8_t>((n >> i) & 3) } } 在 isqrt2 的每个迭代中，我们是通过移位来取得 的两个位，而在 isqrt<T> 中，我们用 extractTwoBitsAt(n， i) 去取得第 i+1 和 第 i 位。这种改动是因为大整数中可直接取得某个位，而不需另外复制一个大整数来做移位操作。这里的 bitCount() 其实可简单返回 sizeof(T) * 8，但这里加上简单优化，循环找出最高的非零两位。然后，我们只需要设计一个支持上述操作的大整数类型，以 std::vector<U> 储存，U 一般可设置为 uint32_t 或 uint64_t，并加入十六进制流输出：template <typename U> class biguint { public: biguint() : v{0} {} biguint(std::initializer_list<U> 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) _back(inBits) return *this } biguint& operator>>=(size_t shift) { assert(shift <= unitBitCount) U inBits = 0 for (auto itr = in() itr != () ++itr) { U outBits = *itr << (unitBitCount - shift) *itr = (*itr >> shift) | inBits inBits = outBits } if (() == 0) _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 < () i++) { U r = i < () ? rhs.v[i] : 0 U previous = v[i] v[i] -= r + inBorrow inBorrow = v[i] > previous ? 1 : 0 } assert(inBorrow == 0) while (() > 1 && () == 0) _back() return *this } biguint& operator++() { for (auto& x : v) if (++x != 0) return *this _back(1) return *this } biguint& operator--() { assert(!(() == 1 && v[0] == 0)) // non-zero for (auto& x : v) if (x-- != 0) return *this return *this } bool operator<=(const biguint& rhs) const { if (() == ()) { for (auto i = () i-- > 0 ) if (v[i] < rhs.v[i]) return true else if (v[i] > rhs.v[i]) return false return true } else return () < () } friend std::ostream& operator<<(std::ostream& os， const biguint& x) { auto f(s()) os << "0x" << std::hex for (auto itr = in() itr != () ++itr) os << *itr s(f) return os } friend struct isqrt_traits<biguint> private: static const size_t unitBitCount = sizeof(U) * 8 std::vector<U> v } 并为 biguint<U> 提供一个 isqrt_traits：template<typename U> struct isqrt_traits<biguint<U>> { static size_t bitCount(const biguint<U>& n) { return biguint<U>::unitBitCount * (() - 1) + isqrt_traits<U>::bitCount(()) } static uint8_t extractTwoBitsAt(const biguint<U>& n， size_t i) { return static_cast<uint8_t>((n.v[i / biguint<U>::unitBitCount] >> (i % biguint<U>::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<unsigned>(45765) << std::endl // 50! = 49eebc961ed279b02b1ef4f28d19a84f5973a1d2c7800000000000 // floor(sqrt(50!)) = 899310e94a8b185249821ebce70 std::cout << isqrt(biguint<uint32_t>{0x00000000， 0xd2c78000， 0x4f5973a1， 0xf28d19a8， 0xb02b1ef4， 0x961ed279， 0x49eebc}) << std::endl } 输出$ g++ -std=c++11 -o isqrt && ./isqrt 213 213 213 0x899310e94a8b185249821ebce7050! 开平方的结果和 (sqrt(50!))+in+hex 吻合（知乎插入 URL 有 bug）。原整代码在 Big integer square root ?· GitHub注意：未经完整测试。---更新1：按@算海无涯 的提示，时间复杂度的次序应为---更新2：isqrt0() 之前有錯，謝 @LOOP 反馈