diff --git a/gmp-impl.h b/gmp-impl.h
index 8f83aa82..ffbc0390 100644
--- a/gmp-impl.h
+++ b/gmp-impl.h
@@ -784,6 +784,8 @@ void __gmp_default_free _PROTO ((void *, size_t));
 
 int is_likely_prime_BPSW(mp_limb_t n);
 
+mp_limb_t n_sqrt(mp_limb_t r);
+
 void __gmpz_aorsmul_1 _PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr u, mp_limb_t v, mp_size_t sub))) REGPARM_ATTR(1);
 #define mpz_aorsmul_1(w,u,v,sub)  __gmpz_aorsmul_1 (REGPARM_3_1 (w, u, v, sub))
 
diff --git a/mpz/likely_prime_p.c b/mpz/likely_prime_p.c
index fdbaf3b5..35db8741 100644
--- a/mpz/likely_prime_p.c
+++ b/mpz/likely_prime_p.c
@@ -1,7 +1,7 @@
 /*
 Copyright 2009 Jason Moxham
-Copyright (C) 2008 Peter Shrimpton
-Copyright (C) 2008, 2009 William Hart
+Copyright (C) 2008 Peter Shrimpton
+Copyright (C) 2008, 2009 William Hart
 
 This file is part of the MPIR Library.
 
@@ -52,16 +52,36 @@ mp_limb_t n_sqrt(mp_limb_t r)
       mp_limb_t l;
 	} temp;
 
+   mp_limb_t bits32 = (r & GMP_LIMB_HIGHBIT);
+   mp_limb_t r2;
+   // algorithm can't handle 32 bits
+   if (bits32) 
+   {
+      r2 = r;
+      r >>= 2;
+   }
+
 	temp.f = (float) r;
 	temp.l = (CNST_LIMB(0xbe6ec85e) - temp.l)>>1; // estimate of 1/sqrt(y) 
 	x =  temp.f;
 	z =  (float) r*0.5;                        
-    x = (1.5*x) - (x*x)*(x*z);
+   x = (1.5*x) - (x*x)*(x*z);
 	x = (1.5*x) - (x*x)*(x*z);
 	x = (1.5*x) - (x*x)*(x*z);
 	x = (1.5*x) - (x*x)*(x*z);
-    res =  is + ((is+1)*(is+1) <= r);
-    return res - (res*res > r);
+   res =  is + ((is+1)*(is+1) <= r);
+   if (!bits32) return res - (res*res > r);
+   else 
+   {
+      mp_limb_t sq;
+      res = res - (res*res > r);
+      res <<= 1;
+      sq = res*res;
+      res = res - ((sq > r2) || ((sq ^ r2) & GMP_LIMB_HIGHBIT));
+      sq = (res + 1)*(res + 1);
+      res = res + ((sq <= r2) && !((sq ^ r2) & GMP_LIMB_HIGHBIT));
+      return res;
+   }
 #else
 	
 	double x, z;
@@ -70,18 +90,38 @@ mp_limb_t n_sqrt(mp_limb_t r)
       mp_limb_t l;
 	} temp;
 
+   mp_limb_t bits64 = (r & GMP_LIMB_HIGHBIT);
+   mp_limb_t r2;
+   // algorithm can't handle 64 bits
+   if (bits64) 
+   {
+      r2 = r;
+      r >>= 2;
+   }
+
 	temp.f = (double) r;
 	temp.l = (CNST_LIMB(0xbfcdd90a00000000) - temp.l)>>1; // estimate of 1/sqrt(y) 
 	x =  temp.f;
 	z =  (double) r*0.5;                        
-    x = (1.5*x) - (x*x)*(x*z);
+   x = (1.5*x) - (x*x)*(x*z);
 	x = (1.5*x) - (x*x)*(x*z);
 	x = (1.5*x) - (x*x)*(x*z);
 	x = (1.5*x) - (x*x)*(x*z);
 	x = (1.5*x) - (x*x)*(x*z);
-    is = (mp_limb_t) (x*(double) r);
-    res =  is + ((is+1)*(is+1) <= r);
-    return res - (res*res > r);
+   is = (mp_limb_t) (x*(double) r);
+   res =  is + ((is+1)*(is+1) <= r);
+   if (!bits64) return res - (res*res > r);
+   else 
+   {
+      mp_limb_t sq;
+      res = res - (res*res > r);
+      res <<= 1;
+      sq = res*res;
+      res = res - ((sq > r2) || ((sq ^ r2) & GMP_LIMB_HIGHBIT));
+      sq = (res + 1)*(res + 1);
+      res = res + ((sq <= r2) && !((sq ^ r2) & GMP_LIMB_HIGHBIT));
+      return res;
+   }
 #endif
 }
 
diff --git a/tests/mpz/t-likely_prime_p.c b/tests/mpz/t-likely_prime_p.c
index 1c469008..4bf7c5f1 100644
--- a/tests/mpz/t-likely_prime_p.c
+++ b/tests/mpz/t-likely_prime_p.c
@@ -26,6 +26,50 @@ MA 02110-1301, USA. */
 #include "gmp-impl.h"
 #include "tests.h"
 
+void
+check_sqrt (void)
+{
+  gmp_randstate_t  rands;
+  int    i;
+  mpz_t  x;
+  unsigned long bits;
+  mp_limb_t p, m, m2, s;
+
+  mpz_init (x);
+  gmp_randinit_default(rands);
+
+  for (i = 0; i < 2000000; i++)
+    {
+      do
+      {
+         mpz_set_ui(x, GMP_NUMB_BITS/2 - 1);
+         mpz_urandomm(x, rands, x);
+         bits = mpz_get_ui(x) + 2;
+         mpz_rrandomb(x, rands, bits);
+         p = mpz_getlimbn(x, 0);
+         m = p*p;
+         mpz_set_ui(x, MIN(p, 1000));
+         mpz_urandomm(x, rands, x);
+         m2 = m + mpz_get_ui(x) - MIN(p, 1000)/2;
+      } while ((m2 == m) || ((mp_limb_signed_t) (m2 ^ m) < (mp_limb_signed_t) 0) || (m2 == 0) || (m2 == 1));
+
+      s = n_sqrt(m2);
+      if (((m2 < m) && (s != p - 1)) || ((m2 > m) && (s != p)))
+      {
+          printf ("mpz_likely_prime_p\n");
+          printf ("n_sqrt is broken\n");
+#if defined( _MSC_VER ) && defined( _WIN64 )
+          printf ("%llu is returned as n_sqrt(%llu)\n", s, m2);
+#else
+          printf ("%lu is returned as n_sqrt(%lu)\n", s, m2);
+#endif
+          abort();
+      }
+    }
+  mpz_clear (x);
+  gmp_randclear(rands);
+}
+
 void
 check_rand (void)
 {
@@ -38,7 +82,7 @@ check_rand (void)
   mpz_init (x);
   gmp_randinit_default(rands);
 
-  for (i = 0; i < 200000; i++)
+  for (i = 0; i < 2000; i++)
     {
       do
       {
@@ -141,6 +185,7 @@ main (void)
 {
   tests_start ();
   
+  check_sqrt ();
   check_rand ();
   check_small ();