generate.sage

back to table · edit · history · where entries came from · files · download

3736 bytes, as of the version from 2026-09-21 12:16 (current). Recorded here, not run.

"""
Reproduce the intervals in numberdb.org/T10.

Run it with SageMath:

    $ sage -pip install numberdb
    $ sage -python generate.sage

The table stores the union of five published measurements, each widened to
five standard uncertainties. CODATA 2022 is included as a maintenance check:
adding it to the union should leave the stored intervals unchanged.
"""

import re

import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField

RBF = RealBallField(200)

MEASUREMENTS_2014_2020 = [
    ('0.0072973525664(17)', '137.035999139(31)', 'CODATA 2014'),
    ('0.0072973525657(18)', '137.035999150(33)', 'Aoyama et al. 2017'),
    ('0.0072973525713(14)', '137.035999046(27)', 'Parker et al. 2018'),
    ('0.0072973525693(11)', '137.035999084(21)', 'CODATA 2018'),
    ('0.0072973525628(6)', '137.035999206(11)', 'Morel et al. 2020'),
]

CODATA_2022 = ('0.0072973525643(11)', '137.035999177(21)', 'CODATA 2022')


def number_with_uncertainty_to_real_ball(text, standard_deviations=5):
    pattern = r'^([+-]?)(\d*)((?:\.\d*))((?:\(\d+\)))((?:[eE]-?\d+)?)$'
    match = re.match(pattern, text)
    if match is None:
        raise ValueError(text)

    sign, whole, fractional, uncertainty, exponent = match.groups()
    if fractional == '':
        fractional = '.'
    exponent_value = ZZ(0 if exponent == '' else exponent[1:])
    places = ZZ(len(fractional) - 1)
    mantissa = ZZ(int(sign + whole + fractional[1:]))
    radius = ZZ(int(uncertainty[1:-1])) * ZZ(standard_deviations)
    return RBF(mantissa, radius) * ZZ(10) ** (exponent_value - places)


def union(intervals):
    iterator = iter(intervals)
    result = next(iterator)
    for interval in iterator:
        result = result.union(interval)
    return result


def int_to_decimal(coefficient, exponent):
    if coefficient == 0:
        return '0'
    text = str(coefficient)
    digits = len(str(abs(coefficient)))
    point = digits + exponent
    if point >= 7 or point <= -5:
        return '%se%s' % (text[0] + '.' + text[1:], exponent + digits - 1)
    if point <= 0:
        return '0.%s%s' % ('0' * (-point), text)
    if point < digits:
        return text[:point] + '.' + text[point:]
    return text + '0' * exponent


def real_ball_to_string(ball, extra_digits=0):
    exponent = (ball.rad_as_ball().log() / RBF(10).log()).upper().ceil()
    exponent -= extra_digits
    scaled = ball * RBF(10) ** (-exponent)
    center = scaled.center()
    coefficient = center.floor()
    ceiling = center.ceil()
    if center - coefficient > ceiling - center:
        coefficient = ceiling
    radius = (scaled - RBF(coefficient)).abs().upper().ceil()
    return '%s +/- %s' % (
        int_to_decimal(coefficient, exponent),
        int_to_decimal(radius, exponent),
    )


def stored_interval(column):
    return union(
        number_with_uncertainty_to_real_ball(row[column])
        for row in MEASUREMENTS_2014_2020
    )


def interval_with_codata_2022(column):
    return union(
        number_with_uncertainty_to_real_ball(row[column])
        for row in MEASUREMENTS_2014_2020 + [CODATA_2022]
    )


alpha = stored_interval(0)
alpha_inv = stored_interval(1)

print('alpha:', real_ball_to_string(alpha, extra_digits=3))
print('alpha_inv:', real_ball_to_string(alpha_inv, extra_digits=3))

assert real_ball_to_string(alpha, extra_digits=3) == '0.0072973525675 +/- 1.09e-11'
assert real_ball_to_string(alpha_inv, extra_digits=3) == '137.035999113 +/- 2.03e-7'
assert real_ball_to_string(interval_with_codata_2022(0), extra_digits=3) == real_ball_to_string(alpha, extra_digits=3)
assert real_ball_to_string(interval_with_codata_2022(1), extra_digits=3) == real_ball_to_string(alpha_inv, extra_digits=3)