Деление нацело двух double-овых чисел

NataNata

Замечена странная ситуация
Если в с++-ной проге строка типа long C = (longA / B где A и B - double-ы. Нездоровая ситуация возникает, когда A нацело делится на B, скажем, для примера, A = 18.0, B = 6.0, A/B=3. Замечено, что на разных процах получаются разные результаты, если попытаться поделить такие вот (не как в примере, но нацело делящиеся) double-ы. Грубо говоря, где-то получается 2.99999999999999, поэтому C=2, где-то 3.00000000000000, поэтому C=3. Отличие на уровне 1e-14.
Тестировалось на quad core q9650, xeon x5667, amd opteron хр какой-то, xeon x3430. Прога компилировалась на msvs без SSE\SSE2. Переключение м\у debug и release- билдами ничего не меняет. На Xeon X3430 результаты получаются иные, нежели на других приведенных процах, по озвученной причине. Оси везде разные, но это влиять не должно, как я понимаю.
Вопрос. Пусть такова особенность процов, так сказать. Но как избежать приведенной ситуации? Достаточно ли будет поменять код на long C = (longA / B + DBL_EPSILON где DBL_EPSILON - некоторая константа, отвечающая за точность работы в dоuble-ах?
З.Ы. да, экспериментально наблюдалось, что программа вида int main { double A = 18.0; double B = 6.0; long C = (longA / B); printf("%d\n", C); return 0; } для разных процессоров потенциально может выдать разные результаты (если правильно подобрать A и B).

Devid

Сей факт весьма широко известен.
(A / B + DBL_EPSILON) должно помочь, только DBL_EPSILON зависит от того, что ты в итоге хочешь получать, а не от точности дабла.
Т.е. взяв A=18 и B=3.000...0001 ты везде получишь что-то типа 2.999...999, при этом когда 2.999...999 должно превращаться в 3 должен будешь решать ты сам.

Inflict84

Могу предложить другой (медленный) вариант:
C = round(N*A)/round(N*B где, как и в твоём варианте, N зависит от того, что ты хочешь получить.

NataNata

меня смущает, что на разных компах - одни и те же exe-шники, одни и те же входные данные. и результаты получаются разными. Как-то это странно.

Devid

Ты сам указал причину, разные процы.

tokuchu

А не лучше ли в таком случае делать floor(result + 0.5)?

smit1

З.Ы. да, экспериментально наблюдалось, что программа вида int main { double A = 18.0; double B = 6.0; long C = (longA / B); printf("%d\n", C); return 0; } для разных процессоров потенциально может выдать разные результаты (если правильно подобрать A и B).
Будь любезен, запости сюда точный код этой программы.
А так же табличку типа "компилятор | ОС | процессор | результат".

evolet

честно говоря непонятно, в чем твоя проблема
при плавающих вычислениях на границах всегда получается неопределенность, которая возникает в конечном счете из-за машинных округлений промежуточных вычислений, причем эти округления могут проистекать не только из-за органиченной разрадности плавающих чисел, но и из-за неточности аппаратной реализации некоторых операций (причем, как ты сам указал, эта реализация может различается даже на процах одного производителя)
имхо вещественные числа надо рассматривать как интервалы (а-E, a+E) в том смысле, что распечатать это а конечно можно, но чему же равно "математическое число" (не бинарное представление числа, а идеалное — то которое должно быть если бы все промежуточные вычисления были бы абсолютно точными) - неизвестно - извесно, что оно где-то в этом интервале.
Таким образом если у тебя в программе стоит уcловие (3.0-E, 3.0+E) >= 3.0 (какое неявно присутствует при преобразовании в твоем случае) - то это условие может быть как ложным, так и истинным, в зависимости от, можно считать, неизвестной случайно величины.
PS есть еще такой момент, что эта E может очень быстро вырасти до значимых величин и полученный плавающий результат перестанет иметь какой-либо смысл. Для простых вычислений можно по абсолютным/относительным погрешностям на пальцах прикидывать, что ошиба не выйдет за допустимые пределы. Но если вычислять что-то нетривиальное, то оценка этих ошибок очень быстро превращается в математическую теорию и как-бы без нее получится только "вроде работает, но это только до первого случая неудобных входных данных", например (если я не путаю с давних времен "численных методов" ) если систему линейных уравнений решать методом Гаусса - у этого способа получается самая большая ошибка (вроде она там даже к бесконечности стремится при детерминанте не стремящемся к нулю)

Serab

например (если я не путаю с давних времен "численных методов" ) если систему линейных уравнений решать методом Гаусса - у этого способа получается самая большая ошибка (вроде она там даже к бесконечности стремится при детерминанте не стремящемся к нулю)
там число обусловленности используется, да, методом гаусса плохо решать в общем случае.

lubanj

ой ну в методе гаусса же можно выбрать другой столбец

Serab

) какое это отношение имеет к теме треда?
2) насколько я помню, есть контрпримеры, когда выбор главного элемента все равно дает неустойчивое решение.

kill-still

ботай
FNCLEX и FLDCW
не слушай их.
одинаковые наборы инструкций должны давать одинаковые результаты, если в одинаковом состоянии сопроцессор.

evolet

IEEE 754 не определяет с точностью до битов результат операций, чем производители процов и пользуются.
Задрачивать плавающий управляющий регистр в большинстве случаев не надо. Надо писать программы так, чтобы учитывать возможную плавающую ошибку.
PS даже при одинаковом состоянии сопроцессора все равно будут различия, в sin\cos etc точно, для +-*/ - хз, для деления может и такая же фигня получается.

Papazyan

Господа, все обдумано до вас - http://download.oracle.com/docs/cd/E19957-01/806-3568/ncg_go...

Maurog

Господа, все обдумано до вас
а теперь кратко по теме на русском просвети нас

Papazyan

Не знаешь английского что-ли?

Maurog

не знаешь русского что ли?

alexkravchuk

Достаточно ли будет поменять код на long C = (longA / B + DBL_EPSILON где DBL_EPSILON - некоторая константа, отвечающая за точность работы в dоuble-ах?
   Интересное начнётся, когда либо A, либо B — отрицательные :)
   А по теме, лучше использовать функцию lround ( long lround(double) )
#include <stdio.h>
#include <math.h>

int main(void)
{
double aa;

aa = 18.1;
printf("%lf %d %d %d \n", aa, (int) aa, (int) (aa+0.2 (int) lround(aa;

aa = 17.9;
printf("%lf %d %d %d \n", aa, (int) aa, (int) (aa+0.2 (int) lround(aa;

aa = -18.1;
printf("%lf %d %d %d \n", aa, (int) aa, (int) (aa+0.2 (int) lround(aa;

aa = -17.9;
printf("%lf %d %d %d \n", aa, (int) aa, (int) (aa+0.2 (int) lround(aa;

return 0;
}

Serab

Не знаешь английского что-ли?
по теме треда напиши. Статью ту я читал, кстати, целиком, правда где-то полгода назад, но уже запамятовал.

Papazyan

Many programmers may not realize that even a program that uses only the numeric formats and operations prescribed by the IEEE standard can compute different results on different systems. In fact, the authors of the standard intended to allow different implementations to obtain different results. Their intent is evident in the definition of the term destination in the IEEE 754 standard: "A destination may be either explicitly designated by the user or implicitly supplied by the system (for example, intermediate results in subexpressions or arguments for procedures). Some languages place the results of intermediate calculations in destinations beyond the user's control. Nonetheless, this standard defines the result of an operation in terms of that destination's format and the operands' values." (IEEE 754-1985, p. 7) In other words, the IEEE standard requires that each result be rounded correctly to the precision of the destination into which it will be placed, but the standard does not require that the precision of that destination be determined by a user's program. Thus, different systems may deliver their results to destinations with different precisions, causing the same program to produce different results (sometimes dramatically so even though those systems all conform to the standard.
Т.е. если один процессор вычисляет дробь даблов в более точном формате, то приведение к лонг не сработает как ожидается. Другой процессор может округлить результат до дабл и получится ожидаемый результат.

bleyman

IEEE 754 не определяет с точностью до битов результат операций, чем производители процов и пользуются.
Да блджад, при чём тут производители процов вообще?
Ви таки не поверите, но IEEE 754 таки определяет с точностью до битов результат операций, когда эти операции выглядят как "взять эти значения, произвести операцию, выдать результат".
Но оно совершенно не запрещает производителям процов не делать именно этого в случаях когда результат не пишется сразу в память, и производителям компиляторов — компилировать программы в код, который не пишет все промежуточные результаты в память.
I'm shaking in fear in disgust читая про "недетерминированность плавающих точек" и прочее. Всё там детерминированно, наверняка поболее чем в других местах твоей программы (где ты надеешься на two-complement result of signed integer overflow, например).
Добавь /fp:strict или стратегическое volatile на временные переменные.

evolet

мда... в стандарте прописанные операции действительно определены с точностью до битов (кроме неких "rare exceptions" :))
http://www.cs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF
"
Algebraic operations covered by IEEE 754, namely + , - , · , / , Ö and Binary <-> Decimal Conversion with
rare exceptions, must be Correctly Rounded to the precision of the operation’s destination unless the programmer
has specified a rounding other than the default. If it does not Overflow, a correctly rounded operation’s error
cannot exceed half the gap between adjacent floating-point numbers astride the operation’s ideal ( unrounded )
result. Half-way cases are rounded to Nearest Even,
"
но в практическом плане все-равно имеется "недетерминированность плавающих точек", пусть только из-за округлений (потому что котролировать эти округления слишком геморно, если возможно)
PS если хочется повторяемости, то с промежуточными 80 битными вычислениями имхо лучше бороться выставлением режима точности в 64 в fpcw (Precision control, которое у интела почму-то называется "Double Precision (53 bits)" )
( http://www.intel.com/Assets/PDF/manual/253665.pdf )
Оставить комментарий
Имя или ник:
Комментарий: