Алгоритм делителей

Почему бы не найти простую факторизацию каждого числа, как предлагает @SamVarshavchik в комментариях? Использование классического ситового подхода определенно будет длиться вечно для таких больших чисел, с которыми имеет дело OP (т.е. до 2⁶⁴), однако эти числа находятся в пределах диапазона алгоритмов, таких как Алгоритм Полларда Ро. От ответа на вопрос Какой самый быстрый алгоритм факторизации, мы можем видеть, что алгоритм Rho Полларда начинает разрушаться вокруг 2⁷⁰ диапазон, поэтому мы должны быть в порядке. Есть несколько опубликованных C++ реализации легко доступны, однако я предлагаю преобразовать Python код ниже, чтобы C++, Это более быстрый вариант исходного алгоритма, опубликованного Ричардом Брентом в 1980 году (код найден Вот)

def brent(N):
if N%2==0:
return 2
y,c,m = random.randint(1, N-1),random.randint(1, N-1),random.randint(1, N-1)
g,r,q = 1,1,1
while g==1:
x = y
for i in range(r):
y = ((y*y)%N+c)%N
k = 0
while (k<r and g==1):
ys = y
for i in range(min(m,r-k)):
y = ((y*y)%N+c)%N
q = q*(abs(x-y))%N
g = gcd(q,N)
k = k + m
r = r*2
if g==N:
while True:
ys = ((ys*ys)%N+c)%N
g = gcd(abs(x-ys),N)
if g>1:
break

return g

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

Чтобы этот алгоритм работал, вам понадобится функция для быстрого тестирования простоты (я предлагаю Тест первичности Миллера-Рабина). Эта функция будет использоваться на переднем крае в качестве фильтра, чтобы не тратить время на функцию Полларда Ро. Я вручную написал все вышеперечисленные алгоритмы на нескольких языках и могу подтвердить, что это очень простая задача.

Чтобы показать вам, что это приведет к желаемым результатам, я придумал сценарий наихудшего случая ниже и предоставил R код (преобразование в C++ должен запускаться мгновенно), что хорошо работает в течение 15 секунд. factorize функция от gmp пакет и реализует стандартный алгоритм Rho Полларда. Когда я говорю «наихудший случай», я имею в виду, что все 1000 чисел в вашем списке являются полупростыми (это занимает больше всего времени для множителя алгоритма Rho Полларда) очень близко к 2⁶⁴ предел. Если данные числа были просто 1000 случайных чисел меньше, чем 2⁶⁴, приведенный ниже алгоритм будет работать менее чем за пару секунд, так как вероятность того, что число имеет небольшой простой множитель, чрезвычайно велика (именно здесь проявляется истинная сила алгоритма Rho Полларда).

KindergartenGraduation <- function(passTime) {
BegTime <- Sys.time()
set.seed(11)
Samp_NoRep <- nextprime(sample(2^32, 2000, replace = FALSE))  ## generate 2000 random primes

## pick 2000 primes from the sample above with replacement
## so as to guarantee a non-trivial answer
Samp_WiRep <- sample(Samp_NoRep, 2000, replace = TRUE)

## Below, we multiply each odd element with its neighbor
## to create large semi-primes
## all odd indices           ## and the even ones
TestNums <- mul.bigz(Samp_WiRep[seq.int(1,2000,2)], Samp_WiRep[seq.int(2,2000,2)])

## finding the prime factorization should be the slowest step
BegFacTime <- Sys.time()
MyPriFacs <- do.call(c, lapply(TestNums, factorize))
EndFacTime <- Sys.time()
TotFacTime <- EndFacTime - BegFacTime
myFacTime <- strsplit(format(TotFacTime), split = " ")[[1]]
timeFacNum <- myFacTime[1]
timeFacUnit <- myFacTime[2]
print(paste(c("Factoring list of large integers took ",timeFacNum," ",timeFacUnit), collapse = ""))

## sorting in order to apply run length encoding
MyPriFacs <- MyPriFacs[order(asNumeric(MyPriFacs))]

MyRle <- function (x1) {
n1 <- length(x1)
y1 <- x1[-1L] != x1[-n1]
i <- c(which(y1), n1)
list(lengths = diff(c(0L, i)), values = x1[i], uni = sum(y1)+1L)
}

BreakDown <- MyRle(MyPriFacs)

## Find maximum power
MyMax <- max(BreakDown$lengths)

## Determine which prime(s) occur(s) the maximum number of times
MaxRep <- which(BreakDown$lengths == MyMax)

## If we find multiple, we combine them to obtain the single divisor with the maximum power
MyAns <- list(primes = BreakDown$values[MaxRep], divisor = prod.bigz(BreakDown$values[MaxRep]), largest_power = MyMax)
EndTime <- Sys.time()
TotTime <- EndTime - BegTime
myTime <- strsplit(format(TotTime), split = " ")[[1]]
timeNum <- myTime[1]
timeUnit <- myTime[2]
print(paste(c("You found the most frequent divisor in ",timeNum," ",timeUnit), collapse = ""))

if (as.numeric(timeNum) < passTime) {
print("Yay!!! I graduated Kindergarten on Earth 2.0 and I'm now ready to build my first quantum computer!!!")
} else {
print("Oh no!!! I failed Kindergarten and have to recite K&R 10^100 more times!!!")
}

return(MyAns)
}

Называя это мы имеем:

KindergartenGraduation(15)
[1] "Factoring list of large integers took 9.417561 secs"[1] "You found the most frequent divisor in 9.492585 secs"[1] "Yay!!! I graduated Kindergarten on Earth 2.0 and I'm now ready to build my first quantum computer!!!"$primes
Big Integer ('bigz') object of length 4:
[1] 243258997  335592419  2615771369 3202455203

$divisor
Big Integer ('bigz') :
[1] 683854798468133589409768665696700901

$largest_power
[1] 5

Как видите, простые числа: 243258997, 335592419, 2615771369, 3202455203 встречается чаще всего в 5 раз за штуку, таким образом, произведение этих чисел (т.е. 683854798468133589409768665696700901) — это делитель с максимальной мощностью для произведения всех 1000 чисел в данном списке.