adeshere
@adeshere
РАН, Фортран, временные ряды

Генератор случайных чисел ИНОГДА (очень редко!) возвращает NaN?

Я занимаюсь научными расчетами, и внезапно обнаружил очень странный эффект. Если вызывать интеловский (встроенный в Intel Fortran) генератор равномерно распределенных случайных чисел (0,1) много раз подряд, то на некоторых компах (пока таких 4 из 5 протестированных) примерно через A=50 млн вызовов (точная цифра меняется, подробнее см. по ссылке ниже) он начинает квазипериодически, с периодом около B=10 млн вызовов, возвращать NaN.

UPD: Благодаря наводящему вопросу Tony, хочу уточнить один важный нюанс: даже при генерации одинаковых последовательностей ПСЧ (инициализация SEED одной и той же константой) NAN-ы при каждом запуске появляются в разных местах.
UPD: вот подробности (картинка с результатами теста)

На этой картинке три ряда.
Первый и второй ряды (зеленый и желтый) - две это последовательности нормально распределенных ПСЧ, полученные при двух запусках тестовой программы с одинаковым SEED(301). Каждое значение сигнала на графике - это сумма 12-ти значений RANDOM минус 6. Тем самым я получаю нормально распределенное случайное число, а заодно уменьшаю длину ряда в 12 раз. Если хотя бы один из 12-ти вызовов RANDOM вернул NaN, то и значение Noise на графике тоже будет NaN. Таким образом, чтобы получить 1 000 000 значений Noise, я вызываю функцию RANDOM 12 млн раз.
Третий ряд (синий) - это разность между первым и вторым рядом.
NaN-ы, как и на всех остальных картинках, закодированы числом 99. Соответственно, разность двух первых рядов дает либо 0 (если значения совпадают), либо +-99 (примерно).

61525d77af16c947014893.jpeg

Из картинки хорошо видно, что:
1) при первом запуске теста первый NaN появляется примерно после 6.2 млн вызовов функции Noise (75 млн вызовов RANDOM). После этого NaN-ы появляются примерно через 10 млн вызовов RANDOM. Но иногда (не очень часто) этот интервал может быть существенно меньше. Например, на отметках 9 и 9.5 млн.

2) А вот при втором запуске теста первый NaN появляется почти сразу же! "Мертвый период" (без NaN-ов) исчезающе мал. В исходном посте я писал, что этот "мертвый период" чаще всего близок к 60 млн вызовов, и лишь иногда он бывает существенно меньше. Но тогда я гонял тесты, инициализируя ГСЧ комбинацией даты и времени. А вот когда я стал подсовывать в seed небольшие натуральные числа (пробовал от 20 до 400), то оказалось, что при такой инициализации цикл с квазипериодическим появлением NaN может запускаться почти сразу же после старта программы!

3) Теперь смотрим на ряд разности (синий). Видно, что он почти везде равен нулю, т.е. последовательности ПСЧ действительно одинаковые.
За исключением NaN-ов...

4) Еще интересен фрагмент от 6.3 до 9.5 млн. На этом участке NaN-ы в двух последовательностях в точности совпадают! Напомню, что между последовательными NaN-ами примерно по 0.8 млн вызовов Noise, или около 10 млн вызовов RANDOM. То есть, вероятность случайного совпадения двух NaN-ов в рядах Noise примерно E-6. А в приведенном примере они совпали трижды подряд! Ясно, что случайностью это быть не может.

5) Однако, после совпадения трех NaN-ов подряд следующие NaN-ы (на отметке 9 500 000) в двух тестах разные! И дальше последовательности расходятся.

Из анализа этой и других подобных картинок я могу сделать такие выводы:
1) Момент появления первого NaN в последовательности как-то зависит от того, чем инициализирован генератор случайных чисел.
2) Появление NaN зависит как от внутренних переменных генератора случайных чисел (иначе мы бы никогда не увидели точное совпадение NaN в двух разных последовательностях), так и еще от каких-то внешних по отношению к ГСЧ факторов (иначе бы при одинаковом внутреннем состоянии ГСЧ NaN всегда появлялись бы в одинаковые моменты времени).


Будет ли аналогичный эффект наблюдаться в других языках/компиляторах, я проверить не могу - не располагаю компилятором/навыком. Но подозреваю, что у фортрана и Си базовые математические библиотеки часто общие. Даже название функции - RANDOM в современном фортране и Си одинаковое. А другие языки могут что-то таскать из Си-шных библиотек.

Для себя я пока сделал подпорку: обернул встроенный генератор случайных чисел в свою функцию, которая в случае isNan(Random) вызывает генератор повторно. Но внутренний дискомфорт остался. И хочется все же понять причину бага, а заодно и предупредить о нем других "заинтересованных лиц".

И еще один важный нюанс. Баг возникает, только если собирать программу c полной оптимизацией. У меня в фортране это ключи /MP /O2 /QaxSSE3 /QxSSE3 /Qparallel. При сборке без оптимизации результат будет NaN-free.
Меня это, мягко говоря, удивило, так как генератор случайных чисел, по идее, должен в обоих случаях вызываться один и тот же. Поскольку исходный текст функции RANDOM к моему компилятору не прилагается - вероятно, она цепляется из какой-то стандартной библиотеки - я полез в дизассемблер, чтобы понять разницу между двумя версиями генератора случайных чисел (сборка с оптимизацией и без нее). Сам генератор, как и следовало ожидать, действительно в обоих случаях одинаковый. Но внутри кода RANDOM есть
один внешний вызов

