Bill Allombert on Mon, 20 Sep 2004 22:58:48 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
[GMP kernel] using mpn_sqrtrem for sqrtr |
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; }