ИЗДАНО / ИСПРАВЛЕНО: Исправлены коды для прохождения 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://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 секунды, чтобы найти следующее регулярное число, большее 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