Показать сообщение отдельно
Старый 15.11.2014, 18:02   #1014
Paul Kellerman
Gold Member
 
Регистрация: 25.06.2005
Адрес: F000:FFF0
Сообщений: 1,826
По умолчанию

Удалось написать и отладить весьма быстрый алгоритм извлечения целочисленного
корня 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

Последний раз редактировалось Paul Kellerman; 15.11.2014 в 20:29.
Paul Kellerman вне форума   Ответить с цитированием
Реклама