// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

// Pin for the negative binomial's far regime: when the summation's seed
// p^r underflows to an exact zero the recurrence can no longer start,
// and the CDF answers through the regularised beta identity instead.
// The law with r = 100 and p = 1e-4 has its mean at 999 900, so the CDF
// at the mean is near one half; the seed 1e-400 underflows outright, and
// the summation that ignored that answered an exact zero at every k.

package stats

import (
	"math"
	"testing"
)

func TestNegativeBinomialFarTail(t *testing.T) {
	atMean, err := NegativeBinomialCDF(999900, 100, 1e-4)
	if err != nil {
		t.Fatalf("NegativeBinomialCDF(999900, 100, 1e-4): %v", err)
	}
	// The mean is r(1−p)/p = 999 900 and the spread is near 100 000, so
	// the true value at the mean sits within a few 1e-3 of one half
	// (exactly 0.5133 by the exact referent). The identity's accuracy in
	// this regime is bounded by the lgamma noise of its front factor, so
	// the pin holds a loose band and a strict ordering.
	if math.Abs(atMean-0.5133007914) > 1e-6 {
		t.Fatalf("NegativeBinomialCDF at the mean = %.17g, want 0.5133007914 to within 1e-6", atMean)
	}
	below, err := NegativeBinomialCDF(950000, 100, 1e-4)
	if err != nil {
		t.Fatalf("NegativeBinomialCDF(950000, 100, 1e-4): %v", err)
	}
	above, err := NegativeBinomialCDF(1050000, 100, 1e-4)
	if err != nil {
		t.Fatalf("NegativeBinomialCDF(1050000, 100, 1e-4): %v", err)
	}
	if !(below < atMean && atMean < above) {
		t.Fatalf("the CDF is not increasing through the mean: %.17g, %.17g, %.17g", below, atMean, above)
	}
	far, err := NegativeBinomialCDF(2000000, 100, 1e-4)
	if err != nil {
		t.Fatalf("NegativeBinomialCDF(2000000, 100, 1e-4): %v", err)
	}
	if far < 0.9999 {
		t.Fatalf("NegativeBinomialCDF two hundred standard deviations out = %.17g, want ~1", far)
	}
	// The quantile route runs on the same CDF: the median must sit near
	// the mean the law's own moments give, a fraction of the spread
	// below it on the skewed side. The exact referent puts the median at
	// 996 569, and the pin holds it within a thousandth of the spread.
	median, err := NegativeBinomialQuantile(0.5, 1e-4, 100)
	if err != nil {
		t.Fatalf("NegativeBinomialQuantile: %v", err)
	}
	if math.Abs(median-996569) > 100 {
		t.Fatalf("NegativeBinomialQuantile(0.5) = %v, want within 100 of the exact median 996569", median)
	}
}