Is Delphi's Skewness correct

Question

In Delphi one can calculate Skewness using System.Math's function MomentSkewKurtosis().

var m1, m2, m3, m4, skew, k: Extended;
System.Math.MomentSkewKurtosis([1.1,
  3.345,
  12.234,
  11.945,
  14.235,
  16.876,
  20.213,
  11.001,
  7.098,
  21.234], m1, m2, m3, m4, skew, k);

• This will prints skew equal to -0.200371489809269.
• Minitab prints the value -0.24
• SigmaXL prints the value -0.23611

The reason is that Delphi does not not perform adjustment.

Here is my implementation which calculates adjustment:

function Skewness(const X: array of Extended; const Adjusted: Boolean): Extended;
begin
  var AMean := Mean(X);
  var xi_minus_mean_power_3 := 0.0;
  var xi_minus_mean_power_2 := 0.0;
  for var i := Low(X) to High(X) do
  begin
    xi_minus_mean_power_3 := xi_minus_mean_power_3 + IntPower((X[i] - AMean), 3);
    xi_minus_mean_power_2 := xi_minus_mean_power_2 + IntPower((X[i] - AMean), 2);
  end;
  // URL : https://www.gnu.org/software/octave/doc/v4.0.1/Descriptive-Statistics.html
  {                mean ((x - mean (x)).^3)
   skewness (X) = ------------------------.
                     std (x).^3
  }
  var N := Length(X);
  Result := xi_minus_mean_power_3 / N /
    IntPower(Sqrt(1 / N * xi_minus_mean_power_2), 3);
  // URL : https://www.gnu.org/software/octave/doc/v4.0.1/Descriptive-Statistics.html
  {                  sqrt (N*(N-1))   mean ((x - mean (x)).^3)
  skewness (X, 0) = -------------- * ------------------------.
                        (N - 2)             std (x).^3
  }
  if Adjusted then
    Result := Result * Sqrt(N * Pred(N)) / (N - 2);
end;

The helper routine IntPower is as follows:

function IntPower(const X: Extended; const N: Integer): Extended;

  /// <remarks>
  /// Calculate any float to non-negative integer power. Developed by Rory Daulton and used with permission. Last modified December 1998.
  /// </remarks>
  function IntPow(const Base: Extended; const Exponent: LongWord): Extended;
    { Heart of Rory Daulton's IntPower: assumes valid parameters &
    non-negative exponent }
  {$IFDEF WIN32}
  asm
      fld1                      // Result := 1
      cmp     eax, 0            // eax := Exponent
      jz      @@3
      fld     Base
      jmp     @@2
    @@1:    fmul    ST, ST      // X := Base * Base
    @@2:    shr     eax,1
      jnc     @@1
      fmul    ST(1),ST          // Result := Result * X
      jnz     @@1
      fstp    st                // pop X from FPU stack
    @@3:
      fwait
  end;
  {$ENDIF}
  {$IFDEF WIN64}
  begin
    Result := Power(Base, Exponent);
  end;
  {$ENDIF}
begin
  if N = 0 then
    Result := 1
  else if (X = 0) then
  begin
    if N < 0 then
      raise EMathError.Create('Zero cannot be raised to a negative power.')
    else
      Result := 0
  end
  else if (X = 1) then
    Result := 1
  else if X = -1 then
  begin
    if Odd (N) then
      Result := -1
    else
      Result := 1
  end
  else if N > 0 then
    Result := IntPow (X, N)
  else
  begin
    var P: LongWord;
    if N <> Low (LongInt) then
      P := Abs(N)
    else
      P := LongWord(High(LongInt)) + 1;
    try
      Result := IntPow(X, P);
    except
      on EMathError do
      begin
        Result := IntPow(1 / X, P); // try again with another method, perhaps less precise
        Exit;
      end;
    end;
    Result := 1 / Result;
  end;
end;

With that function the adjusted skewness becomes the accurate -0.237611357234441 matching Matlab and Minitab.

The only explanation I found is: https://octave.org/doc/v4.0.1/Descriptive-Statistics.html

"The adjusted skewness coefficient is obtained by replacing the sample second and third central moments by their bias-corrected versions."

Same goes with Kurtosis:

function Kurtosis(const X: array of Extended; const Adjusted: Boolean): Extended;
begin
  var AMean := Mean(X);
  var xi_minus_mean_power_4 := 0.0;
  var xi_minus_mean_power_2 := 0.0;
  for var i := Low(X) to High(X) do
  begin
    xi_minus_mean_power_4 := xi_minus_mean_power_4 + IntPower((X[i] - AMean), 4);
    xi_minus_mean_power_2 := xi_minus_mean_power_2 + IntPower((X[i] - AMean), 2);
  end;

  {     mean ((x - mean (x)).^4)
  k1 = ------------------------
              std (x).^4
  }
  var N := Length(X);
  Result := xi_minus_mean_power_4 / N /
    IntPower(1 / N * xi_minus_mean_power_2, 2);
  {             N - 1
  k0 = 3 + -------------- * ((N + 1) * k1 - 3 * (N - 1))
           (N - 2)(N - 3)
  }
  if Adjusted then
    // Mathlab, Minitab and SigmaXL do not add 3 (which is the kurtosis for normal distribution
    Result := {3 + }(N - 1) / ((N - 2) * (N - 3)) * ((N + 1) * Result - 3 * (N - 1));
end;

What is the reason for such adjustments and why Delphi decided not to implement it?