Оптимизация кодового словаря (полиномиальное деление)

Ответ 1

Ниже на 3 раза быстрее, чем pypy на моей машине (0.04s против 0.15s). Использование Cython:

ctypedef unsigned char uint8_t # does not work with Microsoft C Compiler: from libc.stdint cimport uint8_t
cimport cpython.array as array

cdef uint8_t[::1] gf_exp = bytearray([1, 3, 5, 15, 17, 51, 85, 255, 26, 46, 114, 150, 161, 248, 19,
   lots of numbers omitted for space reasons
   ...])

cdef uint8_t[::1] gf_log = bytearray([0, 0, 25, 1, 50, 2, 26, 198, 75, 199, 27, 104, 
    more numbers omitted for space reasons
    ...])

import cython

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.initializedcheck(False)
def rsenc(msg_in_r, nsym, gen_t):
    '''Reed-Solomon encoding using polynomial division, better explained at http://research.swtch.com/field'''

    cdef uint8_t[::1] msg_in = bytearray(msg_in_r) # have to copy, unfortunately - can't make a memory view from a read only object
    cdef int[::1] gen = array.array('i',gen_t) # convert list to array

    cdef uint8_t[::1] msg_out = bytearray(msg_in) + bytearray(len(gen)-1)
    cdef int j
    cdef uint8_t[::1] lgen = bytearray(gen.shape[0])
    for j in xrange(gen.shape[0]):
        lgen[j] = gf_log[gen[j]]

    cdef uint8_t coef,lcoef

    cdef int i
    for i in xrange(msg_in.shape[0]):
        coef = msg_out[i]
        if coef != 0: # coef 0 is normally undefined so we manage it manually here (and it also serves as an optimization btw)
            lcoef = gf_log[coef] # precaching

            for j in xrange(1, gen.shape[0]): # optimization: can skip g0 because the first coefficient of the generator is always 1! (that why we start at position 1)
                msg_out[i + j] ^= gf_exp[lcoef + lgen[j]] # equivalent (in Galois Field 2^8) to msg_out[i+j] -= msg_out[i] * gen[j]

    # Recopy the original message bytes
    msg_out[:msg_in.shape[0]] = msg_in
    return msg_out

Это ваша самая быстрая версия со статическими типами (и проверка html от cython -a до тех пор, пока петли не будут выделены желтым цветом).

Несколько кратких заметок:

• Cython предпочитает x.shape[0] - len(shape)

• Определяя представления памяти как [::1] promises, они непрерывны в памяти, что помогает

• initializedcheck(False) необходимо избегать множества проверок существования на глобально определенных gf_exp и gf_log. (Возможно, вы можете ускорить свой базовый код Python/PyPy, создав для них локальную ссылку на переменную и используя этот istead)

• Мне пришлось скопировать пару входных аргументов. Cython не может создать представление памяти из объекта readonly (в данном случае msg_in). Я мог бы, возможно, просто сделать это char *). Также gen (список) должен быть в чем-то с быстрым доступом к элементу.

Кроме того, все это довольно прямолинейно. (Я не пробовал никаких вариантов его получения быстрее). Я действительно впечатлен тем, насколько хорошо работает PyPy.

Ответ 2

Основываясь на ответе DavidW, вот реализация, которую я сейчас использую, что примерно на 20% быстрее, используя nogil и параллельное вычисление:

from cython.parallel import parallel, prange

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.initializedcheck(False)
cdef rsenc_cython(msg_in_r, nsym, gen_t) :
    '''Reed-Solomon encoding using polynomial division, better explained at http://research.swtch.com/field'''

    cdef uint8_t[::1] msg_in = bytearray(msg_in_r) # have to copy, unfortunately - can't make a memory view from a read only object
    #cdef int[::1] gen = array.array('i',gen_t) # convert list to array
    cdef uint8_t[::1] gen = gen_t

    cdef uint8_t[::1] msg_out = bytearray(msg_in) + bytearray(len(gen)-1)
    cdef int i, j
    cdef uint8_t[::1] lgen = bytearray(gen.shape[0])
    for j in xrange(gen.shape[0]):
        lgen[j] = gf_log_c[gen[j]]

    cdef uint8_t coef,lcoef
    with nogil:
        for i in xrange(msg_in.shape[0]):
            coef = msg_out[i]
            if coef != 0: # coef 0 is normally undefined so we manage it manually here (and it also serves as an optimization btw)
                lcoef = gf_log_c[coef] # precaching

                for j in prange(1, gen.shape[0]): # optimization: can skip g0 because the first coefficient of the generator is always 1! (that why we start at position 1)
                    msg_out[i + j] ^= gf_exp_c[lcoef + lgen[j]] # equivalent (in Galois Field 2^8) to msg_out[i+j] -= msg_out[i] * gen[j]

    # Recopy the original message bytes
    msg_out[:msg_in.shape[0]] = msg_in
    return msg_out

Мне бы хотелось, чтобы он был быстрее (при реальной реализации данные кодируются со скоростью около 6.4 МБ/с с n = 255, n - размером сообщения + кодового слова).

Основной причиной более быстрой реализации, которую я нашел, является использование подхода LUT (LookUp Table) путем предварительного вычисления массивов умножения и сложения. Однако в моих реализациях Python и Cython подход LUT медленнее, чем вычисление операций XOR и добавления.

Существуют и другие подходы к внедрению более быстрого кодировщика RS, но у меня нет возможностей и времени, чтобы опробовать их. Я оставлю их в качестве ссылок для других заинтересованных читателей:

• "Быстрая реализация программного обеспечения конечных полевых операций", Чэн Хуан и Лихао Сюй, Вашингтонский университет в Сент-Луисе, Tech. Rep (2003). ссылка и правильная реализация кода здесь.
• Ло, Цзяньцян и др. "Эффективные программные реализации больших конечных полей GF (2 n) для приложений безопасного хранения". ACM Transactions on Storage (TOS) 8.1 (2012): 2.
• "Оценка эффективности и изучение библиотек кодирования Erasure для хранения с открытым исходным кодом", Plank, J. S. and Luo, J. and Schuman, C.D. и Xu, L. и Wilcox-O'Hearn, Z, FAST. Том 9. 2009. ссылкаИли также не расширенную версию: "Сравнение производительности библиотек кодирования стирания с открытым исходным кодом для приложений хранения", Plank и Schuman.
• Исходный код библиотеки ZFEC с оптимизацией оптимизации LUT ссылка.
• "Оптимизированная арифметика для кодов Рида-Соломона", Кристоф Паар (1997, июнь). В Международном симпозиуме IEEE по теории информации (стр. 250-250). ИНСТИТУТ ЭЛЕКТРОТЕХНИЧЕСКИХ ИНЖЕНЕРОВ INC (IEEE). ссылка
• "Быстрый алгоритм кодирования (255,233) Рид-Соломонова кода над GF (2 ^ 8)", Р.Л. Миллер и Т.К. Чыонг, И.С. Рид. ссылка
• "Оптимизация арифметики поля Галуа для различных процессорных архитектур и приложений", Гринана, Кевина и М., Этана и Л. Миллера и Томаса Дж. Э. Шварца, Моделирование, анализ и моделирование компьютеров и телекоммуникационных систем, 2008. MASCOTS 2008. IEEE Международный симпозиум. IEEE, 2008. ссылка
• Anvin, H. Peter. "Математика RAID-6". (2007). ссылка и ссылка
• Wirehair library, одна из немногих реализаций Коши Рида-Соломона, которая называется очень быстрой.
• "Логарифмический логический алгоритм времени для параллельного полиномиального деления", Бини, Д. и Пан, В.Ю. (1987), Буквы обработки информации, 24 (4), 233-237. См. Также Бини, Д. и В. Пан. "Быстрые параллельные алгоритмы для полиномиального деления над произвольным полем констант". Компьютеры и математика с приложениями 12.11 (1986): 1105-1118. link
• Kung, H.T. "Быстрая оценка и интерполяция". (1973). ссылка
• Цао, Чжэнцзюнь и Хэньюэ Цао. "Обратите внимание на алгоритм быстрого деления для многочленов, использующих итерацию Ньютона". arXiv preprint arXiv: 1112.4014 (2011). ссылка
• "Введение в поля Галуа и кодирование Рида-Соломона", Джеймс Уэстолл и Джеймс Мартин, 2010. ссылка
• Мамиди, Суман и др. "Расширения набора инструкций для кодирования и декодирования Рида-Соломона". Application-Specific Systems, Архитектурные процессоры, 2005. ASAP 2005. 16-я Международная конференция IEEE. IEEE, 2005. ссылка
• Дюма, Жан-Гийом, Лоран Фюсс и Бруно Сальви. "Одновременная модулярная редукция и замена Кронекера для малых конечных полей". Journal of Symbolic Computation 46.7 (2011): 823-840.
• Грина, Кевин М., Этан Л. Миллер и Томас Шварц. Анализ и построение месторождений galois для повышения надежности хранения. Том 9. Технический отчет UCSC-SSRC-07, 2007. ссылка

Однако, я считаю, что лучше всего использовать эффективную многочленную модульную редукцию вместо полиномиального деления:

• "Модульное сокращение в GF (2 n) без предварительной вычислительной фазы". Kneževic, M., et al. Арифметика конечных полей. Springer Berlin Heidelberg, 2008. 77-87.
• "О вычислении полиномиальной модулярной редукции". Ву, Хуапенг. Технический отчет, Univ. Ватерлоо, Центр прикладных криптографических исследований, 2000 г.
• "Быстрая реализация программного обеспечения для арифметических операций в GF (2n)". De Win, E., Bosselaers, A., Vandenberghe, S., De Gersem, P. и Vandewalle, J. (1996, январь). В Advances in Cryptology-Asiacrypt'96 (стр. 65-76). Спрингер Берлин Гейдельберг. ссылка
• Barnett reduction

/EDIT: на самом деле кажется, что "при вычислении полиномиальной модульной редукции" просто использует тот же подход, что и я, с вариантами rsenc_alt1() и rsenc_alt2() (основная идея состоит в том, что мы предкоммутируем пары коэффициентов, которые мы будем необходимо и сократить их все сразу), и, к несчастью, это не быстрее (это на самом деле медленнее, потому что прекомпция не может быть выполнена один раз для всех, поскольку она зависит от ввода сообщения).

/EDIT: я нашел библиотеку с действительно интересными оптимизациями, много, которые даже не найдены ни в каких научных статьях (которые автор заявил, что он прочитал кстати), и которая, вероятно, является самой быстрой программной реализацией Рида-Соломона: проект wirehair и кодек Коши-Рида-Соломона под названием longhair с аналогичными трюками оптимизации.

/FINAL EDIT: кажется, что самая быстрая реализация доступна на основе этой статьи:

Планк, Джеймс С., Кевин М. Гринан и Этан Л. Миллер. "Кричащий быстрая арифметика поля Галуа с использованием инструкций Intel SIMD." FAST. 2013. ссылка

Реализация в чистом Go доступна здесь и является автором Klaus Post. Это самая быстрая реализация, о которой я когда-либо читал, как в одном потоке, так и в параллельном (он поддерживает оба). Он требует более 1 ГБ/с в одном потоке и более 4 ГБ/с с 8 потоками. Однако он опирается на оптимизированные инструкции SIMD и различные низкоуровневые оптимизации для матричных операций (потому что здесь кодек RS ориентирован матрицей вместо полиномиального подхода, который у меня есть в моем вопросе).

Итак, если вы заинтересованный читатель и хотите найти самый быстрый доступный кодек Рида-Соломона, то тот.

Ответ 3

В качестве альтернативы, если вы знаете C, я бы рекомендовал переписать эту функцию Python в plain C и вызвать ее (скажем, с CFFI). По крайней мере, вы знаете, что вы достигаете максимальной производительности во внутренних циклах ваших функций, не зная о том, что вы используете PyPy или Cython.

Смотрите: http://cffi.readthedocs.org/en/latest/overview.html#performance