gaussian_test.go raw

   1  package gnarlring
   2  
   3  import (
   4  	"math"
   5  	"testing"
   6  )
   7  
   8  func TestGaussSampleZCDT(t *testing.T) {
   9  	sigma := math.Sqrt(float64(N)) * 2.0 // ≈ 10.39
  10  	gs := NewGaussSampler(sigma)
  11  
  12  	// Sample many values, check mean is near 0 and variance is near sigma^2.
  13  	const samples = 50000
  14  	var sum int64
  15  	var sumSq float64
  16  	for i := 0; i < samples; i++ {
  17  		z := gs.SampleZ(0)
  18  		sum += z
  19  		sumSq += float64(z) * float64(z)
  20  	}
  21  	mean := float64(sum) / float64(samples)
  22  	variance := sumSq / float64(samples)
  23  
  24  	// Standard error for mean: σ/√n ≈ 10.4/223 ≈ 0.05.
  25  	// 3σ interval: [-0.15, +0.15].
  26  	if mean < -0.15 || mean > 0.15 {
  27  		t.Fatalf("mean %f outside [-0.15, 0.15]", mean)
  28  	}
  29  
  30  	// Variance target: σ²/(2π) ≈ 108/(2π) ≈ 17.2 for discrete Gaussian.
  31  	// For 50k samples, std err ≈ 17.2×√(2/n) ≈ 17.2×√(2/50000) ≈ 0.11.
  32  	// 5σ interval: [17.2-0.55, 17.2+0.55] ≈ [16.65, 17.75].
  33  	if variance < 16.5 || variance > 18.0 {
  34  		t.Fatalf("variance %f outside [106, 110]", variance)
  35  	}
  36  }
  37  
  38  func TestGaussSampleZBounds(t *testing.T) {
  39  	sigma := math.Sqrt(float64(N)) * 2.0
  40  	gs := NewGaussSampler(sigma)
  41  
  42  	for i := 0; i < 10000; i++ {
  43  		z := gs.SampleZ(0)
  44  		// Should never exceed 13*sigma ≈ 135.
  45  		if z < -140 || z > 140 {
  46  			t.Fatalf("sample %d outside bounds", z)
  47  		}
  48  	}
  49  }
  50  
  51  func TestGaussSamplePoly(t *testing.T) {
  52  	sigma := math.Sqrt(float64(N)) * 2.0
  53  	gs := NewGaussSampler(sigma)
  54  
  55  	p := gs.SamplePoly()
  56  
  57  	// Should be non-zero with high probability.
  58  	if IsZero(p) {
  59  		t.Fatal("SamplePoly returned zero polynomial (astronomically unlikely)")
  60  	}
  61  
  62  	// All coefficients should be small (< 135).
  63  	for i, v := range p.Coeffs {
  64  		// Convert to centered.
  65  		half := uint16(Q / 2)
  66  		var absV uint16
  67  		if v > half {
  68  			absV = Q - v
  69  		} else {
  70  			absV = v
  71  		}
  72  		if absV > 140 {
  73  			t.Fatalf("coefficient %d too large: %d", i, absV)
  74  		}
  75  	}
  76  }
  77  
  78  func TestGaussSamplePolyReproducible(t *testing.T) {
  79  	sigma := math.Sqrt(float64(N)) * 2.0
  80  
  81  	// Two samplers with different RNG should produce different outputs.
  82  	gs1 := NewGaussSampler(sigma)
  83  	gs2 := NewGaussSampler(sigma)
  84  
  85  	p1 := gs1.SamplePoly()
  86  	p2 := gs2.SamplePoly()
  87  
  88  	if Equal(p1, p2) {
  89  		t.Log("two SamplePoly outputs identical (possible but unlikely — retry)")
  90  		// Not a hard failure — genuinely possible at 1/2^243 probability.
  91  	}
  92  }
  93  
  94  func TestGaussSampleZWithCenter(t *testing.T) {
  95  	sigma := math.Sqrt(float64(N)) * 2.0
  96  	gs := NewGaussSampler(sigma)
  97  
  98  	const samples = 20000
  99  	var sum int64
 100  	for i := 0; i < samples; i++ {
 101  		z := gs.SampleZ(5.0)
 102  		sum += z
 103  	}
 104  	mean := float64(sum) / float64(samples)
 105  
 106  	if mean < 4.5 || mean > 5.5 {
 107  		t.Fatalf("centered mean %f outside [4.5, 5.5]", mean)
 108  	}
 109  }
 110