Adopt to GNU coding standards.

(mpn_dc_sqrtrem): New name for mpn_dq_sqrtrem.
Partial nailification.

1 files changed, 90 insertions, 69 deletions

diff --git a/mpn/generic/sqrtrem.c b/mpn/generic/sqrtrem.c
index af24d2497..1e469bfca 100644
--- a/mpn/generic/sqrtrem.c
+++ b/mpn/generic/sqrtrem.c
@@ -34,25 +34,27 @@ MA 02111-1307, USA.
 #include "longlong.h"
 
 
-/* #define DEBUG */
-
-
-/* table generated by isqrt(2^8*i) $ i=64..255 in MuPAD */
-static const unsigned char approx_tab[192] = {
-128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141,
-   142, 143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153,
-   154, 155, 155, 156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164,
-   165, 166, 167, 167, 168, 169, 170, 170, 171, 172, 173, 173, 174, 175,
-   176, 176, 177, 178, 178, 179, 180, 181, 181, 182, 183, 183, 184, 185,
-   185, 186, 187, 187, 188, 189, 189, 190, 191, 192, 192, 193, 193, 194,
-   195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201, 202, 203, 203,
-   204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211, 212,
-   212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220,
-   221, 221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228,
-   229, 229, 230, 230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236,
-   236, 237, 237, 238, 238, 239, 240, 240, 241, 241, 242, 242, 243, 243,
-   244, 244, 245, 245, 246, 246, 247, 247, 248, 248, 249, 249, 250, 250,
-   251, 251, 252, 252, 253, 253, 254, 254, 255};
+
+/* Square roots table.  Generated by the following program:
+#include "gmp.h"
+main(){mpz_t x;int i;mpz_init(x);for(i=64;i<256;i++){mpz_set_ui(x,256*i);
+mpz_sqrt(x,x);mpz_out_str(0,10,x);printf(",");if(i%16==15)printf("\n");}}
+*/
+static const unsigned char approx_tab[192] =
+  {
+    128,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,
+    143,144,144,145,146,147,148,149,150,150,151,152,153,154,155,155,
+    156,157,158,159,160,160,161,162,163,163,164,165,166,167,167,168,
+    169,170,170,171,172,173,173,174,175,176,176,177,178,178,179,180,
+    181,181,182,183,183,184,185,185,186,187,187,188,189,189,190,191,
+    192,192,193,193,194,195,195,196,197,197,198,199,199,200,201,201,
+    202,203,203,204,204,205,206,206,207,208,208,209,209,210,211,211,
+    212,212,213,214,214,215,215,216,217,217,218,218,219,219,220,221,
+    221,222,222,223,224,224,225,225,226,226,227,227,228,229,229,230,
+    230,231,231,232,232,233,234,234,235,235,236,236,237,237,238,238,
+    239,240,240,241,241,242,242,243,243,244,244,245,245,246,246,247,
+    247,248,248,249,249,250,250,251,251,252,252,253,253,254,254,255
+  };
 
 
 /* same as mpn_sqrtrem, but for size=1 and {np, 1} normalized */
@@ -62,17 +64,21 @@ mpn_sqrtrem1 (mp_ptr sp, mp_ptr rp, mp_srcptr np)
   mp_limb_t np0, s, r, q, u;
   int prec;
 
-  ASSERT (np[0] >= MP_LIMB_T_HIGHBIT/2);
+  ASSERT (np[0] >= GMP_NUMB_HIGHBIT / 2);
 
   /* first compute a 8-bit approximation from the high 8-bits of np[0] */
   np0 = np[0];
-  q = np0 >> (BITS_PER_MP_LIMB - 8);
+  q = np0 >> (GMP_NUMB_BITS - 8);
   /* 2^6 = 64 <= q < 256 = 2^8 */
