Поиск пар чисел, которые удовлетворяют простоту и делимости условия


Я написал программу на Python 3 для проверки всех кортежей вида (s,t) где s и t положительные целые числа, удовлетворяющие трем условиям:

  1. Неравенство s <= t^2 и t <= s^2 удержание и значение s+t делит st(s+1)(t+1);
  2. Значение s+1 - простое;
  3. Неравенства как s * [math.ceil( math.ceil( (t^2)/(s+1) ) * ((s+1)/t) )] > s(s+t) и s > t держать.

Меня интересует как часто кортеж (s,t) удовлетворяет каждое из этих условий, в свою очередь; то есть, как часто кортеж удовлетворяет условию (1), против встречи (1) и (2), против встречи (1) - (3).

import math

def IsPrime(n):
# This returns True if a number is prime, and False if the number is composite.
    if n % 2 == 0 and n > 2: 
        return False
    return all(n % i for i in range(3, int(math.sqrt(n)) + 1, 2)) 

def IsGQ(s,t):
# This checks the divisibility condition.
    return (s*t*(s+1)*(t+1)) % (s+t) == 0 and s <= t**2 and t <= s**2

def IsNotTransitive(s,t):
    n = math.ceil((t**2)/(s+1))
    k = math.ceil(n*((s+1)/t))
    return (s*k) > t*(t+s) and s > t

rng = 1000 # The upper limit that `t` will iterate up to
quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) and (2) 
prime_quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) - (3) 
intransitive_quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) - (4)

# The next 5 lines place all prime numbers up to rng^2 into a list. 
# Ideally, this cuts down on how many times s+1 has to be checked to be a prime.
primes = [] 
for i in range(1, rng**2):
    if IsPrime(i):
        primes.append(i)

for t in range(4, rng + 1): # Due to project details, I don't care about t<4. 
    if t % 50 == 0: 
        print("We have reached t = " + str(t) + ".")
    for s in range(2, t**2 + 1): # To satisfy condition (1), I don't want s to iterate above t^2
        if IsGQ(s,t): # Condition (1) and (2)?
            quads += 1
            if s+1 in primes:            # Condition (3)?
                prime_quads += 1
                if IsNotTransitive(s,t): # Condition (4)?
                    intransitive_quads += 1

print(str(quads))
print(str(prime_quads))
print(str(intransitive_quads))  

В настоящее время, это работает для завершения в течение 10 минут rng = 1000. В идеале, работает это для rng = 10000 в разумные сроки-это моя нынешняя цель; вычислять эти цифры для значения rng = 100000 это мои самые оптимистические цели.

Есть ли способ, чтобы изменить порядок мое итерации или поменять способ проверки primeness, что может быть эффективнее? Кажется, что время для s+1 in primes для завершения растет квадратично с rng. Так, после определенного момента, s+1 будет меньше, чем все значения в primes это еще не был проверен, было бы лучше, чтобы написать свое собственное "в" функции вдоль линии

def CheckPrime(s, primes):
    for n in primes:
        if s+1 == n:
            return True
        elif s+1 < n:
            return False

или это все равно будет медленнее, чем текущий код?

Примечание: интерфейс isprime исходит от этой переполнение стека ответ.



248
8
задан 10 апреля 2018 в 07:04 Источник Поделиться
Комментарии
4 ответа

Вы правы в том, что наибольшее время тонуть в этой функции ваш IsPrime функция. Ваш CheckPrime функция в теории должна быть быстрее, однако, есть лучший способ.

Во-первых, относительно небольшая оптимизация: вместо использования IsPrime из связанных поэтому вопрос, использовать решето Эратосфена для создания списка из True / False значения в зависимости от индекса. Вот простая реализация:

def sieve(limit):
primes = [True] * limit
primes[0] = primes[1] = False

for i in range(limit):
if primes[i]:
for n in range(i ** 2, limit, i):
primes[n] = False

return primes

Если вы кэшируете результат этой функции, проверить, является ли число простым становится таким же простым, как cache[n], O(1) вместо грубо O(n).

Это изменение, безусловно, помогает, с rng = 1000на моем компьютере это заняло 4 минуты, 3 секунд для запуска. *

Лучшее изменение, которое я вижу-это уменьшение диапазона перебора s. Если условие 1 не проходит, нет необходимости даже прикасаться к этому значению по s. Имея это в виду, я изменил s in range(2, t**2 + 1) для s in range(int(t ** 0.5)) и проверил код снова. Хотя это, безусловно, поможет с больших диапазонов, это не помогло слишком много с rng = 1000в результате время (опять же) 4 минуты и 3 секунды. Я подозреваю, что это потому, что на верхней границы предела, мы только получаем, чтобы пропустить ~100 итераций. Это будет иметь значение, если вы запустите это rng = 100000.

