Найдите наименьшее регулярное число, которое не меньше N

Можно произвольно создать срез последовательности Хэмминга вокруг n-го члена во времени ~ n^(2/3) путем прямого перечисления троек (i,j,k) такое, что N = 2^i * 3^j * 5^k.

Алгоритм работает с log2(N) = i+j*log2(3)+k*log2(5); перечисляет все возможные k и для каждого, все возможные j, находит верхний i и, таким образом, тройку (k,j,i) и сохраняет его в «полосе», если внутри заданной «ширины» ниже заданного высокого логарифмического верхнего значения (когда width ‹1 может быть не более одного такого i), затем сортирует их по логарифму.

WP говорит, что n ~ (log N)^3, то есть время выполнения ~ (log N)^2. Здесь нас не заботит точное положение найденной тройки в последовательности, поэтому все вычисления подсчета из исходный код можно выбросить:

slice hi w = sortBy (compare `on` fst) b where       -- hi>log2(N) is a top value
  lb5=logBase 2 5 ; lb3=logBase 2 3                  -- w<1 (NB!) is log2(width)
  b  = concat                                        -- the slice
      [ [ (r,(i,j,k)) | frac < w ]                   -- store it, if inside width
        | k <- [ 0 .. floor ( hi   /lb5) ],  let p = fromIntegral k*lb5,
          j <- [ 0 .. floor ((hi-p)/lb3) ],  let q = fromIntegral j*lb3 + p,
          let (i,frac)=properFraction(hi-q) ;    r = hi - frac ]   -- r = i + q
                    -- properFraction 12.7 == (12, 0.7)

-- update: in pseudocode:
def slice(hi, w):
    lb5, lb3 = logBase(2, 5), logBase(2, 3)  -- logs base 2 of 5 and 3
    for k from 0 step 1 to floor(hi/lb5) inclusive:
        p = k*lb5
        for j from 0 step 1 to floor((hi-p)/lb3) inclusive:
           q = j*lb3 + p
           i = floor(hi-q)
           frac = hi-q-i                     -- frac < 1 , always
           r = hi - frac                     -- r == i + q
           if frac < w:
              place (r,(i,j,k)) into the output array
   sort the output array's entries by their "r" component
        in ascending order, and return thus sorted array

Перечислив тройки в срезе, можно легко выполнить сортировку и поиск, требуя практически O(1) времени (для произвольно тонкого среза), чтобы найти первую тройку выше N. На самом деле, для постоянной ширины (логарифмической) количество чисел в срезе (элементы «верхней корки» в (i,j,k)-пространстве ниже log(N) плоскости) снова m ~ n^2/3 ~ (log N)^2, а сортировка занимает m log m времени (так что поиск, даже линейный, тогда требуется ~ m времени выполнения). Но ширина может быть уменьшена для большего Ns, следуя некоторым эмпирическим наблюдениям; а постоянные коэффициенты для перебора троек в любом случае намного выше, чем для последующей сортировки.

Даже с постоянной шириной (логарифмической) он выполняется очень быстро, вычисляя 1 000 000-е значение в последовательности Хэмминга мгновенно и миллиардный за 0,05 с.

Первоначальная идея «верхней полосы троек» принадлежит Луи Клаудеру, как цитируется в моем опубликовать в блогах DDJ в 2008 году.

обновление: как указано GordonBGood в комментарии, нет необходимости во всей группе, а скорее одно или два значения выше и ниже целевого. С этой целью алгоритм легко корректируется. Входные данные также должны быть проверены на предмет того, что они являются числом Хэмминга перед продолжением алгоритма, чтобы избежать проблем с округлением с двойной точностью. Нет проблем с округлением при сравнении логарифмов чисел Хэмминга, которые, как было известно заранее, отличаются (хотя приближается к триллионной записи в последовательности используется около 14 значащих цифр в логарифмических значениях, оставляя только 1-2 цифры в запасе, поэтому на самом деле ситуация может быть сомнительной; но для 1-миллиардной нам нужно только 11 значащих цифр).

update2: оказывается, что двойная точность логарифмов ограничивает это числами от 20 000 до 40 000 десятичных цифр (т. е. от 10 до 100 триллионных чисел Хэмминга). Если в этом есть реальная необходимость для таких больших чисел, алгоритм можно переключить обратно на работу с самими целочисленными значениями, а не с их логарифмами, что будет медленнее.

Вот решение на Python, основанное на ответе Несса, но с некоторыми сокращениями и использованием чистой целочисленной математики, чтобы избежать попадания в пространство журнала ошибки числовой точности:

import math

