SilverSight/python/q16_canonical.py

"""Canonical Q16_16 fixed-point arithmetic.

Single source of truth: CoreFormalism/Q16_16_Spec.lean
All operations MUST produce identical results to the Lean implementation.

Q16_16 represents fixed-point numbers with 16 integer bits and 16 fractional bits.
Range: [-32768.0, 32767.9999847412109375]
Resolution: 1/65536 ≈ 0.0000152587890625

CANONICAL ROUNDING MODE: round-half-up (banker's rounding)
- Values exactly at half-LSB round to nearest even
- All other values round to nearest
"""

import math
import struct

# ============================================================
# §1  CONSTANTS
# ============================================================

Q16_SCALE: int = 65536  # 2^16
Q16_MAX_RAW: int = 2147483647   # INT32_MAX
Q16_MIN_RAW: int = -2147483648  # INT32_MIN
Q16_MAX_FLOAT: float = 32767.9999847412109375  # Max representable
Q16_MIN_FLOAT: float = -32768.0                # Min representable
Q16_RESOLUTION: float = 1.0 / Q16_SCALE        # ≈ 0.0000152587890625


# ============================================================
# §2  CONVERSIONS
# ============================================================

def float_to_q16(f: float) -> int:
    """Convert float to Q16_16 raw value with canonical round-half-up.
    
    Uses banker's rounding: round(x * 65536) with ties to nearest even.
    Result is clamped to [INT32_MIN, INT32_MAX].
    
    Args:
        f: Float value in range [-32768.0, 32767.9999847412109375]
    
    Returns:
        32-bit signed integer representing the Q16_16 value
    
    Raises:
        ValueError: If f is NaN or infinite
    """
    if math.isnan(f) or math.isinf(f):
        raise ValueError(f"Cannot convert non-finite float to Q16_16: {f}")
    
    # Python's round() implements banker's rounding (round-half-to-even)
    # round(x) = nearest integer, ties go to nearest even integer
    scaled = f * Q16_SCALE
    rounded = round(scaled)  # Banker's rounding: ties to even
    
    # Clamp to 32-bit signed range with saturation
    return max(Q16_MIN_RAW, min(Q16_MAX_RAW, rounded))


def q16_to_float(q: int) -> float:
    """Convert Q16_16 raw value to float.
    
    This is exact: no rounding occurs.
    
    Args:
        q: 32-bit signed integer Q16_16 raw value
    
    Returns:
        Float value = q / 65536.0
    """
    return q / Q16_SCALE


def int_to_q16(i: int) -> int:
    """Convert integer to Q16_16 raw value (exact, no rounding).
    
    The integer is scaled by 65536. Clamped to valid range.
    
    Args:
        i: Integer in range [-32768, 32767]
    
    Returns:
        Q16_16 raw value = clamp(i * 65536)
    """
    scaled = i * Q16_SCALE
    return max(Q16_MIN_RAW, min(Q16_MAX_RAW, scaled))


def q16_to_int(q: int) -> int:
    """Convert Q16_16 to integer (truncates toward zero).
    
    Args:
        q: Q16_16 raw value
    
    Returns:
        Integer part = q // 65536 (toward zero)
    """
    # Python's // truncates toward negative infinity, so we need
    # to handle negative values correctly for toward-zero truncation
    if q >= 0:
        return q // Q16_SCALE
    else:
        return -(-q // Q16_SCALE)


# ============================================================
# §3  ARITHMETIC OPERATIONS (all with saturation)
# ============================================================

def q16_add(a: int, b: int) -> int:
    """Add two Q16_16 values with saturation.
    
    result = clamp(a + b)
    """
    result = a + b
    return max(Q16_MIN_RAW, min(Q16_MAX_RAW, result))


def q16_sub(a: int, b: int) -> int:
    """Subtract two Q16_16 values with saturation.
    
    result = clamp(a - b)
    """
    result = a - b
    return max(Q16_MIN_RAW, min(Q16_MAX_RAW, result))


def q16_mul(a: int, b: int) -> int:
    """Multiply two Q16_16 values with canonical rounding.
    
    result = canonical_round((a * b) / 65536)
    
    Uses 64-bit intermediate, then applies banker's rounding.
    """
    # Use Python's arbitrary precision integers (no overflow issue)
    prod_64 = a * b
    
    # Divide by scale with banker's rounding
    # prod_64 / 65536 with half-to-even
    scaled = prod_64 / Q16_SCALE  # This is a float division for correct rounding
    rounded = round(scaled)       # Banker's rounding
    
    return max(Q16_MIN_RAW, min(Q16_MAX_RAW, rounded))


def q16_div(a: int, b: int) -> int:
    """Divide two Q16_16 values with canonical rounding.
    
    result = canonical_round((a * 65536) / b)
    
    Args:
        a: Dividend (Q16_16 raw value)
        b: Divisor (Q16_16 raw value), must not be zero
    
    Raises:
        ZeroDivisionError: If b is zero
    """
    if b == 0:
        raise ZeroDivisionError("Q16_16 division by zero")
    
    # (a * 65536) / b with banker's rounding
    num = a * Q16_SCALE
    scaled = num / b  # Float division for correct rounding
    rounded = round(scaled)
    
    return max(Q16_MIN_RAW, min(Q16_MAX_RAW, rounded))


# ============================================================
# §4  COMPARISON OPERATIONS
# ============================================================

def q16_eq(a: int, b: int) -> bool:
    return a == b

def q16_lt(a: int, b: int) -> bool:
    return a < b

def q16_le(a: int, b: int) -> bool:
    return a <= b


# ============================================================
# §5  UTILITY FUNCTIONS
# ============================================================

def q16_from_bytes(raw_bytes: bytes) -> int:
    """Convert 4 bytes (little-endian int32) to Q16_16 raw value."""
    return struct.unpack('<i', raw_bytes)[0]


def q16_to_bytes(q: int) -> bytes:
    """Convert Q16_16 raw value to 4 bytes (little-endian int32)."""
    # Clamp first to ensure valid int32
    clamped = max(Q16_MIN_RAW, min(Q16_MAX_RAW, q))
    return struct.pack('<i', clamped)


def q16_is_valid(q: int) -> bool:
    """Check if a raw value is in the valid Q16_16 range."""
    return Q16_MIN_RAW <= q <= Q16_MAX_RAW


def q16_repr(q: int) -> str:
    """Return human-readable representation of Q16_16 value."""
    return f"Q16_16({q} / 65536 = {q16_to_float(q)})"


# ============================================================
# §6  EXPORTS FOR C INTEROP
# ============================================================

# These functions provide the C-compatible interface for the roundtrip test
def c_float_to_q16(f: float) -> int:
    """C-compatible wrapper for float_to_q16."""
    return float_to_q16(f)


def c_q16_to_float(q: int) -> float:
    """C-compatible wrapper for q16_to_float."""
    return q16_to_float(q)