| Paul Kellerman |
15.11.2014 18:02 |
Удалось написать и отладить весьма быстрый алгоритм извлечения целочисленного
корня N-й степени из целого числа произвольной разрядности. В основе старый доб-
рый метод Ньютона, однако его сходимость напрямую зависит от начального прибли-
жения, и при неудачном выборе - миллиарды итераций, и крайне медленно сходится.
Выкладываю дзен-решение, формирующее начальное приближение, за счет которого
метод Ньютона сходится всего за 2-15 итераций, и, например, корень 50000 степени
из 1000000-разрядного целого числа, вычисляется за 6 итераций, по времени 0,1 сек.
P.S. FirstNonZeroBit - номер старшего ненулевого бита числа в двоичном представлении,
по сути - очень быстрый способ вычисления целочисленного логарифма по основанию 2.
Код:
restart:
with(Bits):
IntRoot:= proc(x,n) local a,b,q,u,v,f,it:
if (n = 1) then
it:= 0:
return(x,it):
end if:
if ((x = 0) and (n > 1)) then
it:= 0:
return(0,it):
end if:
if ((x >= 1) and (n > 1)) then
q:= FirstNonzeroBit(x):
u:= trunc(q/n):
v:= frac(q/n):
f:= ceil(n*(2^v)):
a:= iquo((2^u)*f,n):
b:= iquo(x,a^(n-1)):
if (a < b) then
u:= trunc((q+1)/n):
v:= frac((q+1)/n):
f:= ceil(n*(2^v)):
a:= iquo((2^u)*f,n):
b:= iquo(x,a^(n-1)):
end if:
it:= 1:
while (a > b) do
a:= iquo((n-1)*a+b,n):
b:= iquo(x,a^(n-1)):
it:= it + 1:
end do:
return(a,it):
end if:
end proc:
x:= 7*10^1000000:
n:= 50000:
IntRoot(x,n);
a = 100003891896030419242, it = 6
a, b, x: MPInteger
n, q, u, f: Int64
v: Extended |
