feat: add population and sample standard deviation calculations

README.md

@@ -1,8 +1,8 @@
 # sampdo
 
-`sampdo` is a lightweight statistical aggregation library for Go. It collects a large number of `float64` samples and efficiently computes minimum, maximum, mean, median, and arbitrary percentiles from a sorted view. For non-negative values it uses radix sort, which is faster than a comparison sort.
+`sampdo` is a lightweight statistical aggregation library for Go. It collects a large number of `float64` samples and efficiently computes minimum, maximum, mean, median, population/sample standard deviation, and arbitrary percentiles from a sorted view. For non-negative values it uses radix sort, which is faster than a comparison sort.
 
-`sampdo` は Go 用の軽量な数値集計ライブラリです。大量の `float64` サンプルに対して、ソート済みの状態から最小値・最大値・平均値・中央値・任意のパーセンタイルを効率的に取得できます。非負数に対しては基数ソートを使い、比較ソートより高速に集計できます。
+`sampdo` は Go 用の軽量な数値集計ライブラリです。大量の `float64` サンプルに対して、ソート済みの状態から最小値・最大値・平均値・中央値・母標準偏差・標本標準偏差・任意のパーセンタイルを効率的に取得できます。非負数に対しては基数ソートを使い、比較ソートより高速に集計できます。
 
 ## Features / 特徴
 
@@ -107,9 +107,45 @@ P99    : 9.91
 | `Min()` | Smallest value; returns `0` when empty. | 最小値。空の場合は `0`。 |
 | `Max()` | Largest value; returns `0` when empty. | 最大値。空の場合は `0`。 |
 | `Mean()` | Arithmetic mean; sums in input order. Returns `0` when empty. | 算術平均。入力順に合計し、空の場合は `0`。 |
+| `PopulationStdDev()` | Population standard deviation, divisor n. | 母標準偏差。分母は n。 |
+| `SampleStdDev()` | Sample standard deviation, divisor n-1; requires two points. | 標本標準偏差。分母は n-1。2件以上が必要。 |
 | `Median()` | 50th percentile; averages the two middle values for an even count. | 50 パーセンタイル。要素数が偶数の場合は中央 2 値の平均。 |
 | `Percentile(p)` | Returns the `p`th percentile (0–100) using linear interpolation; returns `0` for out-of-range `p`. | `p` パーセンタイル（0〜100）を線形補間で返す。範囲外は `0`。 |
 
+### Standard deviation / 標準偏差
+
+```go
+population, err := sorted.PopulationStdDev()
+if err != nil {
+    log.Fatal(err)
+}
+sample, err := sorted.SampleStdDev()
+if err != nil {
+    log.Fatal(err)
+}
+fmt.Printf("Population: %.6f, Sample: %.6f\n", population, sample)
+// For samples 1 through 10: Population: 2.872281, Sample: 3.027650
+```
+
+`PopulationStdDev()` computes `sqrt(sum((x - mean)^2) / n)`;
+`SampleStdDev()` computes `sqrt(sum((x - mean)^2) / (n - 1))`.
+Both return `(float64, error)` in the same units as the input, using O(n) time
+and O(1) extra space per call. Calculations center and scale the values to reduce
+rounding errors and avoid intermediate overflow/underflow.
+
+`PopulationStdDev()` は母標準偏差（分母 `n`）、`SampleStdDev()` は標本標準偏差
+（分母 `n-1`、不偏分散の平方根）を返します。戻り値は `(float64, error)` で、
+標準偏差の単位は入力と同じです。呼び出しごとに O(n) 時間・O(1) 追加メモリで計算します。
+値の中心を移動してスケールを調整し、丸め誤差や中間計算のオーバーフロー・アンダーフローを抑えます。
+
+- Empty views return `0, ErrNoPoints` for both methods.
+  空の `Sorted` は、どちらも `0, ErrNoPoints` を返します。
+- A single finite point has population standard deviation `0`; sample standard
+  deviation returns `0, ErrInsufficientPoints` for any single point.
+  有限値1件の母標準偏差は `0`。標本標準偏差は値によらず1件では `0, ErrInsufficientPoints` です。
+- With sufficient points, any infinity yields `NaN, nil`.
+  必要な件数を満たしていて無限大を含む場合は `NaN, nil` を返します。
+
 ### How percentiles are calculated / パーセンタイルの計算方法
 
 `Percentile(p)` computes the position `k` with the formula below and linearly interpolates between the two surrounding points.

sorted.go

@@ -8,6 +8,9 @@ import (
 var ErrNoPoints = errors.New("no points available")
 var ErrInvalidPercentile = errors.New("invalid percentile")
 
+// ErrInsufficientPoints indicates that sample standard deviation needs at least two points.
+var ErrInsufficientPoints = errors.New("at least two points required")
+
 type Sorted struct {
 	points   []float64
 	sum      float64
@@ -46,6 +49,59 @@ func (s *Sorted) Mean() (float64, error) {
 	return sum / float64(len(s.points)), nil
 }
 
+// PopulationStdDev returns the population standard deviation (divisor n).
+// Empty input returns 0 and ErrNoPoints. Input containing infinity returns NaN.
+func (s *Sorted) PopulationStdDev() (float64, error) {
+	return s.stdDev(0)
+}
+
+// SampleStdDev returns the sample standard deviation (divisor n-1).
+// Empty input returns 0 and ErrNoPoints; one point returns 0 and
+// ErrInsufficientPoints. Otherwise, input containing infinity returns NaN.
+func (s *Sorted) SampleStdDev() (float64, error) {
+	return s.stdDev(1)
+}
+
+func (s *Sorted) stdDev(correction int) (float64, error) {
+	n := len(s.points)
+	if n == 0 {
+		return 0, ErrNoPoints
+	}
+	if n <= correction {
+		return 0, ErrInsufficientPoints
+	}
+	low, high := s.points[0], s.points[n-1]
+	if math.IsInf(low, 0) || math.IsInf(high, 0) {
+		return math.NaN(), nil
+	}
+	if low == high {
+		return 0, nil
+	}
+	// Center before scaling to preserve small differences on a large offset.
+	// Halving the endpoints first avoids overflow for extreme finite inputs.
+	center := low/2 + high/2
+	scale := math.Max(math.Abs(low-center), math.Abs(high-center))
+	// Compensate the sum of normalized offsets; the cached input-order sum
+	// can overflow or lose the very differences needed by this calculation.
+	sum, compensation := 0.0, 0.0
+	for _, v := range s.points {
+		y := (v-center)/scale - compensation
+		next := sum + y
+		compensation = (next - sum) - y
+		sum = next
+	}
+	mean := sum / float64(n)
+	squares, compensation := 0.0, 0.0
+	for _, v := range s.points {
+		d := (v-center)/scale - mean
+		y := d*d - compensation
+		next := squares + y
+		compensation = (next - squares) - y
+		squares = next
+	}
+	return scale * math.Sqrt(squares/float64(n-correction)), nil
+}
+
 func (s *Sorted) Median() (float64, error) {
 	if len(s.points) == 0 {
 		return 0, ErrNoPoints

sorted_test.go

@@ -38,6 +38,97 @@ func TestSortedBasic(t *testing.T) {
 
 }
 
+func TestSortedStdDev(t *testing.T) {
+	tests := []struct {
+		name               string
+		points             []float64
+		population, sample float64
+	}{
+		{"textbook", []float64{2, 4, 4, 4, 5, 5, 7, 9}, 2, math.Sqrt(32.0 / 7)},
+		{"negative", []float64{-5, -3, -1}, math.Sqrt(8.0 / 3), 2},
+		{"mixed", []float64{-1, 0, 1}, math.Sqrt(2.0 / 3), 1},
+		{"constant", []float64{42, 42, 42}, 0, 0},
+		{"zeros", []float64{math.Copysign(0, -1), 0}, 0, 0},
+		{"large offset", []float64{1e16, 1e16 + 2, 1e16 + 4}, math.Sqrt(8.0 / 3), 2},
+		{"large magnitude", []float64{-1e200, 1e200}, 1e200, math.Sqrt2 * 1e200},
+		{"small magnitude", []float64{-1e-200, 1e-200}, 1e-200, math.Sqrt2 * 1e-200},
+		{"overflowing sum", []float64{math.MaxFloat64, math.MaxFloat64}, 0, 0},
+		{"extreme range", []float64{-math.MaxFloat64, math.MaxFloat64}, math.MaxFloat64, math.Inf(1)},
+		{"subnormal", []float64{0, 2 * math.SmallestNonzeroFloat64}, math.SmallestNonzeroFloat64, math.SmallestNonzeroFloat64},
+		{"positive infinity", []float64{1, math.Inf(1)}, math.NaN(), math.NaN()},
+		{"negative infinity", []float64{math.Inf(-1), 1}, math.NaN(), math.NaN()},
+		{"both infinities", []float64{math.Inf(-1), math.Inf(1)}, math.NaN(), math.NaN()},
+		{"equal infinities", []float64{math.Inf(1), math.Inf(1)}, math.NaN(), math.NaN()},
+	}
+	for _, tt := range tests {
+		t.Run(tt.name, func(t *testing.T) {
+			data := New()
+			require.NoError(t, data.Append(tt.points...))
+			sorted, err := data.Sorted()
+			require.NoError(t, err)
+			for _, view := range []*Sorted{sorted, {points: slices.Clone(sorted.points)}} {
+				population, err := view.PopulationStdDev()
+				require.NoError(t, err)
+				sample, err := view.SampleStdDev()
+				require.NoError(t, err)
+				for i, got := range []float64{population, sample} {
+					want := []float64{tt.population, tt.sample}[i]
+					switch {
+					case math.IsNaN(want):
+						require.True(t, math.IsNaN(got))
+					case want == 0 || math.IsInf(want, 0) || math.Abs(want) < 1e-300:
+						require.Equal(t, want, got)
+					default:
+						require.InEpsilon(t, want, got, 1e-14)
+					}
+				}
+			}
+		})
+	}
+}
+
+func TestSortedStdDevSmallInputs(t *testing.T) {
+	s := &Sorted{}
+	got, err := s.PopulationStdDev()
+	require.ErrorIs(t, err, ErrNoPoints)
+	require.Zero(t, got)
+	got, err = s.SampleStdDev()
+	require.ErrorIs(t, err, ErrNoPoints)
+	require.Zero(t, got)
+	for _, point := range []float64{42, math.MaxFloat64, math.Inf(1)} {
+		s = &Sorted{points: []float64{point}}
+		got, err = s.PopulationStdDev()
+		require.NoError(t, err)
+		if math.IsInf(point, 0) {
+			require.True(t, math.IsNaN(got))
+		} else {
+			require.Zero(t, got)
+		}
+		got, err = s.SampleStdDev()
+		require.ErrorIs(t, err, ErrInsufficientPoints)
+		require.Zero(t, got)
+	}
+}
+
+func TestSortedStdDevFrozenView(t *testing.T) {
+	data := New()
+	require.NoError(t, data.Append(3, 1))
+	before, err := data.Sorted()
+	require.NoError(t, err)
+	require.NoError(t, data.Append(5))
+	after, err := data.Sorted()
+	require.NoError(t, err)
+	got, err := before.PopulationStdDev()
+	require.NoError(t, err)
+	require.Equal(t, 1.0, got)
+	got, err = before.SampleStdDev()
+	require.NoError(t, err)
+	require.Equal(t, math.Sqrt2, got)
+	got, err = after.SampleStdDev()
+	require.NoError(t, err)
+	require.Equal(t, 2.0, got)
+}
+
 func TestSortedMedian(t *testing.T) {
 	tests := []struct {
 		name string