?

Log in

No account? Create an account
Previous Entry Share Flag Next Entry
Программистское
goliafffff
Далее под катом немного кода на СИ и одна проблема. Если у вас есть немного свободного времени и вы программист, то я буду очень рад помощи.

Сразу скажу: я не профессиональный программист. Большую часть знаний по программированию приобрёл самостоятельно. Программы приходится писать для решения ряда задач нелинейной динамики. В ходе решения очередной задачи возникла проблема: при добавлении синуса (см. пункт 2) программа начинает вести себя по-разному на различных компьютерах. Причём, если в программе синуса нет (см. пункт 1), то на всех доступных компьютерах программа работает абсолютно идентично. Я запускал задачу на трёх компьютерах: MacBook Pro (Intel Core 2 Duo), Desktop (Intel i7, на нём, кстати, самый хороший счёт наилучшим образом совпадающий с теорией) и ещё один Desktop с двуядерным интеловским процессором (точную модификацю забыл).

Поясню немного чем я занимаюсь. В ходе своего исследования я анализирую некоторые статистические характеристики одномерного нелинейного отображения в режиме хаоса. Отображение представляет из себя обыкновенную функцию, которая при каждой итерации вычисляет следующее значение переменной x, зная предыдущее значение переменной. Если отложить по оси X итерации, а по оси Y значения переменной x, то мы получим временной ряд, статистические характеристики которого можно исследовать. Сложность заключается в том, что описанное отображение x_{n+1} = (a*x_{n} - x_{n}^3)*exp(-(x_{n}^2)/b) находится в режиме динамического хаоса. То есть если мы запустим две копии отображения и зададим мельчайшее различие начальных условий переменной x, то уже через несколько итераций значения переменных станут абсолютно разными. Строго говоря, мельчайшее отклонение начальных условий нарастает экспоненциально, и мы получим два различных временных ряда. На основе получившегося временного ряда считается некая статистическая зависимость (назовём её p = f(g) ), которая должна подчинятся определённому математическому закому.

Если в системе нет синуса, то на всех компьютерах значения рассчитанной величины p абсолютно идентичны. Но стоит добавить синус, как три компьютера показывают три различные значения p. Причём Desktop с Core i7 считает просто идеально, зависимость p = f(g) ведёт себя так как надо и подчиняется теоретическому закому, а вот на остальных двух компьютерах какая-то околесица (вроде что-то похоже на правду, но не то).

Я не знаю, как в компьютере реализована функция синуса. Могу только предположить, что синус задан разложением в ряд Тейлора и на каждом из компьютеров различное количество элементов разложения. Либо в процессоре Core i7 реализована какая-то отдельная штука для вычисления математических функций (может какой-то математический сопроцессор). К сожалению, всё это лишь мои догадки.

1

Есть ещё один интересный момент. В ходе поиска ошибки коллега предложил поиграться с аргументами синуса и сделать следующее:

2

И вот что получилось: зависимость статистической величины p = f(g), на всех компьютерах стала одинаковой, но неправильной с точки зрения математической теории! Сами значения величины p не совпали, но стали отличаться всего на несколько процентов. В чём проблема? Как с этим бороться?


  • 1
Обычно так бывает из-за того, что в программе есть неинициализированные переменные, в которых на разных машинах оказывается разный мусор. Тем не менее, the IEEE standard does not guarantee that the same program will deliver identical results on all conforming systems, так что, если вы работаете с хаосом, где это и должно происходить, то, может, вам как раз удалось сконструировать именно тот пример, где результаты должных быть совсем другие. Есть возможность запостить программу целиком?

Спасибо, что откликнулись.

Код программы довольно большой. Там много функций, и чтобы понять принцип их работы нужно иметь представление об одной математической теории. И, да, возможно есть пара неинициализированных переменных, но я ими не пользуюсь.

И спасибо за ссылку.

Edited at 2013-12-14 07:52 pm (UTC)

Попробуйте откомпилировать с выводом абсолютно всех диагностик (в gcc, например, флаги -Wall -Wextra)

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