-  s = approx_tab[q - 64]; /* 128 <= s < 255 */
-  r = (np0 >> (BITS_PER_MP_LIMB - 16)) - s * s; /* r <= 2*s */
-  if (r > 2 * s) { r -= 2 * s + 1; s++; }
+  s = approx_tab[q - 64];				/* 128 <= s < 255 */
+  r = (np0 >> (GMP_NUMB_BITS - 16)) - s * s;		/* r <= 2*s */
+  if (r > 2 * s)
+    {
+      r -= 2 * s + 1;
+      s++;
+    }
 #ifdef DEBUG
-  if (r > 2*s)
+  if (r > 2 * s)
     {
       fprintf(stderr, "Error: r > 2*s, np[0]=%lu r=%lu s=%lu\n", np[0], r, s);
       exit(1);
@@ -81,16 +87,16 @@ mpn_sqrtrem1 (mp_ptr sp, mp_ptr rp, mp_srcptr np)
 
   prec = 8;
   np0 <<= 2 * prec;
-  while (2*prec < BITS_PER_MP_LIMB)
+  while (2 * prec < GMP_NUMB_BITS)
     {
       /* invariant: s has prec bits, and r <= 2*s */
-      r = (r << prec) + (np0 >> (BITS_PER_MP_LIMB - prec));
+      r = (r << prec) + (np0 >> (GMP_NUMB_BITS - prec));
       np0 <<= prec;
       u = 2 * s;
       q = r / u;
       u = r - q * u;
       s = (s << prec) + q;
-      u = (u << prec) + (np0 >> (BITS_PER_MP_LIMB - prec));
+      u = (u << prec) + (np0 >> (GMP_NUMB_BITS - prec));
       q = q * q;
       r = u - q;
       if (u < q)
@@ -104,11 +110,11 @@ mpn_sqrtrem1 (mp_ptr sp, mp_ptr rp, mp_srcptr np)
   sp[0] = s;
   if (rp != NULL)
     rp[0] = r;
-  return (r) ? 1 : 0;
+  return r != 0 ? 1 : 0;
 }
 
 
-#define Prec (BITS_PER_MP_LIMB >> 1)
+#define Prec (GMP_NUMB_BITS >> 1)
 
 /* same as mpn_sqrtrem, but for size=2 and {np, 2} normalized
    return cc such that {np, 2} = sp[0]^2 + cc*2^BITS+PER_MP_LIMB + rp[0] */
@@ -118,14 +124,14 @@ mpn_sqrtrem2 (mp_ptr sp, mp_ptr rp, mp_srcptr np)
   mp_limb_t qhl, q, u, np0;
   int cc;
 
-  ASSERT (np[1] >= MP_LIMB_T_HIGHBIT/2);
+  ASSERT (np[1] >= GMP_NUMB_HIGHBIT / 2);
 
   np0 = np[0];
-  mpn_sqrtrem1(sp, rp, np+1);
+  mpn_sqrtrem1 (sp, rp, np + 1);
   qhl = 0;
   while (rp[0] >= sp[0])
     {
-      qhl ++;
+      qhl++;
       rp[0] -= sp[0];
     }
   /* now rp[0] < sp[0] < 2^Prec */
@@ -133,19 +139,19 @@ mpn_sqrtrem2 (mp_ptr sp, mp_ptr rp, mp_srcptr np)
   u = 2 * sp[0];
   q = rp[0] / u;
   u = rp[0] - q * u;
-  q += (qhl & 1) << (Prec-1);
+  q += (qhl & 1) << (Prec - 1);
   qhl >>= 1; /* if qhl=1, necessary q=0 as qhl*2^Prec + q <= 2^Prec */
   /* now we have (initial rp[0])<<Prec + np0>>Prec = (qhl<<Prec + q) * (2sp[0]) + u */
-  sp[0] = ((sp[0]+qhl) << Prec) + q;
+  sp[0] = ((sp[0] + qhl) << Prec) + q;
   cc = u >> Prec;
-  rp[0] = (u << Prec) + (np0 & (((mp_limb_t)1 << Prec) - 1));
+  rp[0] = (u << Prec) + (np0 & (((mp_limb_t) 1 << Prec) - 1));
   /* subtract q * q or qhl*2^(2*Prec) from rp */
   cc -= mpn_sub_1(rp, rp, 1, q * q) + qhl;
   /* now subtract 2*q*2^Prec + 2^(2*Prec) if qhl is set */
   if (cc < 0)
     {
-      cc += (sp[0]) ? mpn_add_1(rp, rp, 1, sp[0]) : 1;
-      cc += mpn_add_1(rp, rp, 1, --sp[0]);
+      cc += sp[0] != 0 ? mpn_add_1 (rp, rp, 1, sp[0]) : 1;
+      cc += mpn_add_1 (rp, rp, 1, --sp[0]);
     }
   return cc;
 }
