c语言中浮点数4舍5入,与mathematica相比,C中的浮点数学四舍五入

以下帖子解决,问题是因为

http://www.cplusplus.com/reference/random/piecewise_constant_distribution/错过了公式的解释,强烈建议读者考虑页面：

http://en.cppreference.com/w/cpp/numeric/random/piecewise_constant_distribution

我有以下奇怪的现象让我感到困惑！：

我有一个分段常数概率密度

using RandomGenType = std::mt19937_64;

RandomGenType gen(51651651651);

using PREC = long double;

std::array intervals {0.59, 0.7, 0.85, 1, 1.18};

std::array weights {1.36814, 1.99139, 0.29116, 0.039562};

// integral over the pdf to normalize:

PREC normalization =0;

for(unsigned int i=0;i<4;i++){

normalization += weights[i]*(intervals[i+1]-intervals[i]);

}

std::cout << std::setprecision(30) << "Normalization: " << normalization << std::endl;

// normalize all weights (such that the integral gives 1)!

for(auto & w : weights){

w /= normalization;

}

std::piecewise_constant_distribution

distribution (intervals.begin(),intervals.end(),weights.begin());

当我从这个分布中绘制n个随机数(以毫米为单位的球体半径)并计算球体的质量并将它们总结为：

unsigned int n = 1000000;

double density = 2400;

double mass = 0;

for(int i=0;i

auto d = 2* distribution(gen) * 1e-3;

mass += d*d*d/3.0*M_PI_2*density;

}

我得到质量= 4.3283千克(见LIVE here)

在Mathematica中做完全相同的事情,如：

提供4.5287千克的正确值. (见mathematica)

哪个不一样,也有不同的种子,C和Mathematica从不匹配！ ？是数字不准确,我怀疑它是什么……？

问题：C中的采样有什么问题？

简单Mathematica代码：

pdf[r_] = 2*Piecewise[{{0, r < 0.59}, {1.36814, 0.59 <= r <= 0.7},

{1.99139, Inequality[0.7, Less, r, LessEqual, 0.85]},

{0.29116, Inequality[0.85, Less, r, LessEqual, 1]},

{0.039562, Inequality[1, Less, r, LessEqual, 1.18]},

{0, r > 1.18}}];

pdfr[r_] = pdf[r] / Integrate[pdf[r], {r, 0, 3}];(*normalize*)

Plot[pdf[r], {r, 0.4, 1.3}, Filling -> Axis]

PDFr = ProbabilityDistribution[pdfr[r], {r, 0, 1.18}];

(*if you put 1.18=2 then we dont get 4.52??*)

SeedRandom[100, Method -> "MersenneTwister"]

dataR = RandomVariate[PDFr, 1000000, WorkingPrecision -> MachinePrecision];

Fold[#1 + (2*#2*10^-3)^3 Pi/6 2400 &, 0, dataR]

(*Analytical Solution*)

PDFr = ProbabilityDistribution[pdfr[r], {r, 0, 3}];

1000000 Integrate[ 2400 (2 InverseCDF[PDFr, p] 10^-3)^3 Pi/6, {p, 0, 1}]

更新：

我做了一些分析：

>读入Mathematica生成的数字(64位双打)

C – >计算总和,它与Mathematica相同

通过还原计算的质量：4.52528010260687096888432279229

>将从C(64位双倍)生成的数字读入Mathematica – >计算总和,它给出相同的4.32402

>我几乎得出结论,使用std :: piecewise_constant_distribution的采样是不准确的(或者像64bit浮点数那样准确)或者有一个错误……或者我的权重有问题？

>在http://coliru.stacked-crooked.com/a/ca171bf600b5148f中错误地计算密度std :: piecewise_constant_distribution ===>这似乎是一个错误！

CPP的直方图生成值与所需分布相比较：

file = NotebookDirectory[] <> "numbersCpp.bin";

dataCPP = BinaryReadList[file, "Real64"];

Hpdf = HistogramDistribution[dataCPP];

h = DiscretePlot[ PDF[ Hpdf, x], {x, 0.4, 1.2, 0.001},

PlotStyle -> Red];

Show[h, p, PlotRange -> All]