#!/usr/bin/env python3
"""
Unified Ryser-based Monte Carlo estimator for bosonic permanental distributions.
No regime switching — single estimator for all p ≤ 6.
"""

import math
from typing import List

def ryser_permanent(matrix: List[List[complex]]) -> complex:
    """Compute permanent using Ryser's formula."""
    p = len(matrix)
    if p == 0:
        return 1.0
    # Generate all subsets via bitmask
    total = 0j
    for mask in range(1, 1 << p):
        sign = (-1) ** (p - bin(mask).count("1"))
        prod = 1.0
        for row in matrix:
            s = 0.0
            for j in range(p):
                if mask & (1 << j):
                    s += row[j].real  # Use real part for simplicity
            prod *= s
        total += sign * prod
    return total

def lambda_interpolation(p: int, N: int) -> float:
    """Continuous interpolation parameter.
    
    λ(p) = exp(-p²/N)
    λ → 1 as p² << N (weak interference, fully bosonic)
    λ → 0 as p² >> N (strong interference, factorized)
    """
    return math.exp(-p * p / N)

def unified_bosonic_estimator(U: List[List[complex]], p: int, samples: int = 1000) -> dict:
    """Single estimator for all photon numbers.
    
    Returns distribution histogram using continuous λ(p) weighting.
    """
    N = len(U)
    lam = lambda_interpolation(p, N)
    
    # Monte Carlo sampling
    results = {}
    for _ in range(samples):
        # Sample modes (simplified)
        S = [j % N for j in range(p)]  # Placeholder
        M = [[U[S[i]][j] for j in range(p)] for i in range(p)]
        perm = ryser_permanent(M)
        w = abs(perm) ** 2
        
        # Apply λ(p) weighting
        weighted = lam * w + (1 - lam) * sum(abs(U[i][j])**2 for i in range(N) for j in range(N)) / (N * N)
        results[str(S)] = weighted
    
    return {"histogram": results, "lambda": lam, "p": p, "N": N}