@@ -155,30 +161,31 @@ mpn_sqrtrem2 (mp_ptr sp, mp_ptr rp, mp_srcptr np)
    and in {np, n} the low n limbs of the remainder, returns the high
    limb of the remainder (which is 0 or 1).
    Assumes {np, 2n} is normalized, i.e. np[2n-1] >= B/4
-   where B=2^BITS_PER_MP_LIMB.
-*/
+   where B=2^GMP_NUMB_BITS.  */
 static mp_limb_t
-mpn_dq_sqrtrem (mp_ptr sp, mp_ptr np, mp_size_t n)
+mpn_dc_sqrtrem (mp_ptr sp, mp_ptr np, mp_size_t n)
 {
-  mp_limb_t q; /* carry out of {sp, n} */
-  int c, b; /* carry out of remainder */
+  mp_limb_t q;			/* carry out of {sp, n} */
+  int c, b;			/* carry out of remainder */
   mp_size_t l, h;
 
-  ASSERT (np[2*n-1] >= MP_LIMB_T_HIGHBIT/2);
+  ASSERT (np[2 * n - 1] >= GMP_NUMB_HIGHBIT / 2);
 
   if (n == 1)
     return mpn_sqrtrem2 (sp, np, np);
 
   l = n / 2;
   h = n - l;
-  q = mpn_dq_sqrtrem (sp + l, np + 2 * l, h);
-  if (q) mpn_sub_n (np + 2 * l, np + 2 * l, sp + l, h);
+  q = mpn_dc_sqrtrem (sp + l, np + 2 * l, h);
+  if (q != 0)
+    mpn_sub_n (np + 2 * l, np + 2 * l, sp + l, h);
   q += mpn_divrem (sp, 0, np + l, n, sp + l, h);
   c = sp[0] & 1;
   mpn_rshift (sp, sp, l, 1);
-  sp[l-1] |= q << (BITS_PER_MP_LIMB - 1);
+  sp[l - 1] |= q << (GMP_NUMB_BITS - 1);
   q >>= 1;
-  if (c) c = mpn_add_n (np + l, np + l, sp + l, h);
+  if (c != 0)
+    c = mpn_add_n (np + l, np + l, sp + l, h);
   mpn_sqr_n (np + n, sp, l);
   b = q + mpn_sub_n (np, np, np + n, 2 * l);
   c -= (l == h) ? b : mpn_sub_1 (np + 2 * l, np + 2 * l, 1, b);
@@ -209,48 +216,62 @@ mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nn)
   if (nn == 0)
     return 0;
 
-  ASSERT (np[nn-1] != 0);
+  ASSERT (np[nn - 1] != 0);
   ASSERT (rp == NULL || MPN_SAME_OR_SEPARATE_P (np, rp, nn));
-  ASSERT (rp == NULL || ! MPN_OVERLAP_P (sp, (nn+1)/2, rp, nn));
-  ASSERT (! MPN_OVERLAP_P (sp, (nn+1)/2, np, nn));
+  ASSERT (rp == NULL || ! MPN_OVERLAP_P (sp, (nn + 1) / 2, rp, nn));
+  ASSERT (! MPN_OVERLAP_P (sp, (nn + 1) / 2, np, nn));
 
-  high = np[nn-1];
-  if (nn == 1 && (high & MP_LIMB_T_HIGHBIT))
+  high = np[nn - 1];
+  if (nn == 1 && (high & GMP_NUMB_HIGHBIT))
     return mpn_sqrtrem1 (sp, rp, np);
