О Random и RandomRange

[Сквозняк]

В начале работы программы функция Random работает хорошо, как на полигоне. Но когда запущено много потоков, значительно больше чем ядер, и они что-то колбасят и гоняют между собой данные, то Random работает звиздец как плохо. Наверно это связано с особенностями ОС и железа - данные "энтропии" используемые этой функцией поступают в неё в упорядоченном виде. В результате, наблюдал, что выхлоп Random(1000); почти никогда не был меньше 140 и тяготел к верхнему пределу. Но если запустить

Код: Выделить всё

for A:=1 to 10000 do write(random(1000),' ');
writeln;


то выхлоп становится "мягким и пушистым", практически идеальным. Печать в консоль стабилизирует работу функции Random, а если выхлоп не печатать, то он упорядычивается и можно ждать до "морковкиного заговения" когда выпадет результат в мелком диапазоне. Для себя решил, что в многопоточных приложениях, много использующих random(1000); и random(100); , надо запасаться в массивы данными рандома сразу после старта программы, а потом их оттуда использовать. 100К-200К рандомов хватит не на один час без повторов. Ну а если повторы это плохо, то можно данные в массивах перемешать - это всё равно получится рандомнее, чем читать random(1000); из потоков.

Тестировал RandomRange, но оказалось что они с Random глючат одинаково, а значит в подключении ради неё модуля math нет никакого смысла. Только запас данных, сделанный "в благословенные времена" не глючит.

[Vadim]

... то выхлоп становится "мягким и пушистым", практически идеальным.

Вы не поверите, но подобные вещи (я имею в виду явную или скрытую где-то в другом модуле многопоточность) наблюдаются и в других языках тоже. Например, попросили меня посмотреть прогу на VB, которая выдавала эксцепшен. Я её "умаслил" в некоторых местах логами, т.е. сделал лог, чтобы выяснить, что-же проге не нравится. И прога перестала выдавать эксцепшен...
Конкретно по поводу Random(). Несмотря на то, что в современном FP используется более-менее продвинутый метод Мерсен-Твистер (или вихрь Мерсена), трудно ожидать, что он будет генерировать величины равномерно в указанном диапазоне. Как раз для подобной фигни придумали отдельный алгоритм равномерного распределения.

[Сквозняк]

Конкретно по поводу Random(). Несмотря на то, что в современном FP используется более-менее продвинутый метод Мерсен-Твистер (или вихрь Мерсена), трудно ожидать, что он будет генерировать величины равномерно в указанном диапазоне. Как раз для подобной фигни придумали отдельный алгоритм равномерного распределения.

У однопоточного хелловорда имеющийся рандом работает с приемлемым качеством. Хоть неидеально, но с ним можно работать. Но при усложнении программы рандом портится. Пока что других способов поддержать рандом в нормальном состоянии, кроме как запасаться данными заранее, а потом их тасовать, например порченым рандомом, пока не нашёл. Тут не до идеальных алгоритмов, пусть хоть приемлемые будут.

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

[Снег Север]

У многопоточного надо при каждом старте потока вызывать Randomize. А может даже и при каждом обращении к неосновному потоку.

[iskander]

А разве где-то утверждается, что Random() потокобезопасен?

[zub]

рандом емнип не потокобезопасен

>>У многопоточного надо при каждом старте потока вызывать Randomize
нет, это будет уже не рандом

[runewalsh]

Random работает на глобальных переменных и потому непотокобезопасен (и нет, threadvar были бы ещё хуже), я это одной упоротой личности пытался втолковывать на багтрекере, когда предлагал заменить алгоритм на что-нибудь нормальное. Потому что вихрь Мерсенна СЛИШКОМ жирный, 2 Кб состояния — это несерьёзно, разобрано здесь: https://www.pcg-random.org/posts/too-big-to-fail.html.

Сделай свой рандом как объект (type Random = record ... procedure Init(const seed); function NextUint32: uint32; ...), создавай его локально (var rng: Random;) и дёргай методы (rng.NextUint(10)) — это будет потокобезопасно.

В качестве внутреннего алгоритма рекомендую Xoshiro128**/32 или PCG64/32, вот мои реализации:

Код: Выделить всё

{$modeswitch advancedrecords}
type
   Pcg64_32 = record
      procedure Setup(const seed: uint64);
      function Next: uint32;

      // Even though delta is an unsigned integer, we can pass a signed integer to go backwards, it just goes "the long way round".
      procedure Jump(const delta: uint64);
   const
      GoldenJump = uint64(11400714819323198485);

   private
      state: uint64;
   const
      Multiplier = uint64(6364136223846793005);
      Increment = 3; // Может быть любым нечётным числом. В оригинале Increment неконстантен, но это вводит в заблуждение, см. https://pcg.di.unimi.it/pcg.php.
   end;

{$push} {$rangechecks off} {$overflowchecks off}
   procedure Pcg64_32.Setup(const seed: uint64);
   var
      i: integer;
   begin
      self.state := seed;
      for i := 0 to 1 do Next; // По идее достаточно 1 раз, но тогда seed=0 выдаёт первым Next'ом 0, что некрасиво. :(
   end;

   function Pcg64_32.Next: uint32;
   var
      oldState: uint64;
   begin
      oldState := state;
      state := oldState * Multiplier + Increment;
      result := RorDWord(uint32(((oldState shr 18) xor oldState) shr 27), oldState shr 59);
   end;

   procedure Pcg64_32.Jump(const delta: uint64);
   var
      deltaLeft, accMul, accPlus, curMul, curPlus: uint64;
   begin
      // Это общий код для продвижения любого линейного конгруэнтного генератора на основании числа шагов (delta),
      // состояния (state), множителя (Multiplier) и инкремента (Increment).
      deltaLeft := delta;
      curMul := Multiplier;
      curPlus := Increment;
      accMul := 1;
      accPlus := 0;
      while deltaLeft <> 0 do
      begin
         if deltaLeft and 1 <> 0 then
         begin
            accMul *= curMul;
            accPlus := accPlus * curMul + curPlus;
         end;
         curPlus := (curMul + 1) * curPlus;
         curMul *= curMul;
         deltaLeft := deltaLeft shr 1;
      end;
      state := accMul * state + accPlus;
   end;
{$pop}

type
   // http://prng.di.unimi.it/
   // Xoshiro128**: http://prng.di.unimi.it/xoshiro128starstar.c
   // Xoshiro256**: http://prng.di.unimi.it/xoshiro256starstar.c

   // Готовые полиномы прыжков, которых нет в оригинальных реализациях, взяты отсюда:
   // http://peteroupc.github.io/jump.html
   //
   // Для Xoshiro128** — 2^32, 2^48, (2^64), (2^96).
   // Для Xoshiro256** — 2^32, 2^48, 2^64, 2^96, (2^128), 2^160, (2^192), 2^224.
   //
   // Оттуда же полиномы для прыжков на (период_генератора / φ), где φ — золотое сечение.
   // Утверждается, что повторные применения таких прыжков, врапающихся через полный период, дадут состояния,
   // отстоящие друг от друга на расстояние, близкое к максимально возможному.
   //
   // Там же описывается способ выполнять произвольные прыжки, но это очень сложна.

   // The state must be seeded so that it is not everywhere zero.

   Xoshiro128ss_32 = record
   type
      StateVec = array[0 .. 3] of uint32;

      procedure Setup(const seed: uint64);
      function Next: uint32; {$ifdef inline_nexts} inline; {$endif}
      procedure Jump(const poly: StateVec);

   private
      s: StateVec;

   public const
      PolyGolden: StateVec = (uint32($d1cf96b6), uint32($8fda8fd5), uint32($59016992), uint32($338b58d0)); // 0x338b58d0590169928fda8fd5d1cf96b6
      Poly2p32: StateVec = (uint32($f7afe108), uint32($f3be07b8), uint32($730b948d), uint32($0f8aed94)); // 0xf8aed94730b948df3be07b8f7afe108
      Poly2p48: StateVec = (uint32($cb56667c), uint32($87a4583d), uint32($dec5bb9a), uint32($deaa4ca2)); // 0xdeaa4ca2dec5bb9a87a4583dcb56667c
      Poly2p64: StateVec = (uint32($8764000b), uint32($f542d2d3), uint32($6fa035c3), uint32($77f2db5b)); // 0x77f2db5b6fa035c3f542d2d38764000b
      Poly2p96: StateVec = (uint32($b523952e), uint32($0b6f099f), uint32($ccf5a0ef), uint32($1c580662)); // 0x1c580662ccf5a0ef0b6f099fb523952e
   end;

{$push} {$rangechecks off} {$overflowchecks off}
   procedure Xoshiro128ss_32.Setup(const seed: uint64);
   var
      i: integer;
   begin
      if seed = 0 then
      begin
         Setup(High(uint64));
         exit;
      end;

      s[0] := Lo(seed);
      s[1] := Hi(seed);
      s[2] := Lo(seed);
      s[3] := Hi(seed);

      // У Xoshiro проблемы с декорреляцией, затягивающиеся на примерно 30 шагов: https://www.pcg-random.org/posts/xoshiro-repeat-flaws.html
      // Поэтому, чтобы seed=N и seed=N+1 разошлись, желательно промотать хотя бы на 30 состояний вперёд, а лучше больше.
      for i := 0 to 63 do Next;
   end;

   function Xoshiro128ss_32.Next: uint32;
   {var
      t: uint32;
   begin
      result := RolDWord(s[1] * 5, 7) * 9;
      t := s[1] shl 9;

      s[2] := s[2] xor s[0];
      s[3] := s[3] xor s[1];
      s[1] := s[1] xor s[2];
      s[0] := s[0] xor s[3];
      s[2] := s[2] xor t;
      s[3] := RolDWord(s[3], 11);
   end;}
   // Free Pascal очень тупой, вручную расписанный вариант намного быстрее.
   var
      s0, s1, s2, s3: uint32;
   begin
      s0 := s[0];
      s1 := s[1];
      s2 := s[2] xor s0;
      s3 := s[3] xor s1;
      result := RolDWord(s1 * 5, 7) * 9;

      s[0] := s0 xor s3;
      s[1] := s1 xor s2;
      s[2] := s2 xor (s1 shl 9);
      s[3] := RolDWord(s3, 11);
   end;

   procedure Xoshiro128ss_32.Jump(const poly: StateVec);
   var
      ns: StateVec;
      i, ip, b: SizeInt;
   begin
      FillChar(ns, sizeof(ns), 0);
      for ip := 0 to High(poly) do
         for b := 0 to bitsizeof(uint32) - 1 do
         begin
            if poly[ip] and (uint32(1) shl b) <> 0 then
               for i := 0 to High(ns) do
                  ns[i] := ns[i] xor s[i];
            Next;
         end;
      s := ns;
   end;
{$pop}

[Vadim]

Есть такая интересная штука - переменная RandSeed - это начальные данные для генерации Random(). Randomize как раз занимается тем, что получает кол-во тиков у процессора в эту переменную. Это так называемая начальная неопределённость (или источник энтропии). Можно эту переменную заменить на свою и заполнять другим способом, например, пихнуть туда текущую скорость сети (или взять текущий долг США, который растёт с сумашедшей скоростью и вовсе не линейно. ). Тогда это будет потокобезопасно.
Естественно каждый порождаемый поток должен иметь свою переменную RandSeed (во встроенной, скорее всего 0 в начале. Простите, никогда не проверял... ).

[iskander]

Естественно каждый порождаемый поток должен иметь свою переменную RandSeed

Неа, у всех потоков она одна и та же, к сожалению.

Сделай свой рандом как объект (type Random = record ... procedure Init(const seed); function NextUint32: uint32; ...), создавай его локально (var rng: Random;) и дёргай методы (rng.NextUint(10)) — это будет потокобезопасно.

Имхо, тут есть нюанс(с).
Программный генератор последовательностей, похожих на случайные, должен хранить своё состояние между вызовами Random(). Иначе его требуется каждый раз инициализировать чем похожим на элементы случайной последовательности. А зачем тогда эта лишняя сущность - генератор?

[Vadim]

Неа, у всех потоков она одна и та же, к сожалению.

Если делать свой ГСПЧ, а такое для генерации равномерного распределения, да ещё и, вдобавок, потокобезопасного, по-любому придётся делать, то предусмотреть собственную хранилку энтропии - это прямо- таки сомо собой напрашивается. Я как раз такой случай и имел в виду. А название RandSeed - это чтобы все понимали, о чём речь.

[runewalsh]

А зачем тогда эта лишняя сущность - генератор?

Ты же сам и ответил, зачем — хранить состояние. Такие генераторы может создавать, и удерживать нужное время, каждый, кому они нужны — например, если мы выносим работу в поток, её персональный генератор можно держать в поле соответствующего наследника TThread.

>инициализировать чем похожим на элементы случайной последовательности
Не обязательно, принято делать так, чтобы получались независимые генераторы при инициализации любыми разными числами, хоть 0, 1, 2. Поэтому обычно инициализируют наносекундным временем — оно де-факто всегда меняется между 2 вызовами и очень нескоро повторится. Можно xor'нуть его с uint64(ThreadID) shl 32. :)

Есть другой вариант (я прочитал о нём у автора Xo(ro)shiro) — явно разделить весь период генератора на части. Многие генераторы предоставляют возможность прыгнуть вперёд на большое расстояние — у меня это методы .Jump. Можно начать с одного генератора и раздавать его версию, пропрыгнутую на разные расстояния.

