与mathematica相比,C++中的浮点数学变得怪异

Question

以下帖子已经解决,问题是由于http://www.cplusplus.com/reference/random/piecewise_constant_distribution/上对公式的解释错误而引起的. 强烈建议读者考虑以下页面:http://en.cppreference .COM/W/CPP /数字/随机/ piecewise_constant_distribution

我有以下奇怪的现象让我感到困惑!:

给定的分段常数概率密度为

using RandomGenType = std::mt19937_64;
RandomGenType gen(51651651651);

using PREC = long double;
std::array<PREC,5> intervals {0.59, 0.7, 0.85, 1, 1.18};
std::array<PREC,4> 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<PREC>
distribution (intervals.begin(),intervals.end(),weights.begin());

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

unsigned int n = 1000000;
double density = 2400;
double mass = 0;

for(int i=0;i<n;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}]

更新:我做了一些分析:

1. 读入从Mathematica生成的数字(64位双精度)到C++ - >计算总和,它给出与Mathematica Mass相同
   的减少计算:4.52528010260687096888432279229

2. 读入从C++(64位双倍)到Mathematica生成的数字 - >计算总和并给出相同的4.32402

3. 我几乎得出的结论std::piecewise_constant_distribution是不准确(或者与64位浮点数一样准确)或有错误......或者我的权重有问题？

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

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]

该文件在此处生成:数字生成CPP

Answer 1

[为了正确性，对以下段落进行了编辑。——编者注]

Mathematica 可能使用也可能不使用 IEEE 754 浮点数。来自 Wolfram 文档：

Wolfram 语言具有复杂的内置自动数值精度和准确度控制。但对于数值计算的特殊目的优化或研究数值分析，Wolfram 语言还允许对精度和准确度进行详细控制。

和

Wolfram 语言处理任意位数的整数和实数，并在适当时自动标记数值精度. Wolfram 语言内部使用了几种高度优化的数字表示形式，但仍然为数字和精度操作提供了统一的接口，同时允许数值分析师在需要时研究表示细节。