Sexy primes, «медленный питон» или как я бился о стену непонимания

Многие разработчики, особенно принимающие активное участие в проектировании системы, наверняка сталкивались с подобной ситуацией: приходит коллега (разраб, проектлид или продажник не суть важно) с очередной идеей-фикс: давай перепишем все на java, scala и т.д. (любимое подставить).

Вот и мне в очередной раз «спустили» такую идею в немаленьком-таком legacy проекте. Не совсем переписать, не совсем все (ну в перспективе). В общем перейти с питона (а у нас там еще и тикль модульно присутствует) на scala. Речь пока шла о разработке новых модулей и сервисов, т.е. начинать с наименее привязанных к middle-level и front-nearby API's. Как я понял в перспективе возможно совсем.

Человек — не разработчик, типа нач-проекта и немного продажник (для конкретного клиента) в одном лице.

Я не то, чтобы против. И скалу уважаю и по-своему люблю. Обычно я вообще открыт ко всему новому. Так, например, местами кроме тикля и питона у нас появляются сервисы или модули на других языках. Так, например, мы переехали с cvs на svn, а затем на git (а раньше, давно-давно, вообще MS-VSS был). Примеров на самом деле масса, объединяет их один момент — так решили или как минимум одобрили сами разработчики (коллективно ли, или была группа инициаторов — не суть важно). Да и дело в общем в аргументах за и против.

Проблема в том, что иногда для аргументированной дискуссии «Developer vs. Anybody-Else» у последнего не дотягивает уровень знаний «материи» или просто невероятно сложно донести мысль — т.е. как-бы разговор на разных языках. И хорошо если это кто-нибудь типа software architect. Хуже, если имеем «беседу» например с чистым «продажником», огласившим например внезапные «требования» заказчика.

Ну почему никто не предписывает, например, плиточнику — каким шпателем ему работать (типа с зубцами 10мм клея же больше уйдет, давайте может все же 5мм. А то что там полы-стены кривущие никого не волнует). И шуруп теоретически тоже можно «закручивать» молотком, но для этого же есть отвертка, а позже был придуман шуруповёрт. Утрирую конечно, но иногда действительно напоминает такой вот абсурд.

Это я к тому, что инструмент в идеале должен выбирать сам разработчик или иметь в этом как минимум последнее слово — ему с этим инструментом работать или мучиться.

Но что-то я отвлекся. В моей конкретной истории аргументов — за scala, у человека как всегда почти никаких.
Хотя я мог бы долго говорить про вещи, типа наличие разрабов, готовые наработки, отточенную и отлаженную систему и т.д. и т.п. Но зацепился за его «Питон очень медленный». В качестве примера он в меня кинул ссылкой на Interpreting a benchmark in C, Clojure, Python, Ruby, Scala and others — Stack Overflow, которую он даже до конца не прочитал (ибо там почти прямым текстом есть — не так плохо все с питоном).

Имелось ввиду именно вот это (время указано в секундах):

  Sexy primes up to:        10k      20k      30k      100k
  ---------------------------------------------------------
  Python2.7                1.49     5.20    11.00       119     
  Scala2.9.2               0.93     1.41     2.73     20.84
  Scala2.9.2 (optimized)   0.32     0.79     1.46     12.01

Разговаривать с ним на уровне, про чисто техническую сторону, нет ни малейшего желания. Т.е. про выделение памяти/пул объектов, GC (у нас не чистый CPython, скорее похож на RPython или pypy с собственным Jit, MemMgr, GC), про всякий сахар, с которым человек, писавший бенчь, немного переборщил и т.д.

Мой ответ был «любят разрабатывать на питоне не за это» и «на нем так не пишут, по крайней мере критичный к скорости код». Я немного слукавил и естественно понимаю, что этот пример для benchmark искусственный, ну т.е. значит чисто померить скорость. Но и болячки, которые при таком тесте повыскакивают в CPython, мне тоже известны.
Поэтому все же постарался на этом конкретном примере показать почему этот тест не целесообразен, ну или почему не совсем объективный что ли.

Начал с того, что показал ему этот тест и результаты исполнения (помечены звездой) в PyMNg (что есть наша сборка), pypy2, pypy3 и python3 (который то же был по тем же понятным причинам медленный). Железо, конечно, другое, но порядок, думаю, примерно одинаков.

  Sexy primes up to:        10k      20k      30k      100k
  ---------------------------------------------------------
  PyMNg                *   0.10        -        -      2.46
  pypy2                *   0.13        -        -      5.44
  pypy3                *   0.12        -        -      5.69
  ---------------------------------------------------------
  scala 2.10.4 (warmed up) *           -        -      6.57
  scala 2.10.4 (cold)      *           -        -      7.11
  ---------------------------------------------------------
  Python3.4            *   1.41        -        -    117.69
  ---------------------------------------------------------
  Python2.7                1.49     5.20    11.00    119.00
  Scala2.9.2               0.93     1.41     2.73     20.84
  Scala2.9.2 (optimized)   0.32     0.79     1.46     12.01

