algorithm - 一种计算任意大整数的整数平方根 (isqrt) 的有效算法

Question

通知

有关 Erlang 或 C/C++ 的解决方案，请转到下面的试用 4。

---

维基百科文章

Integer square root

• 可以在这里找到“整数平方根”的定义

Methods of computing square roots

• 可在此处找到执行“位魔术”的算法

---

[ 试用 1:使用库函数 ]

代码

isqrt(N) when erlang:is_integer(N), N >= 0 ->
    erlang:trunc(math:sqrt(N)).

问题

此实现使用 C 库中的 sqrt() 函数，因此它不适用于任意大的整数(请注意，返回的结果与输入不匹配。正确答案应为 12345678901234567890):

Erlang R16B03 (erts-5.10.4) [source] [64-bit] [smp:8:8] [async-threads:10] [hipe] [kernel-poll:false]

Eshell V5.10.4  (abort with ^G)
1> erlang:trunc(math:sqrt(12345678901234567890 * 12345678901234567890)).
12345678901234567168
2> 

---

[ 试验 2:仅使用 Bigint + ]

代码

isqrt2(N) when erlang:is_integer(N), N >= 0 ->
    isqrt2(N, 0, 3, 0).

isqrt2(N, I, _, Result) when I >= N ->
    Result;

isqrt2(N, I, Times, Result) ->
    isqrt2(N, I + Times, Times + 2, Result + 1).

描述

此实现基于以下观察:

isqrt(0) = 0   # <--- One 0
isqrt(1) = 1   # <-+
isqrt(2) = 1   #   |- Three 1's
isqrt(3) = 1   # <-+
isqrt(4) = 2   # <-+
isqrt(5) = 2   #   |
isqrt(6) = 2   #   |- Five 2's
isqrt(7) = 2   #   |
isqrt(8) = 2   # <-+
isqrt(9) = 3   # <-+
isqrt(10) = 3  #   |
isqrt(11) = 3  #   |
isqrt(12) = 3  #   |- Seven 3's
isqrt(13) = 3  #   |
isqrt(14) = 3  #   |
isqrt(15) = 3  # <-+
isqrt(16) = 4  # <--- Nine 4's
...

问题

这个实现涉及仅 bigint 添加，所以我希望它运行得很快。但是，当我用 1111111111111111111111111111111111111111 * 11111111111111111111111111111111111111 喂它时，它似乎永远在我的(非常快的)机器上运行。

---

[ 试验 3:仅对 Bigint +1、-1 和 div 2 使用二进制搜索]

代码

变体 1(我的原始实现)

isqrt3(N) when erlang:is_integer(N), N >= 0 ->
    isqrt3(N, 1, N).

isqrt3(_N, Low, High) when High =:= Low + 1 ->
    Low;

isqrt3(N, Low, High) ->
    Mid = (Low + High) div 2,
    MidSqr = Mid * Mid,
    if
        %% This also catches N = 0 or 1
        MidSqr =:= N ->
            Mid;
        MidSqr < N ->
            isqrt3(N, Mid, High);
        MidSqr > N ->
            isqrt3(N, Low, Mid)
    end.

变体 2(修改以上代码，使边界改为 Mid+1 或 Mid-1，引用 answer by Vikram Bhat )

isqrt3a(N) when erlang:is_integer(N), N >= 0 ->
    isqrt3a(N, 1, N).

isqrt3a(N, Low, High) when Low >= High ->
    HighSqr = High * High,
    if
        HighSqr > N ->
            High - 1;
        HighSqr =< N ->
            High
    end;

isqrt3a(N, Low, High) ->
    Mid = (Low + High) div 2,
    MidSqr = Mid * Mid,
    if
        %% This also catches N = 0 or 1
        MidSqr =:= N ->
            Mid;
        MidSqr < N ->
            isqrt3a(N, Mid + 1, High);
        MidSqr > N ->
            isqrt3a(N, Low, Mid - 1)
    end.

问题

