Как наиболее Ъ занулить последние байты мантиссы для double и float?

Меня вообще младшие биты не волнуют. Но Вы тред то почитайте, в итоге совсем другое решение.

Пиши так, чтобы чтобы минимизировать накопление ошибки. Это же самое главное в численных методах.

max(abs(0.00000001*a )+abs(0.00000001*b))

Потому что мне нужно сравнивать double. Есть альтернативы?

да. Не нужно сравнивать на точное равенство double, я уверен, что тебе нужен диапазон. Вот его и ищи. И считать его в double получается обычно быстрее, чем переводить в long. Т.е. вместо x==y, считай (x-y)/((x+y)/2)<eps. Данную формулу конечно можно и оптимизировать, но чаще оно быстрее и так считать. Ты как раз и пытаешься это сделать.

Тред читать уже не модно? Во первых, деление это очень дорогая операция. Во вторых, я уже это сделал.

Я знаю, что у меня последние 16 бит невалидные, а актуальные значение (которое должно сравниваться) лежит в первых скольки-то-там битах. Зачем мне сравнивать числа с учетом заведомо ошибочного хвоста?

потому-что диапазон ошибки вовсе не обязан укладываться в целое число бит, и иметь середину на нуле, как ты почему-то полагаешь. К тому-же, тебе всё равно считать ненужные биты, даже если ты их и обнулишь. Но если ты их обнулишь, то ты во первых сдвинешь ошибку (она не нулевая, а ты приравняешь её к нулю. Это приведёт к накоплению положительной ошибки), а во вторых потеряешь значение ошибки, т.е. фактически увеличишь её диапазон. Например, если у тебя неправильные биты 0 и 1, а правильные 2 и 3, и при этом биты 2 и 3 не равны нулю оба, то сбрасывая битвы 3210, ты сбрасываешь ПРАВИЛЬНЫЕ биты 32, УВЕЛИЧИВАЯ погрешность.

ЗЫЖ если ты хочешь быстрое сравнение, то сравнивай сначала порядки. Очевидно, что числа равны, тогда и только тогда, когда модуль разницы порядков не более 1. Отбрасываешь случаи большего модуля, и получаешь 9 случаев, из которых 2 тривиальных, и получается 7 случаев. Денормализуешь мантиссу отдельно для всех 6и случаев, а для 7го уже денормализовано. Ну и сравниваешь теперь их как целые. Разность это и есть твой нужный тебе eps с учётом порядка.

Во первых, деление это очень дорогая операция.

а кто говорил про деление? В твоём случае деление отлично меняется на вычитание порядков. Ведь тебе известен диапазон мантисс (1<m≤½), потому делить мантисы не нужно (конечно при этом диапазон ошибки расширяется, но это и не страшно, ибо тебе этот диапазон совсем не нужно знать с «как можно меньшей точностью».)

Во вторых, я уже это сделал.

я вижу. В принципе потянет, хотя и не самое оптимальное. Только не хватает проверки на 0, NaN и INF. С ними твоя функция может работать неправильно.

Так да, фиксед пойнт хватит всем же. Ну, почти всем.

-e < (a-b)/(a+b) < e

Разница минимальная

Деление это ЕМНИП порядка 10ти тактов.

ЛОР тебе говорит, что ты делаешь не то и спрашиваешь не так. Какой вопрос - такой ответ

Тебе, ИМХО, лучше вообще не писать программы

Код: http://paste.debian.net/11144/

Результаты: http://paste.debian.net/11145/

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

Не говори мне что мне надо делать, и я не скажу тебе куда идти(с)

#include<omp.h>
#include<stdio.h>
#include<stdlib.h>

const int N = 1000;
double arr[N];

int main(int argc, const char** argv){
	for(int i=0; i<N; i++) arr[i] = i;

	if(argc!=2){ printf("usage: ./a.out CPU_freq[MHz]\n"); return 1; }
	double Nf = 1e6*atof(argv[1])/N;
	double t0 = omp_get_wtime();
	for(int i=0; i<Nf; i++) for(int j=0; j<N-1; j++) arr[j] /= 1.001;
	double t1 = omp_get_wtime();
	for(int i=0; i<Nf; i++) for(int j=0; j<N-1; j++) arr[j] -= .1;
	double t2 = omp_get_wtime();
	printf("x=%g div=%g sum=%g\n", arr[N-1], t1-t0, t2-t1);
}

$ g++ -O3 test.cpp -lgomp
$ ./a.out 1e3
x=999 div=41.5334 sum=1.01054

разница в 40 раз. Чем рассказывать мне про кэш, лучше подумайте почему у Вас такое медленное суммирование.

Хмм. Не думал, что деление настолько медленнее.

Внезапно, умножение тоже занимает один такт. И умножения на фоне сложения (или наоборот - зависит от того чего больше) не видно вообще.

Странно. Почему так? Разве алгоритмическая сложность деления не такая же, как и умножения?

Для умножения есть эффективные алгоритмы. И так уж устроено нынешнее железо...

А у Вас в примере наск я понимаю во первых латентность операций заиграла (результат обработки след элемента зависит от пред.), во вторых данные здоровые, не закэшировались. Ну и собирали мб без оптимизации...

Нашел статью, в которой описывается деление через умножение (Goldschmidt division algorithm). Но в универе нам рассказывали про деление только на сумматорах. Сейчас даже в книге проверил и не нашел умножителей.

UPD; на вики все написано.

Вообще всегда в коде вместо

... x/3. ...

const double _3 = 1./3.;

... x*_3 ...

Хмм. Не думал, что деление настолько медленнее.

Потому его и заменяют умножением.