Дальше можно было в принципе не продолжать, просто хотелось сделать попытку объяснить человеку, в чем он не прав, что выбор языка не оценивается по бенчмаркам типа «Hello, world» и т.д.

Итак, задача — ищем «сексуальные» пары простых чисел на питоне. Много и быстро.

Для тех кому разбор полетов не интересен, ссылка на скрипт на github, если есть желание поиграться.

Результаты исполнения каждого варианта под спойлером соответственно.
100K для всех вариантов.


10M начиная с варианта max (остальные просто медленные на таком массиве).


100M начиная с варианта siev1 (по той же причине).


Кстати, такой разброс между CPython и PyPy и является часто одной из причин, почему люди переходят на PyPy, пишут собственные алокаторы, менеджеры памяти, GC, stackless и иже с ними, используют сторонние модули (например NumPy) и делают собственные сборки. Например, когда важна скорость исполнения и как здесь, имеем «агрессивное» использование пула объектов / множественные вызовы и т.д. Так же было когда-то давно и в нашем случае, когда мы переехали с чистого питона. Там еще было много чего, и тормозящий multithreding, и refcounting и т.д. Но сам переезд был обдуманным решением всей команды, а не «спущенной» сверху причудой. Если есть интерес и найду время, можно будет попробовать тиснуть про это статью.

Для этой же конкретной «задачи» можно было бы написать собственный C-binding, использовать модули типа numpy и т.д. Я же попробовал убедить коллегу, что оно на коленке за незначительное время решается практически «алгоритмически», если знаешь, как питон тикает внутри.

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

Оригинальный скрипт, немного измененный мной для читабельности и поддержки третьего питона, этот вариант у меня в скрипте-примере называется org. _{(Только, плз, не надо здесь про «xrange vs range» — я прекрасно представляю, в чем различие, и здесь конкретно оно не суть важно, тем более в 3-м питоне, кроме того, и итерация как-бы «завершенная»).}

def is_prime(n):
  return all((n % i > 0) for i in range(2, n))
# def sexy_primes_below(n):
#   return [[j-6, j] for j in range(9, n+1) if is_prime(j) and is_prime(j-6)]
def sexy_primes_below(n):
  l = []
  for j in range(9, n+1):
    if is_prime(j-6) and is_prime(j):
      l.append([j-6, j])
  return l

Т.к. даже на 10M имеем всего 100K sexy пар, изменение оригинальной primes_below на мой вариант с прямым циклом не сильно влияет на время исполнения, просто оно наглядней для изменений в последующих вариантах (например сито). Весь цимес в реализации is_prime, во всяком случае пока.

1. Во-первых, использование такого «сахара» как в оригинале (тем более в «сборках» типа PyPy, ну или нашего PyMNg) не поощряется, ну или как минимум, как и в этом случае, больно бьет отдачей в виде снижения скорости. Перепишем это как вариант mod1

def is_prime(n):
  i = 2
  while True:
    if not n % i:
       return False
    i += 1
    if i >= n:
       return True
  return True

Уже быстрее, как минимум в PyPy. Но дело не в этом…

2. Код стал сразу наглядней и видно, что его можно переписать как mod2 в половину быстрее, если не проверять четные номера (которые, кроме двойки, изначально не prime).

def is_prime(n):
  if not n % 2:
    return False
  i = 3
  while True:
    if n % i == 0:
       return False
    i += 2
    if i >= n:
       return True
  return True

Поправим это в оригинале — org2 это то же самое что и mod2, но в одну строчку используя «сахар».

def is_prime(n):
  return n % 2 and all(n % i for i in range(3, n, 2))

3. Если проверять делители до значения квадратного корня (правильно было бы до округленного, но мы сделаем проще — это же просто тест), то все можно сделать еще намного быстрее, получим вариант max:

def is_prime(n):
  if not n & 1:
    return 0
  i = 3
  while 1:
    if not n % i:
       return 0
    if i * i > n:
       return 1
    i += 2
  return 1

Намного быстрее, правда.

Опять правим это в оригинале — orgm.

def is_prime(n):
  #return ((n & 1) and all(n % i for i in range(3, int(math.sqrt(n))+1, 2)))
  return ((n & 1) and all(n % i for i in range(3, int(n**0.5)+1, 2)))

И видим, что как минимум в PyPy оно опять выполняется медленнее (хотя частично, возможно, и из-за прямого подсчета «корня», который в range).

