diff -r 2b662fbd0183 gmp-impl.h
--- a/gmp-impl.h	Tue Jul 14 20:48:57 2026 +0200
+++ b/gmp-impl.h	Fri Jul 17 11:18:34 2026 +0200
@@ -4036,6 +4036,9 @@
 #define PP_FIRST_OMITTED 3
 #endif
 
+#define mpn_sqrt_exact __MPN(sqrt_exact)
+__GMP_DECLSPEC int mpn_sqrt_exact (mp_ptr, mp_srcptr, mp_size_t);
+
 #define mpn_probab_perfect_square_p __MPN(probab_perfect_square_p)
 __GMP_DECLSPEC int mpn_probab_perfect_square_p (mp_srcptr, mp_size_t);
 
diff -r 2b662fbd0183 mpn/generic/perfsqr.c
--- a/mpn/generic/perfsqr.c	Tue Jul 14 20:48:57 2026 +0200
+++ b/mpn/generic/perfsqr.c	Fri Jul 17 11:18:34 2026 +0200
@@ -184,6 +184,8 @@
 {
   ASSERT (usize >= 1);
 
+  TRACE (gmp_printf ("mpn_perfect_square_p %Nd\n", up, usize));
+
   /* The first test excludes 212/256 (82.8%) of the perfect square candidates
      in O(1) time.  */
   {
@@ -242,8 +244,7 @@
     TMP_MARK;
     root_ptr = TMP_ALLOC_LIMBS ((usize + 1) / 2);
 
-    /* Iff mpn_sqrtrem returns zero, the square is perfect.  */
-    res = ! mpn_sqrtrem (root_ptr, NULL, up, usize);
+    res = mpn_sqrt_exact (root_ptr, up, usize);
     TMP_FREE;
 
     return res;
diff -r 2b662fbd0183 mpn/generic/sqrtrem.c
--- a/mpn/generic/sqrtrem.c	Tue Jul 14 20:48:57 2026 +0200
+++ b/mpn/generic/sqrtrem.c	Fri Jul 17 11:18:34 2026 +0200
@@ -304,6 +304,70 @@
 }
 #endif
 
+#define SQRTAPPR_MIN_SUPPORTED 5
+/* Writes in {sp, n} the square root of {np - 2n + 1, 2n}<<s.
+   The result may be up to 4 units too large.
+   Assumes {np}<<s is normalized, i.e. np[0] >= B/4
+   where B=2^GMP_NUMB_BITS.
+   Needs a scratch {tp, n + h + h/2 + 1}, should fit in 2n */
+static void
+mpn_sqrtappr (mp_ptr sp, mp_srcptr np, mp_size_t n, unsigned s, mp_ptr tp)
+{
+  mp_limb_t q, high;			/* carry out of {sp, n} */
+  int c, b;				/* carry out of remainder */
+  mp_size_t l, h;
+
+  ASSERT (n > 0);
+  /* ASSERT (n > 1); */
+  ASSERT (*np != 0);
+
+  l = n >> 1;
+  h = n - l;
+  np -= n + h;
+
+  if (s == 0)
+    {
+      MPN_COPY (tp, np + 1, n + h);
+    }
+  else
+    {
+      mpn_rshift (tp, np, n + h + 1, GMP_NUMB_BITS - 2 * s);
+      ASSERT (tp [n + h] == 0);
+    }
+  ASSERT ((tp [n + h - 1] & (GMP_NUMB_HIGHBIT | GMP_NUMB_HIGHBIT / 2)) != 0);
+
+  ASSERT (n >= SQRTAPPR_MIN_SUPPORTED);
+  if ((SQRTAPPR_MIN_SUPPORTED < 3) && (h == 1)) /* n < 3 */
+    {
+      q = CALL_SQRTREM2_INPLACE (sp + l, tp + l);
+      if ((SQRTAPPR_MIN_SUPPORTED < 2) && (l == 0)) /* n == 1 */
+	return;
+    }
+  else
+    q = mpn_dc_sqrtrem (sp + l, tp + l, h, 0, tp + n + h);
+  if (q != 0)
+    ASSERT_CARRY (mpn_sub_n (tp + l, tp + l, sp + l, h));
+
+  if ((SQRTAPPR_MIN_SUPPORTED >= 5) || (h > 2))
+    q += mpn_divappr_q (sp, tp, n, sp + l, h, tp); /* tp is used for both value and scratch */
+  else /* n < 5 */
+    {
+      mpn_div_q (tp + n + 1, tp, n, sp + l, h, tp);
+      MPN_COPY (sp, tp + n + 1, l);
+      q += tp [n + 1 + l];
+    }
+
+  if (UNLIKELY (q > 1))
+    {
+      MPN_FILL (sp, l, GMP_NUMB_MAX);
+    }
+  else
+    {
+      mpn_rshift (sp, sp, l, 1);
+      sp[l - 1] |= (q << (GMP_NUMB_BITS - 1)) & GMP_NUMB_MASK;
+    }
+}
+
 /* writes in {sp, n} the square root (rounded towards zero) of {np, 2n-odd},
    returns zero if the operand was a perfect square, one otherwise.
    Assumes {np, 2n-odd}*4^nsh is normalized, i.e. B > np[2n-1-odd]*4^nsh >= B/4
@@ -425,6 +489,340 @@
   return c;
 }
 
+#define IS_SQUARE 1
+#define IS_NOT_SQUARE 0
+/* Computes the square root of the input {np, nn} if it exists in the
+   set N of natural numbers.  If the input is a perfect square, it
+   returns IS_SQUARE and writes the root in {sp, (nn+1)/2}.  If it is
+   not, it returns IS_NOT_SQUARE.  It assumes that the limbs pointed
+   by sp are writable, and it uses them as a scratch space, this means
+   that if the input is not a perfect square, {sp, (nn+1)/2} may
+   contain random data.
+
+   If the highest limb sp[(nn+1)/2] gets overwritten, then new value
+   is non-zero.  I.e. if sp was pointed by a valid mpz_t of size
+   (nn+1)/2, after the work of this function it will keep on being a
+   valid mpz_t of the same size, even if the value will be nonsense.
+ */
+
+/* #define SMALLSIZE_OPT 0 /\* For a smaller code, but slower on smal inputs. *\/ */
+#define SMALLSIZE_OPT 1 /* To enable optimizations for small sizes */
+
+int
+mpn_sqrt_exact (mp_ptr sp, mp_srcptr np, mp_size_t nn)
+{
+  mp_limb_t high, lo;
+  unsigned clz, ctz, odd;
+
+  TRACE(printf ("sqrt_exact (,,%u)\n", (unsigned) nn));
+  for (odd = 0; UNLIKELY ((lo = *np) == CNST_LIMB (0)); ++np, --nn)
+    if ((odd ^= 1) == 0)
+      *sp++ = CNST_LIMB (0);
+  ASSERT (lo != 0);
+
+  if (SMALLSIZE_OPT && (nn == 1))
+    {
+      if (lo & (GMP_NUMB_HIGHBIT | (GMP_NUMB_HIGHBIT / 2)))
+	{
+	  if (mpn_sqrtrem1 (sp, lo) != 0)
+	    return IS_NOT_SQUARE;
+	  if (odd != 0)
+	    *sp <<= GMP_NUMB_BITS / 2;
+	  return IS_SQUARE;
+	}
+
+      count_leading_zeros (clz, lo);
+      clz -= GMP_NAIL_BITS;
+      clz >>= 1; /* we have to shift left by 2clz bits to normalize {np, nn} */
+      if (mpn_sqrtrem1 (sp, lo << clz * 2) != 0)
+	return IS_NOT_SQUARE;
+      if (odd != 0)
+	*sp <<= GMP_NUMB_BITS / 2 - clz;
+      else
+	*sp >>= clz;
+      return IS_SQUARE;
+    }
+
+  high = np[nn - 1];
+  ASSERT (high != 0);
+
+  if (high & (GMP_NUMB_HIGHBIT | (GMP_NUMB_HIGHBIT / 2)))
+    clz = 0;
+  else
+    {
+      count_leading_zeros (clz, high);
+      clz -= GMP_NAIL_BITS;
+      clz >>= 1; /* we have to shift left by 2clz bits to normalize {np, nn} */
+      if (SMALLSIZE_OPT && (nn == 2) && ((lo << clz * 2) == 0))
+	{
+	  if (mpn_sqrtrem1 (sp, (high << clz * 2) |
+			    (lo >> (GMP_NUMB_BITS - 2 * clz ))) != 0)
+	    return IS_NOT_SQUARE;
+	  if (odd != 0)
+	    {
+	      mp_limb_t sp0 = sp[0];
+	      sp[1] = sp0 >> clz;
+	      sp[0] = sp0 << GMP_NUMB_BITS - clz;
+	    }
+	  else
+	    sp[0] <<= GMP_NUMB_BITS / 2 - clz;
+	  return IS_SQUARE;
+	}
+    }
+
+  count_trailing_zeros (ctz, lo);
+#if GMP_NUMB_BITS % 2 == 1
+  ctz += ctz+odd+1 & 1;	/* increment iff (ctz+odd) is even */
+#else
+  ctz |= 1;		/* increment only if even */
+#endif
+  lo >>= ctz;		/* shift down to remove lowest non-zero bit */
+  if (((lo & 3) != 0) || (SMALLSIZE_OPT && ((np[1] & (ctz == GMP_NUMB_BITS - 1)) != 0)))
+    return IS_NOT_SQUARE;
+
+  ctz >>= 1; /* shift right by 2ctz bits {np, nn} for the odd part */
+  if (nn < 20) /* Should this be tuned? */
+    {
+      mp_size_t tn;
+      mp_ptr tp;
+      mp_bitcnt_t t;
+      mp_size_t rn;
+      TMP_SDECL;
+      TMP_SMARK;
+
+      tn = nn;
+      if (ctz != 0) {
+	tp = TMP_SALLOC_LIMBS (3 * nn + 3);
+	TRACE(printf("rshift(,,%u,%u)\n", (unsigned) tn, ctz*2));
+	mpn_rshift (tp, np, tn, ctz * 2);
+	tn -= tp [tn - 1] == 0;
+	np = tp;
+	tp += tn;
+      } else
+	tp = TMP_SALLOC_LIMBS (2 * nn + 3);
+
+      t = nn * GMP_NUMB_BITS / 2 - clz - ctz;
+      ASSERT ((t > GMP_NUMB_BITS / 2) || ! SMALLSIZE_OPT );
+      rn = 1 + t / GMP_NUMB_BITS;
+
+      MPN_FILL (tp, rn, CNST_LIMB(0));
+      ASSERT (((*np & 7) == 1) || ! SMALLSIZE_OPT);
+      /* Follows a rewrite of the function mpn_bsqrt.
+	 We should improve that implementation instead,
+	 and simply call it here.*/
+      TRACE(printf("bsqrtinv(,,%u,)\n", (unsigned) t));
+      if (! mpn_bsqrtinv (tp, np, t, tp + rn) && ! SMALLSIZE_OPT)
+	{
+	  TMP_SFREE;
+	  return IS_NOT_SQUARE;
+	};
+
+      --t;
+      t %= GMP_NUMB_BITS;
+      rn -= t == GMP_NUMB_BITS - 1;
+
+      if (SMALLSIZE_OPT && (rn == 1))
+	{
+	  mp_limb_t sp0;
+	  sp0 = *np * *tp;
+	  if (((sp0 >> t) & 1) == 0)
+	    sp0 = -sp0;
+	  sp0 &= GMP_NUMB_MAX >> GMP_NUMB_BITS - t - 1;
+	  *sp = sp0;
+	  umul_ppmm (tp[1], tp[0], sp0, sp0);
+	}
+      else
+	{
+	  TRACE(printf("mullo(,,,%u)\n", (unsigned) rn));
+	  mpn_mullo_n (sp, np, tp, rn);
+	  if (((sp [rn - 1] >> t) & 1) == 0)
+	    {
+	      *sp ^= 1; /* It's mpn_neg */
+	      mpn_com (sp, sp, rn);
+	    }
+	  sp [rn - 1] &= GMP_NUMB_MAX >> GMP_NUMB_BITS - t - 1;
+
+	  TRACE(printf("sqr(,,%u)\n", (unsigned) rn));
+	  mpn_sqr (tp, sp, rn);
+	}
+      int cmp = mpn_cmp (np + rn - 1, tp + rn - 1, tn - rn + 1) != 0;
+      TMP_SFREE;
+      if (cmp != 0)
+	return IS_NOT_SQUARE;
+
+      ctz += odd ? GMP_NUMB_BITS/2 : 0;
+      if (ctz != 0)
+	{
+	  TRACE(printf("lshift(,,%u,%u)\n", (unsigned) tn, ctz));
+	  high = mpn_lshift (sp, sp, rn, ctz);
+	  if (high != 0)
+	    sp [rn] = high;
+	}
+    }
+  else
+    {
+      mp_size_t bpart, hpart;
+      mp_ptr tp;
+      TMP_DECL;
+      TMP_MARK;
+      ASSERT (nn >= 11);
+
+      bpart = nn >> 3;
+      hpart = (nn + 1 >> 1) - bpart;
+      ++bpart;
+      ASSERT (hpart >= 5);
+      /* Scratch needed by the bpart, 3*bpart+3 */
+      /* Scratch needed by the hpart, 2*hpart */
+      ASSERT (2 * hpart >= 3 * bpart + 3);
+      tp = TMP_ALLOC_LIMBS (nn + 2);
+
+      TRACE(printf("split %u : %u [%u],\n", (unsigned) bpart, (unsigned) hpart, odd));
+      /* Bpart */
+      {
+	mp_srcptr lnp;
+	mp_ptr ltp;
+	mp_size_t tn = bpart;
+	mp_bitcnt_t t;
+
+	if (ctz != 0) {
+	  TRACE(printf("rshift(,,%u,%u)\n", (unsigned) tn + 1, ctz*2));
+	  mpn_rshift (tp, np, tn + 1, ctz * 2);
+	  lnp = tp;
+	  ltp = tp + tn;
+	} else {
+	  lnp = np;
+	  ltp = tp;
+	}
+
+	t = tn * GMP_NUMB_BITS - 1 - clz - ctz;
+	MPN_FILL (ltp, tn, CNST_LIMB(0));
+	ASSERT (((*lnp & 7) == 1) || ! SMALLSIZE_OPT);
+	TRACE(printf("bsqrt(,,%u,)\n", (unsigned) t));
+	if (! mpn_bsqrt (sp, lnp, t, tp + tn) && ! SMALLSIZE_OPT)
+	  {
+	    TMP_FREE;
+	    return IS_NOT_SQUARE;
+	  };
+      }
+      mp_limb_t bsafe = sp[--bpart] << 1 + ctz + clz;
+      bsafe |= sp[bpart - 1] >> GMP_NUMB_BITS - 1 - ctz - clz;
+      /* Bit offset of bsafe wrt sp+bpart */
+      int bbitoff = - 1 - clz + (odd ? GMP_NUMB_BITS/2 : 0);
+
+      /* Hpart */
+      ASSERT (bpart + hpart == (nn+1)/2);
+      TRACE(printf("sqrtappr(,,%u,%u,)\n", (unsigned) hpart, (unsigned) clz));
+      mpn_sqrtappr (sp + bpart + (odd & (nn ^ 1)), np + nn - 1, hpart, clz, tp);
+
+      /* Do the two parts agree for the overlapping bits?	*/
+      /* The bottom "b" part should be
+	 lshift-ed ctz + (odd ? GMP_NUMB_BITS/2 : 0) bits.
+	 The high "h" part should be
+	 rshift-ed clz + (odd ^ nn & 1 ? GMP_NUMB_BITS/2:0)
+
+	 Cases n is odd/even, odd is 1/0:
+	 odd/0 no more limbs, >> half;
+	 odd/1 no more limbs, aligned;
+	 even/0 no more limbs, aligned;
+	 even/1 1 more limb, >> half.
+
+	 odd/0, even/0, aligned;
+	 odd/1, even/1, << half.
+	*/
+
+      /* nn = 12, np = [11][10][09][08][07][06][05][04][03][02][01][00]
+	 sn =  6, sp = [05][04][03][02][01]:[00]
+	 bp =  1 {0}, hp = 5 {0}, sn = 6.
+
+	 nn = 11, np = [10][09][08][07][06][05][04][03][02][01][00]
+	 sn =  6, sp = [05][04][03][02][01]:[00]
+	 bp =  1 {0}, hp = 5{-b}, sn = 6.
+
+	 odd
+	 nn = 12, np = [11][10][09][08][07][06][05][04][03][02][01][00]
+	 sn =  6, sp = [**][05][04][03][02]:[01][00]
+	 bp =  1 {b}, hp = 5 {b}, sn = 6+1.
+
+	 nn = 11, np = [10][09][08][07][06][05][04][03][02][01][00]
+	 sn =  6, sp = [05][04][03][02][01]:[00]
+	 bp =  1 {b}, hp = 5 {0}, sn = 6.
+      */
+
+      if ((nn & 1) == 0)
+	{
+	  mp_limb_t sph0 = (sp [bpart + odd] & (GMP_NUMB_MAX >> 1)) >> 2;
+	  if ((sph0 == (~bsafe >>3)) || (sph0-1 == (~bsafe >>3)))
+	    {
+	      *sp ^= 1;
+	      TRACE(printf("com(,,%u) [%i]\n", (unsigned) bpart, (sph0-1 == (~bsafe >>3))));
+	      mpn_com (sp, sp, bpart + odd);
+	      sp [bpart + odd] |= GMP_NUMB_MAX >> 1;
+	    }
+	  else if ((sph0 == (bsafe >>3)) || (sph0-1 == (bsafe >>3)))
+	    {
+	      TRACE(printf("[%i]", (sph0-1 == (bsafe >>3))));
+	      sp [bpart + odd] &= GMP_NUMB_HIGHBIT;
+	    }
+	  else
+	    {
+	      TMP_FREE;
+	      return IS_NOT_SQUARE;
+	    };
+	  sp [bpart + odd] ^= bsafe >> 1;
+	  TRACE(printf("CORR 1, %u %u\n", (unsigned) nn & 1, (unsigned) odd));
+	}
+      else
+	{
+	  mp_limb_t sph0 = (sp [bpart] >> GMP_NUMB_BITS/2 - 1) |
+	    (sp [bpart + 1] << GMP_NUMB_BITS/2 + 1);
+	  if (sph0 == ~bsafe)
+	    {
+	      TRACE(printf("neg(,,%u)\n", (unsigned) bpart));
+	      *sp ^= 1;
+	      mpn_com (sp, sp, bpart);
+	    }
+	  else if (sph0 != bsafe)
+	    {
+	      TMP_FREE;
+	      return IS_NOT_SQUARE;
+	    };
+	}
+      if (odd ^ nn & 1)
+	clz += GMP_NUMB_BITS/2;
+      if (clz != 0)
+	{
+	  TRACE(printf("rshift(,,%u,%u)\n", (unsigned) hpart, clz));
+	  mpn_rshift (sp + bpart + (odd & (nn ^ 1)),
+		      sp + bpart + (odd & (nn ^ 1)), hpart, clz);
+	}
+      if (odd)
+	ctz += GMP_NUMB_BITS/2;
+      if (ctz != 0)
+	{
+	  TRACE(printf("lshift(,,%u,%u)\n", (unsigned) (bpart + (odd & (nn ^ 1))), ctz));
+	  mp_limb_t high = mpn_lshift (sp, sp, bpart + (odd & (nn ^ 1)), ctz);
+	  if (odd & nn)
+	    {
+	      sp[bpart] &= GMP_NUMB_MAX << ctz;
+	      sp[bpart] |= high;
+	      TRACE(printf("*CORR 2* "));
+	    }
+	}
+
+      mp_size_t sn = mpn_sqrmod_bnm1_next_size ((nn + 1 + odd) >> 1);
+      mp_size_t rn = mpn_sqrmod_bnm1_next_size ((nn + 1) >> 1);
+      mp_ptr scratch = TMP_ALLOC_LIMBS (mpn_sqrmod_bnm1_itch (rn, sn));
+      TRACE(printf("sqrmod_bnm1(,,%u)\n", (unsigned) nn));
+      mpn_sqrmod_bnm1 (tp, rn, sp, sn, scratch);
+      mp_limb_t bw = mpn_sub_n (tp + odd, tp + odd, np, rn - odd);
+      MPN_DECR_U (tp, rn, bw);
+      int cmp = mpn_cmp (np + rn - odd, tp, nn - rn + odd) != 0;
+      TMP_FREE;
+      if (cmp)
+       	return IS_NOT_SQUARE;
+    }
+  return IS_SQUARE;
+}
 
 mp_size_t
 mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nn)
