Hello, all

May I ask you if there is any chance of fixing the error in the Math unit that I reported here:

https://bugs.freepascal.org/view.php?id=32804

and provided a patch? Here is a test case:

Program Vartest;
uses Math;
const Size = 1000000;
var dataS: array of Single;
    dataD: array of Double;
    dataE: array of Extended;
    i,n: longint;
begin
  WriteLn('Each run should return a value near unity.');
  WriteLn('Single:');
  SetLength( dataS, Size );
  for n := 1 to 4 do
  begin
    for i := 0 to Size - 1 do
    begin
      dataS[i] := 10000000 + RandG(0,1);
    end;
    WriteLn( Math.Variance( dataS ):5:3 );
  end;
  WriteLn('Double:');
  SetLength( dataD, Size );
  for n := 1 to 4 do
  begin
    for i := 0 to Size - 1 do
    begin
      dataD[i] := 1000000000000000 + RandG(0,1);
    end;
    WriteLn( Math.Variance( dataD ):5:3 );
  end;
  WriteLn('Extended:');
  SetLength( dataE, Size );
  for n := 1 to 4 do
  begin
    for i := 0 to Size - 1 do
    begin
      dataE[i] := 1000000000000000000 + RandG(0,1);
    end;
    WriteLn( Math.Variance( dataE ):5:3 );
  end;
end.

Yes, the formula used is  SS-Sqr(S)/N, and this is the bad one. Explained somewhere in the Numerical Recipes, but also in Widepedia (https://en.wikipedia.org/wiki/Variance#Formulae_for_the_variance): "it is not advised for computer calculations as it suffers from catastrophic cancellation <https://en.wikipedia.org/wiki/Catastrophic_cancellation> if the two components of the equation are similar in magnitude and floating point arithmetic is used."

I have checked and applied the patch.

I will add your example to the documentation, if I have your permission ?

Michael.

Michael Van Canneyt to Anton Shepelev:

Thank you.

I will add your example to the documentation, if I have your permission?

Feel free to do so. The credit is due my friend Pe- ter Karpov, who originally discovered this bug, but you need not mention our names :-)

Good to see that Delphi 7 did not calculate the variance of double values correctly either (it can't do Single and Extended). XE2 returns the correct value, however.

‎Hello,

You may want to subscribe to the mailing list and set your subscription to the non-delivery mode if you want your messages sent via Gmane to arrive to the list faster (as it is now, they all have to wait for the moderator to let them through‎).

Tomas
(one of FPC mailing list moderators)