## Coverage Report

#### Created: 2026-01-25 15:05

next uncovered line (L), next uncovered region (R), next uncovered branch (B)

    /home/noah/src/realizar/src/bench/statistics.rs

    Line
    Count
    Source
    1
    //! Statistical measurement and analysis for benchmarking
    2
    //!
    3
    //! Extracted from bench/mod.rs (PMAT-802) to reduce module size.
    4
    //! Contains:
    5
    //! - BENCH-004: MeasurementProtocol
    6
    //! - BENCH-005: LatencyStatistics
    7
    //! - BENCH-006: Outlier Detection (MAD-based)
    8
    //! - BENCH-007: Regression Detection
    9
    //! - BENCH-008: Welch's t-test for Statistical Significance
    10
    11
    #![allow(clippy::cast_precision_loss)]
    12
    13
    use std::time::Duration;
    14
    15
    // ============================================================================
    16
    // BENCH-004: MeasurementProtocol (following SPEC-BENCH-001)
    17
    // ============================================================================
    18
    19
    /// Complete measurement protocol for benchmarking
    20
    ///
    21
    /// Follows MLPerf™ Inference benchmarking principles for scientific rigor.
    22
    #[derive(Debug, Clone)]
    23
    pub struct MeasurementProtocol {
    24
        /// Number of latency samples to collect
    25
        pub latency_samples: usize,
    26
        /// Percentiles to compute (e.g., 50, 90, 95, 99, 99.9)
    27
        pub latency_percentiles: Vec<f64>,
    28
        /// Duration for throughput measurement
    29
        pub throughput_duration: Duration,
    30
        /// Ramp-up time before throughput measurement
    31
        pub throughput_ramp_up: Duration,
    32
        /// Number of memory samples to collect
    33
        pub memory_samples: usize,
    34
        /// Interval between memory samples
    35
        pub memory_interval: Duration,
    36
    }
    37
    38
    impl Default for MeasurementProtocol {
    39
    2
        fn default() -> Self {
    40
    2
            Self {
    41
    2
                latency_samples: 100,
    42
    2
                latency_percentiles: vec![50.0, 90.0, 95.0, 99.0, 99.9],
    43
    2
                throughput_duration: Duration::from_secs(60),
    44
    2
                throughput_ramp_up: Duration::from_secs(10),
    45
    2
                memory_samples: 10,
    46
    2
                memory_interval: Duration::from_secs(1),
    47
    2
            }
    48
    2
        }
    49
    }
    50
    51
    impl MeasurementProtocol {
    52
        /// Create a new measurement protocol with default values
    53
        #[must_use]
    54
    1
        pub fn new() -> Self {
    55
    1
            Self::default()
    56
    1
        }
    57
    58
        /// Set the number of latency samples
    59
        #[must_use]
    60
    1
        pub fn with_latency_samples(mut self, samples: usize) -> Self {
    61
    1
            self.latency_samples = samples;
    62
    1
            self
    63
    1
        }
    64
    65
        /// Set the percentiles to compute
    66
        #[must_use]
    67
    1
        pub fn with_percentiles(mut self, percentiles: Vec<f64>) -> Self {
    68
    1
            self.latency_percentiles = percentiles;
    69
    1
            self
    70
    1
        }
    71
    72
        /// Set the throughput measurement duration
    73
        #[must_use]
    74
    1
        pub fn with_throughput_duration(mut self, duration: Duration) -> Self {
    75
    1
            self.throughput_duration = duration;
    76
    1
            self
    77
    1
        }
    78
    79
        /// Set the number of memory samples
    80
        #[must_use]
    81
    1
        pub fn with_memory_samples(mut self, samples: usize) -> Self {
    82
    1
            self.memory_samples = samples;
    83
    1
            self
    84
    1
        }
    85
    }
    86
    87
    // ============================================================================
    88
    // BENCH-005: LatencyStatistics (following SPEC-BENCH-001 Section 7.1)
    89
    // ============================================================================
    90
    91
    /// Comprehensive latency statistics following MLPerf™ reporting standards
    92
    #[derive(Debug, Clone)]
    93
    pub struct LatencyStatistics {
    94
        /// Mean latency
    95
        pub mean: Duration,
    96
        /// Standard deviation
    97
        pub std_dev: Duration,
    98
        /// Minimum latency
    99
        pub min: Duration,
    100
        /// Maximum latency
    101
        pub max: Duration,
    102
        /// 50th percentile (median)
    103
        pub p50: Duration,
    104
        /// 90th percentile
    105
        pub p90: Duration,
    106
        /// 95th percentile
    107
        pub p95: Duration,
    108
        /// 99th percentile
    109
        pub p99: Duration,
    110
        /// 99.9th percentile (tail latency)
    111
        pub p999: Duration,
    112
        /// Number of samples
    113
        pub samples: usize,
    114
        /// 95% confidence interval (lower, upper)
    115
        pub confidence_interval_95: (Duration, Duration),
    116
    }
    117
    118
    impl LatencyStatistics {
    119
        /// Compute statistics from a slice of duration samples
    120
        ///
    121
        /// # Panics
    122
        /// Panics if samples is empty
    123
        #[must_use]
    124
    4
        pub fn from_samples(samples: &[Duration]) -> Self {
    125
    4
            assert!(!samples.is_empty(), "samples must not be empty"0);
    126
    127
    4
            let n = samples.len();
    128
    4
            let n_f64 = n as f64;
    129
    130
            // Compute mean
    131
    4
            let sum_nanos: u128 = samples.iter().map(Duration::as_nanos).sum();
    132
    4
            let mean_nanos = sum_nanos / n as u128;
    133
    4
            let mean = Duration::from_nanos(mean_nanos as u64);
    134
    135
            // Compute standard deviation
    136
    4
            let variance: f64 = samples
    137
    4
                .iter()
    138
    215
                .map4(|s| {
    139
    215
                    let diff = s.as_nanos() as f64 - mean_nanos as f64;
    140
    215
                    diff * diff
    141
    215
                })
    142
    4
                .sum::<f64>()
    143
    4
                / (n_f64 - 1.0).max(1.0);
    144
    4
            let std_dev_nanos = variance.sqrt();
    145
    4
            let std_dev = Duration::from_nanos(std_dev_nanos as u64);
    146
    147
            // Sort for percentile computation
    148
    4
            let mut sorted: Vec<Duration> = samples.to_vec();
    149
    4
            sorted.sort();
    150
    151
            // Min/max
    152
    4
            let min = sorted[0];
    153
    4
            let max = sorted[n - 1];
    154
    155
            // Percentiles using nearest-rank method
    156
    20
            let percentile4 = |p: f64| -> Duration {
    157
    20
                let idx = ((p / 100.0) * n_f64).ceil() as usize;
    158
    20
                sorted[idx.saturating_sub(1).min(n - 1)]
    159
    20
            };
    160
    161
    4
            let p50 = percentile(50.0);
    162
    4
            let p90 = percentile(90.0);
    163
    4
            let p95 = percentile(95.0);
    164
    4
            let p99 = percentile(99.0);
    165
    4
            let p999 = percentile(99.9);
    166
    167
            // 95% confidence interval using t-distribution approximation
    168
            // For large n, t ≈ 1.96
    169
    4
            let t_value = if n >= 30 { 1.962 } else { 2.0 + 4.0 / n_f642 };
    170
    4
            let margin = std_dev_nanos * t_value / n_f64.sqrt();
    171
    4
            let lower = Duration::from_nanos((mean_nanos as f64 - margin).max(0.0) as u64);
    172
    4
            let upper = Duration::from_nanos((mean_nanos as f64 + margin) as u64);
    173
    174
    4
            Self {
    175
    4
                mean,
    176
    4
                std_dev,
    177
    4
                min,
    178
    4
                max,
    179
    4
                p50,
    180
    4
                p90,
    181
    4
                p95,
    182
    4
                p99,
    183
    4
                p999,
    184
    4
                samples: n,
    185
    4
                confidence_interval_95: (lower, upper),
    186
    4
            }
    187
    4
        }
    188
    }
    189
    190
    // ============================================================================
    191
    // BENCH-006: Outlier Detection (MAD-based)
    192
    // ============================================================================
    193
    194
    /// Detect outliers using Median Absolute Deviation (MAD) method
    195
    ///
    196
    /// More robust than standard deviation for non-normal distributions.
    197
    /// Uses the modified Z-score method with configurable threshold.
    198
    ///
    199
    /// # Arguments
    200
    /// * `samples` - Slice of f64 samples
    201
    /// * `threshold` - Modified Z-score threshold (typically 3.5 for strict, 2.0 for lenient)
    202
    ///
    203
    /// # Returns
    204
    /// Vector of indices that are considered outliers
    205
    5
    pub fn detect_outliers(samples: &[f64], threshold: f64) -> Vec<usize> {
    206
    5
        if samples.len() < 3 {
    207
    0
            return Vec::new();
    208
    5
        }
    209
    210
        // Calculate median
    211
    5
        let mut sorted = samples.to_vec();
    212
    102
        sorted5.sort_by5(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    213
    5
        let median = if sorted.len().is_multiple_of(2) {
    214
    3
            f64::midpoint(sorted[sorted.len() / 2 - 1], sorted[sorted.len() / 2])
    215
        } else {
    216
    2
            sorted[sorted.len() / 2]
    217
        };
    218
    219
        // Calculate MAD (Median Absolute Deviation)
    220
    42
        let mut deviations5: Vec<f64>5 = samples5.iter5().map5(|x| (x - median).abs()).collect5();
    221
    116
        deviations5.sort_by5(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    222
    5
        let mad = if deviations.len().is_multiple_of(2) {
    223
    3
            f64::midpoint(
    224
    3
                deviations[deviations.len() / 2 - 1],
    225
    3
                deviations[deviations.len() / 2],
    226
            )
    227
        } else {
    228
    2
            deviations[deviations.len() / 2]
    229
        };
    230
    231
        // Avoid division by zero
    232
    5
        if mad < f64::EPSILON {
    233
    0
            return Vec::new();
    234
    5
        }
    235
    236
        // Constant for normal distribution approximation
    237
    5
        let k = 1.4826;
    238
    239
        // Find outliers using modified Z-score
    240
    5
        samples
    241
    5
            .iter()
    242
    5
            .enumerate()
    243
    42
            .filter5(|(_, &x)| {
    244
    42
                let modified_z = (x - median) / (k * mad);
    245
    42
                modified_z.abs() > threshold
    246
    42
            })
    247
    5
            .map(|(i, _)| i)
    248
    5
            .collect()
    249
    5
    }
    250
    251
    // ============================================================================
    252
    // BENCH-007: Regression Detection
    253
    // ============================================================================
    254
    255
    /// Single benchmark metric for comparison
    256
    #[derive(Debug, Clone)]
    257
    pub struct BenchmarkMetrics {
    258
        /// Metric name
    259
        pub name: String,
    260
        /// Mean value
    261
        pub mean: f64,
    262
        /// Standard deviation
    263
        pub std_dev: f64,
    264
        /// Number of samples
    265
        pub samples: usize,
    266
    }
    267
    268
    /// Individual regression item
    269
    #[derive(Debug, Clone)]
    270
    pub struct Regression {
    271
        /// Metric that regressed
    272
        pub metric: String,
    273
        /// Baseline value
    274
        pub baseline: f64,
    275
        /// Current value
    276
        pub current: f64,
    277
        /// Percentage change
    278
        pub change_percent: f64,
    279
    }
    280
    281
    /// Report from regression analysis
    282
    #[derive(Debug, Clone)]
    283
    pub struct RegressionReport {
    284
        /// Metrics that exceeded failure threshold
    285
        pub regressions: Vec<Regression>,
    286
        /// Metrics that exceeded warning threshold
    287
        pub warnings: Vec<Regression>,
    288
        /// Metrics that improved significantly
    289
        pub improvements: Vec<Regression>,
    290
        /// Overall pass/fail (no regressions)
    291
        pub passed: bool,
    292
    }
    293
    294
    /// Performance regression detector
    295
    ///
    296
    /// Compares baseline and current benchmark results to detect
    297
    /// performance regressions, warnings, and improvements.
    298
    #[derive(Debug, Clone)]
    299
    pub struct RegressionDetector {
    300
        /// Warning threshold (default: 2%)
    301
        pub warning_threshold: f64,
    302
        /// Failure threshold (default: 5%)
    303
        pub failure_threshold: f64,
    304
    }
    305
    306
    impl Default for RegressionDetector {
    307
    5
        fn default() -> Self {
    308
    5
            Self {
    309
    5
                warning_threshold: 0.02, // 2%
    310
    5
                failure_threshold: 0.05, // 5%
    311
    5
            }
    312
    5
        }
    313
    }
    314
    315
    impl RegressionDetector {
    316
        /// Compare baseline and current metrics
    317
    4
        pub fn compare(
    318
    4
            &self,
    319
    4
            baseline: &BenchmarkMetrics,
    320
    4
            current: &BenchmarkMetrics,
    321
    4
        ) -> RegressionReport {
    322
    4
            let mut regressions = Vec::new();
    323
    4
            let mut warnings = Vec::new();
    324
    4
            let mut improvements = Vec::new();
    325
    326
            // Calculate percentage change (positive = regression for latency-like metrics)
    327
    4
            let change = (current.mean - baseline.mean) / baseline.mean;
    328
    329
    4
            let item = Regression {
    330
    4
                metric: baseline.name.clone(),
    331
    4
                baseline: baseline.mean,
    332
    4
                current: current.mean,
    333
    4
                change_percent: change * 100.0,
    334
    4
            };
    335
    336
    4
            if change > self.failure_threshold {
    337
    1
                regressions.push(item);
    338
    3
            } else if change > self.warning_threshold {
    339
    1
                warnings.push(item);
    340
    2
            } else if change < -self.warning_threshold {
    341
    1
                improvements.push(item);
    342
    1
            }
    343
    344
    4
            RegressionReport {
    345
    4
                passed: regressions.is_empty(),
    346
    4
                regressions,
    347
    4
                warnings,
    348
    4
                improvements,
    349
    4
            }
    350
    4
        }
    351
    }
    352
    353
    // ============================================================================
    354
    // BENCH-008: Welch's t-test for Statistical Significance
    355
    // Per Hoefler & Belli [17], statistical testing is required for valid comparisons
    356
    // ============================================================================
    357
    358
    /// Result of Welch's t-test for comparing two sample means
    359
    #[derive(Debug, Clone)]
    360
    pub struct WelchTTestResult {
    361
        /// Calculated t-statistic
    362
        pub t_statistic: f64,
    363
        /// Welch-Satterthwaite degrees of freedom
    364
        pub degrees_of_freedom: f64,
    365
        /// Two-tailed p-value
    366
        pub p_value: f64,
    367
        /// Whether the difference is statistically significant at given alpha
    368
        pub significant: bool,
    369
    }
    370
    371
    /// Perform Welch's t-test to compare two sample means
    372
    ///
    373
    /// Welch's t-test is used when samples may have unequal variances.
    374
    /// Returns statistical significance information.
    375
    ///
    376
    /// # Arguments
    377
    /// * `sample_a` - First sample
    378
    /// * `sample_b` - Second sample
    379
    /// * `alpha` - Significance level (e.g., 0.05 for 95% confidence)
    380
    ///
    381
    /// # Example
    382
    /// ```
    383
    /// use realizar::bench::welch_t_test;
    384
    ///
    385
    /// let a = vec![10.0, 11.0, 10.5, 10.2, 10.8];
    386
    /// let b = vec![20.0, 21.0, 20.5, 20.2, 20.8];
    387
    /// let result = welch_t_test(&a, &b, 0.05);
    388
    /// assert!(result.significant); // Clearly different means
    389
    /// ```
    390
    7
    pub fn welch_t_test(sample_a: &[f64], sample_b: &[f64], alpha: f64) -> WelchTTestResult {
    391
    7
        let n1 = sample_a.len() as f64;
    392
    7
        let n2 = sample_b.len() as f64;
    393
    394
        // Calculate means
    395
    7
        let mean1 = sample_a.iter().sum::<f64>() / n1;
    396
    7
        let mean2 = sample_b.iter().sum::<f64>() / n2;
    397
    398
        // Calculate sample variances (using n-1 for unbiased estimator)
    399
    7
        let var1 = if n1 > 1.0 {
    400
    36
            sample_a7.iter7().map7(|x| (x - mean1).powi(2)).sum7::<f64>() / (n1 - 1.0)7
    401
        } else {
    402
    0
            0.0
    403
        };
    404
    7
        let var2 = if n2 > 1.0 {
    405
    36
            sample_b7.iter7().map7(|x| (x - mean2).powi(2)).sum7::<f64>() / (n2 - 1.0)7
    406
        } else {
    407
    0
            0.0
    408
        };
    409
    410
        // Handle zero variance case
    411
    7
        let se1 = var1 / n1;
    412
    7
        let se2 = var2 / n2;
    413
    7
        let se_diff = (se1 + se2).sqrt();
    414
    415
    7
        if se_diff < f64::EPSILON {
    416
            // Both samples have zero variance - cannot compute t-statistic
    417
    1
            return WelchTTestResult {
    418
    1
                t_statistic: 0.0,
    419
    1
                degrees_of_freedom: n1 + n2 - 2.0,
    420
    1
                p_value: 1.0,
    421
    1
                significant: false,
    422
    1
            };
    423
    6
        }
    424
    425
        // Calculate t-statistic
    426
    6
        let t_stat = (mean1 - mean2) / se_diff;
    427
    428
        // Welch-Satterthwaite degrees of freedom
    429
    6
        let df_num = (se1 + se2).powi(2);
    430
    6
        let df_denom = if n1 > 1.0 && se1 > f64::EPSILON {
    431
    6
            se1.powi(2) / (n1 - 1.0)
    432
        } else {
    433
    0
            0.0
    434
    6
        } + if n2 > 1.0 && se2 > f64::EPSILON {
    435
    6
            se2.powi(2) / (n2 - 1.0)
    436
        } else {
    437
    0
            0.0
    438
        };
    439
    440
    6
        let df = if df_denom > f64::EPSILON {
    441
    6
            df_num / df_denom
    442
        } else {
    443
    0
            n1 + n2 - 2.0
    444
        };
    445
    446
        // Approximate p-value using normal distribution for large df
    447
        // For small df, we use a more conservative approximation
    448
    6
        let p_value = approximate_t_pvalue(t_stat.abs(), df);
    449
    450
    6
        WelchTTestResult {
    451
    6
            t_statistic: t_stat,
    452
    6
            degrees_of_freedom: df,
    453
    6
            p_value,
    454
    6
            significant: p_value < alpha,
    455
    6
        }
    456
    7
    }
    457
    458
    /// Approximate two-tailed p-value from t-distribution
    459
    ///
    460
    /// Uses normal approximation for large df, conservative approximation for small df
    461
    6
    fn approximate_t_pvalue(t_abs: f64, df: f64) -> f64 {
    462
        // For very large df, use normal approximation
    463
    6
        if df > 100.0 {
    464
            // Use error function approximation for normal CDF
    465
    0
            let z = t_abs;
    466
    0
            let p = erfc_approx(z / std::f64::consts::SQRT_2);
    467
    0
            return p;
    468
    6
        }
    469
    470
        // For smaller df, use a polynomial approximation of t-distribution CDF
    471
        // Based on Abramowitz and Stegun approximation
    472
    6
        let ratio = df / (df + t_abs * t_abs);
    473
    6
        incomplete_beta_approx(ratio, df / 2.0, 0.5)
    474
    6
    }
    475
    476
    /// Approximate complementary error function
    477
    0
    fn erfc_approx(x: f64) -> f64 {
    478
        // Horner form coefficients for erfc approximation
    479
        // From Abramowitz and Stegun, formula 7.1.26
    480
    0
        let a1 = 0.254_829_592;
    481
    0
        let a2 = -0.284_496_736;
    482
    0
        let a3 = 1.421_413_741;
    483
    0
        let a4 = -1.453_152_027;
    484
    0
        let a5 = 1.061_405_429;
    485
    0
        let p = 0.327_591_1;
    486
    487
    0
        let sign = if x < 0.0 { -1.0 } else { 1.0 };
    488
    0
        let x = x.abs();
    489
    490
    0
        let t = 1.0 / (1.0 + p * x);
    491
    0
        let y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * (-x * x).exp();
    492
    493
    0
        if sign < 0.0 {
    494
    0
            2.0 - y
    495
        } else {
    496
    0
            y
    497
        }
    498
    0
    }
    499
    500
    /// Approximate incomplete beta function (simplified for t-test)
    501
    7
    fn incomplete_beta_approx(x: f64, a: f64, b: f64) -> f64 {
    502
        // Use continued fraction expansion for better accuracy
    503
        // Simplified approximation suitable for t-distribution p-values
    504
    7
        if x < (a + 1.0) / (a + b + 2.0) {
    505
    6
            let beta_factor =
    506
    6
                gamma_ln(a + b) - gamma_ln(a) - gamma_ln(b) + a * x.ln() + b * (1.0 - x).ln();
    507
    6
            let beta_factor = beta_factor.exp();
    508
    6
            beta_factor * cf_beta(x, a, b) / a
    509
        } else {
    510
    1
            1.0 - incomplete_beta_approx(1.0 - x, b, a)
    511
        }
    512
    7
    }
    513
    514
    /// Continued fraction for incomplete beta
    515
    #[allow(clippy::many_single_char_names)] // Standard math notation for beta function
    516
    6
    fn cf_beta(x: f64, a: f64, b: f64) -> f64 {
    517
    6
        let max_iter = 100;
    518
    6
        let eps = 1e-10;
    519
    6
        let tiny = 1e-30;
    520
    521
    6
        let qab = a + b;
    522
    6
        let qap = a + 1.0;
    523
    6
        let qam = a - 1.0;
    524
    525
    6
        let mut c = 1.0;
    526
    6
        let mut d = 1.0 - qab * x / qap;
    527
    6
        if d.abs() < tiny {
    528
    0
            d = tiny;
    529
    6
        }
    530
    6
        d = 1.0 / d;
    531
    6
        let mut h = d;
    532
    533
    19
        for m in 1..=max_iter6 {
    534
    19
            let m_f = m as f64;
    535
    19
            let m2 = 2.0 * m_f;
    536
    537
            // Even step
    538
    19
            let aa = m_f * (b - m_f) * x / ((qam + m2) * (a + m2));
    539
    19
            d = 1.0 + aa * d;
    540
    19
            if d.abs() < tiny {
    541
    0
                d = tiny;
    542
    19
            }
    543
    19
            c = 1.0 + aa / c;
    544
    19
            if c.abs() < tiny {
    545
    0
                c = tiny;
    546
    19
            }
    547
    19
            d = 1.0 / d;
    548
    19
            h *= d * c;
    549
    550
            // Odd step
    551
    19
            let aa = -(a + m_f) * (qab + m_f) * x / ((a + m2) * (qap + m2));
    552
    19
            d = 1.0 + aa * d;
    553
    19
            if d.abs() < tiny {
    554
    0
                d = tiny;
    555
    19
            }
    556
    19
            c = 1.0 + aa / c;
    557
    19
            if c.abs() < tiny {
    558
    0
                c = tiny;
    559
    19
            }
    560
    19
            d = 1.0 / d;
    561
    19
            let del = d * c;
    562
    19
            h *= del;
    563
    564
    19
            if (del - 1.0).abs() < eps {
    565
    6
                break;
    566
    13
            }
    567
        }
    568
    569
    6
        h
    570
    6
    }
    571
    572
    /// Approximate log-gamma function (Stirling's approximation)
    573
    #[allow(clippy::excessive_precision)] // Lanczos coefficients require high precision
    574
    18
    fn gamma_ln(x: f64) -> f64 {
    575
    18
        if x <= 0.0 {
    576
    0
            return f64::INFINITY;
    577
    18
        }
    578
    579
        // Lanczos approximation coefficients
    580
    18
        let g = 7.0;
    581
    18
        let c = [
    582
    18
            0.999_999_999_999_81,
    583
    18
            676.520_368_121_885,
    584
    18
            -1_259.139_216_722_403,
    585
    18
            771.323_428_777_653,
    586
    18
            -176.615_029_162_141,
    587
    18
            12.507_343_278_687,
    588
    18
            -0.138_571_095_265_72,
    589
    18
            9.984_369_578_02e-6,
    590
    18
            1.505_632_735_15e-7,
    591
    18
        ];
    592
    593
    18
        let x = x - 1.0;
    594
    18
        let mut sum = c[0];
    595
    144
        for (i, &coef) in c18.iter18().enumerate18().skip18(1) {
    596
    144
            sum += coef / (x + i as f64);
    597
    144
        }
    598
    599
    18
        let t = x + g + 0.5;
    600
    18
        0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + sum.ln()
    601
    18
    }