diff -r 2b662fbd0183 mpz/perfsqrt.c
--- a/mpz/perfsqrt.c	Tue Jul 14 20:48:57 2026 +0200
+++ b/mpz/perfsqrt.c	Fri Jul 17 11:18:34 2026 +0200
@@ -52,7 +52,7 @@
   if (! mpn_probab_perfect_square_p (PTR (op), u_size))
     return 0;
 
-  /* This is the precise size of sqrt(op) because leading limp is non-zero. */
+  /* This is the precise size of sqrt(op) because leading limb is non-zero. */
   mp_size_t rop_size = (SIZ (op) + 1) / 2;
 
   /* mpn_sqrtrem doesn't allow sp == np */
@@ -62,7 +62,7 @@
       TMP_MARK;
 
       mp_ptr rop_ptr = TMP_ALLOC_LIMBS (rop_size);
-      int res = ! mpn_sqrtrem (rop_ptr, NULL, PTR (op), u_size);
+      int res = mpn_sqrt_exact (rop_ptr, PTR (op), u_size);
 
       /* It might be better to always copy so the value in rop doesn't
          depend on if the arguments are aliased or not. */
@@ -78,5 +78,8 @@
 
   /* Make sure rop has exact space for the square root. */
   SIZ(rop) = rop_size;
-  return ! mpn_sqrtrem (MPZ_REALLOC (rop, rop_size), NULL, PTR (op), u_size);
+
+  mp_ptr rp = MPZ_NEWALLOC (rop, rop_size);
+  rp [rop_size - 1] = 1;
+  return mpn_sqrt_exact (rp, PTR (op), u_size);
 }