對大整數N開平方，求其整數部分，有什麼好的演算法嗎？

12-30

N&>&>2^64，很大的整數。

首先，我們定義整數開平方為非負整數映射至非負整數的函數：

$mathrm{isqrt} : mathbb{Z}^* ightarrow mathbb{Z}^*\ n mapsto lfloor sqrt{n} floor$

可利用乘法線性搜尋或二分搜尋，得到最大而平方不超過 $n$ 的根。通過完全平方數（square number）數列，我們還可以在線性搜尋中只用加法，因為兩個完全平方數的差為奇數數列：

$Delta_{i+1} = (i + 1)^2 - i^2 = i^2 + 2i + 1 - i^2 = 2i + 1 = Delta_i + 2$

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; }

因為問題是關於大整數的，我們要把大整數的位數（ $log n$ ）也考慮在內。線性搜尋需要 $lceil sqrt{n} , ceil$ 次迭代，每次迭代的加法需 $O(log n)$ 時間，合共 $O(sqrt{n}log n)$ 。而二分搜尋最壞情況需要 $log n$ 次迭代，每次的乘法需 $O(log^2 n)$ 時間，合共 $O(log^3 n)$ 。

而一些數值方法（如牛頓迭代）只適合計算近似值，而且當中也涉及除法。

我們換一個思路，參考 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; }

這個方法的迭代次數是 $log n$ 次（整數有多少位），每次迭代的加法、減法、移位、比較都是 $O(log n)$ ，合共 $O(log^2 n)$ 時間，時間複雜度比線性和二分搜尋都要低。

由於 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 的每個迭代中，我們是通過移位來取得 $n$ 的兩個位，而在 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 0x899310e94a8b185249821ebce70

50! 開平方的結果和

https://www.wolframalpha.com/input/?i=floor(sqrt(50!))+in+hex

吻合（知乎插入 URL 有 bug）。

原整代碼在 Big integer square root ?· GitHub

注意：未經完整測試。

---

更新1：按@算海無涯的提示，時間複雜度的次序應為

$log^2 n prec log^3 n prec sqrt{n}log n$

---

更新2：isqrt0() 之前有錯，謝 @LOOP 反饋

repo/gmp-6.1: e225fb40ca63 mpn/generic/rootrem.c

1. 整數的平方間隔越來越大：

1，4，9，16，25，。。。

所以當n比較大的時候，n^2 和（n+1）^2之間間隔應該很大，而這之間的數的平方根的整數部分都是n。

2. 應該可以採用newton 法Newton"s method，不過你可以放寬收斂的條件，因為你找的是整數部分。

跑個題。

int l;int main(int o,char **O, int I){char c,*D=O[1];if(o&>0){ for(l=0;D[l ];D[l ++]-=10){D [l++]-=120;D[l]-= 110;while (!main(0,O,l))D[l] += 20; putchar((D[l]+1032) /20 ) ;}putchar(10);}else{ c=o+ (D[I]+82)%10-(I&>l/2)* (D[I-l+I]+72)/10-9;D[I]+=I&<0?0 :!(o=main(c/10,O,I-1))*((c+999 )%10-(D[I]+92)%10);}return o;}

來自IOCCC 2001年的獲獎作品。

摘下眼鏡，風味更佳。

大整數開方 | Tom smith"s

NOIP2011初賽題，用二分做

看到這題腦袋裡蹦出的是NOIP2011（？可能是，不太確定年份了）提高組初賽試題……

你見過哪些神一樣的優化操作？ - 高天的回答

float fast_sqrt(float x) { uint32_t x_bits = 0;
x_bits = *((uint32_t*) x); //猜猜這是什麼？ x_bits = (x_bits &>&> 1) + 532369198; return *((float*) x_bits); }
作者：高天鏈接：https://www.zhihu.com/question/44213758/answer/97229375 來源：知乎著作權歸作者所有。商業轉載請聯繫作者獲得授權，非商業轉載請註明出處。

各位大豬所說的都是在迭代基礎上，只有使用計算機才能解出的方式，有一種使用二次函數手動計算出解的方法，想必用二次函數原型編程代碼更少更簡單，查遍了百度以及各大博客也沒見到這個表達式子。可惜了了

一般來說，在乘法複雜度為nlogn(n為位長)前提下

開方的複雜度為上限至少可以低至nlogn*logn

方法是牛頓迭代，

具體是先牛頓迭代計算平方根倒數，1/sqrt(N)=r.

有效精度不少於平方根的位長k位，

計算r*N，同樣有效精度不少於平方根位長k位，

則r*N的前k位就是平方根的近似值，誤差只要有效精度不少於k，則誤差至多為1。

其中牛頓迭代時，精度是倍增的，實際計算時可以只保留不少於當前有效精度即可提高速度！

實測C++代碼，開方的速度大約是乘法的1.5至3倍的時間，

平方根倒數的方法其它答主有提供鏈接

牛頓迭代，複雜度可以做到lognloglogn

參考折半查找演算法

我的直觀想法：平方根本質和除法一樣，如果只求整數部分，對於n位數只需要O= $frac{n}{2}+1$ 次運算（一次運算指一個乘法和一個減法）。

實現比二分法複雜，但運算次數永遠是二分法(lg n / lg2)的一半。但二分法也不慢，例如2^1024隻需要300次運算。

a=sqrt(n)

print int(a)

vb編程

牛頓迭代或者2分時間複雜度應該都是logn的

個人感覺2分好寫一點

我來說個C++實現的思路：

1.數據結構使用bitset和string,方便大數的表示和顯示;

2.演算法採用二分搜索;

3.搜索範圍依據。舉個咧子，10 * 10 = 100；100 * 100 = 10000；1000 * 1000 = 1000000...100和10000之間的值在10和100之間搜;10000和1000000之間的值在100和1000之間搜。二進位表示的數，我想也有類似規律。

附：感覺有開源代碼，大神們應該早就已經實現了

int sqrt(int x) { // write your code here if(x &<= 0 ) return 0; int left = 1, right = x; int middle; while(1) { middle = left + (right - left) / 2; if((middle * middle) == x) return middle; if(middle &> x / middle) { // 防止溢出啊！ right = middle; } else if((middle &<= (x / middle)) (right - middle) == 1) { return middle; } else if((middle &<= (x / middle)) (right - middle) &> 1) { left = middle; } } }

不是大整數求根，純粹是插隊雷神之錘裡面用到的快速演算法

https://zh.m.wikipedia.org/wiki/平方根倒數速演算法

https://en.m.wikipedia.org/wiki/Fast_inverse_square_root