def next_regular(target):
    """
    Find the next regular number greater than or equal to target.
    """
    # Check if it's already a power of 2 (or a non-integer)
    try:
        if not (target & (target-1)):
            return target
    except TypeError:
        # Convert floats/decimals for further processing
        target = int(math.ceil(target))

    if target <= 6:
        return target

    match = float('inf') # Anything found will be smaller
    p5 = 1
    while p5 < target:
        p35 = p5
        while p35 < target:
            # Ceiling integer division, avoiding conversion to float
            # (quotient = ceil(target / p35))
            # From https://stackoverflow.com/a/17511341/125507
            quotient = -(-target // p35)

            # Quickly find next power of 2 >= quotient
            # See https://stackoverflow.com/a/19164783/125507
            try:
                p2 = 2**((quotient - 1).bit_length())
            except AttributeError:
                # Fallback for Python <2.7
                p2 = 2**(len(bin(quotient - 1)) - 2)

            N = p2 * p35
            if N == target:
                return N
            elif N < match:
                match = N
            p35 *= 3
            if p35 == target:
                return p35
        if p35 < match:
            match = p35
        p5 *= 5
        if p5 == target:
            return p5
    if p5 < match:
        match = p5
    return match

На английском языке: перебирайте каждую комбинацию 5 и 3, быстро находя следующую степень 2> = цель для каждой пары и сохраняя наименьший результат. (Это пустая трата времени - перебирать все возможные кратные 2, если только одно из них может быть правильным). Он также возвращается раньше, если когда-либо обнаруживает, что цель уже является обычным числом, хотя в этом нет строгой необходимости.

Я протестировал его довольно тщательно, проверяя каждое целое число от 0 до 51200000 и сравнивая со списком в OEIS http://oeis.org/A051037, а также многие большие числа, которые ± 1 от обычных чисел, и т. д. Это теперь доступен в SciPy как fftpack.helper.next_fast_len, чтобы найти оптимальные размеры для БПФ (исходный код).

Я не уверен, что метод журнала работает быстрее, потому что мне не удалось заставить его работать достаточно надежно, чтобы проверить его. Думаю, у него такое же количество операций? Я не уверен, но это достаточно быстро. Требуется ‹3 секунды (или 0,7 секунды с gmpy), чтобы вычислить, что 2 ¹⁴² × 3 ⁸⁰ × 5 ⁴⁴⁴ - следующее обычное число выше 2 ² × 3 ⁴⁵⁴ × 5 ²⁴⁹ +1 (100000000-е регулярное число, состоящее из 392 цифр)

Вы хотите найти наименьшее число m, то есть m >= N и m = 2^i * 3^j * 5^k, где все i,j,k >= 0.

Логарифмируя уравнения можно переписать как:

 log m >= log N
 log m = i*log2 + j*log3 + k*log5

Вы можете вычислить log2, log3, log5 и logN с (достаточно высокой, в зависимости от размера N) точностью. Тогда эта проблема выглядит как проблема целочисленного линейного программирования, и вы можете попытаться решить ее с помощью одного известных алгоритмов для этой NP-трудной задачи.

ИСПРАВЛЕНО / ИСПРАВЛЕНО: исправлены коды для передачи scipy tests:

Вот ответ, основанный на ответе эндолита, но почти исключающий длинные целочисленные вычисления с разной точностью, используя представления логарифма float64 для создания базы сравнение, чтобы найти тройные значения, которые соответствуют критериям, прибегая к сравнениям с полной точностью только тогда, когда есть вероятность, что значение логарифма может быть недостаточно точным, что происходит только тогда, когда цель очень близка либо к предыдущему, либо к следующему регулярному числу:

import math

def next_regulary(target):
    """
    Find the next regular number greater than or equal to target.
    """
    if target < 2: return ( 0, 0, 0 )
    log2hi = 0
    mant = 0
    # Check if it's already a power of 2 (or a non-integer)
    try:
        mant = target & (target - 1)
        target = int(target) # take care of case where not int/float/decimal
    except TypeError:
        # Convert floats/decimals for further processing
        target = int(math.ceil(target))
        mant = target & (target - 1)

    # Quickly find next power of 2 >= target
    # See https://stackoverflow.com/a/19164783/125507
    try:
        log2hi = target.bit_length()
    except AttributeError:
        # Fallback for Python <2.7
        log2hi = len(bin(target)) - 2

    # exit if this is a power of two already...
    if not mant: return ( log2hi - 1, 0, 0 )

    # take care of trivial cases...
    if target < 9:
        if target < 4: return ( 0, 1, 0 )
        elif target < 6: return ( 0, 0, 1 )
        elif target < 7: return ( 1, 1, 0 )
        else: return ( 3, 0, 0 )

    # find log of target, which may exceed the float64 limit...
    if log2hi < 53: mant = target << (53 - log2hi)
    else: mant = target >> (log2hi - 53)
    log2target = log2hi + math.log2(float(mant) / (1 << 53))

    # log2 constants
    log2of2 = 1.0; log2of3 = math.log2(3); log2of5 = math.log2(5)

    # calculate range of log2 values close to target;
    # desired number has a logarithm of log2target <= x <= top...
    fctr = 6 * log2of3 * log2of5
    top = (log2target**3 + 2 * fctr)**(1/3) # for up to 2 numbers higher
    btm = 2 * log2target - top # or up to 2 numbers lower

    match = log2hi # Anything found will be smaller
    result = ( log2hi, 0, 0 ) # placeholder for eventual matches
    count = 0 # only used for debugging counting band
    fives = 0; fiveslmt = int(math.ceil(top / log2of5))
    while fives < fiveslmt:
        log2p = top - fives * log2of5
        threes = 0; threeslmt = int(math.ceil(log2p / log2of3))
        while threes < threeslmt:
            log2q = log2p - threes * log2of3
            twos = int(math.floor(log2q)); log2this = top - log2q + twos

            if log2this >= btm: count += 1 # only used for counting band
            if log2this >= btm and log2this < match:
                # logarithm precision may not be enough to differential between
                # the next lower regular number and the target, so do
                # a full resolution comparison to eliminate this case...
                if (2**twos * 3**threes * 5**fives) >= target:
                    match = log2this; result = ( twos, threes, fives )
            threes += 1
        fives += 1

    return result

print(next_regular(2**2 * 3**454 * 5**249 + 1)) # prints (142, 80, 444)

Поскольку большинство длинных вычислений с высокой точностью было исключено, gmpy не требуется, а в IDEOne приведенный выше код занимает 0,11 секунды вместо 0,48 секунды для решения endolith, чтобы найти следующее регулярное число больше 100-миллионного, как показано; требуется 0,49 секунды вместо 5,48 секунды, чтобы найти следующее регулярное число после миллиардного (следующее будет (761,572,489) после (1334,335,404) + 1), и разница будет становиться еще больше по мере увеличения диапазона, поскольку вычисления с высокой точностью становятся все длиннее для эндолита версия по сравнению с почти никакой здесь. Таким образом, эта версия может вычислить следующее регулярное число от триллионного числа в последовательности примерно за 50 секунд на IDEOne, тогда как с эндолитной версией это, вероятно, займет больше часа.

Описание алгоритма на английском языке почти такое же, как и для эндолитовой версии, но отличается следующим: 1) вычисляет логарифмическую оценку целевого значения аргумента (мы не можем использовать встроенную функцию log напрямую, так как диапазон может быть слишком велик для представления в виде 64-битного числа с плавающей запятой), 2) сравнивает значения представления журнала при определении квалифицирующих значений внутри оценочного диапазона выше и ниже целевого значения, состоящего только из двух или трех чисел (в зависимости от округления), 3 ) сравнивать значения с разной точностью только в том случае, если в пределах указанной выше узкой полосы; 4) выводит тройные индексы, а не полное длинное целое с разной точностью (будет около 840 десятичных цифр для единицы после миллиардной, в десять раз больше, чем для триллионной ), который при необходимости можно легко преобразовать в длинное значение с множественной точностью.

