from sympy.ntheory import sqrt_mod_iter

# https://oeis.org/A358016
# https://docs.sympy.org/latest/modules/ntheory.html#sympy.ntheory.residue_ntheory.sqrt_mod_iter
def A358016(n):
    return max(filter(lambda k: k<=n-2, sqrt_mod_iter(1, n)))

print (A358016(15))
print (A358016(100))

rt=0
for n in range(3, 2*10**7+1):
    rt+=A358016(n)

print (rt)