По поводу неточности двоичного представления в форматах с плавающей точкой.
Округления до целого в FPU (типичное его использование под Windows — при генерации 32-битного кода) выполняется во время сохранения вершины стека в памяти инструкцией FISTP. Режим округления задаётся битами в Control Word. В Borland Delphi 5 в Control Word занесены значения соответствующие режиму округления к ближайшему, а в случае значения посередине — к чётному (\Source\RTL\Sys\system.pas)
[В случае необходимости сменить режим округления в некоторой функции перед её завершением происходит восстановление режима округления.]procedure _ROUND;
asm
{ -> FST(0) Extended argument }
{ <- EDX:EAX Result }
SUB ESP,8
FISTP qword ptr [ESP]
FWAIT
POP EAX
POP EDX
end;
В Borland Delphi 5 значения типа Longint возвращаются через регистры EAX и EAX. «Функция» выделяет дополнительно два двойных слова на стеке (SUB ESP,8), а затем перед «завершением» восстанавливает состояние стека (POP EAX, POP EDX).Следующий текст тестит округление при использовании FPU
Код:
var
i: Integer;
d1, d1m, d2, d2m, d3, d3m : Double;
f: Text;
begin
i: Integer;
d1, d1m, d2, d2m, d3, d3m : Double;
f: Text;
begin
AssignFile(f, 'out.txt');
Rewrite(f);
for i:= 0 to 50
do begin
d1:= 0.05+i/10; d1m:= 10*d1;
d2:= 0.0005+i/1000; d2m:= 1000*d2;
d3:= 0.00005+i/10000; d3m:= 10000*d3;
writeln(f,d1:4:2,' -> ',Round(d1m)/10:3:1,'; ',d2:6:4,' -> ',Round(d2m)/1000:5:3,'; ',d3:7:5,' -> ',Round(d3m)/10000:6:4);
end;
CloseFile(f)
end.
Содержимое файла out.txt (при выполнении на процессоре AMD, на Intel не пробовал):
- 0.05 -> 0.0; 0.0005 -> 0.000; 0.00005 -> 0.0000
- 0.15 -> 0.2; 0.0015 -> 0.002; 0.00015 -> 0.0001
- 0.25 -> 0.2; 0.0025 -> 0.002; 0.00025 -> 0.0002
- 0.35 -> 0.4; 0.0035 -> 0.004; 0.00035 -> 0.0004
- 0.45 -> 0.4; 0.0045 -> 0.004; 0.00045 -> 0.0004
- 0.55 -> 0.6; 0.0055 -> 0.006; 0.00055 -> 0.0006
- 0.65 -> 0.6; 0.0065 -> 0.006; 0.00065 -> 0.0006
- 0.75 -> 0.8; 0.0075 -> 0.008; 0.00075 -> 0.0008
- 0.85 -> 0.8; 0.0085 -> 0.008; 0.00085 -> 0.0008
- 0.95 -> 1.0; 0.0095 -> 0.010; 0.00095 -> 0.0010
- 1.05 -> 1.0; 0.0105 -> 0.010; 0.00105 -> 0.0010
- 1.15 -> 1.2; 0.0115 -> 0.012; 0.00115 -> 0.0012
- 1.25 -> 1.2; 0.0125 -> 0.012; 0.00125 -> 0.0012
- 1.35 -> 1.4; 0.0135 -> 0.014; 0.00135 -> 0.0014
- 1.45 -> 1.4; 0.0145 -> 0.014; 0.00145 -> 0.0014
- 1.55 -> 1.6; 0.0155 -> 0.016; 0.00155 -> 0.0016
- 1.65 -> 1.6; 0.0165 -> 0.016; 0.00165 -> 0.0016
- 1.75 -> 1.8; 0.0175 -> 0.018; 0.00175 -> 0.0018
- 1.85 -> 1.8; 0.0185 -> 0.018; 0.00185 -> 0.0018
- 1.95 -> 2.0; 0.0195 -> 0.020; 0.00195 -> 0.0020
- 2.05 -> 2.0; 0.0205 -> 0.020; 0.00205 -> 0.0020
- 2.15 -> 2.2; 0.0215 -> 0.022; 0.00215 -> 0.0022
- 2.25 -> 2.2; 0.0225 -> 0.022; 0.00225 -> 0.0022
- 2.35 -> 2.4; 0.0235 -> 0.024; 0.00235 -> 0.0024
- 2.45 -> 2.4; 0.0245 -> 0.024; 0.00245 -> 0.0024
- 2.55 -> 2.6; 0.0255 -> 0.026; 0.00255 -> 0.0026
- 2.65 -> 2.6; 0.0265 -> 0.026; 0.00265 -> 0.0026
- 2.75 -> 2.8; 0.0275 -> 0.028; 0.00275 -> 0.0028
- 2.85 -> 2.8; 0.0285 -> 0.028; 0.00285 -> 0.0028
- 2.95 -> 3.0; 0.0295 -> 0.030; 0.00295 -> 0.0030
- 3.05 -> 3.0; 0.0305 -> 0.030; 0.00305 -> 0.0031
- 3.15 -> 3.2; 0.0315 -> 0.032; 0.00315 -> 0.0032
- 3.25 -> 3.2; 0.0325 -> 0.032; 0.00325 -> 0.0032
- 3.35 -> 3.4; 0.0335 -> 0.034; 0.00335 -> 0.0034
- 3.45 -> 3.4; 0.0345 -> 0.034; 0.00345 -> 0.0034
- 3.55 -> 3.6; 0.0355 -> 0.036; 0.00355 -> 0.0036
- 3.65 -> 3.6; 0.0365 -> 0.036; 0.00365 -> 0.0036
- 3.75 -> 3.8; 0.0375 -> 0.038; 0.00375 -> 0.0038
- 3.85 -> 3.8; 0.0385 -> 0.038; 0.00385 -> 0.0038
- 3.95 -> 4.0; 0.0395 -> 0.040; 0.00395 -> 0.0040
- 4.05 -> 4.0; 0.0405 -> 0.040; 0.00405 -> 0.0040
- 4.15 -> 4.2; 0.0415 -> 0.042; 0.00415 -> 0.0042
- 4.25 -> 4.2; 0.0425 -> 0.042; 0.00425 -> 0.0042
- 4.35 -> 4.4; 0.0435 -> 0.044; 0.00435 -> 0.0044
- 4.45 -> 4.4; 0.0445 -> 0.044; 0.00445 -> 0.0044
- 4.55 -> 4.6; 0.0455 -> 0.046; 0.00455 -> 0.0046
- 4.65 -> 4.6; 0.0465 -> 0.046; 0.00465 -> 0.0046
- 4.75 -> 4.8; 0.0475 -> 0.048; 0.00475 -> 0.0048
- 4.85 -> 4.8; 0.0485 -> 0.048; 0.00485 -> 0.0048
- 4.95 -> 5.0; 0.0495 -> 0.050; 0.00495 -> 0.0050
- 5.05 -> 5.0; 0.0505 -> 0.050; 0.00505 -> 0.0050
Видно, что с увеличением позиции вправо от десятичной точки на 3, по сравнению с округлением до целого, начинают проявляться «сбои» округления. В частности, 0.00015 округляется до 0.0001, а не до ожидаемого 0.0002.
Выглядит так, что Extended (в который преобразуется Double при загрузке в FPU и в котором выполняются все операции в FPU, в том числе умножение)
недостаточно для корректного вычисления произведения:
10000*(0.00005+1/1000) в формате Double будет равно 3FF7FFFFFFFFFFFF =
00111111111 10111111111111111111111111111111111111111111111111111
1.5 в формате Double будет равно 3FF8000000000000 =
00111111111 11000000000000000000000000000000000000000000000000000.
Видно, что мантисса первого числа меньше мантиссы второго числа, поэтому округление и приводит к 0.0001.
В 64-битном коде под Windows, который используется в пакетах или генерируется компиляторами, как правило, используется SSE/AVX (нет поддержки Extended). Результат выполнения (на процессоре AMD), скомпилированного под Win64 в Embarcadero Delphi XE7 предыдущего исходного текста (файл out.txt):
- 0.05 -> 0.0; 0.0005 -> 0.000; 0.00005 -> 0.0000
- 0.15 -> 0.2; 0.0015 -> 0.002; 0.00015 -> 0.0002
- 0.25 -> 0.2; 0.0025 -> 0.002; 0.00025 -> 0.0002
- 0.35 -> 0.4; 0.0035 -> 0.004; 0.00035 -> 0.0004
- 0.45 -> 0.4; 0.0045 -> 0.005; 0.00045 -> 0.0004
- 0.55 -> 0.6; 0.0055 -> 0.006; 0.00055 -> 0.0006
- 0.65 -> 0.6; 0.0065 -> 0.007; 0.00065 -> 0.0006
- 0.75 -> 0.8; 0.0075 -> 0.008; 0.00075 -> 0.0008
- 0.85 -> 0.8; 0.0085 -> 0.008; 0.00085 -> 0.0008
- 0.95 -> 1.0; 0.0095 -> 0.010; 0.00095 -> 0.0010
- 1.05 -> 1.0; 0.0105 -> 0.010; 0.00105 -> 0.0010
- 1.15 -> 1.2; 0.0115 -> 0.012; 0.00115 -> 0.0012
- 1.25 -> 1.2; 0.0125 -> 0.012; 0.00125 -> 0.0012
- 1.35 -> 1.4; 0.0135 -> 0.014; 0.00135 -> 0.0013
- 1.45 -> 1.4; 0.0145 -> 0.014; 0.00145 -> 0.0014
- 1.55 -> 1.6; 0.0155 -> 0.016; 0.00155 -> 0.0016
- 1.65 -> 1.6; 0.0165 -> 0.016; 0.00165 -> 0.0016
- 1.75 -> 1.8; 0.0175 -> 0.018; 0.00175 -> 0.0017
- 1.85 -> 1.8; 0.0185 -> 0.018; 0.00185 -> 0.0018
- 1.95 -> 2.0; 0.0195 -> 0.020; 0.00195 -> 0.0020
- 2.05 -> 2.0; 0.0205 -> 0.020; 0.00205 -> 0.0020
- 2.15 -> 2.2; 0.0215 -> 0.022; 0.00215 -> 0.0022
- 2.25 -> 2.2; 0.0225 -> 0.022; 0.00225 -> 0.0023
- 2.35 -> 2.3; 0.0235 -> 0.024; 0.00235 -> 0.0024
- 2.45 -> 2.4; 0.0245 -> 0.024; 0.00245 -> 0.0024
- 2.55 -> 2.6; 0.0255 -> 0.026; 0.00255 -> 0.0026
- 2.65 -> 2.6; 0.0265 -> 0.026; 0.00265 -> 0.0026
- 2.75 -> 2.8; 0.0275 -> 0.028; 0.00275 -> 0.0028
- 2.85 -> 2.8; 0.0285 -> 0.028; 0.00285 -> 0.0028
- 2.95 -> 2.9; 0.0295 -> 0.030; 0.00295 -> 0.0030
- 3.05 -> 3.0; 0.0305 -> 0.030; 0.00305 -> 0.0031
- 3.15 -> 3.2; 0.0315 -> 0.032; 0.00315 -> 0.0032
- 3.25 -> 3.2; 0.0325 -> 0.032; 0.00325 -> 0.0032
- 3.35 -> 3.4; 0.0335 -> 0.034; 0.00335 -> 0.0034
- 3.45 -> 3.4; 0.0345 -> 0.034; 0.00345 -> 0.0034
- 3.55 -> 3.6; 0.0355 -> 0.036; 0.00355 -> 0.0036
- 3.65 -> 3.6; 0.0365 -> 0.036; 0.00365 -> 0.0036
- 3.75 -> 3.8; 0.0375 -> 0.038; 0.00375 -> 0.0038
- 3.85 -> 3.8; 0.0385 -> 0.038; 0.00385 -> 0.0038
- 3.95 -> 4.0; 0.0395 -> 0.040; 0.00395 -> 0.0039
- 4.05 -> 4.0; 0.0405 -> 0.040; 0.00405 -> 0.0040
- 4.15 -> 4.1; 0.0415 -> 0.042; 0.00415 -> 0.0042
- 4.25 -> 4.2; 0.0425 -> 0.042; 0.00425 -> 0.0042
- 4.35 -> 4.4; 0.0435 -> 0.044; 0.00435 -> 0.0044
- 4.45 -> 4.4; 0.0445 -> 0.044; 0.00445 -> 0.0044
- 4.55 -> 4.6; 0.0455 -> 0.046; 0.00455 -> 0.0045
- 4.65 -> 4.6; 0.0465 -> 0.046; 0.00465 -> 0.0046
- 4.75 -> 4.8; 0.0475 -> 0.048; 0.00475 -> 0.0048
- 4.85 -> 4.8; 0.0485 -> 0.048; 0.00485 -> 0.0048
- 4.95 -> 5.0; 0.0495 -> 0.050; 0.00495 -> 0.0049
- 5.05 -> 5.0; 0.0505 -> 0.050; 0.00505 -> 0.0050
Уже при округлении до десятых наблюдается «неожиданные» результаты:
2.35->2.3; 2.95->2.9; 4.15->4.1.
При округлении до сотых:
.0045->.005; .0065->.007.
При округлении до тысячных:
.00135->.0013; .00175->.0017; .00225->.0023; .00305->.0031 и др.
Напрашивается вывод: если бы для умножения был использован binary128 (аппаратный, например, Itanium или эмуляция на регистрах общего назначения), то можно было бы округлять до тысячных или даже до больших знаков после запятой.
Помимо преобразования десятичных к рациональным в Wolfram Mathematica можно задать точность. Не прикидывал, а просто взял 50 (Mathematica 12):
Код:
Table[{.00005 + i/10000, Round[N[10000*(.00005`50 + i/10000)]]/10000.0}, {i, 0, 50}] // Column
{{{0.00005, 0.}},
{{0.00015, 0.0002}},
{{0.00025, 0.0002}},
{{0.00035, 0.0004}},
{{0.00045, 0.0004}},
{{0.00055, 0.0006}},
{{0.00065, 0.0006}},
{{0.00075, 0.0008}},
{{0.00085, 0.0008}},
{{0.00095, 0.001}},
{{0.00105, 0.001}},
{{0.00115, 0.0012}},
{{0.00125, 0.0012}},
{{0.00135, 0.0014}},
{{0.00145, 0.0014}},
{{0.00155, 0.0016}},
{{0.00165, 0.0016}},
{{0.00175, 0.0018}},
{{0.00185, 0.0018}},
{{0.00195, 0.002}},
{{0.00205, 0.002}},
{{0.00215, 0.0022}},
{{0.00225, 0.0022}},
{{0.00235, 0.0024}},
{{0.00245, 0.0024}},
{{0.00255, 0.0026}},
{{0.00265, 0.0026}},
{{0.00275, 0.0028}},
{{0.00285, 0.0028}},
{{0.00295, 0.003}},
{{0.00305, 0.003}},
{{0.00315, 0.0032}},
{{0.00325, 0.0032}},
{{0.00335, 0.0034}},
{{0.00345, 0.0034}},
{{0.00355, 0.0036}},
{{0.00365, 0.0036}},
{{0.00375, 0.0038}},
{{0.00385, 0.0038}},
{{0.00395, 0.004}},
{{0.00405, 0.004}},
{{0.00415, 0.0042}},
{{0.00425, 0.0042}},
{{0.00435, 0.0044}},
{{0.00445, 0.0044}},
{{0.00455, 0.0046}},
{{0.00465, 0.0046}},
{{0.00475, 0.0048}},
{{0.00485, 0.0048}},
{{0.00495, 0.005}},
{{0.00505, 0.005}}}