24.11.2020 08:30
vdsina_m 7 1400 Источник


Pixelmusement создаёт видео об играх для MS-DOS и программном обеспечении. Каждое видео завершается коротким случайно выбранным списком тех, кто поддержал канал финансово. В видео ADG Filler #57 Крис рассказал, как происходит процесс выбора. Оказалось. что он абсолютно вписывается в основную тему канала: генерация производится программой на QBasic. В его программе используется встроенный генератор псевдослучайных чисел QBasic (pseudo random number generator, PRNG). Даже учитывая ограничения платформы, этот PRNG гораздо хуже, чем он мог бы быть. Давайте обсудим его слабые стороны и разберёмся, как сделать выбор более справедливым.

В программе Криса порождающее значение (seed) PRNG связано с системными часами (RANDOMIZE TIMER, идиомой QBasic). Программа заполняет массив спонсорами, представленными в виде целых чисел (индексов), а потом постоянно перемешивает список, пока пользователь не нажмёт клавишу, после чего, наконец, выводит случайную выборку из массива. Вот упрощённая версия программы (примечание: комментарии в QBasic начинаются с апострофа '):

CONST ntickets = 203  ' input parameter
CONST nresults = 12

RANDOMIZE TIMER

DIM tickets(0 TO ntickets - 1) AS LONG
FOR i = 0 TO ntickets - 1
    tickets(i) = i
NEXT

CLS
PRINT "Press any key to stop shuffling..."
DO
    i = INT(RND * ntickets)
    j = INT(RND * ntickets)
    SWAP tickets(i), tickets(j)
LOOP WHILE INKEY$ = ""

FOR i = 0 to nresults - 1
    PRINT tickets(i)
NEXT

Код вполне должен быть понятен, даже если вы не знаете QBasic. Примечание: в реальной программе более щедрые спонсоры получают больше «билетов» (tickets), что повышает их вес в результатах. Это учитывается в последнем цикле, чтобы никто не встречался в списке более одного раза. Всё это почти не относится к обсуждаемому в статье, поэтому я решил игнорировать подобные особенности.

В конечном итоге, готовый результат — это функция от всего трёх множеств входных данных:

  1. Системных часов (TIMER)
  2. Общего количества «билетов»
  3. Количества итераций цикла до нажатия клавиши

Второй пункт имеет удобное свойство — став спонсором, вы влияете на результат.

RND в QBasic


PRNG языка QBasic представляет собой такой 24-битный линейный конгруэнтный метод (Linear Congruential Generator, LCG):

uint32_t
rnd24(uint32_t *s)
{
    *s = (*s*0xfd43fd + 0xc39ec3) & 0xffffff;
    return *s;
}

Результатом является полное 24-битное состояние. RND делит его на 2^24 и возвращает результат как число с плавающей запятой одиночной точности, поэтому вызывающая программа получает значение в интервале от 0 до 1 (не включая 1).

Не стоит говорить, что это очень плохой PRNG. Константы LCG выбраны разумно, но странно, что было решено ограничить состояние 24 битами. Согласно 16-битному ассемблерному коду QBasic (примечание: перечисленные здесь константы LCG неверны), реализация является полным значением, кратным 32 битам и состоящим из 16-битных частей, и при сохранении результата она выделяет память и записывает все 32 бита. Как и можно ожидать для 8086, от использования только нижних 24 бит реализация ничего не выигрывает.

Чтобы показать, насколько он плох, покажу рандограмму этого PRNG, в которой заметна очевидная структура. (Это небольшой срез рандограммы размером 4096x4096, в которой каждый из 24-битных сэмплов 2^23 нанесён на график как две 12-битные координаты.)


Общепризнанно, что это сильно ограничивает PRNG. При 24-битном состоянии он хорош только для 4 096 (2^12) выходных значений, после чего он больше не следует парадоксу дней рождения: никакие выходные значения не повторяются, хотя мы должны увидеть повторения. Однако, как я покажу ниже, на самом деле это неважно.

Вместо того, чтобы отбрасывать верхние 8 бит — самые высококачественные выходные биты — проектировщики QBasic должны были отказаться от нижних 8 бит, превратив реализацию в урезанный 32-битный LCG:

uint32_t
rnd32(uint32_t *s)
{
    *s = *s*0xfd43fd + 0xc39ec3;
    return *s >> 8;
}

Этот LCG имел бы такую же производительность, но со значительно лучшим качеством. Вот рандограмма этого PRNG, и она тоже сильно ограничена (больше 65 536 (2^16) выходных значений).


Качественное улучшение, и получили мы его совершенно бесплатно!

RANDOMIZE в QBasic


Но это ещё не конец наших проблем. Конструкция RANDOMIZE получает seed с двойной точностью (т.е. 64-битной). Верхние 16 бит его двоичного представления по IEEE 754 подвергаются XOR с следующими верхними 16 битами. Верхним 16 битам состояния PRNG присваивается этот результат. Самые нижние 8 бит сохраняются.

Чтобы это стало понятнее, приведу реализацию на C, проверенную по QBasic 7.1:

uint32_t s;

void
randomize(double seed)
{
    uint64_t x;
    memcpy(&x ,&seed, 8);
    s = (x>>24 ^ x>>40) & 0xffff00 | (s & 0xff);
}

Другими словами, RANDOMIZE задаёт генератору PRNG только одно из 65 536 возможных состояний.

В качестве последнего примера: вот как реализован RND, тоже проверенный по QBasic 7.1:

float
rnd(float arg)
{
    if (arg < 0) {
        memcpy(&s, &arg, 4);
        s = (s & 0xffffff) + (s >> 24);
    }
    if (arg != 0.0f) {
        s = (s*0xfd43fd + 0xc39ec3) & 0xffffff;
    }
    return s / (float)0x1000000;
}

Seed по системным часам


Функция TIMER возвращает количество секунд с одиночной точностью, прошедших с полночи, и точность составляет ~55 мс (т.е. равна счётчику прерываний таймера 18,2 Гц). Это только время суток, в результат не входит текущая дата, как это происходит, например в unix-времени.

Это означает, что TIMER может возвращать всего 1 572 480 уникальных значений. Это количество мало даже без учёта того, что конструкцией RANDOMIZE они сопоставляются всего с 65 536 возможными seed. К счастью, все эти значения реализуемы с помощью TIMER.

Из всех трёх множеств входящих данных для случайного выбора очень плохо выглядит первое.

Итерации цикла


Идея Криса по постоянному перемешиванию массива до нажатия клавиши в большой степени компенсирует слабость PRNG языка QBasic. Он позволяет выполнять более 200 000 обменов массива, проходя более 2% периода PRNG, и сам массив используется как расширенное состояние PRNG, дополняя 24-битное состояние RND.

Поскольку итерации происходят так быстро, точное количество итераций становится ещё одним источником энтропии. Результаты при выполнении в течение 214 600 итераций будут сильно отличаться от результатов при 273 500 итераций.

Возможный вариант улучшения: выполнять выход из цикла только при нажатии определённой клавиши. Если нажимается любая другая клавиша, то эти входящие данные и TIMER примешиваются к состоянию PRNG. Долбёжка по клавиатуре во время выполнения цикла добавляет энтропии.

Заменяем PRNG


Поскольку встроенный PRNG настолько плох, мы можем улучшить ситуацию, реализовав новый в самом QBasic. Сложность в том, что в QBasic нет беззнаковых integer, нет даже операторов с беззнаковыми integer (например, >>> из Java и JavaScript), а переполнение значений со знаком является ошибкой времени выполнения. Мы даже не можем заново реализовать сам LCG языка QBasic без выполнения программного умножения long, потому что промежуточный результат переполняет его 32-битный LONG.

При таких ограничениях популярными вариантами выбора являются генератор Парка-Миллера (как мы видели это в Bash) или запаздывающий генератор Фибоначчи (он используется в Emacs, который долгое время был ограничен 29-битными integer).

Однако у меня есть идея получше: PRNG, основанный на RC4. Точнее, мой собственный проект называется Sponge4, функцией губки поверх RC4. Если вкратце: для примешивания дополнительных входящих данных достаточно снова выполнить развёртывание ключа. Для реализации этого PRNG требуются всего две простые операции: сложение по модулю 2^8 и обмен массивов. В QBasic есть конструкция SWAP, поэтому такая система очень ему подходит!

Sponge4 (RC4) обладает гораздо более качественными результатами, чем 24-битный LCG, и к нему можно примешать больше источников энтропии. При своём 1700-битном состоянии он без потерь может впитать довольно много энтропии.

Изучаем QBasic


До предыдущих выходных я не касался QBasic около 23 лет и мне, по сути, пришлось изучать его с нуля. Однако всего за пару часов я уже знал его лучше, чем в прошлом. В основном так получилось, потому что у меня теперь гораздо больше опыта, но ещё и потому, что практически все туториалы по QBasic ужасны. Не удивительно, что они написаны для начинающих, но такое ощущение, что и писали их начинающие. Вскоре у меня сложилось впечатление, что сообщество QBasic в целом было ещё одним примером ситуации «слепой ведёт слепого».

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

  • Основные типы: INTEGER (int16), LONG (int32), SINGLE (float32), DOUBLE (float64) и два вида STRING — с постоянной и переменной длиной. Последние версии также имели неполную поддержку 64-битного типа CURRENCY 10,000x с фиксированной запятой.
  • Переменные по умолчанию имеют тип SINGLE и их не нужно объявлять заранее. Массивы по умолчанию содержат 11 элементов.
  • Переменные, константы и функции могут иметь суффикс, если их тип не равен SINGLE: INTEGER %, LONG &, SINGLE !, DOUBLE #, STRING $ и CURRENCY @. Для функций это тип возвращаемого значения.
  • Каждый тип переменных имеет собственное пространство имён, т.е. i% отличается от i&. Массивы также имеют собственное пространство имён, т.е. i% отличается от i%(0) и от i&(0).
  • Переменные можно объявлять явным образом при помощи DIM. При объявлении переменной с помощью DIM можно не использовать суффикс. Также такое объявление исключает это имя из других пространств имён, т.е. при объявлении DIM i AS LONG делает невозможным использование i% в той же области видимости. Хотя массивы и скаляры всё равно могут иметь такое же имя даже при объявлениях через DIM.
  • Численные операции со смешанными типами косвенно меняют тип, как в C.
  • Функции и подпроцедуры имеют единое, общее пространство имён, вне зависимости от суффикса функции. Из-за этого в местах вызова функций суффикс (обычно) можно пропускать. Особыми в этом случае являются встроенные фукнции.
  • Несмотря на первое впечатление, QBasic имеет статическую типизацию.
  • По умолчанию передача выполняется по ссылке. Для передачи по значению следует использовать BYVAL.
  • В объявлениях массивов параметр — это не размер, а наибольший индекс. Есть поддержка многомерных массивов. Массивы не обязаны индексироваться с нуля (например, есть (x TO y)), хотя это используется по умолчанию.
  • Строки являются не массивами, а самостоятельными типами со специальными конструкциями и функциями.
  • Областями видимости являются модуль, подпроцедура и функция. «Глобальные» переменные нужно объявлять с помощью SHARED.
  • Пользователи могут задавать собственные структуры с помощью TYPE. Функции не могут возвращать заданные пользователями типы и вместо этого используют передачу по ссылке.
  • Грубое подобие динамического распределения поддерживается с помощью REDIM, изменяющего размер массивов $DYNAMIC во время выполнения. ERASE освобождает выделенную память.

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

Реализуем Sponge4


Как и в RC4, мне нужен 256-элементный массив байтов и два 1-байтных индекса, i и j. Sponge4 также хранит третий 1-байтный счётчик k для подсчёта входящих данных.

TYPE sponge4
    i AS INTEGER
    j AS INTEGER
    k AS INTEGER
    s(0 TO 255) AS INTEGER
END TYPE

QBasic не имеет типа byte. Обычно в таком случае подходит 256-байтная строка фиксированной длины, но поскольку здесь они не являются массивами, строки несовместимы со SWAP и не индексируются эффективно. Поэтому я решил смириться с неиспользуемым пространством и использовать для всего 16-битные integer.

Для этой структуры есть четыре «метода». Три из них являются подпроцедурами, потому что не возвращают значение, а модифицируют губку. Последний, squeeze, возвращает следующий байт как INTEGER (%).

DECLARE SUB init (r AS sponge4)
DECLARE SUB absorb (r AS sponge4, b AS INTEGER)
DECLARE SUB absorbstop (r AS sponge4)
DECLARE FUNCTION squeeze% (r AS sponge4)

Инициализация выполняется в соответствии с RC4:

SUB init (r AS sponge4)
    r.i = 0
    r.j = 0
    r.k = 0
    FOR i% = 0 TO 255
        r.s(i%) = i%
    NEXT
END SUB

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

SUB absorb (r AS sponge4, b AS INTEGER)
    r.j = (r.j + r.s(r.i) + b) MOD 256
    SWAP r.s(r.i), r.s(r.j)
    r.i = (r.i + 1) MOD 256
    r.k = (r.k + 1) MOD 256
END SUB

SUB absorbstop (r AS sponge4)
    r.j = (r.j + 1) MOD 256
END SUB

Для сжатия байта может сначала потребоваться примешивание состояния, после чего производится обычное выполнение генератора RC4.

FUNCTION squeeze% (r AS sponge4)
    IF r.k > 0 THEN
        absorbstop r
        DO WHILE r.k > 0
            absorb r, r.k
        LOOP
    END IF

    r.j = (r.j + r.i) MOD 256
    r.i = (r.i + 1) MOD 256
    SWAP r.s(r.i), r.s(r.j)
    squeeze% = r.s((r.s(r.i) + r.s(r.j)) MOD 256)
END FUNCTION

Вот и весь генератор на QBasic! Однако нам будет полезна пара вспомогательных функций. Одна поглощает целые строки, а вторая возвращает 24-битные результаты.

SUB absorbstr (r AS sponge4, s AS STRING)
    FOR i% = 1 TO LEN(s)
        absorb r, ASC(MID$(s, i%))
    NEXT
END SUB

FUNCTION squeeze24& (r AS sponge4)
    b0& = squeeze%(r)
    b1& = squeeze%(r)
    b2& = squeeze%(r)
    squeeze24& = b2& * &H10000 + b1& * &H100 + b0&
END FUNCTION

QBasic не имеет операций битового сдвига, поэтому приходится пользоваться умножением. &H — это шестнадцатеричная запись.

Используем губку в деле


Одна из проблем исходной программы заключается в том, что в seed используется только время суток. Даже если смешивать его получше, то при запуске ровно в один и тот же момент в два разных дня мы получим одинаковый seed. Функция DATE$ возвращает текущую дату, поэтому мы можем поглотить её в губку, сделав частью входящих данных дату целиком.

DIM sponge AS sponge4
init sponge
absorbstr sponge, DATE$
absorbstr sponge, MKS$(TIMER)
absorbstr sponge, MKI$(ntickets)

За этим следует таймер. Он преобразуется в строку при помощи MKS$, возвращающей двоичное представление с одиночной точностью в little-endian в виде 4-байтной строки. MKI$ делает то же самое для INTEGER в виде 2-байтной строки.

Одной из проблем исходной программы было смещение: умножение RND на константу с последующим усечением результата до integer в большинстве случаев неоднородно. Некоторые числа выбираются чуть чаще, чем остальные, потому что входящие данные 2^24 не могут равномерно быть сопоставлены, допустим, с 10 выходящими данными. Учитывая, что в исходной программе используется перемешивание массива, это, вероятно, не имеет на практике разницы, но я бы хотел этого избежать.

В своей программе я учитываю это, генерируя ещё одно число, если оно попадает в дополнительный «хвост» распределения входящих данных (это очень маловероятно для малого значения ntickets). Функция squeezen выполняет равномерную генерацию числа от 0 до N (не включая N).

FUNCTION squeezen% (r AS sponge4, n AS INTEGER)
    DO
       x& = squeeze24&(r) - &H1000000 MOD n
    LOOP WHILE x& < 0
    squeezen% = x& MOD n
END FUNCTION

Затем выполняем тасование Фишера-Йетса и вывод первых N элементов:

FOR i% = ntickets - 1 TO 1 STEP -1
    j% = squeezen%(sponge, i% + 1)
    SWAP tickets(i%), tickets(j%)
NEXT

FOR i% = 1 TO nresults
    PRINT tickets(i%)
NEXT

Но если вам очень понравилась идея Криса, то можно сделать так:

PRINT "Press Esc to finish, any other key for entropy..."
DO
    c& = c& + 1
    LOCATE 2, 1
    PRINT "cycles ="; c&; "; keys ="; k%

    FOR i% = ntickets - 1 TO 1 STEP -1
        j% = squeezen%(sponge, i% + 1)
        SWAP tickets(i%), tickets(j%)
    NEXT

    k$ = INKEY$
    IF k$ = CHR$(27) THEN
        EXIT DO
    ELSEIF k$ <> "" THEN
        k% = k% + 1
        absorbstr sponge, k$
    END IF
    absorbstr sponge, MKS$(TIMER)
LOOP

Если вы хотите проверить программу самостоятельно, допустим, в DOSBox, вот её полные исходники: sponge4.bas.


На правах рекламы


Эпичные серверы — это виртуальные серверы для любых задач. Вы можете создать собственный тарифный план, максимальная конфигурация — 128 ядер CPU, 512 ГБ RAM, 4000 ГБ NVMe!