[GMP kernel] using mpn_sqrtrem for sqrtr

[Bill Allombert <allomber@math.u-bordeaux.fr>]

Hello PARI-dev,

Here a patch that implement sqrtr using the GMP function
mpn_sqrtrem() (when the GMP kernel is used).

In fact it provide two different implementations.

sqrtr_with_gmp() add two extra words for rounding whereas
sqrtr_with_gmp2() only add one.
Somehow sqrtr_with_gmp() seems to be faster with a slightly
more accurate rounding.

Any comment/tuning welcome.

Cheers,
Bill

Index: src/kernel/gmp/mp.c
===================================================================
RCS file: /home/cvs/pari/src/kernel/gmp/mp.c,v
retrieving revision 1.53
diff -u -r1.53 mp.c
--- src/kernel/gmp/mp.c	18 Sep 2004 22:12:47 -0000	1.53
+++ src/kernel/gmp/mp.c	20 Sep 2004 19:30:45 -0000
@@ -1187,6 +1187,68 @@
   res= cgeti(l);
   res[1] = evalsigne(1) | evallgefint(l);
   mpn_sqrtrem(LIMBS(res),NULL,LIMBS(a),NLIMBS(a));
+  return res;
+}
+
+GEN sqrtr_with_gmp(GEN a)
+{
+  GEN res, b, c;
+  long pr=RNLIMBS(a);
+  long e=expo(a),er=e>>1;
+  if (!signe(a)) return realzero_bit(er);
+  res = cgetr(2 + pr);
+  if (e&1)
+  {
+    b = new_chunk((pr << 1)+2);
+    xmpn_mirrorcopy(b+pr+2, RLIMBS(a), pr);
+  }
+  else
+  {
+    c = ishiftr_spec(a,pr+2,-1);
+    b = new_chunk(pr);
+    /*xmpn_zero below will overwrite the code word of c, Yay!*/
+  }
+  xmpn_zero(b,pr+2);
+  c = new_chunk(pr+1);
+  mpn_sqrtrem(c,NULL,b,(pr<<1)+2);
+  if (((ulong)(c[0]))&HIGHBIT)
+    if (mpn_add_1(c+1,c+1,pr,1))
+      err(bugparier,"sqrtr_with_gmp");
+  xmpn_mirrorcopy(RLIMBS(res),c+1,pr);
+  res[1] = evalsigne(1) | evalexpo(er);
+  avma=(pari_sp)res;
+  return res;
+}
+
+GEN sqrtr_with_gmp2(GEN a)
+{
+  GEN res, b, c;
+  long pr=RNLIMBS(a);
+  long e=expo(a),er=e>>1;
+  if (!signe(a)) return realzero_bit(er);
+  res = cgetr(2 + pr);
+  if (e&1)
+  {
+    b = new_chunk((pr << 1)+1);
+    xmpn_mirrorcopy(b+pr+1, RLIMBS(a), pr);
+  }
+  else
+  {
+    c = ishiftr_spec(a,pr+2,-1);
+    b = pr<=1?c+pr-1:new_chunk(pr-1);
+    /*xmpn_zero below will overwrite the code word of c, Yay!*/
+  }
+  xmpn_zero(b,pr+1);
+  c = new_chunk(pr+1);
+  mpn_sqrtrem(c,NULL,b,(pr<<1)+1);
+  if (mpn_lshift(c,c,pr+1,BITS_IN_HALFULONG))
+    err(talker,"sqrtr_with_gmp22(%Z)",a);
+  if (((ulong)(c[0]))&HIGHBIT)
+    if (mpn_add_1(c+1,c+1,pr,1))
+      err(bugparier,"sqrtr_with_gmp");
+  xmpn_mirrorcopy(RLIMBS(res),c+1,pr);
+  res[1] = evalsigne(1) | evalexpo(er);
+  avma=(pari_sp)res;
   return res;
 }