Last active
September 23, 2017 18:26
-
-
Save zed/f43537c9fc9889b02382f3aaccb2d876 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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" |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment