Комментарии 83
Да, с FMA тоже бывают приколы. Суть в том, что это (float)(a * b + c), а не (float)((float)(a * b) + c). То есть не происходит промежуточного округления после умножения. Поэтому результат может отличаться.
Ведь, грубо говоря, если у нас есть последовательность знаков 0.10100000, которая вычислена с погрешностью < 0.5 последнего бита, то дальше в ней либо все (до бесконечности) нули, либо где-то встретится единичка и результат будет ближе к верхнему значению округления (0.11). Вероятность второго случая кажется сильно больше первого.
Ровно так же последовательность 0.1101111111 может быть просто приближением снизу к точному значению 0.111(0), но тут мы всегда округляем вниз, считая (вполне обоснованно) что минимум один из неизвестных битов равен нулю. Почему бы не делать то же — округлять вверх хвост из нулей, считая, что где-то дальше за ними появится неизвестная сейчас единичка?
Или я где-то ошибаюсь?
Благодарю за вопрос! Вы ошибаетесь в том, что для округления вверх может быть другой пример, когда у нас, например, число 1.11000000001, и вы округлите вверх, получая 1.111, но вы заранее не знаете, этот паровоз из нулей — это уже правильный паровоз, или при последующих итерациях алгоритма на самом деле будет 1.1011111111? И тогда округление вверх должно дать 1.110. Опять ошибка в последнем бите. Короче говоря, для любого варианта округления есть своя ситуация, когда мы находимся очень близко к границе между двух направлений округления, и чтобы точно понять, с какой стороны от этой границы мы находимся, нужно чтобы паровоз из нулей или единиц точно закончился.
Ещё добавлю к цитате:
то дальше в ней либо все (до бесконечности) нули, либо где-то встретится единичка и результат будет ближе к верхнему значению округления
Нет, есть ещё третий вариант: нули превратятся в единички по ходу дальнейших расчётов. Потому что эти нули за пределами 0.5ulp. И вот если это произойдёт, тогда округление вверх даст другой результат.
UPD: Даже могу ещё точнее пример привести. Допустим, дробная часть мантиссы занимает два бита. Берём два числа, которые отличаются друг от друга меньше чем на 0.5ulp:
0.11000001
0.10111111
Наши 0.5ulp равны 0.125 (ценность последнего — второго — бита равна 0.25). Разница между числами значительно меньше. Но если округлять вверх оба, то будут разные результаты. Поэтому чтобы понять, куда же округлять на самом деле, нам нужно добраться в точности до 8-го бита. Аналогичные примеры можно придумать для любых вариантов округления.
Спасибо за напоминание, полезно.
ps:
> #define FAU(x) (*(unsigned int*)(&x))
> #define DAU(x) (*(unsigned long long*)(&x))
Не думаю, что это повлияло на результат, но вообще-то в С++ так делать нельзя, это UB (только в С можно). Доступ к памяти некастуемого типа разрешен только через char* (::std::memcpy() в частности или std::bit_cast (since C++20)).
Я с вами полностью согласен, это UB, но поясню вам свою позицию. Дело в том, что в моей работе код пишется только под одну архитектуру и если нужно написать под другую, пишется другой код, такое правило. Причём пишется заново. Поэтому те вещи, которые в обычной жизни считаются UB, у меня считаются нормой и всегда работают в точности так, как я себе это представляю. То же касается битовых трюков, которые также нельзя было бы использовать. Ещё я пишу на этакой смеси C и C++, это тоже обусловлено постоянной работой на низком уровне. В общем, не судите строго, благодарю за замечание!
Здесь мне важно, чтобы читатели поняли суть, а как они это реализуют для себя — это уже пусть смотрят сами.
Есть МК51, имеет 100% совпадение с проверкой на бумаге, но иногда не совпадает с GCC, и для этого есть простое объяснение. Этот древний мамонт считает в десятичной системе!!!
В двоичной системе для округления достаточно отбросить лишнее — без анализа содержимого. Вы уже выбрали необходимую вам точность, остальное туда просто не влезет.
В десятичной системе так не получится — один младший бит может влиять на старший десятичный разряд, даже без округления (разговор про одинарную и двойную точность).
Именно по этой причине округление чисел для вывода в десятичном варианте — нужно проводить уже после перевода из двоичного в десятичный. А все предварительные манипуляции с последним битом признать ересью, и
В десятичной системе так не получится — один младший бит может влиять на старший десятичный разряд, даже без округления (разговор про одинарную и двойную точность).
Да ладно? Такое возможно только в "пограничных" случаях, и такие пограничные случаи в десятичной системе встречаются в 5 раз реже чем в двоичной!
А вот для одинарной и двойной точности в двоичном варианте — огромное количество. Тут я лично споткнулся, когда писал свой вариант printf. Либо получается быстро с уникальными ошибками перевода, либо точно- но ещё медленные чем стандартная функция.
Кстати, printf для устранения ошибки перевода — почти двух кратно повышает точность. И его перевод считается образцовым, но не идеально точным. Потому как обратное преобразование достаточно часто даёт ошибку.
При чём тут вообще битовое представление? Преобразование между системами счисления — это далеко не единственное место, где может образовываться погрешность.
Да, если у нас нужно сложить или перемножить два числа — это моно сделать в десятичной системе счисления и получить точный результат. Но уже с делением всё сложнее, а начиная с извлечения квадратного корня разница в какой системе счисления работать пропадает.
При чём тут вообще битовое представление? Преобразование между системами счисления — это далеко не единственное место, где может образовываться погрешность.
При том что это очень эффектный способ демонстрации погрешности. Именно демонстрации!!! Когда младший бит влияет на старшую десятичную цифру.
Но вы-то изначально говорили не про демонстрацию проблемы, а про различие между десятичной и двоичной системами счисления. Которых на самом деле нет.
Десятичная арифметика 0,1 + 0,2 = 0,3. Число 0,2 представлено как 2e-1 в BCD формате, после двойки все остальные числа равны нулю. Как учили в школе.
Нужно просто перестать натягивать сову на глобус. Методы округления чисел в двоичном и десятичном варианте — должны выполняться в пределах формата числа, без принудительного перевода формата.
Благодарю! Вы правильно понимаете, если имеете в виду, что не само решение должно быть индивидуальным, а трудоёмкость округления (величина m) для каждой функции и для каждого режима округления подбирается индивидуально. Если мы точно знаем величину m, то есть то, до скольки битов после запятой нужно рассчитать промежуточный результат (до округления), то мы гарантированно получаем правильный ответ после округления, при этом верифицировать ничего не нужно, ответ будет гарантированно правильным сам по себе. Просто по определению величины m.
Сложность в том, чтобы эту величину m отыскать. Я привёл вариант решения, когда её ищут полным перебором с отсечениями. Для binary128 полный перебор как-то не хочет работать даже с отсечениями. Поэтому ждать вторую статью на эту тему от меня пока не надо, эту открытую научную проблему ещё никто не решил, и я пока тоже за неё не взялся :)
Еще полезно знать, что для чисел с плавающей точкой данные утверждения бывают истины:
- a + b — b != a
- a * b / b != a
- (a + b) + c != a + (b + c)
И т.д.
Благодарю, но это далеко не все особенности, о которых полезно знать. У меня на канале есть по этому поводу информация, если интересно. Да и вообще в интернете много про это пишут.
Есть ещё полезная статья David Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, 1991
Денежные расчёты принято делать в десятичной арифметике с плавающей запятой. Если будет интересно, объясню почему.
Ну можно и с фиксированной, однако есть же Стандарт IEEE-754, который описывает типы данных Decimal64, например как раз для этого дела. Причина проста. В финансовой математике действуют законы государства, прописанные в десятичной системе счисления. Например, если написано по закону, что, скажем, налог составляет 12.4 рубля, значит должно быть ровно столько, а когда мы воспользуемся типом данных binary64, будет (ближайшее число):
12.4000000000000003552713678800500929355621337890625
А это уже прямое нарушение закона. Когда выполняется множество действий с такими вот числами (видели, наверное, эти ипотечные формулы с возведением в 365-ю степень и прочее), тогда погрешность будет нарастать, и нарушение закона будет нарастать ещё больше. Как-то так, если кратко.
Ещё лучше их делать в десятичной системе с фиксированной запятой. Но это уже не любой ЯП поддерживает...
Верно, но бывают сложные формулы, где в промежутке могут возникать большие числа (возведение процентов в 365-ю степень, например) и запятая будет плавать. Для этого существует тип данных Decimal64, прописанный в Стандарте IEEE-754.
В таких формулах десятичная система — что мёртвому припарки. Либо мы считаем точно (и тогда нам нужна фиксированная точки либо вовсе рациональные числа) — либо мы считаем с точностью до половины младшего разряда, но тогда нам и двоичная система норм.
Достоинство Decimal64 — в том, что большую часть времени он работает как будто имеет фиксированную точку, а когда мы хотим возвести проценты в 365-ю степень, что точно сделать невозможно по вами же перечисленным причинам — Decimal64 автоматически переходит в "примерный" режим. Ну и в итоге всё как бы всех устраивает, пока какая-нибудь ерунда вроде переполнения мантиссы при сложении не произошла.
Благодарю за комментарий! Вы говорите о совершенно другой проблеме. Вы говорите о том, насколько точным будет итоговый результат после множества операций с числами с плавающей запятой. Эти вещи просчитываются аналитически и описаны в книге [3] (там есть оценки на погрешность при сложениях, умножениях и т. д.). Я же говорю о том, что вычислять значения элементарных функций нужно точно до последнего бита и причина необходимости делать это описана в последних абзацах.
А вот для обычных прикладных программистов эта проблема (описанная в статье) не важна. Им без разницы, что там с битами, им важно, чтобы конечный ответ был бы с определённой точностью. Между 51 и 52 битами точности для их целей может не быть разницы вообще.
им важно, чтобы конечный ответ был бы с определённой точностьюА разве бывает по-другому? У нас так и так точное значение недоступно, так и так значения округляются — причём в логарифмическом, а не линейном масштабе.
Опять мы говорим о разном. Я же приводил пример с числами 1,999999 и 2,0. Прикладному программисту не важно, что у него получилось "почти два, но не совсем". Он округлит — и дело с концом. Разработчику библиотеки ВАЖНО, чтобы результат был корректно округлён до последнего бита, ему ВАЖНО получить именно 2,0 (в данном примере).
Насколько принципиальна разница между 51 и 52 битами точности?
Дело касается переносимости результатов вычислений.
Стандарт есть, но в нём нет жёсткого требования о корректном округлении последнего бита (об этом сказано в статье), потому что обеспечить такое округление наука пока не в силах для всех случаев.
Основные сложности у программистов возникают с округлением «к ближайшему, но в случае равного удаления от ближайших — к тому, у которого последняя цифра чётная». Да, именно так переводится этот режим округления, который западной литературе называют короче: «Round nearest ties to even».
Сверился с википедией на которую вы даёте ссылку, что бы проверить ваше «неарест тиес то евен», так вот, ссылка по википедии под таким названием ведёт на «Round half to even», что является банковским округлением, в котором чётная цифра смотрится не после, а до округляемого разряда и никаких сложностей не представляет. И мне понятно о чём вы говорите: если для математического округления ">=" это единый случай и полностью покрыается одним разрядом, то для банковского, где разделяют ">" и "=", и в случае равенства делаются лишние телодвижения, то для трансцендентных чисел само понятия равно в дробной части не имеет смысла, потому что числа не кончаются и после любого количестава нулей может пойти «нормальная» серия. И всё что делается, это долбёжка чисел до упора в поисках любой еденички, что бы мифическое равно отбросить и перейти к чёткому ">". Соотвественно, банковское округление для транцендентых чисел бессмысленно, потому что равенства не бывает. Но использование для банкинга трансцендентных чисел, это вообще говоря, сомнительное занятие, тем более, что банковского округление это частный случай и не решает оно проблему полностью, а является лишь костылями подходящими лишь для определённых операций… По сути, проблема решается тем, что если и появляются где-то транцендентные числа, то для них надо использовать математическое округление и уже только дальше банковское/случайное, если есть какие-то статистически значимые операции.
А разве я говорил где-то, что смотрится чётная цифра ПОСЛЕ округляемой? Вы что-то не то говорите. Цифра, которая в моих примерах подчёркивается — это так называемый round-bit (назовём его b), по его значению проверяем направление округления, есть ещё sticky-bit (s), который является логическим или над всеми остальными битами, идущими после round-bit. Далее смотрится та самая цифра, которая остаётся последней после округления это least-bit (l), именно она должна стать чётной в случае, если b=1 и s=0. Таким образом, имеем три бита lrs информации по которым однозначно определяем округление. Но я намеренно не хотел как-то путать читателей этими тремя битами, а постарался объяснить на словах.
Для трансцендентных функций никогда не бывает s=0, вы правы, там всегда будет какая-то единичка дальше. Проблема лишь в том, что мы можем не знать b=0 или b=1. Алгоритм может плясать от одного к другому и обратно, пока не стабилизируется. Думаю, вы правильно меня поняли, но мне просто слово банковское непонятно, не знаком с такой терминологией.
Округление до ближайшего чётного (в английском языке известно под названием англ. banker's rounding — «округление банкира») — округление для этого случая происходит к ближайшему чётному числу, то есть 2,5 → 2; 3,5 → 4.
Спасибо, буду знать, это это банковское округление :) А бывает ещё nearest-ties-to-away, это когда величина посередине округляется дальше от нуля. То есть это то, что у нас в математике принято: 12.5→13; 13.5→14. Именно этот режим прописан в Стандарте для десятичных чисел с плавающей запятой, и я думал, что банковская система должна работать по этому варианту округления, раз она работает с десятичными числами.
Как это нет вопросов? :) Я же столько написал о том, что мы НЕ ЗНАЕМ, с какой стороны мы от середины :)
Если хотите, я вам сейчас объясню подробнее. Допустим, нам нужно округлить корректно всего до двух битов после запятой. Перед нами два числа:
1,01011111111111111111111111111111111111111 → 1,01
1,01100000000000000000000000000000000000001 → 1,10
Между ними разница 2 в степени -40. Оба эти числа были получены в ходе вычислений пошаговым алгоритмом (одно на 100-м шаге, другое — на 101-м), которое из них по-вашему нужно взять в качестве достаточно-точного для того, чтобы сформировать правильное округление?
Как видим, 40 битов точности недостаточно, чтобы понять, с какой стороны мы от серединки. А сколько достаточно? Не знаем! Вот в этом и состоит суть проблемы TMD.
Что-то я не понимаю. Какие такие гарантии точного знака? Не слышал. Метод Ньютона даёт, грубо говоря, удвоение знаков, но это жаргон такой, так просто принято говорить. В действительности он даёт возведение точности в квадрат по сравнению с предыдущим шагом. То есть, скажем, была точность 2 в минус 2, стала 2 в минус 4, потом 2 в минус 16.
Мне кажется, вы не понимаете суть. В моём примере чтобы получить точными ВСЕГО ДВА БИТА, неизвестно как долго ещё считать! Уже получили 41 бит после запятой, а всё ещё не знаем, какими точно будут ВСЕГО ДВА бита. А сколько ещё считать? Никто не знает! Это это незнание характерно для ВСЕХ алгоритмов всех трансцендентных функций, никаких гарантий никто никогда не давал.
Вот вы сами пишите: "соотвественно точность не достаточна".
Вопрос в том КОГДА она будет достаточной? Этого никто не знает.
0,12345X, округляя до 0,1235 мы внесём погрешность в 0,00005 максимум и разряд X отличный от нуля, может только уменьшить нашу погрешность. Для итогового же числа 0,1235, было бы уместно записать погрешность в +-0,00005, хотя в реальности она была меньше, просто это покажет, что число округлённое, и 4 знак после запятой не является точным и могло быть как 0,1234x так и 0,1235x.
Собственно говоря, с какой бы стороны мы не приближались к числу, математическое округление внесёт наименьшую погрешность величиной в половину округляемого разряда. В вашем примере, округлив до 2-39 мы полчим одинаковые числа. Понимаю, о чём вы говорите, о том, что сколь угодно малое число может превращать 1цу в набор девяток. Но если мы берём еденицу и +- 0,005 погрешность, то это будут числа от 0,995 до 1,004 и эти числа всегда будут округляться до 1цы, фаткически, пока погрешность не вырастет до 0,5 (а она будет только падать 0,9995 1,0004 и т.п) мы можем быть уверены в этой еденице, никакие младшие разряды 0.5 не переплюнут по величине. Поэтому при повышении точности старшие разряды, после округления в младшем разряде, не будут никак изменяться. Но при этом, если вы не уверены в фактических значениях разрядов вы не можете колдовать с 0.5 значением для десятичной, 0.4 для восмеричной, 0.3 для шестиричной, 0.2 для четырёхричной и 0.1 для двоичной, как в банковском округлении, потому что фаткически этот разряд у вас не настолько точный в любой позиции.
0,12345X, округляя до 0,1235 мы внесём погрешность в 0,00005 максимум и разряд X отличный от нуля, может только уменьшить нашу погрешность.
Ну да, если у нас есть число 0,12345... с гарантированно верными 5 цифрами — оно может быть округлено только до 0,1235. Вопрос в том, как получить эти самые 5 верных цифр; всегда может оказаться что на самом деле у нас 0,123449876... и округлять надо до 0,1234.
Да нету такой гарантии! Вот я взял число 0.0199999994**2, и попробовал найти квадратный корень. Первые 10 итераций получились такими:
0: 0.500199999988
1: 0.2504998400339936
2: 0.12604832367255547
3: 0.06461085479793548
4: 0.03540088231324618
5: 0.023350017521141124
6: 0.020240312312009873
7: 0.020001426015725795
8: 0.01999999945087718
9: 0.0199999994Третий разряд после запятой "стабилизировался" только на 9й итерации! Ну и где тут удвоение числа верных разрядов-то?
Квадратный корень, кстати, с точки зрения поставленной автором задачи вообще чуть ли не самая "вредная" функция: он не является трансцендентной функцией, а потому для него многоуровневая стратегия Зива может зациклиться на огромном множестве точек...
Откройте, пожалуйста, хотя бы википедию и посмотрите там раздел Proof of quadratic convergence for Newton's iterative method.
Квадратичная сходимость — это когда точность возводится в квадрат на каждой итерации. Про стабилизацию битов это не говорит. Про десятичные цифры тоже. Когда мы говорим про то, что количество битов удваивается, это не совсем правильные слова, а устоявшееся название, но все вроде понимают их смысл, кто в теме. Это как в математическом анализе, мы говорим "почти всюду", а имеем в виду "за исключением множества меры нуль". Вроде все понимают.
В вашем примере с корнем вы взяли очень далекое приближение, поэтому точность всё равно будет возрастать квадратично, но слишком долго приходится идти до первого правильного бита, то удвоение почти ничего не даст. Вопрос первого приближения очень важен, если взять не то число, алгоритм будет долго-долго тормозить сначала, что у вас и происходит.
Там уже на 4-5й итерациях число довольно близкое к правильному, но число верных разрядов, в отличии от точности, удваиваться всё равно не хочет.
Вот в таком виде хорошо видно как удвоение точности, так и то что происходит с "верными" разрядами:
4: 0.035400882313246180 ± 0.01540088291324
5: 0.023350017521141124 ± 0.00335001812114
6: 0.020240312312009873 ± 0.00024031291200
7: 0.020001426015725795 ± 0.00000142661572
8: 0.019999999450877180 ± 5.08771809404e-11Я дал вам источник, с которого можно начать разбираться. Дальше я рекомендую взять формулы, которые там участвуют, и проверить, что ваш пример в точности в них укладывается, либо убедиться, что вы сделали неправильный пример. Поскольку мы отклонились от темы статьи, я предлагаю на этом закончить.
Благодарю за дискуссию!
Один вопрос — зачем вы дали мне источник? Я и так знаю, что мои формулы в него укладываются, спорю-то я не с вами...
Потому что вы пытаетесь спорить с жаргонизмом. Все люди (из моего окружения), кто произносит фразу про удвоение правильных битов имеют в виду совершенно определённую вещь, то есть возведение погрешности в квадрат (при соблюдении ряда условий, и я об этом дважды предупредил). При этом все они понимают, что может не быть ни одной правильной цифры. Просто так говорят условно, и в большинстве случаев это и правда так. Это был мой последний ответ по данной теме в этой ветке комментариев.
Я спорю не с жарнонизмом, а с вот с этой цитатой:
Вы сами приводили как пример, вычисление квадратного корня, где количество ТОЧНЫХ знаков удваивается после каждой итерации. Точные — означает, что они не меняются при дальнейших вычислениях.
Вот и я привёл пример., где вроде как "точный" разряд ещё как меняется при дальнейших вычислениях.
Ох, так и быть, отвечу, нарушу своё обещание :)
Вы цитируете человека, который приписал мне то, чего я не говорил. Вот с ним, пожалуйста, и дискутируйте.
Хотя может это я не прав, не вижу, кому отвечаете. Тогда прошу прощения, не освоился, видимо :)
Да, всё, я увидел, что вы спорите не со мной, ещё раз прошу прощения :) У меня возникли основания полагать, что вы поверили, что я сказал глупость про точные числа в расчёте квадратного корня, хотя я просто не мог такое сказать.
Да, так пишут потому что это просто удобно так думать, это условно, при этом обычно подразумевается, что есть два обстоятельства:
- Бывают пограничные случаи, где ни одной правильной цифры не будет, под удвоением подразумевают возведение точности в квадрат.
- Бывают условия, при которых квадратичная сходимость скатывается до линейной. И тогда одна итерация, грубо говоря (условно), добавляет только одну цифру.
Нет, двумя округлениями проблема не решается, это лишь усугубляет проблему. И проблема возникает не только для банковского округления, а для любого. При любом способе округления есть определённые пограничные значения, в которых можно ошибиться.
Для примера рассмотрю простейшее округление — усечение. Вот мы считали-считали какую-то формулу, и получили в итоге число 1,00 с абсолютной погрешностью не более 0,005, сокращенно 1,00 ± 0,005. Как его усекать — до 0,99 или до 1,00? Хорошо если ответ "1,00 ± 0,005" нас устраивает, но что делать если у нас задача — обязательно получить два корректных разряда?
Как правило, разница в скорости исполнения между 1ulp функцией и корректно округлённой функцией — 2-4 раза. А ведь тут мы говорим, про [почти] один распоследний бит ответа, который никто не заметит.
Для реализации надо иметь только округление вверх и округление вниз.
Будет медленно, но математически точно (диапазон в котором лежит ответ).
Есть принятый стандарт: IEEE 1788-2015.
Отмечу, что для интервальной арифметики и анализа Уильямом Кэхэн так же сделал очень много.
Для интервалов гарантии одинаковости нет, но есть гарантия того, что на любой платформе точный результат f(x)=0 будет внутри вычисленного интервала. Кроме того, пересечение интервалов, полученных от всех платформ, так же будет содержать точный результат. Но как обычно у всего есть цена.
Ну, для начала, интересно, на какой платформе вы делали описанные тесты. Версия компилятора может ничего не значить, особенно в Unix мире.
Вот на Ubuntu 18.04/amd64 мне одинаково GCC и Clang выдали два раза 1.00206053 (хотя я мог не проверять на разных — очевидно, что эти функции сидят в glibc — хотя и выставлены отдельно в libm). А вот FreeBSD выдала такие же значения, как у вас в тесте — с разницей в 1 ULP, причём одинаково на 32 и 64 битах. Но что FreeBSDʼшная libm достаточно часто даёт ошибку в 1 ULP, а glibc — нет, известно, наверно, каждому, кто этим занимается серьёзно.
Теперь к этому:
Причина ошибок не всегда кроется в лени или низкой квалификации разработчиков, а в одном фундаментальном обстоятельстве, преодолеть которое современная наука пока не смогла.
А суть проблемы TMD в том, что мы заранее не знаем, до скольких битов рассчитать величину z, чтобы получить правильными столько битов после округления, сколько нам требуется.
Вообще-то известно, но решается чуть не так: для всех IEEEʼшных методов округления, и даже для расширений (типа 8 методов в SystemZ) это максимум 3 дополнительных бита (у них даже свои названия есть: guard, round и sticky — вы их знаете, но оставим для остальных) плюс факт ненулевого значения дальше этих битов (точное значение уже не нужно, нужно просто знать, получилось там твёрдое 0 или нет). Всё. Причём сам факт ненулевого хвоста превращается в установку sticky bit в 1, после чего эти три бита вместе с последним значащим подаются в логику округления.
В отдельных реализациях могут делать больше 3 бит, но умеренно — посмотрите например в Berkeley softfloat library — для float32 это 6 бит, но потому, что вообще у них работа идёт с 30-битовыми значениями в 32-битном поле, и удобно делать именно так.
Да, тут решение достигается за счёт того факта "хвост не ноль", который при формальном точном следовании подходу "сколько бит?" мог бы требовать неограниченной точности. Но в том и соль, что это ограничение такими формальными рамками совершенно искусственное. И для большинства этих функций для большинства аргументов, более того, не надо сложных вычислений для "хвост не ноль": есть очень мало аргументов, когда он может стать нулём (типа случай sin 30° == 0.5, тут можно эти 30° в радианах вынести в отдельную проверку), в остальных случаях можно "не ноль" сразу постулировать.
А разрешение ошибки в 1 ULP для трансцедентных функций вызвано, насколько мне известно, не проблемой округления, а тем, что решили, что делать тут полную точность при том, что промежуточные значения могут иметь ту же точность, оказывается слишком дорого. Напомню, что IEEE754-1985 писался по аппаратным реализациям Intel и Motorola, где вся логика вшивалась в сам процессор. Это уже потом поняли, что лучше в процессоре оставить только базовую арифметику и служебные операции, а трансцедентные вынести в софт (как оно сейчас и работает с SSE, Neon и всеми ровесниками и более поздними реализациями). Рядом говорят, что этот один ULP может требовать в 2-4 раза больше времени. Да, возможно. Но сама реализация в софте возможна за подъёмные деньги, в отличие от железа.
UPD: прочёл ваши замечания в другом комментарии про "пограничные" случаи, где микроошибка может перевернуть длинную цепочку бит. Но таких случаев от всех, где реализации дают ошибку в 1 ULP, дай бог чтобы 1 из 1000. Остальные — это просто экономия на последних шагах, которые не просто дороги, а крайне сложны в реализации. Если посмотреть на софтовые реализации трансцедентных функций — видно, какие тонкие игры там применяются, чтобы не потерять этот самый один бит, оставаясь в поле той же точности. Вот это и есть, мне кажется, основной источник проблемы. Если бы можно было легко нарастить точность промежуточных значений при расчёте какого-нибудь log10 рядами — мы бы об этой проблеме ничего не слышали.
(Собственно, 80-битный формат с 64-битной мантиссой FPU Intel и Motorola — как раз попытка решения этой проблемы, насколько известно из исторических документов.)
причём разные компиляторы ошибаются в разных местах, но ни один не будет работать правильно.
Проверьте на линуксе со свежим glibc. Уже интересно, что там найдёте. Они достаточно точно реализовали операции, насколько мне известно, причём достаточно давно. Ну а проверить можно через mpfr, у неё гарантированно своя реализация с управляемо неограниченной точностью.
Во-первых, задача носит фундаментальный характер.
… и давно таки решена? — посмотрите на тот же MPFR. Ну или покажите, что она тоже ошибается в 1 ULP… что-то я такого про неё не слышал.
Замечание по форме: есть такой формат у printf — %a. Лучше пользоваться им, чем переводить состояние в памяти в hex-форму и печатать из неё; и лучше задавать константы в таком контексте тоже лучше в виде 0x1.921fb54442d18p+1, чем в длинной десятичной константе.
Благодарю за совет по поводу %a, но в остальном я не могу прокомментировать ваши замечания, потому что из них очевидно, что вы не поняли содержание статьи. Возможно, вы захотите мне возразить, но не стоит, а лучше напишите сразу тем авторам из списка источников, которые более 15 лет уже пытаются решить эту задачу и всё ещё далеки от решения. Возможно, они удивятся, что оказывается всё уже давно решено.
Другая часть ваших замечаний в точности повторяют то, что я и так написал в статьи (например, про то, что мало всяких критических пограничных случаев и что для программиста важнее точность, а не правильность последних битов), поэтому неясно с кем вы спорите. Ясно что не со мной. Даже про числа с бесконечной точностью я говорил (но не называл конкретного пакета). В общем, давайте уважать друг друга. Я внимательно писал статью, а ваша задача — внимательно её прочитать, если намереваетесь дискутировать.
но в остальном я не могу прокомментировать ваши замечания, потому что из них очевидно, что вы не поняли содержание статьи.
Как раз из вашего комментария очевидно, что я его таки понял, но не считаю настолько важной акцентуацию именно на тех пограничных случаях с "нечёткими битами" (буду пользоваться таким описательным выражением), о которых вы говорили (причём только в комментариях подчеркнули это явно — до того было непонятно некоторым другим комментаторам).
Вы почему-то сочли, что основная или даже единственная причина (тут уточните сами) того самого разрешения погрешности в 1 ULP — это нечёткие биты на пограничных значениях. Так вот — это не просто сомнительно, это скорее всего неправда, и я встречал объяснение (в документах по истории IEEE754), почему это неправда.
Я внимательно писал статью, а ваша задача — внимательно её прочитать, если намереваетесь дискутировать.
Мой комментарий действительно построен из изначального расчёта на то, что вы просто не учли разницу между качеством реализаций (которое в некоторых намеренно сокращено для экономии вычислительных затрат, раз уж стандарт позволяет) и с учётом того, что вы явно неверно описали цели стандарта. Но, если вы внимательно прочтёте его, увидите, что я внёс поправку на проблемы пограничных значений. Если вы намерены осмысленно дискутировать, то уже ваша задача внимательно прочесть замечания оппонента. И если бы вас эта тема всерьёз интересовала на предмет "возможно ли это?", вы бы попросили уточнить источники и тесты моего утверждения про точность glibc, а не начали обвинять в нечтении.
Если хотите тут действительно осмысленного рассмотрения и статьи, которая имела бы реально значение повыше простого копирования известных истин — я считаю, что вам следует:
Опереться на opensource — его тут достаточно (в смысле количества вариантов), а не на неизвестные "компиляторы" мира Windows, когда вы даже не определили, от какого компонента зависит реализация, и как связаны версионности компилятора и этого компонента (hint: связь неоднозначна, а версия компилятора недостаточна).
Провести сравнение типовых реализаций: MSun исходная; MSun осовремененная (в рамках FreeBSD, например); IBM Accurate Mathematical Library исходная; она же доработанная; MPFR. Насколько есть данные, и про x87 FPU (и про его точность — там есть несколько парадоксальных моментов).
Взять хотя бы пару функций для примера и проверить утверждение про бо́льшую (или полную) точность glibc. Опровержение меня тут (на какой-то функции, которую я не проверял) было бы на порядки полезнее всего это "вы меня не читаете" ;)
Объяснить, с помощью каких же приёмов код той же самой IBM AML ухитряется обеспечить точность до последней цифры (всегда или почти всегда) без расширения точности используемых чисел.
Вот это было бы действительно полезно и потянуло бы на полноценную статью (возможно, в исходно более научном издании;)) Сейчас же у вас только описательная затравка, за счёт неточного описания сбивающая с толку.
Table-Maker's Dilemma, или почему почти все трансцендентные элементарные функции округляются неправильно