Далее, я вспомнил, что читал, что стоимость вызова функции в Python является достаточно высокой [источник] - 150ns на вызов. Это может показаться не так много, но когда вы вызываете функцию около 300 сто миллионов раз (IsGQ функция, rng = 1000), оно может быть значительным. Имея это в виду, я встроен в IsGQ функции и бросил ряд проверок, так как они сейчас обрабатываются range функции и перенесли s > t проверить условие 3 (4?) в случае, если заявление. Это хоть немного помогало, меньше, чем я ожидал, в результате чего время выполнения 3 минуты и 54 секунды.

Там, наверное, больше оптимизации, что может быть сделано здесь, но я не эксперт в Python, и ничего не выпрыгивает на меня.

Вот код, который я использовал для окончательного запуска:

import math

rng = 1000

def sieve(limit):
primes = [True] * limit
primes[0] = primes[1] = False

for i in range(limit):
if primes[i]:
for n in range(i ** 2, limit, i):
primes[n] = False

return primes

print('Populating prime cache')
primeCache = sieve(rng ** 2 + 2)
print('Prime cache populated')

def IsNotTransitive(s,t):
n = math.ceil((t**2)/(s+1))
k = math.ceil(n*((s+1)/t))
return (s*k) > t*(t+s)

quads = 0
prime_quads = 0
intransitive_quads = 0

for t in range(4, rng + 1):
if t % 50 == 0:
print("We have reached t = " + str(t))
for s in range(int(t ** 0.5), t**2 + 1):
if s * t * (s + 1) * (t + 1) % (s + t) == 0:
quads += 1
if primeCache[s + 1]:
prime_quads += 1
if s > t and IsNotTransitive(s,t):
intransitive_quads += 1

print(str(quads))
print(str(prime_quads))
print(str(intransitive_quads))

* Все сроки указаны процессорного времени, не в реальном времени. В реальном времени было только несколько секунд выше. 4m7s для первого запуска.

2
ответ дан 10 апреля 2018 в 09:04 Источник Поделиться

Вы работаете на несколько проблем с производительностью.

Во-первых: если вы хотите создать список простых чисел, вы должны использовать сито, которое является гораздо более эффективным.

Во-вторых: вы проверяете s+1 in primesнеэффективных список. Вместо этого вы должны сделать primes будет набор, так что ваш поиск времени является постоянной.
В-третьих: вместо for s in range(2, t**2 + 1):вы должны использовать for s in range(floor(t**.5), t**2 + 1):, так как это мешает тратить время на проверку значений, где s^2<t.

Между этими двумя изменениями, которые вы должны увидеть довольно большой прирост скорости.

2
ответ дан 10 апреля 2018 в 09:04 Источник Поделиться

Если вы хотите улучшить производительность, вы должны знать, где вы проводите свое время. Вы можете использовать профайлер или грех в этом случае достаточно, если вы вставляете некоторые заявления, которые получают фактическое время, как timer.timer() разница функции возвращают значения между эта функция истекшее время в секундах.

percentages
t time %pr %lp cnt %quad quad
50 0.05 20 80 42864 4.1 1766
100 0.32 4.1 96 338239 1.7 5769
150 1.26 2.7 97 1136114 0.94 10675
200 2.84 2.4 98 2686489 0.63 16994
250 6.02 2.3 98 5239364 0.46 24122
300 10.43 2.1 98 9044739 0.36 32359

Так для t=300 прошедшее время был 10.43 секунды, 2.1% этого времени, где используется для вычисления простых чисел, 98% времени программа потратила в курсе. cntэто количество, как часто функция IsGQ называлась и quadsпоказывает, как часто IsGQ был правда и quad был увеличен. %quad показывает вам, сколько из этих cnt номера, где на самом деле проходят IsGQ тест.

в %пр/%числа ЛП показать вам, что оптимизация первой части вашего algorithn, были простыми числами обнаружил, ничего не меняет. Даже если создание primt занимает 0 секунд, ваш алгоритм не будет получать быстрее.

Я сравнил ваш главный алгоритм создания алгоритма просеивания, предложенные
@OscarSmith или @Gerrit0 и получили следующие результаты

prime creation
t t*t org imp
50 2500 0.01 0.00
100 10000 0.01 0.00
150 22500 0.03 0.01
200 40000 0.07 0.01
250 62500 0.14 0.02
300 90000 0.22 0.02