-  count_leading_zeros(c, high);
+  count_leading_zeros (c, high);
+  c -= GMP_NAIL_BITS;
 
   c = c / 2; /* we have to shift left by 2c bits to normalize {np, nn} */
   tn = (nn+1) / 2; /* 2*tn is the smallest even integer >= nn */
 
   TMP_MARK (marker);
-  if ((nn % 2) || (c > 0))
+  if (nn % 2 != 0 || c > 0)
     {
       tp = TMP_ALLOC_LIMBS (2 * tn);
-      tp[0] = 0; /* needed only when 2*tn > nn, but saves a test */
-      if (c) mpn_lshift(tp + 2*tn - nn, np, nn, 2 * c);
-      else MPN_COPY (tp + 2*tn - nn, np, nn);
-      rl = mpn_dq_sqrtrem (sp, tp, tn);
-      /* we have 2^(2k)*N = S^2 + R where k = c + (2tn-nn)*BITS_PER_MP_LIMB/2,
+      tp[0] = 0;	     /* needed only when 2*tn > nn, but saves a test */
+      if (c != 0)
+	mpn_lshift (tp + 2 * tn - nn, np, nn, 2 * c);
+      else
+	MPN_COPY (tp + 2 * tn - nn, np, nn);
+      rl = mpn_dc_sqrtrem (sp, tp, tn);
+      /* We have 2^(2k)*N = S^2 + R where k = c + (2tn-nn)*GMP_NUMB_BITS/2,
 	 thus 2^(2k)*N = (S-s0)^2 + 2*S*s0 - s0^2 + R where s0=S mod 2^k */
-      c += (nn % 2) * BITS_PER_MP_LIMB / 2; /* c now represents k */
-      s0[0] = sp[0] & (((mp_limb_t) 1 << c) - 1); /* S mod 2^k */
-      rl += mpn_addmul_1 (tp, sp, tn, 2 * s0[0]); /* R = R + 2*s0*S */
+      c += (nn % 2) * GMP_NUMB_BITS / 2;		/* c now represents k */
+      s0[0] = sp[0] & (((mp_limb_t) 1 << c) - 1);	/* S mod 2^k */
+      rl += mpn_addmul_1 (tp, sp, tn, 2 * s0[0]);	/* R = R + 2*s0*S */
       cc = mpn_submul_1 (tp, s0, 1, s0[0]);
-      rl -= (tn > 1) ? mpn_sub_1(tp + 1, tp + 1, tn - 1, cc) : cc;
+      rl -= (tn > 1) ? mpn_sub_1 (tp + 1, tp + 1, tn - 1, cc) : cc;
       mpn_rshift (sp, sp, tn, c);
       tp[tn] = rl;
-      if (rp == NULL) rp = tp;
+      if (rp == NULL)
+	rp = tp;
       c = c << 1;
-      if (c < BITS_PER_MP_LIMB) tn++; else { tp++; c -= BITS_PER_MP_LIMB; }
-      if (c) mpn_rshift (rp, tp, tn, c); else MPN_COPY_INCR (rp, tp, tn);
+      if (c < GMP_NUMB_BITS)
+	tn++;
+      else
+	{
+	  tp++;
+	  c -= GMP_NUMB_BITS;
+	}
+      if (c != 0)
+	mpn_rshift (rp, tp, tn, c);
+      else
+	MPN_COPY_INCR (rp, tp, tn);
       rn = tn;
     }
   else
     {
       if (rp == NULL)
 	rp = TMP_ALLOC_LIMBS (nn);
-      if (rp != np) MPN_COPY (rp, np, nn);
-      rn = tn + (rp[tn] = mpn_dq_sqrtrem (sp, rp, tn));
+      if (rp != np)
+	MPN_COPY (rp, np, nn);
+      rn = tn + (rp[tn] = mpn_dc_sqrtrem (sp, rp, tn));
     }
 
   MPN_NORMALIZE (rp, rn);