整數的計算方法？

首先，我們定義整數開平方為非負整數映射至非負整數的函數：可利用乘法線性搜尋或二分搜尋，得到最大而平方不超過 的根。通過完全平方數（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) 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<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 * (n.v.size() - 1) + isqrt_traits<U>::bitCount(n.v.back()); } 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.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 反饋