BC3 aserti9-2d reference implementation
aserti9-2d is a DAA that uses the exact same mathematical formula as BCH, but approximates f(x) = 2^x over [0,1] using a 9th-degree polynomial with 31-bit fractional arithmetic, instead of BCH's cubic approximation with 16-bit fractional arithmetic.
Using reference nBits of 0x1c0ffff0, a 172800-second half-life and a powLimit of 2^224 - 1, a sweep of all integer schedule errors from -172800 to +172800 seconds found that aserti9-2d matched the ideal ASERT nBits from the mathematical formula in 99.9809% of cases, compared with 0.3061% for aserti3-2d. After compact encoding, its decoded targets had approximately 126.4 times lower mean absolute relative error to the ideal mathematical target.
https://bc3.network/
vout=5 contains BC3Dev's signature over the payload in vout=7 | vout=6 contains BC3Dev's signature over the payload in vout=8 | vout=7 contains the aserti9-2d constant-generation Python script | vout=8 contains the aserti9-2d mainnet arithmetic reference implementation, as an algorithmically equivalent port of the C++ implementation in BitcoinIII-Core
-----BEGIN PGP SIGNATURE-----

iKcEABMJAC8WIQTZaKbwh8YXCwkrkAIPDcoDEzprhgUCasmz9xEcYmMzZGV2QHBy
b3Rvbi5tZQAKCRAPDcoDEzprhtUIAX9/sS/yWVvogXdjy7orrfQyfUloCGGszA3N
ms7vcUKMx6z3DKqFonm7gY3sFf1j+RMBf09RGeewnoQfHF2xWCFpk+ShzRU4ZV3k
JWyqRcta4DI44xw3TlhvIUJ2GIDFRICNoQ==
=H1k2
-----END PGP SIGNATURE-----

-----BEGIN PGP SIGNATURE-----

iKcEABMJAC8WIQTZaKbwh8YXCwkrkAIPDcoDEzprhgUCasm0IBEcYmMzZGV2QHBy
b3Rvbi5tZQAKCRAPDcoDEzprhjTzAX95vBR+TSbSUb3nSLWh5YnkIXJ4El8cBxh/
E6oIxTWi7mpqdvCbG/GUjPQoIlPbIkoBf29D+a/aIECJMRa7GyV1y0Roaa0oDNbE
dNU1Dt6q3wusIPsnHACCbJVB8inEx6nxog==
=rPd6
-----END PGP SIGNATURE-----

#!/usr/bin/env python3

from mpmath import mp

PLACES = 60
DECIMAL_SCALE = 10**PLACES
SCALE = 1 << 31


def check(condition: bool, message: str) -> None:
    if not condition:
        raise RuntimeError(message)


def decimal_coefficients(dps: int) -> tuple:
    with mp.workdps(dps):
        ln2 = mp.log(2)
        # P(x) = 1 + x + x*(1-x)*Q(x) -> Q itself has degree 7
        nodes = [(1 - mp.cos((mp.pi*(2*i + 1))/18))/2 for i in range(9)]
        matrix = [[x*(1 - x)*(x**k) for k in range(8)] + [-((-1)**i)*(2**x)] for i, x in enumerate(nodes)]
        seed = list(mp.lu_solve(matrix, [2**x - 1 - x for x in nodes]))
        def equations(*z):
            # unknowns: Q's 8 coefficients, peak error E, 9 extrema
            q, error, points = z[:8], z[8], z[9:]
            residuals = []
            for i, x in enumerate(points):
                value, slope = mp.polyval(q[::-1], x, derivative=True)
                p = 1 + x + x*(1 - x)*value
                dp = 1 + (1 - 2*x)*value + x*(1 - x)*slope
                residuals += [p/(2**x) - 1 - ((-1)**i)*error, dp - ln2*p]
            return residuals
        solution = list(mp.findroot(equations, seed + nodes, tol=mp.mpf(10)**(20 - dps), maxsteps=50))
        points = [mp.mpf(0)] + solution[9:] + [mp.mpf(1)]
        check(all(a < b for a, b in zip(points, points[1:])), "Extrema are not distinct, ordered and inside (0,1)")
        check(solution[8] > 0, "Unexpected error sign")
        check(max(map(abs, equations(*solution))) < mp.mpf(10)**(30 - dps), "Insufficient numerical convergence")
        q = solution[:8]
        c = [mp.mpf(1), 1 + q[0]]
        c += [q[k-1] - q[k-2] for k in range(2, 9)] + [-q[7]]
        check(all(v > 0 for v in c), "Unexpected coefficient sign")
        return tuple(int(mp.floor(v*DECIMAL_SCALE + mp.mpf("0.5"))) for v in c)


