"""Bures volumes of quantum states of fixed rank -- numberdb.org/T572.

Run with SageMath (or compatible passagemath):
    sage -pip install numberdb
    sage -python generate.py
    sage -python generate.py --fill  # fill the draft, NUMBERDB_API_KEY required

The default command verifies all stored entries. --fill writes entries to the
existing draft; it does not publish or offer the table. For private-draft reads,
the owner must set NUMBERDB_API_KEY. Reproduction is also available without
network access: import BuresVolumes and call enumerate() and value().

Sommers and Zyczkowski, arXiv:quant-ph/0304041v3 (17 February 2005),
equation (5.15), with nullity d=N-r. This is the Bures metric with the qubit
hemisphere radius 1/2, not the four-times-larger statistical metric.
All gamma factors at integers or half-integers are reduced over QQ first.
A 64-bit guard beyond the requested decimal precision leaves more than 115
significant digits in all default balls (see numerical-evidence.json).
Only the final rational multiple of an integer power of pi uses Arb.
The argument selection and independent checks accompany this file.
"""
import os
import sys
from math import factorial
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

WORKING_GUARD = 64


def gamma_half(twice_argument):
    """Return (q,e) with Gamma(twice_argument/2) = q*sqrt(pi)^e."""
    m = int(twice_argument)
    if m < 1:
        raise ValueError('gamma argument must be positive')
    if m % 2 == 0:
        return QQ(factorial(m // 2 - 1)), 0
    k = (m - 1) // 2
    return QQ(factorial(2*k)) / (QQ(4)**k * QQ(factorial(k))), 1


def coefficient(field, N, rank):
    """Return exact q,k with V=q*pi^k; rank always means actual rank."""
    beta = {'real': 1, 'complex': 2}[field]
    N, r = int(N), int(rank)
    if not (N >= 2 and 1 <= r <= N):
        raise ValueError('need N >= 2 and 1 <= rank <= N')
    d = N-r
    D = r-1 + beta*r*(2*N-r-1)//2
    q = QQ(2)**(-D)
    p = D+1
    gamma_q, gamma_p = gamma_half(D+1)
    q /= gamma_q
    p -= gamma_p
    for j in range(1, r+1):
        for m, sign in ((j*beta, 1), (2+(2*d+j-1)*beta, 1),
                        ((d+j)*beta, -1), (2+(d+j-1)*beta, -1)):
            gq, gp = gamma_half(m)
            q *= gq**sign
            p += sign*gp
    assert p % 2 == 0 and p > 0
    return q, p//2


class BuresVolumes(numberdb.Generator):
    table = 'T572'
    parameters = ('field', 'N', 'rank')
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'argument-selection.md',
             'numerical-evidence.md', 'numerical-evidence.json',
             'skill-manifest.json', 'runtime-versions.json', 'program.py')

    def enumerate(self, max_N=16):
        for field in ('real', 'complex'):
            for N in range(2, int(max_N)+1):
                for r in range(1, N+1):
                    yield {'field': field, 'N': N, 'rank': r}

    def value(self, params, digits):
        q, k = coefficient(params['field'], params['N'], params['rank'])
        R = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
        value = R(q)*R.pi()**k
        notes = {
            ('real', 2, 1): r'The pure rebit states form a circle of radius $1/2$, with length $\pi$.',
            ('real', 2, 2): r'The full-rank rebit states form an open hemisphere of a sphere of radius $1/2$, with area $\pi/2$.',
            ('complex', 2, 1): r'The pure qubit states form a sphere of radius $1/2$, with area $\pi$.',
            ('complex', 2, 2): r'The full-rank qubit states form an open three-dimensional hemisphere of radius $1/2$, with volume $\pi^2/8$.'}
        key = (params['field'], int(params['N']), int(params['rank']))
        return {'number': value, 'comment': notes[key]} if key in notes else value


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=', 1)[0].replace('_', '').isalnum():
            key = key.split('=', 1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY'] = key
    g = BuresVolumes()
    if '--fill' in sys.argv:
        print(g.publish(message='Fill the Bures rank-volume draft with certified 100-digit values',
                        assisted_by='codex-cli'))
    else:
        report = g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)