Вроде как аппаратной реализации тригонометрических функций в современных процессорах для персоналок нет.

Немного работал с Фортраном. Там возможны сюрпризы с ключами компилятора при оптимизации под конкретный процессор, есть ключи при которых скорость выполнения арифметических операций увеличивается за счет нарушения стандарта IEEE.

Есть, в FPP. Через ряд. 32 бита точность почему-то.

Точность вещественных арифметических вычислений задается в коде. Наиболее распространенными являются одинарная точность (float, 32 двоичных разряда) и двойная точность (double, 64 разряда). Есть еще увеличенная (80 разрядов) и половинная (16). С заданием типов переменных, аргументов, констант из-за этого приходится быть аккуратным. Вроде здесь немного описано, касательно тригонометрических функций http://msdn.microsoft.com/ru-ru/library/wkbss70y.aspx

Есть вопрос как именно в коде или в конкретном компиляторе определяется константа M_PI.

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

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

Например, 3.14/10.0 уже достаточно, чтобы внести в вычисления произвольную чушь.
Потому что не делится на двоичных машинах.

Чего вы хотите добиться-то?

Edited at 2013-12-14 09:40 pm (UTC)

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

3.14/10.0 написал в целях упрощения, для текущей задачи этот параметер не принципиален, могу заменить обратно на M_PI. Почему 3.14/10.0 не делиться на двоичных машинах? Что будет в результате такого деления?

Если FPP не даёт хорошей точности, то что можно использовать вместо него? Что такое четыре режима округления?
Буду признателен, если вы ответите на мои вопросы или хотя бы дадите ссылки, где про это можно прочитать.

Пока что я только нашёл два обсуждения на stackoverflow:
http://stackoverflow.com/questions/2284860/how-does-c-compute-sin-and-other-math-functions
http://stackoverflow.com/questions/345085/how-do-trigonometric-functions-work

При вычислении 3.14/10.0*ds.t возникают ошибки округления (число 3.14 непредставимо в виде конечной двоичной дроби), и результат "прыгает" в зависимости от того, в какую сторону будет округляться итог. А Вы в результате изучаете особенности не отображения, а архитектуры компьютера. :)

Собственно, это предположение можно легко проверить. Я соорудил программку, которая считает первое отображение (в котором никаких синусов нет), но с разной точностью (двойной и четверной). Разница результатов начинается, естественно, с 10^{-16}, но уже за несколько десятков итераций доходит до 10^0.

Так что Вам просто повезло - первое отображение все используемые компьютеры считают одинаково (и, по-видимому, одинаково неверно). Во втором же случае появляются какие-то мелкие различия, которые приводят к разным результатам.

FPP, последний раз, когда я это исследовал, давал всего 32 бита мантиссы на тригонометрических функциях.
Четыре режима округления (ими можно управлять) - это к плюс бесконечности, к минус бесконечности, к нулю, к ближайшему чётному.

Вы делите на 10, т.е. на 2 и на 5. На 2 делить ошибки нет, а на 5 деление получается неточное, потому что это не степерь двойки, дробь бесконечная. и там остаток отбрасывается.

Да, в целях сохранения душевного здоровья - лучше считать что в разных "системах" - разные синусы. Под "системами" я понимаю связку не просто процессор, а связку процессор+математическая библиотека+компилятор.

Это можно победить, но лучше считать что они разные.

Ну да, динамические ситесмы они такие, но как бы и не страшно - погрешность измерительных приборов гораздо больше, чем погрешность double. Есть ещё более страшная проблема - погрешность модели. Если вы эмулируете физический процесс, то вы не учли ветер, изменение гравитации от высоты, или ещё что-нибудь.

У меня тут больше математическая задача. Если в аргументе синуса стоит M_PI/30.0*ds.t, то я получаю различные результаты. Если же заменить M_PI/30.0 на приблизительное 0.10471975511966, то всё встаёт на свои места. На сколько я понимаю корень зла кроется в ошибки округления при делении на 30.0.

