Портал аспирантов - Вычисление целочисленного логарифма

Пусть заданы произвольные целые числа A и B, причем оба могут быть
сверхбольшими (миллионы или миллиарды битов), и попытки свести их
к известным 32/64-битным целым типам, или 80-битным вещественным,
поддерживаемых аппаратно, априори отпадают. Необходимо вычислить
целочисленный логарифм N = ilog[B](A), так что B^N <= A < B^(N+1).

Особо отмечу, что A, B принадлежат бесконечному множеству целых
чисел, и поэтому целочисленный логарифм не имеет никакого отноше-
ния к дискретному логарифму конечных групп, колец или полей целых
чисел, где все операции выполняются по модулю некоторого целого p.

В математических пакетах есть специальная функция ilog[b](a), однако,
быстро она работает только для оснований B = 2 и 10. При других осно-
ваниях (даже относительно небольших, например, 3 или 11, очень долго
считает, если A - достаточное большое целое число, 1024000-битное).
Но может кто-то знает более быстрые и элегантные алгоритмы ilog[b](a).

Я долго ломал голову и навскидку набросал следующий алгоритм, кото-
рый быстро считает при любых основаниях B, даже сверхбольших целых.

Код:

FastIntLog:= proc(a,b) local n, u, v, w, c:

  u:= ilog[2](a):

  if (b = 2) then

    n:= u:

    return(n):

  end if:

  v:= ilog[2](b) + 1:

  n:= trunc(u/v): 

  c:= b^n: 

  while (a > c) do

    w:= trunc((u-ilog[2](c))/v):

    if (w > 0) then

      n:= n + w:

    else 

      n:= n + 1: 

    end if:

    c:= b^n:

  end do:

  if (a < c) then 

    n:= n - 1: 

  end if:

  return(n):

end proc:



a:= 2^1024000:

b:= 19^5000:



FastIntLog(a,b);



48

ilog[2] - это по сути номер старшего бита числа в двоичном представлении.
Вычисляется достаточно быстро при поддержке команд сканирования битов.