Mathematical floating point rounding is weird in C ++ compared to math

Question

The following article is resolved, the problem arose due to an error in interpreting the formula at http://www.cplusplus.com/reference/random/piecewise_constant_distribution/ The reader is strongly recommended to consider the page: http://en.cppreference.com/w/cpp/numeric / random / piecewise_constant_distribution

I have the following strange phenomenon that puzzles me !:

I have a piecewise constant probability density indicated as

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()); 

When I draw n random numbers (the radius of the sphere in millimeters) from this distribution and calculate the mass of the sphere and sum them up like this:

 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; } 

I get mass = 4.3283 kg (see LIVE here )

Doing the EXACT identical thing in Mathematica:

It gives an estimated correct value of 4.5287 kg . (see mathematica )

What's not the same, also with different seeds, C ++ and Mathematica never match !? Is this a numerical inaccuracy in which I doubt that it is ...? Question: what is wrong to crack using fetch in C ++?

Simple math code:

 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}] 

Update : I did some analysis:

• Reading the numbers (64-bit doubles) generated from Mathematica in C ++ → calculated the sum and gives the same as Mathematica
  Mass calculated by reduction: 4.52528010260687096888432279229

• Read in the numbers generated from C ++ (64-bit double) in Mathematica -> calculate the sum, and it gives the same 4.32402

• I almost complete the selection with std::piecewise_constant_distribution inaccurately (or exactly the same as with 64-bit floats) or have an error ... OR is there something wrong with my weights?

• std::piecewise_constant_distribution are calculated erroneously by std::piecewise_constant_distribution at http://coliru.stacked-crooked.com/a/ca171bf600b5148f ===> This seems to be a mistake!

Histogram Plot CPP Created values ​​compared to desired Distribution:

 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] 

File is created here: Pixel pricing

Answer 1

[The following paragraph has been edited for correctness. - Editor’s Note]

Mathematica may or may not use IEEE 754 floating-point numbers. From the Wolfram documentation:

The Wolfram language has sophisticated built-in automatic numerical accuracy and accuracy control. But for special optimization of numerical calculations or for the study of numerical analysis, the Tungsten language also allows detailed control of accuracy and accuracy.

and

Wolfram language processes both integers and real numbers with any number of digits, automatically putting numerical precision when necessary. Wolfram internally uses several highly optimized representations of numbers, but nonetheless provides a single interface for digital and precise manipulation, allowing numerical analysts to study the details of the representation as desired.

Answer 2

The probability formula seems to be spelled incorrectly for std::piecewise_constant_distribution on http://www.cplusplus.com/reference/random/piecewise_constant_distribution/

Summing weights is done without multiplying interval lengths!

The correct formula is: http://en.cppreference.com/w/cpp/numeric/random/piecewise_constant_distribution

This solves all the silly quirks previously discovered as an error with an error / floating point, etc.