Спасибо. Именно об этом и был мой вопрос. Попробую посмотреть, может чего и пойму :)
А то жалко иметь хороший аппаратный double, а расчеты производить на счетах, потому что точнее выходит :)
1. я как раз про свой конкретный случай :) в моем процессоре можно 2 double умножения за такт, но ни одного int64*int64. Эксперимент с sin_e7 примером я провел. Реализовал команду умножения int64*int64 через последовательность умножений int32*int32 и сложений.
2. насколько я помню, FMA избегает ненужного округления перед сложением. Результат будет немного точнее. Но у меня нет и FMA :) Только умножить и сложить.
3. Да, я использовал такой подход в синусе. Тот еще геморрой :)
«Всё-таки long double достаточно экзотический вид данных и его поддержка в современных процессорах не является приоритетной».
А если в процессоре нет целочисленного умножения 64*64? :) Можно конечно через 32*32,
но это в 10 раз медленнее чем double*double.
Правильно ли я понимаю, что если в сложении-вычитании double можно производить вычисления большей точности за счет double.double, то с умножением такой фокус не проходит? Умножили и хвост безнадежно потерян. Кстати, в процессоре можно легко сделать вычисления double+double c сохранением обоих частей результата. Аналогично можно было бы сделать и с умножением. т.е. по аналогии с целыми числами int*int = long long.
Был бы какой-то толк от подобных инструкций?
Просмотрел код. Оказалось, что константы полинома я взял свои, а в s_sin они почему-то другие. Подставил правильные и ничего не поменялось :) сложно это — точные вычисления :)
Спасибо. Ваш ответ я понимаю как то, что результат 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) константы из таблицы
Запустил 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 на данном значении.
Спасибо. Я вот за утро поделил свой полином на две части (чет и нечет) и немного ускорил вычисления. Полином у меня не очень длинный (до x**15), да еще с нечетным количеством членов. Но все равно на самом полиноме я выиграл 25% в скорости. Все ваши слова о точности пока подтверждаются и в моих вычислениях. Но пока это было несколько тысяч когда-то заготовленных значений. Буду сейчас внимательно читать, что написал dorex и попробую сравнить результаты на «плохих» аргументах. :)
Не знаю как в крутых процессорах, но вот на моем VLIW DSP in-oder, запросто можно считать одновременно синус и косинус. Также можно еще что-нибудь посчитать (за те же такты), т.к. эти полиномы очень неудобны для конвейера — текущий результат нужен для следующей операции. А из-за этого пузыри на конвейере. В эти пузыри можно впихнуть что-то полезное. Своеобразный ручной гипертрейдинг :) Главное чтобы эти параллельные вычисления не нуждались в переходах, ну и вычислительная сложность была приблизительно одинаковой.
При написании библиотеки для своего VLIW DSP процессора первая проблема с синусом была в том, как точнее привести аргумент к нужному диапазону. Иначе дальше уже бессмысленно что-то точно считать. А еще при реализации функции одновременного вычисления синуса и косинуса (получались те же такты что и у одиночной функции) у меня почему-то синус получался точнее чем косинус. Надеюсь, что благодаря последним статьям о точных вычислениях с double мне удастся улучшить свои функции :)
Вы говорите про общий случай. Но для своего теста мне был удобнее именно такой подход. Спасибо за ссылку о сортировке, а то я еще не дочитал, но уже пишу :)
Этот алгоритм проще было бы показать на числах float, а точный результат посчитать с помощью double. Еще суммирование таких чисел зависит от порядка их следования. Так что если не отсортировать, то точность не будет идеальной. После первой статьи я написал пару тестов на своем процессоре и случайно обнаружил, что «старое» суммирование, но с SIMD подходом, дало лучший результат чем алгоритм Кахана. Просто порядок следования чисел был для SIMD удачнее :)
Преимущество R-V над PIC16F18326, в том, что сейчас это модно и привлекательно. Также в статье было указано, что на базе R-V будет целая линейка процессорных ядер, а это уже какой-то системный подход к проектированию. Это и можно назвать прорывом :) Ну а для конкурентности нужно еще поработать.
Я попробовал этот алгоритм на сумме float массива из 1К элементов и действительно этот алгоритм дал в 10 раз меньшую ошибку по сравнению с ошибкой обычного суммирования. За эталон брал сумму в double формате. Возможно удачные попались массивы чисел :)
Я не типичный юзер :) у меня свой VLIW in-order DSP с SIMD возможностями. Никаких ОС. Ряд библиотечных функций приходится писать самому. При реализации различных стандартных функций есть возможность за те же такты вычислить точнее. Также порой нужен какой-то более точный (пусть и медленный) алгоритм, чтобы оценить погрешность более скоростного вычисления. Отсюда и интерес к подобным статьям :)
Интересная для меня тема. Недавно пришлось писать DSP библиотеку для процессора, так совсем замучился с float/double числами. Как только используешь SIMD вычисления, сразу результат отличается от результата не SIMD алгоритма. И порой непонятно то-ли баг в программе, то-ли это такая разбежка в результатах :) И какой результат все-таки точнее. Боюсь, что такие точные вычисления сильно убавят производительность. Жду следующую статью :)
Огромное спасибо за ссылки. Особенно за «Стековые процессоры, или новое — это хорошо забытое новое». Даже о дофине-1630 упомянуто. Правда не совсем точно. В свое время NC4016 поразил меня своей простотой. Я тогда не мог понять почему не имея никакой технологии у нас пытаются переделать 80х86 или PDP. А здесь на 2 микронах куча мегагерц :) Для 1630 я когда-то даже написал С-компилятор, отладчик и всю библиотеку функций. Куча трудов коту под хвост :)
Был на Интеграле еще и 8-ми разрядный стековый процессор. Но автор статьи видимо об этом не знал.
PS. Для любителей экзотики :).Asm код С-компилятора 1630.
/* j = --i + a++; */
>at=sp-29;
$_RD(>t,at);
t-0x1;
*at>=t;
>at=sp-27;
$_RD(>t,at);
>t+0x1;
*at>=t>;
t+n; n>;
>at=sp-28;
*at>=t;
А то жалко иметь хороший аппаратный double, а расчеты производить на счетах, потому что точнее выходит :)
2. насколько я помню, FMA избегает ненужного округления перед сложением. Результат будет немного точнее. Но у меня нет и FMA :) Только умножить и сложить.
3. Да, я использовал такой подход в синусе. Тот еще геморрой :)
А если в процессоре нет целочисленного умножения 64*64? :) Можно конечно через 32*32,
но это в 10 раз медленнее чем double*double.
Правильно ли я понимаю, что если в сложении-вычитании double можно производить вычисления большей точности за счет double.double, то с умножением такой фокус не проходит? Умножили и хвост безнадежно потерян. Кстати, в процессоре можно легко сделать вычисления double+double c сохранением обоих частей результата. Аналогично можно было бы сделать и с умножением. т.е. по аналогии с целыми числами int*int = long long.
Был бы какой-то толк от подобных инструкций?
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) константы из таблицы
Есть немного удивительные для меня результаты.
У меня есть таблица эталонов из 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 на данном значении.
Был на Интеграле еще и 8-ми разрядный стековый процессор. Но автор статьи видимо об этом не знал.
PS. Для любителей экзотики :).Asm код С-компилятора 1630.
/* j = --i + a++; */
>at=sp-29;
$_RD(>t,at);
t-0x1;
*at>=t;
>at=sp-27;
$_RD(>t,at);
>t+0x1;
*at>=t>;
t+n; n>;
>at=sp-28;
*at>=t;