Finding the biggest $k\leq k_0$ such that the interval $k[A, B]$ contains an integer

Question

Question

Let $0<A<B$ be rational numbers (with possibly large denominators) and $k_0$ a natural number. What is the biggest $k \in \{0,1,\dots, k_0\}$ such that the interval $[kA, kB]$ contains an integer? (Notice that the maximum exists, since at least $k=0$ works.) Is there a better algorithm than just testing every value from $k_0$ downwards?

Background

If it's any help, the context in which this arises is finding the "next smaller $c$" for which the Diophantine equation

$$2ay + 2bx = c \tag{1}$$

has a non-negative solution ($y$ comes first for another reason, I'll keep them like that so I don't mess up with notation). I'm taking a particular solution $(y_0, x_0)$ to $2ay+2bx=g_0:=\gcd(2a, 2b)$. The general solution to $(1)$ with $c=kg_0$ is then

$$(ky_0-\frac{2b}{g_0}t, kx_0+\frac{2a}{g_0}t), \space t\in\mathbb{Z}.$$

In order for it to be non-negative, I need (assuming $x_0<0$ and $y_0>0$) $$\frac{k(-x_0)g_0}{2a} \leq t \leq \frac{ky_0g_0}{2b}.$$ (We can see they are in that order by putting $-x_0 = \frac{2ay_0-g_0}{2b}$ so in fact the interval is $\frac{kg_0}{2b}[y_0-\frac{g_0}{2a}, y_0]$)

So I arrive to the question of this post with $A=\frac{-x_0g_0}{2a}$ and $B=\frac{y_0g_0}{2b}$ and $k_0 = \left \lfloor \frac{c}{g_0} \right \rfloor$.

My tries

I have been thinking that $[A, B]$ is a $1$-dimensional rational polytope, so an Ehrhart quasi-polynomial $q(k)$ would count the lattice points inside $k[A,B]$. But I don't know how that helps. How first of all to find $q$ (the coefficient of the constant term could have period up to $\frac{4ab}{g_0^2}$) and then how to find the value of $k\leq k_0$ for which $q(k)<q(k+1)$.

I also thought about bisection search, but I don't see how that could work, since the number of integers on the interval is erratic. Here's an example at Desmos. And here's what I've been doodling about the original problem .

Answer 1

Let $A=\frac{p_A}{q_A}$ and $B=\frac{p_B}{q_B}$.

The number of integers on the interval $[kA, kB]$ is $\lfloor kB \rfloor - \lfloor kA \rfloor$. Let's define the function

$$S(m) = \sum_{k=0}^m \lfloor kB \rfloor - \lfloor kA \rfloor$$

Using Hang Wu's method to calculate $\sum_{k=0}^m \left \lfloor \frac{kp}{q} \right\rfloor$, we can calculate $S(m)$ in $O(\log(\max(q_A, q_B)))$-time. Now do a bisection search for when $S(k-1)<S(k_0)$ (or if $k_0$ works, return it). So we get a $\log ^2$ algorithm.

Here's the recursive version of the function $f(a,b,c,n) = \sum_{k=0}^n \left \lfloor \frac{ak+b}{c} \right\rfloor$ and the bisection search for $k$ coded in Sage (hope I don't have bugs)

#f(a,b,c,n) = sum(floor((ak+b)/c) for k in range(n+1))
#assumes a,b,n >= 0 and c>0
def f(a,b,c,n):
    if n<0: return 0
    if n==0: return floor(b/c)
    if a>=c or b>=c:
        q1, r1 = floor(a/c), a%c
        q2, r2 = floor(b/c), b%c
        return n*(n+1)/2*q1 + (n+1)*q2 + f(r1, r2, c, n)
    else:
        v = floor((a*n+b)/c)
        return n*v - f(c, c-b-1, a, v-1)
    

#S(m) = sum(floor(k*B)-floor(k*A) for k in range(m+1))
#findK(A,B,k0) = max(k for k in range(k0+1) if S(k-1)<S(k))
def findK(A,B,k0):
    hasInts = lambda k: floor(k*B)-floor(k*A)>0
    S = lambda m: f(B.numerator(), 0, B.denominator(), m) - f(A.numerator(), 0, A.denominator(), m)
    S0 = S(k0)
    k1 = 0
    k2 = k0
    while True:
        if hasInts(k2): return k2
        if k2-k1<2: return k1
        mid = floor((k1+k2)/2)
        if S(mid-1)<S0:
            k1 = mid
        else:
            k2 = mid
    return 0

I will also store here the $T$-function involved in the original question

#T(a,b,c) = sum(1 for x in range(0,a+1,2) for y in range(0,b+1,2) if a*y+b*x<=c)
def T(a,b,c):
    n = floor(c/(2*b))
    return (n+1)*(1-1/2*(floor(b/a)+1)*n) + f(2*(a-(b%a)), c, 2*a, n)

And here's a test I did (I printed them to see that I'm testing some interesting cases)

def findKBrute(A,B,k0):
    for k in range(k0, 0, -1):
        if floor(k*B)-floor(k*A)>0: return k
    return 0

import random
def randRat(minQ, maxQ, minVal=0, maxVal=1):
    q = random.randint(minQ, maxQ)
    pMin = floor(minVal*q)+1
    pMax = max(pMin, floor(maxVal*q)-1)
    p = random.randint(pMin, pMax)
    return Integer(p)/q

for _ in range(10):
    A = randRat(2, 1000)
    B = randRat(2, 1000, A, A+0.01)
    #print ("A,B = ", A, B)
    for k0 in range(10, 90):
        solBrute = findKBrute(A,B,k0)
        sol = findK(A,B,k0)
        #print ("k0 = ",k0, "sols = ", solBrute, sol)
        if solBrute != sol:
            print("wrong for A=%s, B=%s, k0=%d, sols = %d, %d" %(A,B,k0,solBrute,sol))
 
print ("Test over")