Комментарии 39
По моему мнению не вполне удачно сочетать в себе две совершенно разные темы. Первая — базовые принципы арифметики с плавающей запятой вроде закона дистрибутивности (почему он нарушается, также как закон ассоциативности). Вторая — сложные элементы из неочевидной и нетривиальной темы, касающейся вычисления функций. Целевая аудитория, которой нужно первое, не будет вникать во второе, потому что не сможет, а вторым совершенно не нужно объяснять то, что для них уже очевидно. Поэтому для первых я рекомендую сначала пройти единственный из известных мне видео-курсов (на русском языке), где всё подробно с самого начала. Прямо вообще с нуля. Будет время, может быть сделаю и текстовый курс. Ну а вторым рекомендую помогать автору в его нелёгком начинании.
необходимо вычисление синуса (s) и косинуса одновременно, что требует в два раза больше вычислений ряда ТейлораSSE решает эту проблему, позволяя одновременно вычислять 2 ряда Тейлора.
Не совсем. Во-первых, тут либо ты находишь 2 синуса одновременно, либо синус и косинус. Забегу вперед, код быстрого, но неточного синуса svml_d_sin4_core_avx2.S. Во-вторых, присмотритесь к коду вычисления обоих функций. Он существенно разный. Значит полной парализации не получится.
Не совсемНе вижу никаких проблем, по схеме Горнера распараллеливаются идеально, в конце только синус домножается на x и всё:
Это как раз способ, который лучше избегать. Посмотрите, если не сложно, параграф «Распределительное свойство умножения в мире чисел с плавающей точкой» и в прошлой статье сравнение sin_e4 и sin_e5.
Этот способ я привёл лишь в качестве наглядной демонстрации. У меня тоже есть опыт реализации быстрого вычисления синуса/косинуса — только не в формате double, а double-double, c 32 (примерно) значащими цифрами после запятой. Как раз то, о чём ArtemKaravaev и писал.
Это очень интересно. Есть где почитать, посмотреть?
Пока нет — это же большая ответственность, и соответственно большой риск — использовать математическую библиотеку, написанную неизвестно кем. К тому же для ускорения вычислений помимо SSE я использовал FPU — что не одобряется сообществом. О матчасти уже писал, использовал 4 таблицы по 64 значения.
Ну матчасть у вас абсолютно такая же, как в libm, которую я описал в этой статье. Но в любом случае, даже если библиотека не готова, интересно узнать, как вы использовали FPU для ускорения вычислений, какие у Вас получились точности, как вы поженили SSE и «double double». Сравнение со стандартной библиотекой и пр. Почему 4 таблицы по 64 значения. Как математически вы подобрали размер таблицы. Ну и, естественно, производительность.
Со стандартной библиотекой сравнивать нет смысла — там же значащих цифр в два раза меньше. По скорости сравнивал только базовые функции с библиотекой отсюда — сложение в 2 раза быстрее, умножение в 10. Зная о способности компиляторов c/c++ к экстремальной оптимизации решил основные операции написать на ассемблере, чтобы совершенно точно быть уверенным в порядке вычислений. О точности вычисления последнего бита не заморачивался вообще. Смысл использования FPU — в 80-битной точности промежуточных вычислений.
Размер таблиц подбирал интуитивно, перебирать все варианты с замером скорости и точности было лень. Две таблицы с синусами, две с косинусами, в итоге два комплексных умножения. Пробовал по 3 таблицы — точность уже падала.
код для сложения двух double-double
ddasm_add_dd_dd PROC uses esi edi dd1:PTR DWORD, dd2:PTR DWORD, dd3:PTR DWORD
local s1s:XMMWORD, s2s:XMMWORD
mov esi, dd1
mov edi, dd2
mov eax, dd3
movupd xmm0, XMMWORD PTR [esi];=a
movupd xmm1, XMMWORD PTR [edi];=b
movapd xmm3, xmm0;=sum
addpd xmm3, xmm1
movapd xmm4, xmm3;=bv
subpd xmm4, xmm0
movapd xmm5, xmm3;=av
subpd xmm5, xmm4
subpd xmm1,xmm4;=br
subpd xmm0,xmm5;=ar
addpd xmm0, xmm1;=err
movupd XMMWORD PTR s1s, xmm3
movupd XMMWORD PTR s2s, xmm0
fld REAL8 PTR s1s
fadd REAL8 PTR s1s+8
fadd REAL8 PTR s2s
fadd REAL8 PTR s2s+8
fstp REAL8 PTR [eax]
fld REAL8 PTR [eax]
fsubr REAL8 PTR s1s
fadd REAL8 PTR s1s+8
fadd REAL8 PTR s2s
fadd REAL8 PTR s2s+8
fstp REAL8 PTR [eax+8]
ret
ddasm_add_dd_dd ENDPРазмер таблиц подбирал интуитивно, перебирать все варианты с замером скорости и точности было лень. Две таблицы с синусами, две с косинусами, в итоге два комплексных умножения. Пробовал по 3 таблицы — точность уже падала.
Со стандартной я имел ввиду конечно, libquadmath (которая является частью стандартной библиотеки). s_sinl.c Там немного другой подход, чем у вас.
Вообще синус — это слишком просто и банально. Интереснее посчитать например логарифм, где разложение в ряд Тейлора совершенно никак не поможет.
Не хочу устраивать беспредметную дискуссию, но логарифм кажется проще. По крайней мере в glibc. После преобразований он просто считается «полиномом, очень похожим на ряд Тейлора.» (не хочу снова вступать в дискуссию о том, можно ли называть рядом Тейлора полином у которого коэффициенты различаются на несколько ULP).
Любое решение кажется простым, когда его можно посмотреть в готовом виде в коде. Вы посчитайте логарифм, не подсматривая в исходники glibc и опираясь лишь на школьные знания. Синус проще, потому что его достаточно на интервале 0-пи/2 посчитать, при этом он хорошо аппроксимируется рядом Тейлора, в отличие от многих других.
Я бы не согласился, что посчитать его на интервале от 0 до Pi/2 достаточно. Перед этим есть целая проблема — выполнить редукцию аргумента x до этого интервала; есть алгоритмы, которые призваны делать это с максимально возможной точностью. Просто выполнить операцию типа x ← x-2×Pi×k не получится, будут большие потери. Но мне кажется вы об этом знаете, только забыли упомянуть.
Если считать синус с точностью single используя double — то никакая это не проблема. Для double соответственно есть extended double с 80-битной точностью, придуманный как раз для этих целей. Или можно просто искусственно округлять результат до 13-и знаков. В любом случае, все эти проблемы — алгоритмического, а не математического характера. Во многих других функциях редукция аргумента невозможна в принципе, и считать их на (практически) бесконечном интервале — такое себе удовольствие.
Refridgerator
Какую функцию Вам кажется было бы посчитать сложнее всего. Напишите, я посмотрю, и может быть сделаю статью)))
Какую функцию Вам кажется было бы посчитать сложнее всего. Напишите, я посмотрю, и может быть сделаю статью)))
С логарифмом я нашёл хитрое решение, в которой первые 15 цифр считаются через FPU, а остальные 15 считаются через экспоненту используя свойства логарифма. А вот быстро вычислить (дополненную) функцию ошибок с 30-ю знаками оказалось посложнее. Ещё сложнее оказалось вычислить интегральные синус и косинус. Гамма-функцию в принципе тоже можно через редукцию аргумента посчитать — но там точность быстро падает из-за многократного умножения. А посчитать логарифм от гамма-функции — вот это уже интересно.
Не знаю как в крутых процессорах, но вот на моем VLIW DSP in-oder, запросто можно считать одновременно синус и косинус. Также можно еще что-нибудь посчитать (за те же такты), т.к. эти полиномы очень неудобны для конвейера — текущий результат нужен для следующей операции. А из-за этого пузыри на конвейере. В эти пузыри можно впихнуть что-то полезное. Своеобразный ручной гипертрейдинг :) Главное чтобы эти параллельные вычисления не нуждались в переходах, ну и вычислительная сложность была приблизительно одинаковой.
эти полиномы очень неудобны для конвейера — текущий результат нужен для следующей операцииТак полином можно же элементарно разбить на сумму из двух и считать их поочерёдно — никакой зависимости не будет.
Изменение порядка вычислений на точности никак не скажется? или отличие несущественное?
В своих экспериментах я не увидел качественного увеличения точности при перетасовке аргументов при сложении. Характер искажений — меняется, но незначительно, и чуть большая точность на одних входных данных компенсируется чуть меньшей на других. Вот в этой статье я считал синус схемой Горнера и получил вполне неплохие результаты — при том, что решал задачу «в лоб» и оптимальные коэффициенты никак не подбирал.
Я посмотрел внимательно Вашу статью. Теперь я понял, откуда берутся Ваши заблуждения о простоте синуса, независимости точности при перетасовке аргументов и пр. Я бы не хотел никого обидеть, но ваши эксперименты сделаны очень поверхностно, к сожалению, я вынужден сказать, что безграмотно. На примере вашей же статьи я поясню своё утверждение.
Мне не удалось в точности воспроизвести ваши условия (VC++). Много лет уже не пользовался Windows.))) Компилировал GNU и clang для x86 и ARM в режиме -O0 и -O2. Результаты идентичные между системами и Вашими результатами в таблице в статье.
Что же плохо в коде Вашей программы:
1) Использование чисел вида 1.0 и 0.5 подобных сильно уменьшает ошибку вычисления. Добавил немного «сложных» чисел, которые заполняют всю мантиссу.
2) Делать «позитивные» выводы по 5-7 точкам просто неправильно. В моих статьях я использую 16 млн точек и считаю, что надо больше.
Чтобы показать, что ваш подход даёт плохую точность надо немного изменить вашу программу program.cpp. А именно поменять входные значения для расчётов и сделать вывод немного более информативным:
И чтобы не быть голословным — сразу привожу результаты выполнения программы.
Посмотрите, например на 3-ий x = 0.047619… Результаты отличаются на 2 ULP даже при таком реально малом значении x при котором классический ряд Тейлора сходится на 9 степени.
Посмотрите на x=0.99999914192480920949. Погрешность в 4 ULP!!! При значении близком к единице. Реально такой подход нуждается в глубокой переработке. Ваш расчёт за такое же количество операций умножения (а ещё плюс деление), как и ряд Тейлора даёт гораздо худшую точность.
Чтобы не было сомнения, вы можете проверить результаты в Mathematica. Например, sin(0.99999914192480920949). Видно, что это значение гораздо ближе к «внутреннему» синусу, чем к Вашему.
Уважаемый Refridgerator, я призываю Вас перепроверить Ваши результаты, например, используя framework моей программы. Я так же думаю, что было бы очень честно, если прежде чем писать комментарии Вы бы разобрались бы так же глубоко в моих статьях, как я в Вашей. А иначе Ваши комментарии просто сбивают людей с толку и направляют на ложный путь.
Мне не удалось в точности воспроизвести ваши условия (VC++). Много лет уже не пользовался Windows.))) Компилировал GNU и clang для x86 и ARM в режиме -O0 и -O2. Результаты идентичные между системами и Вашими результатами в таблице в статье.
Вывод на экран исходной программы
**********************
* rational poly test *
**********************
x = 0.5
-------------------------
sin(x): 0.479425538604203
double: 0.479425538604203
x = 1
-------------------------
sin(x): 0.8414709848078965
double: 0.8414709848078965
x = 2
-------------------------
sin(x): 0.9092974268256817
double: 0.9092974268256816
x = 3
-------------------------
sin(x): 0.1411200080598672
double: 0.1411200080585958
x = 4
-------------------------
sin(x): -0.7568024953079282
double: -0.7568024960833886
x = 5
-------------------------
sin(x): -0.9589242746631385
double: -0.9589243758030122
Что же плохо в коде Вашей программы:
1) Использование чисел вида 1.0 и 0.5 подобных сильно уменьшает ошибку вычисления. Добавил немного «сложных» чисел, которые заполняют всю мантиссу.
2) Делать «позитивные» выводы по 5-7 точкам просто неправильно. В моих статьях я использую 16 млн точек и считаю, что надо больше.
Чтобы показать, что ваш подход даёт плохую точность надо немного изменить вашу программу program.cpp. А именно поменять входные значения для расчётов и сделать вывод немного более информативным:
typedef union {std::uint64_t i; double x;} num_t; // Compare binary represenation.
...
cout << std::setprecision(19) << std::defaultfloat;
for (double x : {1E-4/3.0, 1.0/163.0, 1.0/21.0, 2.0/ 3.0,
0.99999914192480920949, 1. + 1.0/3.0})
{
cout << " x = " << x << endl;
cout << " -------------------------" << endl;
num_t sn, si;
sn.x = sin(x);
si.x = padesin<double>(x);
cout << "sin(x), hex, dx: " << sn.x << ", " << std::hex
<< setfill('0') << setw(19) << right
<< sn.i << ", " << abs(sinl(x) - sn.x) << endl;
cout << "double, hex, dx: " << si.x << ", " << std::hex
<< setfill('0') << setw(19) << right
<< si.i << ", " << abs(sinl(x) - si.x) << endl;
cout << endl;
}
...
И чтобы не быть голословным — сразу привожу результаты выполнения программы.
**********************
* rational poly test *
**********************
x = 3.3333333333333334931e-05
-------------------------
sin(x), hex, dx: 3.3333333327160497624e-05, 0003f0179ec9caf9bb8, 2.2003004293910536152e-21
double, hex, dx: 3.3333333327160490848e-05, 0003f0179ec9caf9bb7, 4.5759631486433490974e-21
x = 0.0061349693251533743768
-------------------------
sin(x), hex, dx: 0.0061349308407180058733, 0003f7920f0f52f8960, 1.4568966692773965832e-19
double, hex, dx: 0.0061349308407180067407, 0003f7920f0f52f8961, 7.2167207106066388889e-19
x = 0.047619047619047616404
-------------------------
sin(x), hex, dx: 0.047601053042734112197, 0003fa85f2a4c25bc92, 2.215838190017249687e-18
double, hex, dx: 0.047601053042734098319, 0003fa85f2a4c25bc90, 1.1661949617797207068e-17
x = 0.66666666666666662966
-------------------------
sin(x), hex, dx: 0.61836980306973698962, 0003fe3c9af78209765, 1.0950441942103594783e-17
double, hex, dx: 0.6183698030697368786, 0003fe3c9af78209764, 1.0007186052041205926e-16
x = 0.99999914192480920949
-------------------------
sin(x), hex, dx: 0.84147052118758247641, 0003feaed5396218f4b, 5.3342746886286818153e-17
double, hex, dx: 0.84147052118758280947, 0003feaed5396218f4e, 2.7972416050126014397e-16
x = 1.3333333333333332593
-------------------------
sin(x), hex, dx: 0.97193790136331270624, 0003fef1a1d8383254a, 5.2638015474171240271e-17
double, hex, dx: 0.97193790136331281726, 0003fef1a1d8383254b, 5.8384286988344413771e-17
3.3306690738754696213e-16 0.99999914192480920949
Посмотрите, например на 3-ий x = 0.047619… Результаты отличаются на 2 ULP даже при таком реально малом значении x при котором классический ряд Тейлора сходится на 9 степени.
Посмотрите на x=0.99999914192480920949. Погрешность в 4 ULP!!! При значении близком к единице. Реально такой подход нуждается в глубокой переработке. Ваш расчёт за такое же количество операций умножения (а ещё плюс деление), как и ряд Тейлора даёт гораздо худшую точность.
Чтобы не было сомнения, вы можете проверить результаты в Mathematica. Например, sin(0.99999914192480920949). Видно, что это значение гораздо ближе к «внутреннему» синусу, чем к Вашему.
Уважаемый Refridgerator, я призываю Вас перепроверить Ваши результаты, например, используя framework моей программы. Я так же думаю, что было бы очень честно, если прежде чем писать комментарии Вы бы разобрались бы так же глубоко в моих статьях, как я в Вашей. А иначе Ваши комментарии просто сбивают людей с толку и направляют на ложный путь.
Если бы вы действительно внимательно прочитали мою статью, то обратили бы внимание, что там речь идёт о вычислении полинома, а не о вычислении синуса, и рассматривается не абсолютная, а накопленная погрешность вычислений, полученная при использовании различных типов данных. Вы же обратили внимание, что там есть не только double, но и single? Кроме того, раз он получен через аппроксимацию Паде, то, как и в случае с рядом Тейлора, погрешность в сравнении с реальным синусом растёт вместе с аргументом, что также наглядно продемонстрировано, а внимание на этом не заостряется очевидно потому, что статья вообще не об этом.
было бы очень честно, если прежде чем писать комментарии Вы бы разобрались бы так же глубоко в моих статьях, как я в ВашейЕсли вы помните, то мои комментарии были посвящены вычислению синуса с 30-ю знаками — то есть в 2 раза больше, чем у Вас. Поэтому, если говорить о честности, то и сравнивать нужно на равных — адаптируйте сначала свой алгоритм до формата double-double, а затем можно и посравнивать.
В своих экспериментах я не увидел качественного увеличения точности при перетасовке аргументов при сложении. Характер искажений — меняется, но незначительно, и чуть большая точность на одних входных данных компенсируется чуть меньшей на других.
Данное утверждение идёт в разрез не только с моими данными, но и данными этой статьи. Поэтому оно требует большего обоснования с Вашей стороны.
Вот в этой статье я считал синус схемой Горнера и получил вполне неплохие результаты — при том, что решал задачу «в лоб» и оптимальные коэффициенты никак не подбирал.
В чём «неплохость» Вашего результата? Точность, как видно нет, скорость тоже нет.
Если бы вы действительно внимательно прочитали мою статью, то обратили бы внимание, что там речь идёт о вычислении полинома, а не о вычислении синуса, и рассматривается не абсолютная, а накопленная погрешность вычислений, полученная при использовании различных типов данных.
Цитирую Вашу статью:
Используя Wolfram Mathematica и команду PadeApproximant[Sin[x], {x, 0, {11, 11}}]
получим вот такой вот рациональный полином для аппроксимации синуса, обеспечивающий точность double в диапазоне примерно от -2 до 2:
Т.е. сейчас вы говорите, что там просто какой-то полином. А зачем тогда вы его сравниваете с синусом? Что тогда с выводом статьи? Как абстрактный полином является реальным вычислением? Цитирую Вашу статью:
Преимущество Posit, демонстрируемое на единичных вычислениях — не более чем фокус, расплатой за которое является катастрофическое снижение точности в «тяжёлых» реальных вычислениях.
Я очень внимательно прочитал Вашу статью. Вывод из неё один. Вы нашли полином, который плохо аппроксимирует синус в double и очень плохо в posit. Никаких реальных вычислений, о которых Вы говорите в выводе статьи в самой статье не показано. Реальные вычисления приведены ниже.
Если вы помните, то мои комментарии были посвящены вычислению синуса с 30-ю знаками — то есть в 2 раза больше, чем у Вас.
Если я правильно помню, то Вы не хотите поделится с нами вашими результатами? Я не собираюсь считать в double-double. У меня нет такой цели. Всё, что я увидел из Вашей статьи, что вы не умете грамотно считать ни в double ни в posit. Добавьте в свой код вот эту функцию и вы увидите, что на некоторых значениях (например на 0.5) posit на 1 десятичный знак(!) точнее.
template< typename T >
T sin_smp(T x) {
T xx = x * x;
T res = 1.0;
for(int i = 31; i >= 3; i -= 2) {
res = T(1.0) - res * xx / T(float(i * (i - 1)));
}
return x * res;
}
Данное утверждение идёт в разрез не только с моими данными, но и данными этой статьи. Поэтому оно требует большего обоснования с Вашей стороны.Я же написал — в своих экспериментах — и никакого обобщения на все случаи жизни не делал. А в своих экспериментах я оценивал накопленную погрешность, в частности, через двойное преобразование Фурье. Вполне возможно, что эти эксперименты были проведены некорректно.
В чём «неплохость» Вашего результата? Точность, как видно нет, скорость тоже нет.Что значит «нет точности» и «нет скорости»? 0 правильных знаков и бесконечное время выполнения? А «неплохость» не моего результата — а аппрксимации Паде. Он вполне достаточен для того, чтобы нарисовать график синуса на экране.
Т.е. сейчас вы говорите, что там просто какой-то полином. А зачем тогда вы его сравниваете с синусом? Что тогда с выводом статьи? Как абстрактный полином является реальным вычислением?Потому что полином может аппроксимировать не только аналитические функции, но и произвольный набор экспериментальных данных — и это не делает его «абстрактным». Я его сравниваю с синусом потому так интереснее, чем ни к чему не привязанный случайный набор коэффициентов. Я использовал Паде, а не Тейлора, чтобы «утяжелить» вычисления и максимизировать погрешности — поскольку при тестировании более важен наихудший, а не наилучший результат.
Если я правильно помню, то Вы не хотите поделится с нами вашими результатами?Нет, я говорил что это ответственность, которую я брать на себя пока не готов. Но если Вы тоже поработаете в этом направлении — у меня будет больше мотивации довести библиотеку до ума и выложить в свободный доступ. Ну или наоборот — станет понятно, что она никуда не годится.
Всё, что я увидел из Вашей статьи, что вы не умете грамотно считать ни в double ни в posit.Жаль, что вы не увидели заголовок и не прочитали введение — цели статьи там вполне конкретно обозначены и «грамотные» вычисления в них не входят.
Спасибо. Я вот за утро поделил свой полином на две части (чет и нечет) и немного ускорил вычисления. Полином у меня не очень длинный (до x**15), да еще с нечетным количеством членов. Но все равно на самом полиноме я выиграл 25% в скорости. Все ваши слова о точности пока подтверждаются и в моих вычислениях. Но пока это было несколько тысяч когда-то заготовленных значений. Буду сейчас внимательно читать, что написал dorex и попробую сравнить результаты на «плохих» аргументах. :)
Еще было бы интересно услышать про x87 FSIN/FCOS — как они реализованы у Intel/AMD и почему deprecat'нулись в эпоху SSE.
Спасибо за вопрос. Под прошлой статьей обсудили подробно здесь. Если у Вас будут ещё вопросы, лучше, наверное, продолжить обсуждение в том потоке.
При написании библиотеки для своего VLIW DSP процессора первая проблема с синусом была в том, как точнее привести аргумент к нужному диапазону. Иначе дальше уже бессмысленно что-то точно считать. А еще при реализации функции одновременного вычисления синуса и косинуса (получались те же такты что и у одиночной функции) у меня почему-то синус получался точнее чем косинус. Надеюсь, что благодаря последним статьям о точных вычислениях с double мне удастся улучшить свои функции :)
Не понимаю как удалить комментарий :)
Запустил s_sin.c на своем процессоре.
Есть немного удивительные для меня результаты.
У меня есть таблица эталонов из 1024 элементов, которые я когда-то подготовил в MVS
для интервала от 0 до 4.69398121 с некоторым шагом.
Функция s_sin.c дала 100% совпадение с эталоном.
При этом моя функция на половине значений дает ошибку в младшем бите.
И есть несколько значений где разница мантисс равна 2, а не 1.
И это при том, что длина полинома у меня чуть больше чем у s_sin.c :)
Затем решил я проверить одно число, с которым у меня были проблемы. Входной аргумент
val = 0х40F7A79D c28f5b10 довольно большое число (96889.859999995911) и требуется
серьезная редукция. Результаты:
в MVS = 0xBF4FCE42 611d2228
моя функция = 0хBF4FCE42 611d1da9
функция s_sin.c = 0хBF4FCE42 611d1da9
Результат вольфрама 0хBF4FCE42 60D79B24
Почему у меня плохой результат я понимаю, но почему такой же и у s_sin.c :)
Может кто перепроверит результат s_sin.c на данном значении.
Есть немного удивительные для меня результаты.
У меня есть таблица эталонов из 1024 элементов, которые я когда-то подготовил в MVS
для интервала от 0 до 4.69398121 с некоторым шагом.
Функция s_sin.c дала 100% совпадение с эталоном.
При этом моя функция на половине значений дает ошибку в младшем бите.
И есть несколько значений где разница мантисс равна 2, а не 1.
И это при том, что длина полинома у меня чуть больше чем у s_sin.c :)
Затем решил я проверить одно число, с которым у меня были проблемы. Входной аргумент
val = 0х40F7A79D c28f5b10 довольно большое число (96889.859999995911) и требуется
серьезная редукция. Результаты:
в MVS = 0xBF4FCE42 611d2228
моя функция = 0хBF4FCE42 611d1da9
функция s_sin.c = 0хBF4FCE42 611d1da9
Результат вольфрама 0хBF4FCE42 60D79B24
Почему у меня плохой результат я понимаю, но почему такой же и у s_sin.c :)
Может кто перепроверит результат s_sin.c на данном значении.
Мой результат WolframAlpha в точности совпадает с s_sin
sin(96889.85999999591149389743804931640625)
sin(96889.85999999591149389743804931640625)
Спасибо. Ваш ответ я понимаю как то, что результат s_sin для
sin(0х40F7A79D c28f5b10) есть 0хBF4FCE42 60D79B24. К сожалению, не знаю как в вольфраме посмотреть результат в hex формате :)
Буду смотреть в чем у меня причина иного поведения.
Я провел сравнение результатов своей функции и функции s_sin на интервале от 0 до
пи\2 (чтобы не задействовать случай с редукцией) с шагом равным делению интервала на количество итераций.
Для миллиона итераций отличие не больше 2-х младших бит мантиссы. Для 100 млн итераций плохими оказались 3 бита. При этом одинаково ведут себя как функция со строгим порядком вычисления полинома, так и в случае деления полинома на две части.
Я при вычислении использую таблицу из 5-ти значений, т.е. интервал [0-пи\2] у меня разбит на 5 секторов. Отличия от s_sin появляются когда входное значение попадает в тот случай где s_sin использует таблицу и более точное вычисление.
Сам синус я вычисляю через sin(a+b) = sin(a)*cos(b) + cos(a)*sin(b)
a есть значение сектора, т.е. sin(a) cos(a) константы из таблицы
sin(0х40F7A79D c28f5b10) есть 0хBF4FCE42 60D79B24. К сожалению, не знаю как в вольфраме посмотреть результат в hex формате :)
Буду смотреть в чем у меня причина иного поведения.
Я провел сравнение результатов своей функции и функции s_sin на интервале от 0 до
пи\2 (чтобы не задействовать случай с редукцией) с шагом равным делению интервала на количество итераций.
Для миллиона итераций отличие не больше 2-х младших бит мантиссы. Для 100 млн итераций плохими оказались 3 бита. При этом одинаково ведут себя как функция со строгим порядком вычисления полинома, так и в случае деления полинома на две части.
Я при вычислении использую таблицу из 5-ти значений, т.е. интервал [0-пи\2] у меня разбит на 5 секторов. Отличия от s_sin появляются когда входное значение попадает в тот случай где s_sin использует таблицу и более точное вычисление.
Сам синус я вычисляю через sin(a+b) = sin(a)*cos(b) + cos(a)*sin(b)
a есть значение сектора, т.е. sin(a) cos(a) константы из таблицы
Просмотрел код. Оказалось, что константы полинома я взял свои, а в s_sin они почему-то другие. Подставил правильные и ничего не поменялось :) сложно это — точные вычисления :)
О правильном порядке вычислений :)
Переписал s_sin на ассемблер, чтобы оценить скорость вычисления. При переборе синуса на 100 млн чисел с некоторым шагом, на 2-х числах ассемблерная реализация не совпала с Си-шной (отличался младший бит). Оказалось, что в ассемблере выражение x * xx * (sn3 + xx * sn5) я сделал так: xx * sn5, +sn3, * xx, *x :)
Переписал s_sin на ассемблер, чтобы оценить скорость вычисления. При переборе синуса на 100 млн чисел с некоторым шагом, на 2-х числах ассемблерная реализация не совпала с Си-шной (отличался младший бит). Оказалось, что в ассемблере выражение x * xx * (sn3 + xx * sn5) я сделал так: xx * sn5, +sn3, * xx, *x :)
Зарегистрируйтесь на Хабре , чтобы оставить комментарий
Точные и быстрые вычисления для чисел с плавающей точкой на примере функции синуса. Часть 2: libm