I was looking inside the basics of ecc and found the examples from Internet either uses continuous domain curve or use a very small prime number p like 17 in a discrete domain to show the points.
I am really curious that if I can find a point with a really big p in practice. For example, secp256k1 is using a really big p=2^256−2^32−977 in domain (p,a,b,G,n,h).
Below is the python code I use to deduct possible integer of y from solving equation with range of integer x. To my surprise that there is no finding even in 1 million range!
So my question is, is below code right? And secondly, if it is right or corrected by some real expert, which value range should I try?
P.S.
I am wondering how generator point G is selected, too. But that might need deeper understanding to the topic.
import math
# secp256k1
# y**2 = x**3 + 7 (mod p)
P = 2**256 - 2**32 - 977
A = 0
B = 7
# nist P256
P = 115792089210356248762697446949407573530086143415290314195533631308867097853951
A = -3
B = 41058363725152142129326129780047268409114441015993725554835256314039467401291
def in_curve(x):
    curve = x**3 + A*x + B
    y_float = math.sqrt(curve)
    if abs(math.modf(y_float)[0]) < 0.0001 or \
          (1 - abs(math.modf(y_float)[0]) < 0.0001):
        # print(y_float)
        # bug: y_int = int(math.modf(y_float)[1])
        y_int = int(round(y_float))
        if y_int * y_int == (curve):
            print(y_int)
            return y_int
    return None
 
for x in range(1, 1000000):
    y = in_curve(x)
    if y is not None:
        print(x, y)
Update 1
The previous code is wrong, since floating point sqrt() will cause unacceptable error when converting it back to integer.
But, after replacing math.sqrt() to math.isqrt(), it still doesn't make things reasonable.
Update 2
Thanks to the tips from all replied in the thread. Using generator point to verify my algorithm, I now clearly know why I was failing.
The point is, besides using % for all multiplication and additions, I should also use modular square root to find the solution, instead of integer square root along with %. That is totally a misuse of % for sure.
The modified code pass the test with some test vectors.
import modular_sqrt
#  e.g. I was using code from https://eli.thegreenplace.net/2009/03/07/computing-modular-square-roots-in-python
# please ask permission if usage is beyond educational hobbyist area and always list the credit being a decent human ;-)
# nist P256, taken from https://nvlpubs.nist.gov/nistpubs/SpecialPublications/NIST.SP.800-186-draft.pdf
P = 115792089210356248762697446949407573530086143415290314195533631308867097853951
A = -3
B = 41058363725152142129326129780047268409114441015993725554835256314039467401291
Gx = 48439561293906451759052585252797914202762949526041747995844080717082404635286
Gy = 36134250956749795798585127919587881956611106672985015071877198253568414405109
def get_y_in_curve(x):
    y2 = x**3 + A*x + B
    y_int = modular_sqrt(y2, P)
    if y_int and ((y_int * y_int) % P) == (y2 % P):
        return y_int
    return None
assert get_y_in_curve(Gx) == Gy