generate.py

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

3517 bytes, as of the version from 2026-10-31 20:22 (current). Recorded here, not run.

"""Bettin–Conrey cotangent sums -- numberdb.org/T639.

Run with SageMath:
    sage -pip install numberdb
    sage -python generate.py                 # verify all stored entries
    sage -python generate.py --fill          # write entries with NUMBERDB_API_KEY

The --fill operation sends entries to the existing table. It does not publish
or offer the draft. The default selection is every reduced residue h modulo k,
2 <= k <= 30, ordered by (k, h); see argument-selection.md for its rationale.

The signed convention is c0(h/k) = -sum(m/k*cot(pi*m*h/k), m=1..k-1).
Exact arithmetic in Q(zeta_lcm(4,k)) recognizes rational answers from the
algebraic definition. All other answers are computed with Arb balls at the
requested precision plus 96 bits. The enclosure, not precision agreement,
certifies the digits. Independent checks are in checks.py.
"""
import sys
from functools import lru_cache
from math import gcd, lcm
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.number_field.number_field import CyclotomicField

MAX_DENOMINATOR = 30
WORKING_GUARD = 96

@lru_cache(maxsize=None)
def cyclotomic_cotangents(k):
    n = lcm(4, k)
    field = CyclotomicField(n)
    z = field.gen()
    imaginary_unit = z ** (n // 4)
    root = z ** (n // k)
    return tuple(imaginary_unit * (root**r + 1) / (root**r - 1)
                 for r in range(1, k))

def algebraic_sum(k, h):
    cotangents = cyclotomic_cotangents(k)
    return -sum((QQ(m) / k) * cotangents[(m * h) % k - 1]
                for m in range(1, k))

def cotangent_ball(k, h, digits):
    ball = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    pi = ball.pi()
    total = ball(0)
    for m in range(1, k):
        residue = (m * h) % k
        # cot(pi/2) is exactly zero, including at larger even moduli.
        if 2 * residue == k:
            continue
        angle = pi * ball(QQ(residue) / k)
        total -= ball(QQ(m) / k) * angle.cot()
    assert total.is_finite()
    return total

class BettinConrey(numberdb.Generator):
    table = 'T639'
    parameters = ('k', 'h')
    type = 'R'
    digits = 100
    rigour = 'proven'

    def enumerate(self, maximum=MAX_DENOMINATOR):
        for k in range(2, maximum + 1):
            for h in range(1, k):
                if gcd(h, k) == 1:
                    yield {'k': k, 'h': h}

    def value(self, params, digits):
        k, h = int(params['k']), int(params['h'])
        if k < 2 or not 1 <= h < k or gcd(h, k) != 1:
            raise ValueError('Expected k >= 2 and a reduced residue 1 <= h < k')
        exact = algebraic_sum(k, h)
        if exact in QQ:
            value = QQ(exact)
            comments = {(2, 1): r'$c_0(1/2)=0$ because $\cot(\pi/2)=0$.',
                        (4, 1): r'$c_0(1/4)=1/2$ because $\cot(\pi/4)=1$ and $\cot(3\pi/4)=-1$.',
                        (4, 3): r'$c_0(3/4)=-1/2$ by reflection.'}
            return {'number': value, **({'comment': comments[k, h]} if (k, h) in comments else {})}
        value = cotangent_ball(k, h, digits)
        assert value.rad() < abs(value).lower() * (QQ(1) / 10**(digits + 15))
        return value

if __name__ == '__main__':
    generator = BettinConrey()
    if '--fill' in sys.argv:
        print(generator.publish(message='Fill signed c0 sums with certified balls and exact rational cases'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)