4. Тут у коллеги загораются глаза, и он как в том мультфильме (по-моему, «Жадный богач») про скорняка и 7 шапок выдает: «А можешь еще быстрее?». Т.к. по памяти ограничения нет (не emdedded и т.д.) решил ему быстро переделать, используя «половинчатое» сито — half sieve, что есть подготовленный массив флажков по смещению для нечетных чисел, т.е. разделенных на 2. Тем более, что на питоне организовать такое сито на раз-два.
Ну и сразу видоизменяем sexy_primes_below, вызвав в нем генерацию сита ну и чтобы не править is_prime и не вызывать его в цикле, спрашиваем сразу sieve.
Получаем вариант siev1.

def primes_sieve(n):
  """ temporary "half" mask sieve for primes < n (using bool) """
  sieve = [True] * (n//2)
  for i in range(3, int(n**0.5)+1, 2):
    if sieve[i//2]:
      sieve[i*i//2::i] = [False] * ((n-i*i-1)//(2*i)+1)
  return sieve
def sexy_primes_below(n):
  l = []
  sieve = primes_sieve(n+1)
  #is_prime = lambda j: (j & 1) and sieve[int(j/2)]
  for j in range(9, n+1):
    i = j-6
    #if (i & 1) and is_prime(i) and is_prime(j):
    if (i & 1) and sieve[int(i/2)] and sieve[int(j/2)]:
      l.append([i, j])
  return l

Этот вариант настолько быстрый, что в PyPy, например, он на 100M выдает практически то же время, что оригинал на 100K. Т.е. проверяя в 2000 раз больше чисел и генерируя на несколько порядков больший список сексуально-простых пар.

Сразу переписал в вариант siev2, потому что вспомнил о несколько «туповатой» обработке bool в PyPy. Т.е. заменив булевы флажки на 0/1. Этот пример отрабатывает на 100M уже вдвое-трое быстрее оригинала на 100K!

5. Варианты osiev1 и osiev2 написал, чтобы в дальнейшем можно было заменить сито для всех чисел, на множество более коротких, т.е. чтобы иметь возможность осуществлять поиск пар инкрементально или блочно.

Для этого изменил сито смещений на прямое сито хранящее не флаги, а уже сами значения:

sieve = [1, 1, 1, 0, 1, 1 ...]; # для 3, 5, 7, 9, 11, 13 ...
osieve = [3, 5, 7, 0, 11, 13, ...]; # для 3, 5, 7, 9, 11, 13 ...

Вариант osiev1.

def primes_sieve(n):
  """ temporary odd direct sieve for primes < n """
  sieve = list(range(3, n, 2))
  l = len(sieve)
  for i in sieve:
    if i:
      f = (i*i-3) // 2
      if f >= l:
        break
      sieve[f::i] = [0] * -((f-l) // i)
  return sieve
def sexy_primes_below(n):
  l = []
  sieve = primes_sieve(n+1)
  #is_prime = lambda j: (j & 1) and sieve[int((j-2)/2)]
  for j in range(9, n+1):
    i = j-6
    #if (i & 1) and is_prime(i) and is_prime(j):
    if (i & 1) and sieve[int((i-2)/2)] and sieve[int((j-2)/2)]:
      l.append([i, j])
  return l

Ну и второй вариант osiev2 просто «чуть» быстрее, т.к. инициирует сито гораздо оптимальнее.

def primes_sieve(n):
  """ temporary odd direct sieve for primes < n """
  #sieve = list(range(3, n, 2))
  l = ((n-3)//2)
  sieve = [-1] * l
  for k in range(3, n, 2):
    o = (k-3)//2
    i = sieve[o]
    if i == -1:
      i = sieve[(k-3)//2] = k
    if i:
      f = (i*i-3) // 2
      if f >= l:
        break
      sieve[f::i] = [0] * -((f-l) // i)
  return sieve

Эти два метода можно было переделать на итеративное сито (например, искать пары блочно по 10K или 100K). Это позволило бы значительно сэкономить память при подсчете. Например, если сейчас попробовать оба osieve варианта с 1G или 10G, мы практически гарантировано сразу получим MemoryError исключение (ну или вы богатый буратино — и у вас очень много памяти и 64-х битный питон.

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

А на исходных 100K время исполнения последнего уже невозможно было подсчитать — 0.00 mils (я подсчитал его, только увеличив количество итераций исполнения times до 10).

В чем я практически уверен — это то, что увеличить скорость еще на один порядок вряд ли удастся не только на scala, но и на чистом C. Если только снова алгоритмически…

Сам скрипт, буде кто собирается поэкспериментировать, можно спросить помощь, например, так:

sexy-primes-test.py --help

Если что, про простые числа довольно неплохо и очень подробно написано в wikihow.

По просьбам трудящихся добавил опрос…