Например, для выполнения 1000 параллельных единиц работы можно взять 1 состояние генератора Xoshiro128**/32 (период = 2^128) за начальное, и раздавать его версии, пропрыгнутые на 2^96:

Код: Выделить всё

var
　　state: Xoshiro128ss_32;
　　workerStates: array[0 .. 999] of Xoshiro128ss_32;
　　i: SizeInt;
begin
　　state.Setup(0);
　　for i := 0 to High(workerStates) do
　　begin
　　　　state.Jump(state.Poly2p96);
　　　　workerStates[i] := state;
　　end;
end;

Эти генераторы гарантированно выдадут последовательности 2^96 чисел, непересекающиеся с другими. Таких генераторов можно получить 2^32 штук (4 миллиарда), прежде чем пропрыгается весь период (2^96 × 2^32 = 2^128). По аналогии можно прыгать на 2^64 и получить до 2^64 гарантированно непересекающихся последовательностей этой длины.

Если количество требуемых состояний заранее неизвестно, можно прыгать на золотое сечение от периода — у меня это GoldenJump или PolyGolden. Не знаю математики, по которой это работает (там что-то с иррациональностью и рациональными приближениями), но оно реально будет делить весь период генератора на примерно равные, всё уменьшающиеся части.

[Сквозняк]

Random работает на глобальных переменных и потому непотокобезопасен (и нет, threadvar были бы ещё хуже), я это одной упоротой личности пытался втолковывать на багтрекере, когда предлагал заменить алгоритм на что-нибудь нормальное. Потому что вихрь Мерсенна СЛИШКОМ жирный, 2 Кб состояния — это несерьёзно,

Да как бы покласть на два килобайта. "Не жирные" С++ программы настолько деградировали, что их авторы всё чаще забивают на их нативную сборку в линуксе и переходят на сборку для флатпака, как это произошло с PPSSPP https://www.ppsspp.org/downloads.html Просто установить для С++ программы правильный линукс уже недостаточно, теперь ей и правильный флатпак подавай, в котором куча фич программ может сдохнуть. Если бы починили обычный Random, то и покласть на эти 2Кб и можно ещё добавить, для инициализации не штатным Randomize, а гипотетическим RandomizeNishtjak. Понадобился потокобезопасный Random, запустил RandomizeNishtjak и после этого всем потокам стало раздаваться из нового, правильного корыта и не нужно объекты над капотом плодить.

[runewalsh]

Да как бы покласть на два килобайта.

При генерации сразу пачки значений это потенциально два килобайта L1-кэша, потраченные на ерунду, скажем так.
Плюс, я, кажется, привёл статью, которая объясняет, чем плохи 2 Кб именно применительно к ГПСЧ, а не «а вы на C++ посмотрите!! и вообще у меня гигабайты памяти!!!».

не нужно объекты над капотом плодить

А, ну да, забыл, с кем общаюсь. Обрати внимание, что речь о record'ах, а не class'ах, так что если это «плождение объектов», то объявить 1 переменную типа uint64 (состояние PCG64/32) или 4 переменных типа uint32 (состояние Xoshiro128**/32) — это «плождение переменных».

[Сквозняк]

При генерации сразу пачки значений это потенциально два килобайта L1-кэша, потраченные на ерунду, скажем так.

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

А, ну да, забыл, с кем общаюсь. Обрати внимание, что речь о record'ах, а не class'ах,

Сделай свой рандом как объект

Я вижу что в реализации объектов нет, но это слово произнёс здесь не я. И желания их впихнуть в базовую функцию и исковеркать {$mode fpc} приветствовать не могу. Это решение не идеально, а лишь отпихивает проблему на некое отдалённое будущее, как бревно подвешенное на верёвке, потому что всё равно предлагается брать данные как бы из массива, только кладёт их туда формула, а значит, если её угадать, то по фрагменту данных, есть вероятность угадать из какого места массива берутся данные и какие будут следующие. От такой экономии килобайт, если штатный рандом заменить на спрессованный в формуле массив, компилятор потеряет широту специализации.

Добавлено спустя 17 минут 34 секунды:
В принципе, штатный рандом нужно так выдрессировать, что-бы он выжирал нужные ресурсы, но генерировал данные в массив обычной памяти качественно, а потом их оттуда доставал, а по их окончании повторял цикл. Во что необходимо FPC - специальная процедура для инициализации рандома, после которой он станет работать по такому сценарию. Если не надо - так не инициализируй.

[iskander]

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