首先，我們定義整數開平方為非負整數對映至非負整數的函式：可利用乘法線性搜尋或二分搜尋，得到最大而平方不超過 的根。通過完全平方數（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 反饋
