Как стать автором
Обновить

Комментарии 33

с интересом прочитал, но по коду чувства какие то смешанные…
не люблю почему то записи вида

c += pi+t

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

Понимаю, что для большинства эта запись понятна, но все же я поклонник ясного кода вообще, а не для большинства. Да, повторю, многим этот код ясен не менее, чем его более развернутая запись и тем не менее — не люблю такое в коде. Или старая школа или что… не знаю, скорей всего предпочтение ясности перед соглашениями ))))

ПС: все сказанное — ИМХО и не умаляет интересности самого материала

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

да это не замечание, не поймите меня правильно )))
Это просто некое ворчание, что-ли… даже брюзжание )))

Вы — автор кода, вы и пишете так, как считаете нужным.
В своем коде я такими конструкциями предпочитаю не оперировать, вот и ворчу )))

мне проще написать:

a = a + delta

и это я считаю более янсым, чем подобные сокращенные варианты записи… тем более, если подсчитать символы, то «перерасход» — копеечный, а ясность — выше

Но это ИМХО

Лет 20 назад один из преподавателей привёл нам пример, когда запись вроде вашей выглядит нелепо:


Matrix[NUM_OF_ROWS-1][NUM_OF_COLUMNS-1] = Matrix[NUM_OF_ROWS-1][NUM_OF_COLUMNS-1] + 1

Это просто к слову.


Но вообще я за то, чтобы комментарии к статье были бы только по теме статьи. Считаю, что нужно ценить время и внимание читателей, не привлекая их к себе отстранёнными от темы соображениями. Поговорить можно в ЛС, можно с друзьями, но не с десятками тысяч людей (именно столько откроют статью), которые не для этого сюда приходят. Это моё субъективное мнение :)

Выходит вы пишите
for (int i=0; i < n; i=i+1)
вместо
for (int i=0; i < n; i++)
?
Я свидетель префиксного инкремента и мне больно смотреть на оба примера :)
А можно сделать больно всем
for (int i=0; i < n; (i % 2) ? i++ : ++i)
Простите!
Я оптимизировал немного
for (int i=0; i < n; (i & 1) ? i++ : ++i)

Явная запись намерений программиста мне всегда казалась более читабельной

хм, а что тут неявного? «увеличиваем c на pi+t»

Не могу похвастаться тем, что все понял, но интересно.
если принять во внимание научно-популярный характер изложения и бытовую сферу применения описанных знаний.
С научно-популярным характером согласиться можно. А вот насчет бытовой сферы применения… Это в какой именно части бытовой сферы обычно применяют?
Это в какой именно части бытовой сферы обычно применяют?

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


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

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

А если серьезно… Наверное, очень научная должна быть наука, чтобы обычных вычислений с double не хватало. Это же, как я понимаю, точность каких-то исходных данных должна быть сопоставима с ошибкой ограниченности размера double. Это что-то вроде 10^(-14) или 10^(-15). Не говорю, что вопрос не актуальный, просто интересно, откуда такими интегралами ключи достают.

Здесь всё дело в том, что у нас с вами нет точного определения что такое научная и бытовая сферы. В моём миропонимании бытовое — это всё что делается просто для удовлетворения минимального набора потребностей в той или иной сфере, возможно, на скорую руку, по-быстрому, на коленках, назовите как угодно.


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


А точность может понадобиться наибольшая в том случае, когда полученный результат затем подвергается какой-то обработке, после которой погрешность возрастает ещё больше, затем ещё — и так много миллионов раз. Когда речь идёт об устойчивости какого-то решения, ошибка даже в 2-3 последних битах (7 ULP) в исходном числе может после миллиона последующих операций перевести систему из состояния сходится к состоянию расходится.


Короче говоря, для меня быт — это не только ключи доставать, это ещё и нудная тягомотина, которая предшествует реальному научному поиску. Вся эта тягомотина — это быт учёного и занимает массу времени.


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

Лет 20 почти не пользовался математикой, но попробую привести пару возможных примеров.

Если хочется по-исследовать какую-то модель, заданную дифференциальными/интегральными уравнениями и «погулять» по ней по различным траекториям, то для траекторий, проходящих близко от точек бифуркации, нехватка точности может радикально исказить результат. А процессы, происходящие в районе резонансов могут быть интересны и в акустике и в радиоэлектронике и в механизмах.

Или «бытовой» пример. Сейчас играю в игрушку Kerbal Space Program с модом Interstellar Extended.
Космическая баллистика в игре упрощенная — по Кеплеру — аппарат считается находящимся в сфере влияния одного небесного тела и летает по эллипсу/параболе/гиперболе. Для обычных двигателей можно считать, что они работают только в точке маневра, и там всё достаточно просто. Но. В этом дополнении есть двигатели микро-тяги (солнечный парус, ионные двигатели). И, для определения, когда и на сколько месяцев я должен включить эти двигатели, чтобы попасть, например, с Кербина (Земли) на Джул (Юпитер) — будет много «крючёчков». Причем для обеспечения торможения гравитационным маневром об какой-нибудь из спутников, с довольно высокой точностью нужно вовремя попасть в определенное место.

В этой игре баллистика упрощенная — решение можно частично искать аналитически и точности double, скорее всего, хватит. А вот если бы я подобное попытался сделать в Orbiter-е, то, как и в реальном мире, мне пришлось бы «пошагово» считать траекторию аппарата, на движение которого влияет добрая половина планет Солнечной системы, с накоплением ошибок либо в результате слишком крупного шага, либо в результате количества шагов.

К слову, когда читал об игре, встретился любопытный термин «Межпланетная транспортная сеть»(Википедия).
Наверное, очень научная должна быть наука, чтобы обычных вычислений с double не хватало. Это же, как я понимаю, точность каких-то исходных данных должна быть сопоставима с ошибкой ограниченности размера double.
Дело не точности исходных данных, дело в накоплении погрешности при многократных вычислениях. Например, при вычислении FFT алгоритмом Блюстейна на достаточно больших размерах (>50К) точности double вам не хватит даже для исходных данных в single и шумовая полка на высоких частотах может превысить 0 дБ.
Дело не точности исходных данных, дело в накоплении погрешности при многократных вычислениях.
Именно об этом и не подумал, согласен.
Osnovjansky в общем об этом же сказал.
с накоплением ошибок либо в результате слишком крупного шага, либо в результате количества шагов

Да и я почти тоже самое написал:


А точность может понадобиться наибольшая в том случае, когда полученный результат затем подвергается какой-то обработке, после которой погрешность возрастает ещё больше, затем ещё — и так много миллионов раз. Когда речь идёт об устойчивости какого-то решения, ошибка даже в 2-3 последних битах (7 ULP) в исходном числе может после миллиона последующих операций перевести систему из состояния сходится к состоянию расходится.
С формальной стороны да, почти то же самое или совсем то же самое.
Но с другой стороны… Наверное, с педагогической, хотя более правильно будет сказать «с моей личной», Ваше объяснение более «общее»:
полученный результат затем подвергается какой-то обработке, после которой погрешность возрастает ещё больше
С позиций того, кто сталкивался с подобными вещами в обозримом прошлом, все объяснения примерно об одном. Но этот человек не стал бы задавать детские вопросы)
Да, Refridgerator написал примерно про то же самое, но со своим примером. Ну и его сообщение уже читалось как более понятное, поскольку выше уже сказали про шаги интегрирования (баланс между ошибкой от мелкого шага и крупного, Osnovjansky), а это, наверное самый очевидный пример, который здесь возможен.
После серии хороших уроков, а также после портирования на свой процессор пакета для вычисления с quad-double точностью, уже совсем иначе читаешь следующий урок :)
Может кто встречал в сети какие-то отчеты по быстродействию различных процессоров на тестах с quad-double точностью? Было бы интересно посмотреть.
Спасибо, автор! Отличный материал.