现在它以闪电般的速度解出 79 位数字(即 111111111111111111111111111111111111111 * 11111111111111111111111111111111111111)，结果立即显示。 However, it takes 60 seconds (+- 2 seconds) on my machine to solve one million (1,000,000) 61-digit numbers (namely, from 1000000000000000000000000000000000000000000000000000000000000 to 1000000000000000000000000000000000000000000000000000001000000).我想做得更快。

---

[ 试验 4:将牛顿法与 Bigint + 和 div 结合使用]

代码

isqrt4(0) -> 0;

isqrt4(N) when erlang:is_integer(N), N >= 0 ->
    isqrt4(N, N).

isqrt4(N, Xk) ->
    Xk1 = (Xk + N div Xk) div 2,
    if
        Xk1 >= Xk ->
            Xk;
        Xk1 < Xk ->
            isqrt4(N, Xk1)
    end.

C/C++ 代码(供您感兴趣)

递归变体

#include <stdint.h>

uint32_t isqrt_impl(
    uint64_t const n,
    uint64_t const xk)
{
    uint64_t const xk1 = (xk + n / xk) / 2;
    return (xk1 >= xk) ? xk : isqrt_impl(n, xk1);
}

uint32_t isqrt(uint64_t const n)
{
    if (n == 0) return 0;
    if (n == 18446744073709551615ULL) return 4294967295U;
    return isqrt_impl(n, n);
}

迭代变体

#include <stdint.h>

uint32_t isqrt_iterative(uint64_t const n)
{
    uint64_t xk = n;
    if (n == 0) return 0;
    if (n == 18446744073709551615ULL) return 4294967295U;
    do
    {
        uint64_t const xk1 = (xk + n / xk) / 2;
        if (xk1 >= xk)
        {
            return xk;
        }
        else
        {
            xk = xk1;
        }
    } while (1);
}

问题

Erlang 代码在我的机器上用 40 秒(± 1 秒)解决了一百万 (1,000,000) 个 61 位数字，所以这比 Trial 3 快。它能走得更快吗？

---

关于我的机器

处理器: 3.4 GHz Intel Core i7

内存:32 GB 1600 MHz DDR3

操作系统:Mac OS X 10.9.1 版

---

相关问题

Integer square root in python

• answer by user448810 使用“牛顿法”。我不确定使用“整数除法”进行除法是否可以。我稍后会尝试这个作为更新。 [更新 (2015-01-11):这样做是可以的]

• answer by math 涉及使用第 3 方 Python 包 gmpy，这对我来说不是很有利，因为我主要对仅使用内置工具在 Erlang 中解决它感兴趣。

• answer by DSM 看起来很有趣。我不是很明白这是怎么回事，不过好像里面有“位魔法”，所以也不太适合我。

Infinite Recursion in Meta Integer Square Root

• 这道题是针对 C++ 的，AraK(提问者)的算法看起来与上面的Trial 2 出自相同的想法。

Answer 1

像下面这样的二进制搜索如何不需要浮点除法只需要整数乘法(比牛顿方法慢):-

low = 1;

/* More efficient bound

high = pow(10,log10(target)/2+1);

*/


high = target


while(low<high) {

 mid = (low+high)/2;
 currsq = mid*mid;

 if(currsq==target) {
    return(mid);
 }

 if(currsq<target) {

      if((mid+1)*(mid+1)>target) {
             return(mid);
      }    
      low =  mid+1;
  }

 else {

     high = mid-1;
 }

}

这适用于 O(logN) 迭代，因此即使对于非常大的数字也不应该永远运行

Log10(target) 计算(如果需要):-

acc = target

log10 = 0;

while(acc>0) {

  log10 = log10 + 1;
  acc = acc/10;
}

注意: acc/10是整数除法

编辑:-

有效边界:- sqrt(n) 的位数大约是 n 的一半，因此您可以通过 high = 10^(log10(N)/2+1) && low = 10^(log10(N)/2-1) 以获得更严格的约束，它应该提供 2 倍的加速。

评估界限:-

bound = 1;
acc = N;
count = 0;
while(acc>0) {

 acc = acc/10;

 if(count%2==0) {

    bound = bound*10;
 }

 count++;

}

high = bound*10;
low = bound/10;
isqrt(N,low,high);