After removal of all deterministic components via high-fidelity pulsar timing models, residual pulsar timing noise (dominated by spin noise and dispersion measure variations) constitutes a physically irreducible stochastic process — providing a universally accessible, publicly verifiable source of certified true randomness for cryptographic applications that requires only a radio telescope and an open pulsar ephemeris, with no proprietary hardware or institutional trust required.

Residual timing noise from millisecond pulsars (MSPs), after subtraction of a complete deterministic timing model (including spin-down, proper motion, parallax, binary orbital parameters, and dispersion measure polynomial fits), contains a stochastic component that passes all standard randomness tests (NIST SP 800-22, TestU01 BigCrush, Dieharder) at p > 0.01 for each test, exhibits no autocorrelation at lags > 1 epoch at r < 0.05 significance, and cannot be predicted beyond the noise floor by any known physical model — thereby qualifying as a cryptographically suitable true random number generator (TRNG) source. Specifically: for at least 5 MSPs with timing baselines > 5 years and residual RMS < 1 μs, the extracted bit sequences will be statistically indistinguishable from ideal uniform random sequences at α = 0.01 after von Neumann debiasing.

PHASE 1 — Data Acquisition and Baseline Characterization (Days 1–60): Select 10 MSPs from public archives (IPTA DR2, NANOGrav 15-year, PPTA DR3) with known timing solutions. Download all available TOA (time-of-arrival) data and noise budgets. Fit complete timing models using TEMPO2 and PINT independently.

PHASE 2 — Noise Decomposition (Days 30–90): Apply Bayesian noise analysis (enterprise/PTMCMCSampler) to decompose residuals into: (a) white noise (EFAC, EQUAD, ECORR), (b) red spin noise (power-law spectrum), (c) DM variations (chromatic red noise), (d) GWB contribution (Hellings-Downs correlated). Extract the white-noise-subtracted, GWB-subtracted residual as the candidate random source.

PHASE 3 — Randomness Testing (Days 60–120): Convert residuals to bit strings via von Neumann debiasing and hash-based extraction. Apply NIST SP 800-22, TestU01 BigCrush, and Dieharder test suites. Compute autocorrelation functions, spectral analysis, and entropy estimates.

PHASE 4 — Cryptographic Qualification (Days 90–150): Test against AIS 31 (BSI) Class P2 requirements. Measure min-entropy H_min. Compare against hardware TRNGs (Intel RDRAND, Quantis QRNG) as benchmarks. Assess practical bit generation rate.

PHASE 5 — Independent Replication (Days 120–180): Repeat full pipeline on new observations from a second telescope (e.g., MeerKAT + Parkes cross-validation). Verify reproducibility of residual extraction.

# PULSAR TRNG VALIDATION PIPELINE
# Architecture: Data Ingestion → Timing Model Fitting → Noise Decomposition
#               → Bit Extraction → Randomness Testing → Cryptographic Qualification

import numpy as np
import subprocess
from scipy import stats
from hashlib import sha3_256

# ============================================================
# MODULE 1: DATA INGESTION
# ============================================================
class PulsarDataLoader:
    def __init__(self, pulsar_name, data_source="IPTA_DR2"):
        self.pulsar = pulsar_name
        self.source = data_source
        self.toa_file = f"data/{pulsar_name}.tim"
        self.par_file = f"data/{pulsar_name}.par"

    def load_toas(self):
        # Load TOA data: columns = [MJD, TOA_error_us, frequency_MHz, backend]
        # Returns numpy array shape (N_obs, 4)
        toas = np.loadtxt(self.toa_file, comments='#')
        return toas

    def validate_multifrequency(self, toas, min_bands=2, min_separation_mhz=200):
        frequencies = toas[:, 2]
        unique_bands = self._cluster_frequencies(frequencies, separation=min_separation_mhz)
        assert len(unique_bands) >= min_bands, f"Insufficient frequency coverage: {unique_bands}"
        return unique_bands