Этот алгоритм почти не использует память, кроме потенциально очень большого целого целого числа с множественной точностью, промежуточных сравнительных значений примерно того же размера и выходного расширения троек, если это необходимо. Этот алгоритм является усовершенствованием по сравнению с эндолитной версией, поскольку он успешно использует значения логарифма для большинства сравнений, несмотря на их недостаточную точность, и сужает диапазон сравниваемых чисел до нескольких.

Этот алгоритм будет работать для диапазонов аргументов несколько выше десяти триллионов (время вычисления несколько минут при скорости IDEOne), когда он больше не будет правильным из-за отсутствия точности в значениях представления журнала согласно обсуждению @ WillNess; чтобы исправить это, мы можем изменить представление журнала на представление логарифма «сверните свой собственный», состоящее из целого числа фиксированной длины (124 бита для примерно двойного диапазона экспоненты, хорошо для целей, содержащих более ста тысяч цифр, если один готов подождать); это будет немного медленнее из-за того, что небольшие целочисленные операции с множественной точностью будут медленнее, чем операции float64, но не намного медленнее, поскольку размер ограничен (может быть, в три раза или около того).

Теперь ни одна из этих реализаций Python (без использования C, Cython или PyPy или чего-то еще) не является особенно быстрой, поскольку они примерно в сто раз медленнее, чем реализованные на скомпилированном языке. Для справки, вот версия Haskell:

{-# OPTIONS_GHC -O3 #-}

import Data.Word
import Data.Bits

nextRegular :: Integer -> ( Word32, Word32, Word32 )
nextRegular target
  | target < 2                   = ( 0, 0, 0 )
  | target .&. (target - 1) == 0 = ( fromIntegral lg2hi - 1, 0, 0 )
  | target < 9                   = case target of
                                     3 -> ( 0, 1, 0 )
                                     5 -> ( 0, 0, 1 )
                                     6 -> ( 1, 1, 0 )
                                     _ -> ( 3, 0, 0 )
  | otherwise                    = match
 where
  lg3 = logBase 2 3 :: Double; lg5 = logBase 2 5 :: Double
  lg2hi = let cntplcs v cnt =
                let nv = v `shiftR` 31 in
                if nv <= 0 then
                  let cntbts x c =
                        if x <= 0 then c else
                        case c + 1 of
                          nc -> nc `seq` cntbts (x `shiftR` 1) nc in
                  cntbts (fromIntegral v :: Word32) cnt
                else case cnt + 31 of ncnt -> ncnt `seq` cntplcs nv ncnt
          in cntplcs target 0
  lg2tgt = let mant = if lg2hi <= 53 then target `shiftL` (53 - lg2hi)
                      else target `shiftR` (lg2hi - 53)
           in fromIntegral lg2hi +
                logBase 2 (fromIntegral mant / 2^53 :: Double)
  lg2top = (lg2tgt^3 + 2 * 6 * lg3 * lg5)**(1/3) -- for 2 numbers or so higher
  lg2btm = 2* lg2tgt - lg2top -- or two numbers or so lower
  match =
    let klmt = floor (lg2top / lg5)
        loopk k mtchlgk mtchtplk =
          if k > klmt then mtchtplk else
          let p = lg2top - fromIntegral k * lg5
              jlmt = fromIntegral $ floor (p / lg3)
              loopj j mtchlgj mtchtplj =
                if j > jlmt then loopk (k + 1) mtchlgj mtchtplj else
                let q = p - fromIntegral j * lg3
                    ( i, frac ) = properFraction q; r = lg2top - frac
                    ( nmtchlg, nmtchtpl ) =
                      if r < lg2btm || r >= mtchlgj then
                        ( mtchlgj, mtchtplj ) else
                      if 2^i * 3^j * 5^k >= target then
                        ( r, ( i, j, k ) ) else ( mtchlgj, mtchtplj )
                in nmtchlg `seq` nmtchtpl `seq` loopj (j + 1) nmtchlg nmtchtpl
          in loopj 0 mtchlgk mtchtplk
    in loopk 0 (fromIntegral lg2hi) ( fromIntegral lg2hi, 0, 0 )


trival :: ( Word32, Word32, Word32 ) -> Integer
trival (i,j,k) = 2^i * 3^j * 5^k

main = putStrLn $ show $ nextRegular $ (trival (1334,335,404)) + 1 -- (1126,16930,40)

Этот код вычисляет следующее регулярное число после миллиардного за слишком малое время, чтобы его можно было измерить, и следующее за триллионным за 0,69 секунды в IDEOne (и потенциально может работать даже быстрее, за исключением того, что IDEOne не поддерживает LLVM). Даже Джулия будет работать примерно на этой скорости Haskell после "разминки" перед JIT-компиляцией.

EDIT_ADD: код Джулии следующий:

function nextregular(target :: BigInt) :: Tuple{ UInt32, UInt32, UInt32 }
    # trivial case of first value or anything less...
    target < 2 && return ( 0, 0, 0 )

    # Check if it's already a power of 2 (or a non-integer)
    mant = target & (target - 1)

    # Quickly find next power of 2 >= target
    log2hi :: UInt32 = 0
    test = target
    while true
        next = test & 0x7FFFFFFF
        test >>>= 31; log2hi += 31
        test <= 0 && (log2hi -= leading_zeros(UInt32(next)) - 1; break)
    end

    # exit if this is a power of two already...
    mant == 0 && return ( log2hi - 1, 0, 0 )

    # take care of trivial cases...
    if target < 9
        target < 4 && return ( 0, 1, 0 )
        target < 6 && return ( 0, 0, 1 )
        target < 7 && return ( 1, 1, 0 )
        return ( 3, 0, 0 )
    end

    # find log of target, which may exceed the Float64 limit...
    if log2hi < 53 mant = target << (53 - log2hi)
    else mant = target >>> (log2hi - 53) end
    log2target = log2hi + log(2, Float64(mant) / (1 << 53))

    # log2 constants
    log2of2 = 1.0; log2of3 = log(2, 3); log2of5 = log(2, 5)

    # calculate range of log2 values close to target;
    # desired number has a logarithm of log2target <= x <= top...
    fctr = 6 * log2of3 * log2of5
    top = (log2target^3 + 2 * fctr)^(1/3) # for 2 numbers or so higher
    btm = 2 * log2target - top # or 2 numbers or so lower

    # scan for values in the given narrow range that satisfy the criteria...
    match = log2hi # Anything found will be smaller
    result :: Tuple{UInt32,UInt32,UInt32} = ( log2hi, 0, 0 ) # placeholder for eventual matches
    fives :: UInt32 = 0; fiveslmt = UInt32(ceil(top / log2of5))
    while fives < fiveslmt
        log2p = top - fives * log2of5
        threes :: UInt32 = 0; threeslmt = UInt32(ceil(log2p / log2of3))
        while threes < threeslmt
            log2q = log2p - threes * log2of3
            twos = UInt32(floor(log2q)); log2this = top - log2q + twos

            if log2this >= btm && log2this < match
                # logarithm precision may not be enough to differential between
                # the next lower regular number and the target, so do
                # a full resolution comparison to eliminate this case...
                if (big(2)^twos * big(3)^threes * big(5)^fives) >= target
                    match = log2this; result = ( twos, threes, fives )
                end
            end
            threes += 1
        end
        fives += 1
    end
    result
end

Вот еще одна возможность, о которой я только что подумал:

Если N имеет длину X бит, то наименьшее регулярное число R ≥ N будет в диапазоне
[2X-1, 2X]

например если N = 257 (двоичное 100000001), то мы знаем, что R равно 1xxxxxxxx, если R не равно в точности следующей степени 2 (512)

Чтобы сгенерировать все регулярные числа в этом диапазоне, мы можем сначала сгенерировать нечетные регулярные числа (т. Е. Кратные степени 3 и 5), затем взять каждое значение и умножить на 2 (путем битового сдвига) столько раз, сколько необходимо, чтобы получить это в этот диапазон.

В Python:

from itertools import ifilter, takewhile
from Queue import PriorityQueue

def nextPowerOf2(n):
    p = max(1, n)
    while p != (p & -p):
        p += p & -p
    return p

# Generate multiples of powers of 3, 5
def oddRegulars():
    q = PriorityQueue()
    q.put(1)
    prev = None
    while not q.empty():
        n = q.get()
        if n != prev:
            prev = n
            yield n
            if n % 3 == 0:
                q.put(n // 3 * 5)
            q.put(n * 3)

# Generate regular numbers with the same number of bits as n
def regularsCloseTo(n):
    p = nextPowerOf2(n)
    numBits = len(bin(n))
    for i in takewhile(lambda x: x <= p, oddRegulars()):
        yield i << max(0, numBits - len(bin(i)))

def nextRegular(n):
    bigEnough = ifilter(lambda x: x >= n, regularsCloseTo(n))
    return min(bigEnough)

Знаешь что? Я положу деньги на предположение, что на самом деле «тупой» алгоритм самый быстрый. Это основано на наблюдении, что следующее обычное число, как правило, не кажется намного большим, чем данный ввод. Так что просто начните считать, и после каждого приращения рефакторинг и посмотрите, нашли ли вы обычное число. Но создайте по одному потоку обработки для каждого доступного ядра, а для N ядер пусть каждый поток проверяет каждое N-е число. Когда каждый поток нашел число или пересек пороговое значение степени двойки, сравните результаты (сохраните лучшее число) и готово.

Я написал небольшую программу на C # для решения этой проблемы. Он не очень оптимизирован, но это только начало. Это решение довольно быстро работает с числами размером до 11 цифр.

private long GetRegularNumber(long n)
{
    long result = n - 1;
    long quotient = result;

    while (quotient > 1)
    {
        result++;
        quotient = result;

        quotient = RemoveFactor(quotient, 2);
        quotient = RemoveFactor(quotient, 3);
        quotient = RemoveFactor(quotient, 5);
    }

    return result;
}

private static long RemoveFactor(long dividend, long divisor)
{
    long remainder = 0;
    long quotient = dividend;
    while (remainder == 0)
    {
        dividend = quotient;
        quotient = Math.DivRem(dividend, divisor, out remainder);
    }
    return dividend;
}