Найдите наименьшее регулярное число, которое не меньше N
Обычные числа - это числа, которые равномерно делят степени 60. Например, 602 = 3600 = 48 × 75, поэтому и 48, и 75 являются делителями степени 60. Таким образом, они также являются обычными числами.
Это продолжение округления до следующей степени двух.
У меня есть целочисленное значение N, которое может содержать большие простые факторы, и я хочу округлить его до числа, состоящего только из небольших простых факторов (2, 3 и 5)
Примеры:
f(18) == 18 == 21 * 32f(19) == 20 == 22 * 51f(257) == 270 == 21 * 33 * 51
Какой эффективный способ найти наименьшее число, удовлетворяющее этому требованию?
Используемые значения могут быть большими, поэтому я бы хотел избежать перечисления всех обычных чисел, начиная с 1, или сохранения массива всех возможных значений.
8 ответов
Хорошо, надеюсь, третий раз здесь очарование. Рекурсивный алгоритм ветвления для начального ввода p, где N - это число, которое "строится" в каждом потоке. NB 3a-c здесь запускаются в виде отдельных потоков или иным образом (квази) асинхронно.
Рассчитайте следующую наибольшую степень 2 после p, назовите это R. N = p.
N > R? Выйти из этой темы. Состоит ли р только из небольших простых факторов? Вы сделали В противном случае перейдите к шагу 3.
После любого из 3a-c перейдите к шагу 4.
a) Округлите p до ближайшего кратного 2. Это число может быть выражено как m * 2.
б) Округлите p до ближайшего кратного 3. Это число может быть выражено как m * 3.
c) Округлите p до ближайшего числа, кратного 5. Это число может быть выражено как m * 5.Переходите к шагу 2 с p = m.
Я пропустил бухгалтерию, чтобы отслеживать N, но это довольно просто, я так понимаю.
Изменить: Забыл 6, спасибо ypercube.
Изменить 2: Если бы это до 30, (5, 6, 10, 15, 30) понял, что это было ненужным, убрал это.
Редактировать 3: (Последнее, что я обещаю!) Добавлена проверка степени 30, которая помогает предотвратить использование этого алгоритма всей вашей оперативной памятью.
Редактировать 4: Изменено power-of-30 на power-2, по наблюдениям finnw.
Можно создать сколь угодно тонкий фрагмент последовательности Хемминга вокруг 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 время выполнения тогда). Но ширина может быть уменьшена для большего N s, после некоторых эмпирических наблюдений; и постоянные факторы для перечисления троек намного выше, чем для последующей сортировки в любом случае.
Даже с постоянной шириной (логарифмической) он работает очень быстро, мгновенно вычисляя 1000 000-е значение в последовательности Хемминга и миллиардное за 0,05 с.
Первоначальная идея "верхней группы троек" принадлежит Луи Клаудеру, о чем говорилось в моем посте на дискуссии в блогах DDJ в 2008 году.
обновление: как отметил GordonBGood в комментариях, нет необходимости для всей полосы, а только для одного или двух значений выше и ниже цели. Алгоритм легко изменить с этой целью. Перед тем, как приступить к алгоритму, входные данные также должны быть проверены на наличие числа Хэмминга, чтобы избежать проблем округления с двойной точностью. Нет проблем округления при сравнении логарифмов чисел Хэмминга, о которых заранее известно, что они различны (хотя при увеличении до триллионной записи в последовательности используются около 14 значащих цифр в значениях логарифма, оставляя только 1-2 цифры, чтобы сэкономить, поэтому на самом деле ситуация там может оказаться ненадежной, но для одной миллиардной суммы нам нужно только 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://stackru.com/a/17511341/125507
quotient = -(-target // p35)
# Quickly find next power of 2 >= quotient
# See https://stackru.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 142 × 380 × 5444 - это следующее регулярное число выше 22 × 3454 × 5249+1 (100 000 000-ое регулярное число, которое имеет 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:
Вот ответ, основанный на ответе эндолита, но почти исключающий длинные целочисленные вычисления с множественной точностью, используя логарифмические представления 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://stackru.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 секунды, чтобы найти следующее регулярное число, большее 100-миллионного, как показано; для поиска следующего обычного числа, превышающего миллиардную, требуется 0,49 секунды вместо 5,48 секунды. (761,572,489) прошлое (1334,335,404) + 1), и разница будет увеличиваться по мере увеличения диапазона, так как вычисления с высокой точностью становятся все длиннее для версии с эндолитом по сравнению с почти отсутствующими здесь. Таким образом, эта версия может вычислить следующее регулярное число из триллионной в последовательности примерно за 50 секунд в IDEOne, где, вероятно, потребуется более часа с версией эндолита.
Английское описание алгоритма почти такое же, как и для версии эндолита, отличающееся следующим:
1) вычисляет логарифмическую оценку целевого значения аргумента (мы не можем использовать встроенный log работать напрямую, так как диапазон может быть слишком большим для представления в виде 64-разрядного числа с плавающей запятой),
2) сравнивает значения представления журнала в определении квалификационных значений в пределах оценочного диапазона выше и ниже целевого значения только примерно на два или три числа (в зависимости от при округлении),
3) сравнивать значения с множественной точностью, только если в пределах определенной выше узкой полосы,
4) выводить тройные индексы, а не полное длинное целое число с множественной точностью (было бы около 840 десятичных цифр для одного после миллиардной, в десять раз больше, чем на триллион), который может быть легко преобразован в длинное значение мультиточности при необходимости.
Этот алгоритм почти не использует память, кроме как для потенциально очень большого целочисленного целевого значения с множественной точностью, промежуточных значений сравнения оценки примерно одинакового размера и выходного расширения троек, если это необходимо. Этот алгоритм является улучшением по сравнению с версией эндолита в том, что он успешно использует значения логарифма для большинства сравнений, несмотря на их недостаточную точность, и что он сужает диапазон сравниваемых чисел до нескольких.
Этот алгоритм будет работать для диапазонов аргументов, несколько превышающих десять триллионов (время вычисления в несколько минут при скоростях IDEOne), когда он больше не будет корректным из-за отсутствия точности в значениях представления журнала в соответствии с обсуждением @WillNess; чтобы исправить это, мы можем изменить логарифмическое представление на логарифмическое представление roll-your-own, состоящее из целого числа фиксированной длины (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). Даже Джулия будет бегать с такой скоростью на Хаскеле после "разминки" для компиляции 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)
Я написал небольшую программу на 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;
}
Знаешь что? Я положу деньги на предположение, что на самом деле "тупой" алгоритм самый быстрый. Это основано на наблюдении, что следующее регулярное число, в общем, кажется, не намного больше, чем данный ввод. Так что просто начните считать, и после каждого приращения рефакторируйте и посмотрите, нашли ли вы обычный номер. Но создайте один поток обработки для каждого имеющегося у вас ядра, и для N ядер каждый поток проверяет каждое N-е число. Когда каждый поток найдет число или превысит пороговое значение степени 2, сравните результаты (оставьте наилучшее число) и получите результат.