haskell/src/Verified/Entropy.hs

-- True Entropy — Haskell verification layer.
-- H(P) = log2(N) - (1/N) * sum(c_i * log2(c_i))
-- Symbolic representation: exact rational coefficients * log2(integer base).
-- No Double ever in the core computation.

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Verified.Entropy where

import Data.Ratio ((%))
import qualified Data.Map.Strict as Map
import Data.List (foldl', toList)

-- ─── Symbolic entropy ────────────────────────────────────────────────────────

-- | coeff * log2(base)
data SymLog2 = SymLog2
    { coeff :: Rational
    , base  :: Integer
    } deriving (Eq, Show)

-- | Symbolic entropy: sum of SymLog2 terms
-- H = log2(N) + sum(-c_i/N * log2(c_i))
newtype SymbolicEntropy = SymbolicEntropy { terms :: [SymLog2] }
    deriving Show

-- | Exact symbolic entropy from integer counts.
-- No evaluation, no rounding, no Double.
entropy :: (Foldable f, Integral a) => f a -> SymbolicEntropy
entropy counts = SymbolicEntropy (positiveTerm : negativeTerms)
  where
    cs          = map fromIntegral $ filter (> 0) $ toList counts
    n           = sum cs :: Integer
    -- H = log2(N) term
    positiveTerm = SymLog2 { coeff = 1 % n, base = n }
    -- -c_i/N * log2(c_i) terms
    negativeTerms = map (\c -> SymLog2 { coeff = -(c % n), base = c }) cs

-- ─── Interval evaluation ─────────────────────────────────────────────────────

-- | Rational approximation of log2 for interval bounds.
-- Uses convergents of continued fraction for log2(n).
-- Returns (lower, upper) as Rational pair at given precision.
log2Interval :: Integer -> Int -> (Rational, Rational)
log2Interval n prec
    | n <= 0    = error "log2Interval: non-positive argument"
    | n == 1    = (0, 0)
    | otherwise =
        -- log2(n) = log(n) / log(2)
        -- Use rational approximation: ln(2) ≈ 6931471805599453/10000000000000000
        let ln2_lo = 6931471805599453 % 10000000000000000
            ln2_hi = 6931471805599454 % 10000000000000000
            -- ln(n) approximated via Taylor series for small n, else recursion
            lnn = lnRational n prec
            lo  = fst lnn / ln2_hi  -- divide by larger denominator → smaller result
            hi  = snd lnn / ln2_lo
        in (lo, hi)

-- | Rational interval for ln(n), precision as number of terms
lnRational :: Integer -> Int -> (Rational, Rational)
lnRational n prec
    | n == 1    = (0, 0)
    | n == 2    = (6931471805599453 % 10000000000000000,
                   6931471805599454 % 10000000000000000)
    | even n    = let (l, h) = lnRational (n `div` 2) prec
                      (l2, h2) = lnRational 2 prec
                  in (l + l2, h + h2)
    | otherwise = -- ln(n) ≈ ln(n-1) + 2/(2n-1) + ... (first term bound)
                  let (l, h) = lnRational (n - 1) prec
                      delta_lo = 2 % (2 * n - 1 + 1)
                      delta_hi = 2 % (2 * n - 1 - 0)
                  in (l + delta_lo, h + delta_hi)

-- ─── Connection to Bifrost WORM ──────────────────────────────────────────────

-- | Serialize entropy for WORM audit log
serializeEntropy :: SymbolicEntropy -> String
serializeEntropy (SymbolicEntropy ts) =
    "H = " ++ unwords (map showTerm ts)
  where
    showTerm (SymLog2 c b) =
        "(" ++ show (numerator c) ++ "/" ++ show (denominator c) ++
        ")*log2(" ++ show b ++ ")"
    numerator r   = let (n, _) = (floor (r * 10^15), ()) in n
    denominator _ = 10^15 :: Integer

-- ─── Ramanujan connection ─────────────────────────────────────────────────────

-- | Partition entropy: entropy of the partition number sequence p(0)..p(n)
-- This is the entropy of the complexity invariant distribution.
partitionEntropy :: Int -> SymbolicEntropy
partitionEntropy n = entropy (partitionNumbers n)

-- | Euler pentagonal partition numbers (exact)
partitionNumbers :: Int -> [Integer]
partitionNumbers n = take (n + 1) ps
  where
    ps = 1 : [compute k | k <- [1..n]]
    compute i = sum
        [ sign k * safeIdx (i - penta k)
        | k <- [1..i]
        , penta k <= i
        ]
        +
        sum
        [ sign (-k) * safeIdx (i - penta (-k))
        | k <- [1..i]
        , penta (-k) <= i
        ]
    penta k   = k * (3 * k - 1) `div` 2
    sign k    = if odd k then 1 else -1
    safeIdx j = if j < 0 then 0 else ps !! j