def isqrt(n):
if n < 0:
raise ValueError
("x must be non-negative") if n == 0:
return 0
x = 2 ** ((n.bit_length() - 1) // 2 + 1) # x > sqrt(n)
y = x - 1
while y < x:
x = y
y = (x + n//x) // 2
return x
# Quadratic residues mod
d64 = set((0, 1, 4, 9, 16, 17, 25, 33, 36, 41, 49, 57))
d63 = set((0, 1, 4, 7, 9, 16, 18, 22, 25, 28, 36, 37, 43, 46, 49, 58))
d25 = set((0, 1, 4, 6, 9, 11, 14, 16, 19, 21, 24))
d31 = set((0, 1, 2, 4, 5, 7, 8, 9, 10, 14, 16, 18, 19, 20, 25, 28))
d23 = set((0, 1, 2, 3, 4, 6, 8, 9, 12, 13, 16, 18))
d19 = set((0, 1, 4, 5, 6, 7, 9, 11, 16, 17))
d17 = set((0, 1, 2, 4, 8, 9, 13, 15, 16))
d11 = set((0, 1, 3, 4, 5, 9))
def is_square(n):
return (n % 64 in d64 and n % 63 in d63 and n % 25 in d25 and
n % 31 in d31 and n % 23 in d23 and n % 19 in d19 and
n % 17 in d17 and n % 11 in d11 and isqrt(n)**2 == n)
a64 = [1 if i in d64 else 0 for i in range(64)]
a63 = [1 if i in d63 else 0 for i in range(63)]
a25 = [1 if i in d25 else 0 for i in range(25)]
a31 = [1 if i in d31 else 0 for i in range(31)]
a23 = [1 if i in d23 else 0 for i in range(23)]
a19 = [1 if i in d19 else 0 for i in range(19)]
a17 = [1 if i in d17 else 0 for i in range(17)]
a11 = [1 if i in d11 else 0 for i in range(11)]
def is_square2(n):
return (a64[n % 64] and a63[n % 63] and a25[n % 25] and a31[n % 31] and
a23[n % 23] and a19[n % 19] and a17[n % 17] and a11[n % 11] and
isqrt(n)**2 == n)
def isSquare(n):
if 144680414395695635 >> (n%64) & 1 == 0: return False
if 288872697407275667 >> (n%63) & 1 == 0: return False
if 845593731871251 >> (n%52) & 1 == 0: return False
if 615814660213299 >> (n%55) & 1 == 0: return False
if 54655143861879507 >> (n%57) & 1 == 0: return False
if 576239692522003 >> (n%51) & 1 == 0: return False
s = isqrt(n)
return s*s == n
def isSquare(n):
return (144680414395695635 >> (n%64) & 1 and
288872697407275667 >> (n%63) & 1 and
845593731871251 >> (n%52) & 1 and
615814660213299 >> (n%55) & 1 and
54655143861879507 >> (n%57) & 1 and
576239692522003 >> (n%51) & 1 and
isqrt(n) ** 2 == n)
M = (63*25*11*17*19*23*31)
def isSquare2(n): # fenderbender
m = n & 127
if ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a): return False
largeMod = n % M
m = largeMod % 63
if ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008): return False
m = largeMod % 25
if ((m*0x1929fc1b) & (m*0x4c9ea3b2) & 0x51001005): return False
m = 0xd10d829a * (largeMod % 31)
if (m & (m+0x672a5354) & 0x21025115): return False
m = largeMod % 23
if ((m*0x7bd28629) & (m*0xe7180889) & 0xf8300): return False
m = largeMod % 19
if ((m*0x1b8bead3) & (m*0x4d75a124) & 0x4280082b): return False
m = largeMod % 17
if ((m*0x6736f323) & (m*0x9b1d499) & 0xc0000300): return False
m = largeMod % 11
if ((m*0xabf1a3a7) & (m*0x2612bf93) & 0x45854000): return False
s = isqrt(n); return s*s == n
if __name__ == '__main__':
def bm(meth):
cnt = 0
N = 10 ** 6
start = 10 ** 10
for n in range(start, start + N):
if meth(n):
cnt += 1
print
(meth.__name__
, clock() - t0
, cnt
)
bm(isSquare)
bm(isSquare2)
bm(is_square)
bm(is_square2)
"""
N = 10^7
CPython3.6
isSquare 7.860324823604874 50
isSquare2 11.401496343810134 50
is_square 5.880065683130212 50
is_square2 6.948259140474661 50
PyPy3.6
isSquare 2.482781270044901 50
isSquare2 4.008553688261706 50
is_square 2.052393072482361 50
is_square2 1.6373788325849894 50
"""
ZGVmIGlzcXJ0KG4pOgogICAgaWYgbiA8IDA6CiAgICAgICAgcmFpc2UgVmFsdWVFcnJvcigieCBtdXN0IGJlIG5vbi1uZWdhdGl2ZSIpCiAgICBpZiBuID09IDA6CiAgICAgICAgcmV0dXJuIDAKICAgIHggPSAyICoqICgobi5iaXRfbGVuZ3RoKCkgLSAxKSAvLyAyICsgMSkgICAgIyB4ID4gc3FydChuKQogICAgeSA9IHggLSAxCiAgICB3aGlsZSB5IDwgeDoKICAgICAgICB4ID0geQogICAgICAgIHkgPSAoeCArIG4vL3gpIC8vIDIKICAgIHJldHVybiB4CgojIFF1YWRyYXRpYyByZXNpZHVlcyBtb2QKZDY0ID0gc2V0KCgwLCAxLCA0LCA5LCAxNiwgMTcsIDI1LCAzMywgMzYsIDQxLCA0OSwgNTcpKQpkNjMgPSBzZXQoKDAsIDEsIDQsIDcsIDksIDE2LCAxOCwgMjIsIDI1LCAyOCwgMzYsIDM3LCA0MywgNDYsIDQ5LCA1OCkpCmQyNSA9IHNldCgoMCwgMSwgNCwgNiwgOSwgMTEsIDE0LCAxNiwgMTksIDIxLCAyNCkpCmQzMSA9IHNldCgoMCwgMSwgMiwgNCwgNSwgNywgOCwgOSwgMTAsIDE0LCAxNiwgMTgsIDE5LCAyMCwgMjUsIDI4KSkKZDIzID0gc2V0KCgwLCAxLCAyLCAzLCA0LCA2LCA4LCA5LCAxMiwgMTMsIDE2LCAxOCkpCmQxOSA9IHNldCgoMCwgMSwgNCwgNSwgNiwgNywgOSwgMTEsIDE2LCAxNykpCmQxNyA9IHNldCgoMCwgMSwgMiwgNCwgOCwgOSwgMTMsIDE1LCAxNikpCmQxMSA9IHNldCgoMCwgMSwgMywgNCwgNSwgOSkpCgpkZWYgaXNfc3F1YXJlKG4pOgogICAgcmV0dXJuIChuICUgNjQgaW4gZDY0IGFuZCBuICUgNjMgaW4gZDYzIGFuZCBuICUgMjUgaW4gZDI1IGFuZAogICAgICAgICAgICBuICUgMzEgaW4gZDMxIGFuZCBuICUgMjMgaW4gZDIzIGFuZCBuICUgMTkgaW4gZDE5IGFuZAogICAgICAgICAgICBuICUgMTcgaW4gZDE3IGFuZCBuICUgMTEgaW4gZDExIGFuZCBpc3FydChuKSoqMiA9PSBuKQoKYTY0ID0gWzEgaWYgaSBpbiBkNjQgZWxzZSAwIGZvciBpIGluIHJhbmdlKDY0KV0KYTYzID0gWzEgaWYgaSBpbiBkNjMgZWxzZSAwIGZvciBpIGluIHJhbmdlKDYzKV0KYTI1ID0gWzEgaWYgaSBpbiBkMjUgZWxzZSAwIGZvciBpIGluIHJhbmdlKDI1KV0KYTMxID0gWzEgaWYgaSBpbiBkMzEgZWxzZSAwIGZvciBpIGluIHJhbmdlKDMxKV0KYTIzID0gWzEgaWYgaSBpbiBkMjMgZWxzZSAwIGZvciBpIGluIHJhbmdlKDIzKV0KYTE5ID0gWzEgaWYgaSBpbiBkMTkgZWxzZSAwIGZvciBpIGluIHJhbmdlKDE5KV0KYTE3ID0gWzEgaWYgaSBpbiBkMTcgZWxzZSAwIGZvciBpIGluIHJhbmdlKDE3KV0KYTExID0gWzEgaWYgaSBpbiBkMTEgZWxzZSAwIGZvciBpIGluIHJhbmdlKDExKV0KCmRlZiBpc19zcXVhcmUyKG4pOgogICAgcmV0dXJuIChhNjRbbiAlIDY0XSBhbmQgYTYzW24gJSA2M10gYW5kIGEyNVtuICUgMjVdIGFuZCBhMzFbbiAlIDMxXSBhbmQKICAgICAgICAgICAgYTIzW24gJSAyM10gYW5kIGExOVtuICUgMTldIGFuZCBhMTdbbiAlIDE3XSBhbmQgYTExW24gJSAxMV0gYW5kCiAgICAgICAgICAgIGlzcXJ0KG4pKioyID09IG4pCgpkZWYgaXNTcXVhcmUobik6CiAgICBpZiAxNDQ2ODA0MTQzOTU2OTU2MzUgPj4gKG4lNjQpICYgMSA9PSAwOiByZXR1cm4gRmFsc2UKICAgIGlmIDI4ODg3MjY5NzQwNzI3NTY2NyA+PiAobiU2MykgJiAxID09IDA6IHJldHVybiBGYWxzZQogICAgaWYgICAgODQ1NTkzNzMxODcxMjUxID4+IChuJTUyKSAmIDEgPT0gMDogcmV0dXJuIEZhbHNlCiAgICBpZiAgICA2MTU4MTQ2NjAyMTMyOTkgPj4gKG4lNTUpICYgMSA9PSAwOiByZXR1cm4gRmFsc2UKICAgIGlmICA1NDY1NTE0Mzg2MTg3OTUwNyA+PiAobiU1NykgJiAxID09IDA6IHJldHVybiBGYWxzZQogICAgaWYgICAgNTc2MjM5NjkyNTIyMDAzID4+IChuJTUxKSAmIDEgPT0gMDogcmV0dXJuIEZhbHNlCiAgICBzID0gaXNxcnQobikKICAgIHJldHVybiBzKnMgPT0gbgoKZGVmIGlzU3F1YXJlKG4pOgogICAgcmV0dXJuICgxNDQ2ODA0MTQzOTU2OTU2MzUgPj4gKG4lNjQpICYgMSBhbmQKICAgICAgICAgICAgMjg4ODcyNjk3NDA3Mjc1NjY3ID4+IChuJTYzKSAmIDEgYW5kCiAgICAgICAgICAgIDg0NTU5MzczMTg3MTI1MSA+PiAobiU1MikgJiAxIGFuZAogICAgICAgICAgICA2MTU4MTQ2NjAyMTMyOTkgPj4gKG4lNTUpICYgMSBhbmQKICAgICAgICAgICAgNTQ2NTUxNDM4NjE4Nzk1MDcgPj4gKG4lNTcpICYgMSBhbmQKICAgICAgICAgICAgNTc2MjM5NjkyNTIyMDAzID4+IChuJTUxKSAmIDEgYW5kCiAgICAgICAgICAgIGlzcXJ0KG4pICoqIDIgPT0gbikKCk0gPSAoNjMqMjUqMTEqMTcqMTkqMjMqMzEpCgpkZWYgaXNTcXVhcmUyKG4pOiAjIGZlbmRlcmJlbmRlcgogIG0gPSBuICYgMTI3CiAgaWYgKChtKjB4OGJjNDBkN2QpICYgKG0qMHhhMWUyZjVkMSkgJiAweDE0MDIwYSk6IHJldHVybiBGYWxzZQogIGxhcmdlTW9kID0gbiAlIE0KICBtID0gbGFyZ2VNb2QgJSA2MwogIGlmICgobSoweDNkNDkxZGY3KSAmIChtKjB4YzgyNGE5ZjkpICYgMHgxMGYxNDAwOCk6IHJldHVybiBGYWxzZQogIG0gPSBsYXJnZU1vZCAlIDI1CiAgaWYgKChtKjB4MTkyOWZjMWIpICYgKG0qMHg0YzllYTNiMikgJiAweDUxMDAxMDA1KTogcmV0dXJuIEZhbHNlCiAgbSA9IDB4ZDEwZDgyOWEgKiAobGFyZ2VNb2QgJSAzMSkgCiAgaWYgKG0gJiAobSsweDY3MmE1MzU0KSAmIDB4MjEwMjUxMTUpOiByZXR1cm4gRmFsc2UKICBtID0gbGFyZ2VNb2QgJSAyMwogIGlmICgobSoweDdiZDI4NjI5KSAmIChtKjB4ZTcxODA4ODkpICYgMHhmODMwMCk6IHJldHVybiBGYWxzZQogIG0gPSBsYXJnZU1vZCAlIDE5CiAgaWYgKChtKjB4MWI4YmVhZDMpICYgKG0qMHg0ZDc1YTEyNCkgJiAweDQyODAwODJiKTogcmV0dXJuIEZhbHNlCiAgbSA9IGxhcmdlTW9kICUgMTcgCiAgaWYgKChtKjB4NjczNmYzMjMpICYgKG0qMHg5YjFkNDk5KSAmIDB4YzAwMDAzMDApOiByZXR1cm4gRmFsc2UKICBtID0gbGFyZ2VNb2QgJSAxMSAKICBpZiAoKG0qMHhhYmYxYTNhNykgJiAobSoweDI2MTJiZjkzKSAmIDB4NDU4NTQwMDApOiByZXR1cm4gRmFsc2UKICBzID0gaXNxcnQobik7IHJldHVybiBzKnMgPT0gbgoKCmlmIF9fbmFtZV9fID09ICdfX21haW5fXyc6CiAgICBmcm9tIHRpbWUgaW1wb3J0IGNsb2NrCgogICAgZGVmIGJtKG1ldGgpOgogICAgICAgIHQwID0gY2xvY2soKQogICAgICAgIGNudCA9IDAKICAgICAgICBOID0gMTAgKiogNgogICAgICAgIHN0YXJ0ID0gMTAgKiogMTAKICAgICAgICBmb3IgbiBpbiByYW5nZShzdGFydCwgc3RhcnQgKyBOKToKICAgICAgICAgICAgaWYgbWV0aChuKToKICAgICAgICAgICAgICAgIGNudCArPSAxCiAgICAgICAgcHJpbnQobWV0aC5fX25hbWVfXywgY2xvY2soKSAtIHQwLCBjbnQpCgogICAgYm0oaXNTcXVhcmUpCiAgICBibShpc1NxdWFyZTIpCiAgICBibShpc19zcXVhcmUpCiAgICBibShpc19zcXVhcmUyKQoKIiIiCk4gPSAxMF43CkNQeXRob24zLjYKaXNTcXVhcmUgNy44NjAzMjQ4MjM2MDQ4NzQgNTAKaXNTcXVhcmUyIDExLjQwMTQ5NjM0MzgxMDEzNCA1MAppc19zcXVhcmUgNS44ODAwNjU2ODMxMzAyMTIgNTAKaXNfc3F1YXJlMiA2Ljk0ODI1OTE0MDQ3NDY2MSA1MAoKUHlQeTMuNgppc1NxdWFyZSAyLjQ4Mjc4MTI3MDA0NDkwMSA1MAppc1NxdWFyZTIgNC4wMDg1NTM2ODgyNjE3MDYgNTAKaXNfc3F1YXJlIDIuMDUyMzkzMDcyNDgyMzYxIDUwCmlzX3NxdWFyZTIgMS42MzczNzg4MzI1ODQ5ODk0IDUwCiIiIgogICAg