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



通知

有关ErlangC / C++中的解决方案,请转到下面的试用版4


维基百科文章

整数平方根

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

计算平方根的方法

  • 可以在这里找到一种实现"比特魔术">的算法

[试用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 * 1111111111111111111111111111111111111111时,它似乎在我(非常快)的机器上永远运行。


[试用3:仅使用Bigint+1-1div 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(参考Vikram Bhat的答案,修改了上述代码,使边界改为Mid+1或Mid-1)

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位数字(即1111111111111111111111111111111111111111 * 1111111111111111111111111111111111111111),结果立即显示出来。然而,在我的机器上求解一百万(1000000)个61位数字(即从10000000000000000000000000000000000000000000000000000000000001000000000000000000000000000000000000000000000000000001000000)需要60秒(+-2秒)。我想做得更快。


[试验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秒)解决了一百万(1000000)个61位数字,所以这比试用版3更快。它能更快吗?


关于我的机器

处理器:3.4 GHz Intel Core i7

内存:32 GB 1600 MHz DDR3

操作系统:Mac OS X版本10.9.1


相关问题

python 中的整数平方根

  • 用户448810的答案使用"牛顿法">。我不确定用"整数除法"做除法是否可以。我稍后会尝试此更新。【更新(2015-01-11):这样做是可以的】

  • 数学的答案是使用第三方Python包gmpy,这对我来说不是很好,因为我主要感兴趣的是在Erlang中只使用内置的工具来解决它。

  • DSM的答案似乎很有趣。我真的不明白发生了什么,但似乎有"比特魔法">参与其中,所以它也不太适合我。

元整数平方根中的无限递归

  • 这个问题是针对C++的,AraK(提问者)的算法看起来与上面的试验2的想法相同

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

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(目标)计算(如果需要):-

acc = target
log10 = 0;
while(acc>0) {
log10 = log10 + 1;
acc = acc/10;
}

注意:acc/10是整数除法

编辑:-

有效界:-sqrt(n)的位数大约是n的一半,因此您可以通过high = 10^(log10(N)/2+1)&amp;CCD_ 19得到更严格的约束,并且它应该提供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);

相关内容

  • 没有找到相关文章

最新更新