import numpy as np
from sympy import nextprime


c = nextprime(np.random.randint(1000,2000)) * nextprime(np.random.randint(1000,2000))

a,b = 11,13

while True:
    err = c - ((a+1)//1)*((b+1)//1)
    if np.isclose(err,0, rtol=0.01):
        a = (a+1)//1
        b = (b+1)//1
        break
    
    b += 0.001 * err
    s1 = np.mean([c//i for i in range(10,10000)])
    s2 = np.mean([(a*b)//i for i in range(10,1000)])
    err = s1 - s2
    a += 0.001 * err

    print(np.sum(err)**2 ,a, b)