Так что ваш алгоритм достаточно плохо по сравнению с алгоритма решето, но тем не менее Ист влияние на время исполнения программы незначительна.

Я сравнивал свой алгоритм (орг) для усовершенствованного метода для проверки простых чисел, предложенный @Gerrit0 и способ, в котором числа хранятся как Python setи в цикл, в котором условие 2 и 3 не тестируются вообще, и, наконец, в пустой цикл, в котором условие проверяется на всех.

loop time 
t cnt org imp set nochk empty
50 42864 0.04 0.03 0.03 0.03 0.013
100 338239 0.31 0.22 0.21 0.20 0.028
150 1136114 1.23 0.83 0.72 0.71 0.098
200 2686489 2.77 1.62 1.63 1.66 0.22
250 5239364 5.88 3.50 3.57 3.41 0.44
300 9044739 10.21 6.27 6.08 6.10 0.75

Вы видите, что ваш способ найти в крайней мере свой яркий берет почти вдвое время метод @Gerrit0 по. Используя набор не улучшает скорость и нет другого метода, который может улучшить вашу производительность программы, так как даже если вы пропустите премьер-тестирования программы не будут работать быстрее, так как если вы используете @или мой метод Gerrit0 по. Я думаю, есть одна проблема с методом @Gerrit0, что он должен более существенную память, чем мой набор. Для больших t, что может сделать разницу. Если вы 'поставил' на 'пустой', вы потеряете большую часть своего времени в IsGQ которая большую часть времени возвращает Falseв %квад показывает.

На основании этих данных можно попробовать в модели время, необходимое для алгоритм и у меня

$$ \log (в\текст{время})=3,107 \журнала(Т) -6,691$$

"журнал" - это логарифм по основанию 10. Я использовал Excel для построения линейной регрессионной модели. Из этого я оценивал время выполнения цикла часть вашей программы, как

loop  
t est. time
50 0,04
100 0,33
150 1,18
200 2,87
250 5,75
300 10,13
400 24,76
500 49,52
1000 426,71
2000 3676,76
10000 546057,82
100000 698785591,19

Даже если вы можете половину времени с помощью лучшего способа, чтобы проверить, если число является ярким вам понадобится более 3 дней, чтобы запустить программу для T=1000 и около 11 лет, чтобы запустить его для Т=100000.

Время выполнения программы пропорционально T^3. Это похоже на правду, так как часто IsGQ называется? Это всего

$$\sum_{х=1}^ТС^2=\фрац{(2Т+1)(Т+1)т}{6}\приблизительно\фрац{Т^3}{3}$$

(см. A000330)

Так, кажется, вы не можете сделать свой алгоритм значительный быстрее с medhods предложил в постах @OscarSmith или @Gerrit0.

Время выполнения пустого цикла по сравнению со временем исполнения nochk взгляд показывает, что вы можете избежать IsGQ-тестирование чисел.


Может быть, следующая идея может улучшить производительность: в IsGQ вы проверьте, если S+T делит ст(з+1)(Т+1). Если вы разделите ст(з+1)(Т+1) С+Т остаток

$с $^2(х+1)(х-1)$$
или, из-за симметрии S и T

$$Т^2(Т+1)(Т-1)$$