Для начала предупреждаю, что в ассемблере я ни бум-бум, и пишу только на фортране! Поэтому все написанное под этим спойлером может оказаться полным бредом и воспаленным воображением. Я всего лишь хотел попытаться, насколько сумею, собрать хоть какую-то информацию про обнаруженный баг для тех, кто проявит интерес к ситуации, но не имеет возможности открыть интел-фортран и проделать все это лично.

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

005C052B call _for__acquire_semaphore_threaded (596C80h)

Внутри которого (после череды других call) происходит что-то уже совсем-совсем для меня непонятное:

75C34043 call 75C443C5
(...)
75C443C5 jmp dword ptr ds:[75BA007Ch]

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

Но на всякий случай добавлю, что программа у меня однопоточная и, по идее, ничего согласовывать с другими потоками не должна. С другой стороны, в языке есть массивные операции, которые могут распараллеливаться автоматически, и, теоретически, генератор случайных чисел может в этом случае как-то донастраиваться для параллельных вычислений. В справке я об этом ничего не нашел, но это может быть не проблемой справки, а моей личной проблемой (справка написана по-английски, а у меня отношения с ним не сложились :-(.

Более полная информация, включая ассемблерный код в двух вариантах, собрана вот тут:
FORTRAN_URAND_BUG.doc
Это Word-файл, в который я собираю всю информацию про этот баг. Там же приведен исходный фортрановский код, а также дана ссылка на скачивание моего exe-шника с тестом. Если будет не влом - запустите его и напишите мне, что получилось. Пока из 5 протестированных компов Nan-ы появляются на 4-х, а на одном - нет. (Подробнее в том же Word-файле)


И еще одна версия

Глядя на битовые маски:

Random= 11111111110000000000000000000000 NaN
R=1.00 = 00111111100000000000000000000000 1.00000

возникает подозрение, что проблема может возникать на этапе преобразования int -->> real. Насколько я понимаю, исходный алгоритм генератора псевдослучайных чисел обычно работает с integer. Более того, он никогда не дает числа вида 00000000 или FFFFFFFF. Поэтому на выходе фортран-генератора не должно быть ни точных 0, ни точных 1.000. Но как известно, точность представления чисел с плавающей запятой в компе ограничена. Около 0 это не влияет, а вот вблизи 1 некоторые integer-ы могут преобразовываться в R=1.000
С другой стороны, если бы это было именно так, то код этого преобразования работал бы одинаково с оптимизацией и без нее (или нет?).
В общем, если вы заинтересуетесь и попробуете вызвать RANDOM из Си, то включайте максимальную оптимизацию. У меня бажит только оптимизированная программа. Правда, я не очень понимаю, можно ли включить в Си многопоточность одним ключом /MP, или надо специально цеплять для этого всякие либы. В фортране-то параллельность заложена в языке (массивные операторы) и от пользователя не требуется большого ума, чтобы ее подключить...

И еще для справки: вот что написано про ключ /О2 в справке фортрана:
UPD: к сожалению, я столкнулся с ограничением на длину вопроса и эту информацию вынужден отсюда убрать. Она теперь лежит все в том же злосчастном Word-файле, который доступен про ссылке.


Вообще, история уже получается очень длинная, и она постепенно обрастает подробностями. В один вопрос не уместишь - тут надо целую статью написать... Но уместно ли это делать, имея только вопрос без ответа?

Поэтому для тех, кому вопрос действительно актуален, прикладываю ссылку на Word-файл, в котором собраны ВСЕ подробности ситуации, которые я смог пока накопать: FORTRAN_URAND_BUG.doc.

Также я написал некоторые подробности вот в этом комментте
Он менее полный по сравнению с Word-файлом, но его проще посмотреть - не надо ничего скачивать.

Если Вам будет не влом написать свои тесты RANDOM, то было бы очень интересно узнать результаты. Только не забудьте про оптимизацию - без нее никаких глюков не будет (но это не точно ;-)
  • Вопрос задан
  • 1288 просмотров
Пригласить эксперта
Ответы на вопрос 1
@thecove
Я ветку почитал и могу сделать ввод: или магия случается, тоесть списываем все на магнитные бури на Юпитере или таки надо копать.
Если баг воспроизводится абсолютно хаотично то могу придположить ( из своего C++ опыта ) что кто то где то портит кучу и на выходе у вас фигня.
Кстати я так понимаю что Rand в фортране выдает число типа float / double потому что я не понимаю как NaN возможен в целых числах.

Вопрос: если написать простую голую программу в которой будет только один вызов функции main ( ну или что там в фортране, простите я не знаю ) и только в этой функции в цикле вызывать Rand и больше ничего.
NaN - появляется???

PS можете приложить асм листинг функции Rand ?
Ответ написан
Комментировать
Ваш ответ на вопрос

Войдите, чтобы написать ответ

Войти через центр авторизации
Похожие вопросы