Как это можно победить?

В общем случае множество действительных чисел в один double не засунешь.
Нужно вдумываться в конкретную задачу.

Вычислениями на компьютере можно проверить, что два выражения очень близки, но невозможно доказать что два выражения в точности равны.

Например, с точки зрения double ( cos(0.00000001) - 1 ) в точности равно нулю, но это не значит что 0.00000001 -- неподвижная точка для (cos(x)-1+x)

Я согласен с аргументом про общий случай, но данный случай очень частный. Если я правильно понял автора, он сводится к разному поведению деления на двух разных ieee-совместимых машинах. Я бы безусловно ожидал побитового совпадения. А сейчас бы, соответственно, распечатывал двоичное представление М_PI / 30.0 в обоих случаях.

P.S. а до этого посмотрел бы на результат работы препроцессора, во что это M_PI превращается в коде.

В сгенерированный машинный код надо смотреть ещё, M_PI/30 может вычисляться на этапе компиляции, в зависимости от флажков оптимизации компилятора

Надеюсь, что вычисляется на этапе компиляции! Кстати, будет смешно, если окажется, что деление в собственном коде работает правильно, а в кодогенераторе сломано.

Логика рассуждений мне понятна, но, честно говоря, я не знаю как посмотреть сгенерированный машинный код или результат работы препроцессора. Вообще очень странно, что на разных компьютерах разный результат. Просто если дело в округлении результата деления M_PI на 30, то в рамках моей задачи это будет означать немного другую частоту в системе, и статистические характеристики, которые я считаю могут, измениться (и то не сильно) количественно, но должны подчиняться теоретической закономерности. Не понятно почему при задании частоты через M_PI/30.0 появляется систематическая погрешность.

Edited at 2013-12-28 02:22 pm (UTC)

Пост-препроцессорный код увидеть очень легко, см. например
http://stackoverflow.com/questions/5034209/how-do-i-see-what-a-file-looks-like-after-preprocessing

Это будет обычный C, в котором развернуты все #include, #define, #if и проч. Вместо M_PI будет стоять число в десятичной записи, надо убедиться, что оно одинаковое на двух машинах.

На последний вопрос, ничего не зная про задачу, ничего путного не отвечу.

На всякий случай:

1) Есть библиотеки и программы для расчетов с произвольной точностью. "произвольная" не значит "бесконечная". Т.е. вам допустим нужен результат с точностью до 22 знаков после запятой. Вы опытным путём или на бужаке понимаете, что для этого достаточно 1000 знаков (или 10000) после запятой, ну и считаете с такой точностью.

2) Есть библиотеки для символьных вычислений, которые преобразуют формулы, а не числа.

Подходит это или нет - зависит от задачи. Может реорганизацией вычислений хватит и double, а может у вас хаотическая система, которой и "произвольной" точности не хватит (ну невозможно иррациональное число приблизить какой-то точностью), или мощности символьных вычислений не хватит (например, попробуйте доказать что последовательность n = (n mod 2 == 0) ? n/2 : 3n*1 сходится к 1 для всех целых n).

Можете дать ссылки на информацию о библиотеках с произвольной точностью?

Я сам ими не пользовался, так что буду вбивать в гугл arbitrary precision ( https://www.google.com/search?q=arbitrary+precision ) искать столько же времени, сколько и вы

Я бы почитал такое -
1) http://stackoverflow.com/questions/4486882/whats-the-best-for-speed-arbitrary-precision-library-for-c

И возможно, использовал бы не C, а привычную нотацию в
2) калькуляторе bc
3)
python и
а) его модуль decimal (только sin/cos нужно будет считать самому)
б) библиотека sympy (помимо расчетов произвольной точности, она умеет символьные вычисления, типа sympy.summation(1/(n*n), (n, 1, oo)) равно pi**2/6 )
в) библиотека mpmath


Могу переписать ваш сишный код на python по доллару за строчку :)


Для C есть практически стандартная библиотека MPFR.

  • 1