test_josephus3.py

#!/usr/bin/env python3
"""Closed formula (safe_position3(n)) for Josephus3 problem for 1<=n<2**31"""
from decimal import Decimal, getcontext, localcontext


def safe_position3(n, *, precision=getcontext().prec):
    # [jos]: http://user.math.uzh.ch/halbeisen/publications/pdf/jos.pdf
    # j(n, k, n - l) = (n - c_m) * k + d_m;
    # c_m <= n < c_m_plus_1
    # return j(n, 3, n) + 1;
    if n == 1 or n == 4:
        return 1
    elif n == 2 or n == 3:
        return 2
    elif n < 5:
        raise ValueError(n)
    assert n >= 5

    with localcontext() as ctx:
        ctx.prec = precision
        k = 3
        #  theorem 2 [jos] k = 3; l = 0;
        alpha = Decimal('0.8111352513650005314663662293')  # for 2147483647
        q = Decimal(k) / (k - 1)
        m = round((n / alpha).ln() / q.ln())
        for m in range(m, m - 2, -1):
            e_m = alpha * q**m
            c_m = round(e_m)
            if c_m <= n:
                break
        else:  # (c_m <= n) not found
            assert 0, "never happens"

        d_m = (e_m - c_m) < 0
        return (n - c_m) * k + d_m + 1


def safe_position(n, k):
    # https://en.wikipedia.org/wiki/Josephus_problem#The_general_case
    # O(n): g(n, k) = (g(n-1, k) + k) % n; g(1, k) = 0
    g = 0
    for i in range(1, n + 1):
        g = (g + k) % i
    return g + 1


def jos(n, k, i):
    # http://user.math.uzh.ch/halbeisen/publications/pdf/jos.pdf
    assert n >= 1 and k >= 1 and 1 <= i <= n
    if i == 1:
        return (k - 1) % n
    return (k + jos(n - 1, k, i - 1)) % n


def joseph1(n, k, i):
    # http://geofhagopian.net/images/Josephus.htm
    x = k * i
    while x > n:
        x = (k * (x - n) - 1) // (k - 1)
    return x


def test_safe_position3(precision):
    k = 3
    for n in range(1, 2**31):
        expected = joseph1(n, k, n)
        assert (jos(n, k, n) + 1) == expected   # O(n)
        assert safe_position(n, k) == expected  # O(n)
        got = safe_position3(n, precision=precision)
        assert got == expected, (n, precision, expected, got)
        if n == 9103:
            print("n={n} works with precision={precision}".format(**vars()))


if __name__ == '__main__':
    import sys
    sys.setrecursionlimit(10**6)
    for precision in [1, 2, 3, 4, 5, 6, 7, 10, 100, 1000, 10000, 100000, 10**6]:
        try:
            test_safe_position3(precision=precision)
        except AssertionError as e:
            print(e)
        else:
            assert 0, "never happens"