Вместо прохода через все пары (Т,s) (есть о(Т^3)) Вы можете перебрать все возможные модули М=Т+С (есть только o(Т^2) и вычислить всех возможных с (не должно быть слишком много, это зависит от primefactorization м). Может быть, вы получите алгоритм с продолжительностью пропорционально T^2 вместо T^3. Но это достаточно сложная задача, и еще много где IsGQ результаты Falseпотому что это большие.


Наконец, осуществление IsNotTransitive странно. Почему вы используете расчет поплавка. Теперь вы должны рассмотреть, если ошибки плавающей точки оказывают влияние на вашу логику. ceil могут быть выполнены в целочисленной арифметике:

ceil(a/b) = -((-a)//b)

Это, кажется, даже немного быстрее.

2
ответ дан 16 апреля 2018 в 05:04 Источник Поделиться

Перво-наперво: есть из-за одной ошибки в коде, который заставил меня подозревать, что из-за одной ошибки в моей, пока я на него наткнулся:


for i in range(1, rng**2):
if IsPrime(i):
primes.append(i)

for t in range(4, rng + 1):
for s in range(2, t**2 + 1): # To satisfy condition (1), I don't want s to iterate above t^2
...
if s+1 in primes: # Condition (3)?


Когда \$ы = т^2\$, \$с + т = т(т+1)\$ Итак, первое условие всегда выполняется успешно. Но когда \$Т\ долл rng это означает, что мы тестируем rng**2+1 in primesи что может дать ложные негативы.



Есть ли способ, чтобы изменить порядок мое итерации или поменять способ проверки primeness, что может быть эффективнее?

Неправильный вопрос. Это не вопрос о повторном заказе итерации: дело в переборе что-то гораздо меньше.

Вы ищете \$з\$ такие, что$(с+т) \средний ст(з+1)(Т+1)\$. Но \$Х(Х+1) = (с+т)(ь+1-т) + (т-1)т\$, так что это эквивалентно ищу \$з\$ такие, что$(с+т) \средних (Т-1)Т^2(Т+1)\$. Это довольно легко-фактор \$(Т-1)Т^2(Т+1)\$, так как это уже частично учтены. Чтобы сделать его действительно эффективным, вы можете использовать сито (которое вы уже в любом случае понадобится Для теста делением).

Я адаптировал некоторые код, который я написал ранее, чтобы перечислить факторы дали премьер-факторизации, и был способен работать с rng=1000 за 1,4 секунды. Онлайн демо.

Узкое место (около 55% для rng=10000) был порожнее поколение, и это будет серьезной проблемой для памяти rng=100000. Поскольку мы используем его для факторизации до t+1и только для primality испытаний сверх того, я изменил сито, чтобы работать только до rng+2 и IsPrime в случае раскола. Это дало 70% ускорение (предположительно, избиение очевидный предел 55% за счет улучшения кэша), всего за 30 секунд rng=10000. Онлайн демо.

Тогда есть несколько небольших ускорений доступен нажатием серии тестов на s в генератор и не толкает в кучу значений s+t которые являются слишком большими. Я думаю, что было бы ускорение, заметив, что, когда вы инкремент t вы можете повторно использовать factorisations из t и t+1но что на самом деле не похоже на работу. Моя быстрая версия занимает 20 секунд, чтобы подойти к rng=10000:

import math
import heapq

rng = 10000 # The upper limit that `t` will iterate up to
quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) and (2)
prime_quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) - (3)
intransitive_quads = 0 # Counter for tuples (s,t) that satisfy conditions (1) - (4)

# Pre-calculate sieve
rr = rng + 2
max_prime = [0] * rr
for p in range(2, rr):
if max_prime[p] == 0:
max_prime[p] = p
for q in range(p * p, rr, p):
max_prime[q] = p

small_primes = [p for p in range(2, rr) if max_prime[p] == p]

def IsPrime(n):
if n < rng:
return max_prime[n] == n
for p in small_primes:
if n % p == 0:
return False
if p * p > n:
break
return True

def CandidateS(t):
# Factor (t-1)t^2(t+1)
primeFactors = dict()
for q in [t-1, t, t, t+1]:
r = q
while r > 1:
p = max_prime[r]
primeFactors[p] = primeFactors.get(p, 0) + 1
r //= p

primes = [(1, 1)] + sorted(primeFactors.items())
q = [(1, 0, 1)]
lower = t + int(math.ceil(math.sqrt(t)))
cutoff = t * (t + 1)
while len(q) > 0:
# d is the divisor
# i is the index of its largest "prime" in primes
# a is the exponent of that "prime"
(d, i, a) = heapq.heappop(q)
if d >= lower:
yield d - t
if a < primes[i][1]:
next = d * primes[i][0]
if next <= cutoff:
heapq.heappush(q, (next, i, a + 1))
if i + 1 < len(primes):
next = d * primes[i + 1][0]
if next <= cutoff:
heapq.heappush(q, (next, i + 1, 1))
# The condition i > 0 is to avoid duplicates arising because
# d == d // primes[0][0]
if i > 0 and a == 1:
next = d // primes[i][0] * primes[i + 1][0]
if next <= cutoff:
heapq.heappush(q, (next, i + 1, 1))

def IsNotTransitive(s,t):
n = math.ceil((t**2)/(s+1))
k = math.ceil(n*((s+1)/t))
return (s*k) > t*(t+s) and s > t

for t in range(4, rng + 1): # Due to project details, I don't care about t<4.
if t % 500 == 0:
print("We have reached t = " + str(t) + ".")
for s in CandidateS(t):
quads += 1
if IsPrime(s+1): # Condition (3)?
prime_quads += 1
if IsNotTransitive(s,t): # Condition (4)?
intransitive_quads += 1

print(str(quads))
print(str(prime_quads))
print(str(intransitive_quads))

Онлайн демо

1
ответ дан 17 апреля 2018 в 10:04 Источник Поделиться