def main() -> None:
    decimals = decimal_coefficients(100)
    check(decimals == decimal_coefficients(140), "Decimal coefficients changed when precision was increased")
    check(decimals[0] == DECIMAL_SCALE and sum(decimals) == 2*DECIMAL_SCALE, "The 60-place decimal export does not preserve the endpoints")
    # exact Python integer arithmetic onwards
    nearest = [(v*SCALE + DECIMAL_SCALE//2)//DECIMAL_SCALE for v in decimals]
    check(2*SCALE - sum(nearest) == 1, "Unexpected Q31 endpoint deficit")
    index = max(range(1, 10), key=lambda k: decimals[k]*SCALE - nearest[k]*DECIMAL_SCALE)
    check(index == 7, "The endpoint correction no longer selects x^7")
    coefficients = nearest.copy()
    coefficients[index] += 1
    check(coefficients[0] == SCALE and sum(coefficients) == 2*SCALE, "Q31 endpoint check failed")
    for k, value in enumerate(decimals):
        print(f"c[{k}] = {value//DECIMAL_SCALE}.{value % DECIMAL_SCALE:0{PLACES}d}")
    print(f"================\nNearest Q31:            {tuple(nearest)}")
    print(f"Endpoint-corrected Q31: {tuple(coefficients)}")
    print(f"Correction: x^{index}: {nearest[index]} -> {coefficients[index]}")
    print(f"Sum: {sum(coefficients)} = 2 * 2^31")


if __name__ == '__main__':
    main()

#!/usr/bin/env python3

SCALE       =  2**31
UINT32_MAX  =  0xffffffff
UINT64_MAX  =  0xffffffffffffffff
UINT256_MAX =  0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff
INT64_MAX   =  2**63     - 1
INT64_MIN   = -INT64_MAX - 1
POW_LIMIT   =  2**224    - 1 # actual limit, not the decoded 0x1d00ffff target

ASERT_EXPONENT_LIMIT = 256 << 31

COEFFICIENTS = (
    2147483648, 1488522236, 515882495, 119194168, 20654757,
    2863430, 330624, 32990, 2640, 308,
) # derived from `aserti9-constants.py`


def bits_to_target(bits: int) -> int:
    if not 0 <= bits <= UINT32_MAX:
        raise ValueError("nBits must fit into uint32_t")
    size = bits >> 24
    word = bits & 0x007fffff
    # reject overflow before constructing the shifted target
    if (bits & 0x00800000 or size > 34 or (word > 255 and size > 33) or (word > 0xffff and size > 32)):
        raise ValueError("Invalid compact target!")
    target = word >> (8*(3 - size)) if size <= 3 else word << (8*(size - 3))
    if target == 0:
        raise ValueError("Invalid compact target!")
    return target


def target_to_bits(target: int) -> int:
    if not 0 < target <= UINT256_MAX:
        raise ValueError("Target must be a positive uint256")
    size = (target.bit_length() + 7) // 8
    word = target << (8*(3 - size)) if size <= 3 else target >> (8*(size - 3))
    if word & 0x00800000: # reserve the compact format's sign bit
        word >>= 8
        size += 1
    return ((size << 24) | word) & UINT32_MAX


def fractional_factor(frac: int) -> int:
    if not 0 <= frac < SCALE:
        raise ValueError("Fraction must lie in [0, 2**31)")
    p = COEFFICIENTS[-1]
    for coefficient in reversed(COEFFICIENTS[:-1]):
        p = coefficient + ((p*frac + SCALE//2) >> 31) # fits in uint64_t
    return p & UINT32_MAX # static_cast<uint32_t>(p)


def fractional_quotient(rem: int, denom: int) -> int:
    if not 0 <= rem < denom <= INT64_MAX:
        raise ValueError("Invalid fractional quotient parameters")
    if rem <= (UINT64_MAX >> 31):
        return ((rem << 31) // denom) & UINT32_MAX
    res = 0
    for _ in range(31): # binary long division, using only uint64_t
        rem  *= 2
        res <<= 1
        if rem  >= denom:
            rem -= denom
            res |= 1
    return res


def asert_exponent(spacing:     int,
                   time_diff:   int,
                   height_diff: int,
                   half_life:   int) -> int:
    if not (0 < spacing <= INT64_MAX and 0 < half_life <= INT64_MAX
            and 0 <= height_diff < INT64_MAX and INT64_MIN <= time_diff <= INT64_MAX):
        raise ValueError("Invalid ASERT parameters")
    if height_diff + 1 > INT64_MAX // spacing:
        raise ValueError("ASERT scheduled interval would overflow int64_t")
    scheduled = spacing*(height_diff + 1)
    if time_diff < INT64_MIN + scheduled:
        raise ValueError("ASERT signed time error would overflow int64_t")
    delta = time_diff - scheduled
    whole, remainder = divmod(delta, half_life)
    if delta < 0 and remainder != 0: # C++ / and % operators truncate towards zero
        whole += 1
        remainder -= half_life
    if whole >=  256:
        return ASERT_EXPONENT_LIMIT
    if whole <= -256:
        return -ASERT_EXPONENT_LIMIT
    magnitude = -remainder if remainder < 0 else remainder
    fraction = fractional_quotient(magnitude, half_life)
    return whole*SCALE + (-fraction if remainder < 0 else fraction)


def calculate_asert(reference:   int,
                    spacing:     int,
                    time_diff:   int,
                    height_diff: int,
                    pow_limit:   int,
                    half_life:   int) -> int:
    if not 0 < reference <= pow_limit <= POW_LIMIT: # narrow CalculateASERT
        raise ValueError("Need 0 < reference <= pow_limit <= 2**224 - 1")
    exponent = asert_exponent(spacing, time_diff, height_diff, half_life)
    if exponent >=  ASERT_EXPONENT_LIMIT:
        return pow_limit
    if exponent <= -ASERT_EXPONENT_LIMIT:
        return 1
    shifts   = (exponent >> 31) - 31 # signed right shift floors, as in C++20
    fraction = (exponent & UINT32_MAX) & (SCALE - 1)
    target   =  reference*fractional_factor(fraction) # fits uint256
    if shifts <= 0:
        target >>= -shifts
    elif target > (pow_limit >> shifts): # check before shifting
        return pow_limit
    else:
        target <<= shifts
    return max(1, min(target, pow_limit))


def aserti9_nbits(ref_bits:    int, # anchor block nBits
                  time_diff:   int, # parent timestamp minus anchor's PARENT's timestamp!
                  height_diff: int, # parent height minus anchor's height
                  *,
                  spacing:   int = 600,
                  half_life: int = 172800,
                  pow_limit: int = POW_LIMIT) -> int:
    reference = bits_to_target(ref_bits)
    target = calculate_asert(reference, spacing, time_diff, height_diff, pow_limit, half_life)
    return target_to_bits(target)


def main():
    assert COEFFICIENTS[0] == SCALE and sum(COEFFICIENTS) == 2*SCALE
    assert fractional_factor(0)         == SCALE
    assert fractional_factor(SCALE//2)  == 3037000500
    assert fractional_factor(SCALE - 1) == UINT32_MAX
    assert target_to_bits(0x80)         == 0x02008000
    assert target_to_bits(0x12345678)   == 0x04123456 # low-byte truncation
    assert bits_to_target(0x04123456)   == 0x12345600
    assert bits_to_target(0x1d00ffff)    < POW_LIMIT
    assert target_to_bits(POW_LIMIT)    == 0x1d00ffff
    ref_bits    = 0x1c0ffff0
    height_diff = 1000
    scheduled   = 600*(height_diff + 1)
    # signed schedule error in seconds, expected compact target
    vectors = (
        (        0, 0x1c0ffff0), # on schedule
        (       -1, 0x1c0fffeb),
        (        1, 0x1c0ffff4),
        (   -86400, 0x1c0b5043), # target/sqrt(2) (approx)
        (    86400, 0x1c16a087), # target*sqrt(2) (approx)
        (  -172800, 0x1c07fff8), # target halved
        (   172800, 0x1c1fffe0), # target doubled
        (   864000, 0x1d00ffff), # clamped to powLimit
        ( 44236800, 0x1d00ffff), # positive exponent saturation
        (-44236800, 0x01010000), # minimum target = 1
    )
    for delta, expected in vectors:
        actual = aserti9_nbits(ref_bits, scheduled + delta, height_diff)
        assert actual == expected
        print(f"delta={delta:+10d} s  nBits=0x{actual:08x}")
    assert aserti9_nbits(ref_bits, 599,             0, half_life=INT64_MAX) == ref_bits
    assert aserti9_nbits(ref_bits, 600 - (1 << 39), 0, half_life = 1 << 40) == 0x1c0b5043
    assert asert_exponent(1, INT64_MIN + 1, 0, INT64_MAX) == -SCALE
    assert asert_exponent(1, INT64_MIN + 1, 0,         1) == -ASERT_EXPONENT_LIMIT
    assert calculate_asert(POW_LIMIT, 600, 600, 0, POW_LIMIT, 172800) == POW_LIMIT
    for remainder in (1, UINT64_MAX >> 31, (UINT64_MAX >> 31) + 1, INT64_MAX - 1):
        assert fractional_quotient(remainder, remainder + 1) == (remainder*SCALE)//(remainder + 1)
    for delta in (0, 44236800, -44236800):
        try:
            aserti9_nbits(ref_bits, 600 + delta, 0, pow_limit=POW_LIMIT + 1)
        except ValueError:
            continue
        raise AssertionError("Accepted unsupported powLimit")
    for invalid in (0, 0x1d80ffff, 0x23000001, 0xff7fffff):
        try:
            bits_to_target(invalid)
        except ValueError:
            continue
        raise AssertionError(f"Accepted invalid nBits: 0x{invalid:08x}")
    print("All tests passed.")


if __name__ == '__main__':
    main()