Хочу заметить по поводу этого:
Поскольку значение полинома на случайных тестах очень быстро возрастает

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

Кстати, для вычисления полиномов Чебышева, кажется, есть специальные, оптимизированные методы. Может быть, к ним тоже можно применить описываемые вами трюки?

Далее, интерес представляет задача вычисления интерполяционного полинома (полином Лагранжа). Там тоже вычисления подвержены ошибкам округления, а в некоторых задачах нужно вычислять интерполяционные полиномы высокой степени (порядка 10000). Есть специальные трюки для повышения точности (барицентрическая формула), но и они упираются в пределы точности double. Может быть, описываемыми вами методами можно и здесь что-то улучшить?

Задача перемножения полиномов (свёртка). Уже при степени полиномов порядка десятков, «обычные» алгоритмы теряют точность безнадёжно. Я использовал трюки с FFT, помогают. Но, может, можно сделать ещё лучше.

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


Описываемые алгоритмы можно применить везде, где есть необходимость перемножать пару чисел и складывать цепочки чисел любого размера. Но даже если где-то есть деление, оно может быть заменено серией сложений и умножений. Так что можно всё, было бы желание и необходимость.


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

Посит стандарт обновился, его выложили в маил лист. Теперь es=2 для всех битностей. В итоге 16-битный посит имеет quare аккумулятор, который сможет вместить ±10^305 диапазон — примерно такой диапазон у типа double. Что вы думаете об этом?

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

Спасибо что делитесь своим опытом)
Может быть подскажете что-нибудь подобное для разложения Холецкого? Работаю с матрицами ковариации 100*100, и точности double не хватает. Приходится пользоваться длинной арифметикой, что сильно замедляет код

Пожалуйста.


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


Если нужен конкретный ответ, то мои консультации платные, обращайтесь в ЛС. Но я вас уверяю, дешевле будет вам попробовать сделать самостоятельно, это нетрудно.

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

В случае Холецкого, у вас всё равно вычислительная сложность N^3. Попробуйте LU-разложение, SVD или хотя бы метод Гаусса.

Трюки симпатичные, но наверное стоит указать что есть стандартные реализации этого подхода (long double, double-double, float128 в зависимости от реализации) https://en.m.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic

Стандартных реализаций нет, есть чьи-то реализации давно известных приёмов. И моя задача всего лишь в их популяризации. При этом long double к данной теме не относится из-за невозможности применить его так же просто как остальные типы. Более того, сам по себе double-double или float128 не позволяет создать целый каскад сложений с повышением точности до любых значений, а описанные в статьях приёмы — позволяют.

Еще раз спасибо за статью. Я как-то пропустил в свое время существование таких несложных, но полезных, трюков.
Скажите, у Вы не пробовали точное скалярное произведение в методе сопряженных градиентов? Дают ли такие алгоритмы улучшение сходимости в практических задачах?
Вы не пробовали точное скалярное произведение

Поправлю: не точное, а наиболее точное, это разные вещи. Да и то, всегда можно подобрать тест, когда даже эти трюки спасать не будут.


в методе сопряженных градиентов

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

Придется попробовать :) Странно, но в кодах тех открытых библиотек линейной алгебры, в которых я копался не было попыток использовать «наиболее точное» произведение в итерационных методах. Хотя известно, что это один из основных источников проблем.
Трюки действительно интересные, но очень затратные в плане тактов процессора. Без FMA очень проблематично использовать two_prod (); На своем процессоре я влез только в 17 тактов. А если все это еще делать через вызов-возврат функции, экономя код, то уже задумаешься. С two_sum () тоже 8 тактов (с учетом пузырей). Так что если в итоге все потом посчитать, то уже и другие алгоритмы могут быть привлекательными.
Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.