diff --git a/applications/lazstats/source/forms/analysis/statistical_process_control/basicspcunit.pas b/applications/lazstats/source/forms/analysis/statistical_process_control/basicspcunit.pas
index 64f156176..753fdbc5a 100644
--- a/applications/lazstats/source/forms/analysis/statistical_process_control/basicspcunit.pas
+++ b/applications/lazstats/source/forms/analysis/statistical_process_control/basicspcunit.pas
@@ -93,24 +93,6 @@ type
 var
   BasicSPCForm: TBasicSPCForm;
 
-// Constants for correction of standard deviation, needed by some charts.
-const
-  C4: array[2..25] of double = (
-    0.7979, 0.8862, 0.9213, 0.9400, 0.9515, 0.9594, 0.9650, 0.9693,
-    0.9727, 0.9754, 0.9776, 0.9794, 0.9810, 0.9823, 0.9835, 0.9845, 0.9854, 0.9862,
-    0.9869, 0.9876, 0.9882, 0.9887, 0.9892, 0.9896);
-
-  D3: array[1..24] of double = (
-    0, 0, 0, 0, 0, 0.076, 0.136, 0.184, 0.223, 0.256, 0.283, 0.307, 0.328,
-    0.347, 0.363, 0.378, 0.391, 0.403, 0.415, 0.425, 0.434, 0.443,
-    0.451, 0.459
-  );
-  D4: array[1..24] of double = (
-    3.267, 2.574, 2.282, 2.114, 2.004, 1.924, 1.864, 1.816, 1.777,
-    1.744, 1.717, 1.693, 1.672, 1.653, 1.637, 1.622, 1.608, 1.597,
-    1.585, 1.575, 1.566, 1.557, 1.548, 1.541
-  );
-
 
 implementation
 
@@ -139,6 +121,7 @@ begin
 end;
 
 
+// To be overridden by descendant forms
 procedure TBasicSPCForm.Compute;
 begin
   //
diff --git a/applications/lazstats/source/forms/analysis/statistical_process_control/rchartunit.pas b/applications/lazstats/source/forms/analysis/statistical_process_control/rchartunit.pas
index 8992d180f..59078d956 100644
--- a/applications/lazstats/source/forms/analysis/statistical_process_control/rchartunit.pas
+++ b/applications/lazstats/source/forms/analysis/statistical_process_control/rchartunit.pas
@@ -25,6 +25,19 @@ uses
   Math,
   Globals, Utils, MainUnit, DataProcs;
 
+// Constants for correction of standard deviation
+const
+  D3: array[2..25] of double = (
+    0, 0, 0, 0, 0, 0.076, 0.136, 0.184, 0.223, 0.256, 0.283, 0.307, 0.328,
+    0.347, 0.363, 0.378, 0.391, 0.403, 0.415, 0.425, 0.434, 0.443,
+    0.451, 0.459
+  );
+  D4: array[2..25] of double = (
+    3.267, 2.574, 2.282, 2.114, 2.004, 1.924, 1.864, 1.816, 1.777,
+    1.744, 1.717, 1.693, 1.672, 1.653, 1.637, 1.622, 1.608, 1.597,
+    1.585, 1.575, 1.566, 1.557, 1.548, 1.541
+  );
+
 procedure TRChartForm.Compute;
 var
   i, j: Integer;
@@ -116,8 +129,8 @@ begin
   seMean := seMean / sqrt(numValues);
   grandMean := grandMean / numValues;
   grandRange := grandRange / numGrps;
-  D3Value := D3[grpSize-1];
-  D4Value := D4[grpSize-1];
+  D3Value := D3[grpSize];
+  D4Value := D4[grpSize];
   UCL := D4Value * grandRange;
   LCL := D3Value * grandRange;
 
diff --git a/applications/lazstats/source/forms/analysis/statistical_process_control/schartunit.pas b/applications/lazstats/source/forms/analysis/statistical_process_control/schartunit.pas
index 5e5beca15..19039b703 100644
--- a/applications/lazstats/source/forms/analysis/statistical_process_control/schartunit.pas
+++ b/applications/lazstats/source/forms/analysis/statistical_process_control/schartunit.pas
@@ -1,3 +1,11 @@
+{ This unit was checked against a commercial statistical software  and
+  creates correct results.
+
+  Data file for testing: "boltsize.laz"
+  Group variable:         LotNo
+  Selected variable:      BoltLngth
+}
+
 unit SChartUnit;
 
 {$mode objfpc}{$H+}
@@ -40,11 +48,9 @@ var
   numValues: Integer = 0;
   grp: String;
   grpIndex: Integer;
-  X, Xsq: Double;
+  X, Xsq, Xmin, Xmax: Double;
   seMean, grandMean, grandSigma, grandSD: Double;
-  B, C4, gamma, D3Value, D4Value: Double;
-  xmin, xmax: Double;
-  sizeError: Boolean;
+  B, C4: Double;
   i, j: Integer;
   lReport: TStrings;
 begin
@@ -63,7 +69,6 @@ begin
   seMean := 0.0;
   grandMean := 0.0;
   grandSigma := 0.0;
-  sizeError := false;
 
   // Count "good" data points
   numValues := 0;
@@ -101,12 +106,15 @@ begin
     grpSize := count[j];
     if j = 0 then oldGrpSize := grpSize;
     if oldGrpSize <> grpSize then
-      sizeError := true;
+    begin
+      ErrorMsg('All groups must have the same size.');
+      exit;
+    end;
   end;
 
-  if (grpSize < 2) or (grpSize > 25) or sizeError then
+  if (grpSize < 2) then
   begin
-    ErrorMsg('Group size error.');
+    ErrorMsg('Groups with at least two values required.');
     exit;
   end;
 
@@ -117,14 +125,9 @@ begin
   seMean := seMean / sqrt(numValues);
   grandMean := grandMean / numValues;
   grandSigma := grandSigma / numGrps;
-  D3Value := D3[grpSize-1];
-  D4Value := D4[grpSize-1];
-  C4 := sqrt(2.0 / (grpSize - 1));
-  gamma := exp(GammaLn(grpSize / 2.0));
-  C4 := C4 * gamma;
-  gamma := exp(GammaLn((grpSize-1) / 2.0));
-  C4 := C4 / gamma;
-  B := grandSigma * sqrt(1.0 - (C4 * C4)) / C4;
+
+  C4 := CalcC4(grpSize);
+  B := grandSigma * sqrt(1.0 - sqr(C4)) / C4;
   UCL := grandSigma + 3.0 * B;
   LCL := grandSigma - 3.0 * B;
   if (LCL < 0.0) then LCL := 0.0;
@@ -135,6 +138,9 @@ begin
     lReport.Add('Sigma Chart Results');
     lReport.Add('');
     lReport.Add('Number of values:       %8d',   [numValues]);
+    lReport.Add('Number of groups:       %8d',   [numGrps]);
+    lReport.Add('Group size:             %8d',   [grpSize]);
+    lReport.Add('');
     lReport.Add('Grand Mean:             %8.3f', [grandMean]);
     lReport.Add('Standard Deviation:     %8.3f', [grandSD]);
     lReport.Add('Standard Error of Mean: %8.3f', [seMean]);
@@ -157,7 +163,7 @@ begin
   PlotMeans(
     Format('Sigma Chart for "%s"', [GetFileName]),
     GroupEdit.Text, Format('StdDev(%s)', [MeasEdit.Text]),
-    'Group Sigma', 'Mean',
+    'Group Sigma', 'Avg',
     groups, stdDev,
     UCL, LCL, grandSigma,
     NaN, NaN, NaN
diff --git a/applications/lazstats/source/units/mathunit.pas b/applications/lazstats/source/units/mathunit.pas
index a0312a4d2..a22fb2c79 100644
--- a/applications/lazstats/source/units/mathunit.pas
+++ b/applications/lazstats/source/units/mathunit.pas
@@ -25,6 +25,8 @@ function FDensity(x: Double; DF1, DF2: Integer): Double;
 
 function Chi2Density(x: Double; N: Integer): Double;
 
+function CalcC4(n: Integer): Double;
+
 implementation
 
 uses
@@ -286,5 +288,24 @@ begin
   Result := power(x, (N-2.0) * 0.5) / (factor * exp(x * 0.5));
 end;
 
+
+// Returns the value of the C4 factor needed for correction of standard deviations
+//   C4 = sqrt(2 / (n-1)) (n/2-1)! / ((n-1)/2-1)!
+function CalcC4(n: Integer): Double;
+var
+  gamma1, gamma2: Double;
+begin
+  gamma1 := GammaLn(n/2);
+  gamma2 := GammaLn((n-1)/2);
+  Result := sqrt(2.0 / (n-1)) * exp(gamma1 - gamma2);
+  {
+  C4 := sqrt(2.0 / (grpSize - 1));
+  gamma := exp(GammaLn(grpSize / 2.0));
+  C4 := C4 * gamma;
+  gamma := exp(GammaLn((grpSize-1) / 2.0));
+  C4 := C4 / gamma;
+  }
+end;
+
 end.