is_square gives incorrect answers

[charlesgriffiths]

from sympy.ntheory.primetest import is_square
from math import sqrt

print( 216000, "is square:", is_square( 216000 ) )
print( "sqrt( 216000 ) is:", sqrt( 216000 ))

216000 is square: True
sqrt( 216000 ) is: 464.75800154489

related:
from sympy.ntheory import perfect_power

print( perfect_power( 216000, [2] ))
print( perfect_power( 216000, [2], factor = False ))

(60, 3)
False

[sylee957]

I see perfect power does not give False even if the candidates are given.
So the workaround for testing e is even may work.
I don't know the reason the full algorithm had been implemented, but I have the following diff, and see if the performance improves well with python.

diff --git a/sympy/ntheory/primetest.py b/sympy/ntheory/primetest.py
index 2c808ee0ea..8a9ec2dde3 100644
--- a/sympy/ntheory/primetest.py
+++ b/sympy/ntheory/primetest.py
@@ -87,13 +87,54 @@ def is_square(n, prep=True):
             return False
         if n in [0, 1]:
             return True
+
+    # start with mod 128 rejection
     m = n & 127
-    if not ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a):
-        m = n % 63
-        if not ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008):
-            from sympy.ntheory import perfect_power
-            if perfect_power(n, [2]):
-                return True
+    if ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a):
+        return False
+
+    # Other modulii share one big int modulus
+    largeMod = n % (63L*25*11*17*19*23*31)
+
+    # residues mod 63. 75% rejection
+    m = largeMod % 63
+    if ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008):
+        return False
+
+    # residues mod 25. 56% rejection
+    m = largeMod % 25
+    if ((m*0x1929fc1b) & (m*0x4c9ea3b2) & 0x51001005):
+        return False
+
+    # residues mod 31. 48.4% rejection
+    m = 0xd10d829a*(largeMod % 31)
+    if (m & (m+0x672a5354) & 0x21025115):
+        return False
+
+    # residues mod 23. 47.8% rejection
+    m = largeMod % 23
+    if ((m*0x7bd28629) & (m*0xe7180889) & 0xf8300):
+        return False
+
+    # residues mod 19. 47.3% rejection
+    m = largeMod % 19
+    if ((m*0x1b8bead3) & (m*0x4d75a124) & 0x4280082b):
+        return False
+
+    # residues mod 17. 47.1% rejection
+    m = largeMod % 17
+    if ((m*0x6736f323) & (m*0x9b1d499) & 0xc0000300):
+        return False
+
+    # residues mod 11. 45.5% rejection
+    m = largeMod % 11
+    if ((m*0xabf1a3a7) & (m*0x2612bf93) & 0x45854000):
+        return False
+
+    from sympy.ntheory import perfect_power
+    perfect_pwr = perfect_power(n, [2])
+    if perfect_pwr and perfect_pwr[-1] % 2 == 0:
+        return True
     return False

[charlesgriffiths]

or return integer_nthroot( n, 2 )[1]

[smichr]

or return integer_nthroot( n, 2 )[1]

is_square is supposed to give a quick confirmation that something is not a square, so that might be a reason to avoid a more direct calculation of the integer square root.

[charlesgriffiths]

or return integer_nthroot( n, 2 )[1]

is_square is supposed to give a quick confirmation that something is not a square, so that might be a reason to avoid a more direct calculation of the integer square root.

You have a very good point. I ran a test of three different solutions, and which is the fastest depends on the size of n.

In place of "if perfect_power(n, [2]): return True" in the current version

Solution a: if perfect_power(n, [2], factor=False): return True
Solution b: pp = perfect_power(n, [2], big=False)
return pp and 2 == pp[1]
Solution c: return integer_nthroot( n, 2 )[1]

Offset 1
Average time taken by a: 0.08400182247161865
Average time taken by b: 0.10072704315185547
Average time taken by c: 0.07971318721771241

Offset 2^100
Average time taken by a: 0.08720964431762696
Average time taken by b: 0.11071744918823243
Average time taken by c: 0.08440198898315429

Offset 2^10000
Average time taken by a: 0.4392880868911743
Average time taken by b: 0.9468696451187134
Average time taken by c: 0.53567138671875

Offset 2^1000000 (done with a smaller sample size)
Average time taken by a: 536.8559837341309
Average time taken by b: 12332.743930816651
Average time taken by c: 680.6759834289551

Of course very large n would probably benefit significantly from trailing(n) as "if 1 == trailing(n)%2: return False" and tests involving largeMod as posted earlier.

Here is the program I used for timing.
from sympy.ntheory import perfect_power
from sympy.core.power import integer_nthroot

def is_squareA( n, prep=True ):
. if prep: n = int(n)
.. if n < 0: return False
.. if n in [0, 1]: return True

. m = n & 127
. if not ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a):
.. m = n % 63
.. if not ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008):
... if perfect_power(n, [2], factor=False): return True
. return False

def is_squareB( n, prep=True ):
. if prep: n = int(n)
.. if n < 0: return False
.. if n in [0, 1]: return True

. m = n & 127
. if not ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a):
.. m = n % 63
.. if not ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008):
... pp = perfect_power(n, [2], big=False)
... return pp and 2 == pp[1]
. return False

def is_squareC( n, prep=True ):
. if prep: n = int(n)
.. if n < 0: return False
.. if n in [0, 1]: return True

. m = n & 127
. if not ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a):
.. m = n % 63
.. if not ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008):
... return integer_nthroot( n, 2 )[1]
. return False

# import time module to calculate times
import time

# initialise lists to save the times
aTime = []
bTime = []
cTime = []

offset = 1
#offset = pow( 2, 100 )
#offset = pow( 2, 10000 )
#offset = pow( 2, 1000000 )
trials = 100000

for n in range( trials ):
. a = is_squareA( offset + n )
. b = is_squareB( offset + n )
. c = is_squareC( offset + n )
. if a != b or a != c:
.. print( "failed to be equal", a, b, c, n, (offset+n))

for k in range( 50 ):

. start = time.time()
. for n in range( trials ): is_squareA( offset + n )
. stop = time.time()
. aTime.append(stop-start)

. start = time.time()
. for n in range( trials ): is_squareB( offset + n )
. stop = time.time()
. bTime.append(stop-start)

. start = time.time()
. for n in range( trials ): is_squareC( offset + n )
. stop = time.time()
. cTime.append(stop-start)

print( "Offset", offset )
print( "Average time taken by a:", str(sum(aTime)/50))
print( "Average time taken by b:", str(sum(bTime)/50))
print( "Average time taken by c:", str(sum(cTime)/50))