# ============================================================
# MODULE 2: TIMING MODEL FITTING
# ============================================================
class TimingModelFitter:
    def __init__(self, pulsar_name, software="TEMPO2"):
        self.pulsar = pulsar_name
        self.software = software  # "TEMPO2" or "PINT"

    def fit_model(self):
        if self.software == "TEMPO2":
            cmd = f"tempo2 -f {self.pulsar}.par {self.pulsar}.tim -nofit -residuals"
            result = subprocess.run(cmd, shell=True, capture_output=True)
            residuals = self._parse_tempo2_output(result.stdout)
        elif self.software == "PINT":
            # Use PINT Python API
            from pint.models import get_model
            from pint.toa import get_TOAs
            model = get_model(f"{self.pulsar}.par")
            toas = get_TOAs(f"{self.pulsar}.tim")
            residuals = model.residuals(toas).time_resids.to('us').value
        return residuals  # shape: (N_obs,) in microseconds

    def cross_validate_software(self, res_tempo2, res_pint):
        r, p = stats.pearsonr(res_tempo2, res_pint)
        assert abs(r) > 0.95, f"Software disagreement: r={r:.3f}"
        return r

    def check_model_completeness(self, par_file):
        required_params = ['F0','F1','RAJ','DECJ','PMRA','PMDEC','PX','DM','DM1']
        with open(par_file) as f:
            content = f.read()
        missing = [p for p in required_params if p not in content]
        if missing:
            raise ValueError(f"Incomplete timing model, missing: {missing}")

# ============================================================
# MODULE 3: BAYESIAN NOISE DECOMPOSITION
# ============================================================
class BayesianNoiseDecomposer:
    """Wraps enterprise/PTMCMCSampler for noise analysis"""

    def __init__(self, pulsar_name, n_freqs=30, n_iter=2_000_000):
        self.pulsar = pulsar_name
        self.n_freqs = n_freqs
        self.n_iter = n_iter

    def build_noise_model(self):
        # enterprise model: white noise + red spin noise + DM noise + GWB
        from enterprise.pulsar import Pulsar
        from enterprise.signals import (white_signals, red_signals,
                                         gp_signals, parameter)
        # [enterprise model construction — see enterprise documentation]
        # Returns: PTA object with full noise model
        pass

    def run_mcmc(self, pta_model):
        from PTMCMCSampler.PTMCMCSampler import PTSampler
        # Run sampler, return chain
        # chain shape: (n_iter, n_params)
        pass

    def extract_stochastic_residuals(self, residuals, chain, burn_in=100_000):
        # Subtract ML white noise and GWB realizations
        # Returns: stochastic_residuals (spin noise + DM noise)
        ml_params = chain[burn_in:].mean(axis=0)
        gwb_realization = self._compute_gwb_realization(ml_params)
        white_noise_realization = self._compute_white_noise(ml_params)
        stochastic = residuals - gwb_realization - white_noise_realization
        return stochastic

# ============================================================
# MODULE 4: BIT EXTRACTION
# ============================================================
class BitExtractor:
    def __init__(self, method="von_neumann_sha3"):
        self.method = method

    def normalize_residuals(self, residuals):
        mu, sigma = residuals.mean(), residuals.std()
        return (residuals - mu) / sigma

    def von_neumann_debias(self, normalized_residuals):
        """Convert continuous residuals to bits via threshold + von Neumann"""
        # Step 1: Threshold at median → binary sequence
        median = np.median(normalized_residuals)
        bits = (normalized_residuals > median).astype(int)

        # Step 2: Von Neumann debiasing on consecutive pairs
        debiased = []
        for i in range(0, len(bits) - 1, 2):
            if bits[i] == 0 and bits[i+1] == 1:
                debiased.append(0)
            elif bits[i] == 1 and bits[i+1] == 0:
                debiased.append(1)
            # discard 00 and 11 pairs
        return np.array(debiased)

    def hash_extract(self, debiased_bits, block_size=256):
        """Apply SHA-3-256 to blocks for final output"""
        output_bits = []
        for i in range(0, len(debiased_bits) - block_size, block_size):
            block = debiased_bits[i:i+block_size]
            block_bytes = np.packbits(block).tobytes()
            hash_bytes = sha3_256(block_bytes).digest()
            hash_bits = np.unpackbits(np.frombuffer(hash_bytes, dtype=np.uint8))
            output_bits.extend(hash_bits.tolist())
        return np.array(output_bits)

    def compute_bit_rate(self, n_bits, observation_hours):
        return n_bits / observation_hours  # bits/hour

# ============================================================
# MODULE 5: RANDOMNESS TESTING
# ============================================================
class RandomnessTester:
    def __init__(self, bit_sequence):
        self.bits = bit_sequence
        self.n = len(bit_sequence)

    def run_nist_sp800_22(self, output_dir="nist_results"):
        # Write bits to file, call NIST test suite binary
        bit_file = f"{output_dir}/bits.bin"
        np.packbits(self.bits).tofile(bit_file)
        cmd = f"./sts -b {self.n} -i 1 -s 1 -F 0 {bit_file}"
        subprocess.run(cmd, shell=True)
        results = self._parse_nist_results(output_dir)
        pass_rate = sum(results.values()) / len(results)
        return pass_rate, results

    def run_dieharder(self):
        # Pipe bits to dieharder
        cmd = f"dieharder -a -g 201 -f bits.bin"
        result = subprocess.run(cmd, shell=True, capture_output=True)
        return self._parse_dieharder_output(result.stdout)

    def compute_autocorrelation(self, max_lag=100):
        acf = [np.corrcoef(self.bits[:-lag], self.bits[lag:])[0,1]
               for lag in range(1, max_lag+1)]
        # Ljung-Box test
        from statsmodels.stats.diagnostic import acorr_ljungbox
        lb_result = acorr_ljungbox(self.bits, lags=max_lag, return_df=True)
        return np.array(acf), lb_result

    def compute_min_entropy(self):
        # Estimate min-entropy from empirical distribution
        # Use sliding window of length 8 for byte-level analysis
        bytes_arr = np.packbits(self.bits)
        counts = np.bincount(bytes_arr, minlength=256)
        probs = counts / counts.sum()
        max_prob = probs.max()
        h_min = -np.log2(max_prob) / 8  # normalize to bits/bit
        return h_min

# ============================================================
# MODULE 6: PREDICTABILITY TEST
# ============================================================
class PredictabilityTester:
    def __init__(self, residuals, train_frac=0.8):
        n = len(residuals)
        self.train = residuals[:int(n*train_frac)]
        self.test = residuals[int(n*train_frac):]

    def lstm_predictor(self, seq_len=10, epochs=100):
        import tensorflow as tf
        # Build LSTM: 2 layers × 128 units, dropout 0.2
        model = tf.keras.Sequential([
            tf.keras.layers.LSTM(128, return_sequences=True,
                                  input_shape=(seq_len, 1), dropout=0.2),
            tf.keras.layers.LSTM(128, dropout=0.2),
            tf.keras.layers.Dense(1)
        ])
        model.compile(optimizer='adam', loss='mse')
        X_train, y_train = self._create_sequences(self.train, seq_len)
        model.fit(X_train, y_train, epochs=epochs, verbose=0)
        X_test, y_test = self._create_sequences(self.test, seq_len)
        predictions = model.predict(X_test)
        nrmse = np.sqrt(np.mean((predictions.flatten() - y_test)**2)) / y_test.std()
        return nrmse  # > 0.85 = unpredictable

    def gaussian_process_predictor(self):
        from sklearn.gaussian_process import GaussianProcessRegressor
        from sklearn.gaussian_process.kernels import Matern
        gp = GaussianProcessRegressor(kernel=