Two 10-lines functions

[gianrico.fini at terra.es]

Dear developers,

I was challenged... the challenge was: write a 10-lines substitute for the 
function mpn_mulmod_2expp1.
It computes the product of xp and yp modulo 2^b+1. It assumes that xp and yp 
have exactly b bits.
The parameter c is 1 if yp is all zero but its value modulo 2^b+1 is -1; it is 
2 if xp is -1.
This function is needed for a benchmark checking if some Fermat's numbers are 
prime or not.

I send it here, so that you can use it if you need it.

int mpn_mulmod_2expp1(mp_ptr rp,mp_srcptr xp,mp_srcptr yp,int c,unsigned long 
b,mp_ptr tp)
{mp_limb_t c1, c2; mp_size_t l, d;
l = b/GMP_NUMB_BITS; b -= l*GMP_NUMB_BITS; l += (b!=0);
mpn_mul_n(tp, xp, yp, l); d = (GMP_NUMB_BITS-b)*(b!=0);
c1 = mpn_submul_1(tp, xp, l, c&1)+mpn_submul_1(tp, yp, l, c>>1)+(c==3);
c2 = ((c==0)&&(b!=0))?tp[l-1]>>b:0; tp[l-1]&= GMP_NUMB_MAX>>(d*(c==0));
c1+= mpn_sub_1(rp, tp, l, c2);
c1+= mpn_submul_1(rp, tp+l, l, CNST_LIMB(1)<<d);
c2 = mpn_add_1(rp, rp, l, c1); rp[l-1]&= GMP_NUMB_MAX>>d;
return c2;}

I also add the 2^b-1 version, which is simpler.

void mpn_mulmod_2expm1(mp_ptr rp,mp_ptr xp,mp_ptr yp,unsigned long b,mp_ptr 
tp)
{mp_limb_t c; mp_size_t l, d;
l = b/GMP_NUMB_BITS; b -= l*GMP_NUMB_BITS; l += (b!=0);
mpn_mul_n(tp, xp, yp, l); d = (GMP_NUMB_BITS-b)*(b!=0);
c = (b!=0)?tp[l-1]>>b:0; tp[l-1]&= GMP_NUMB_MAX>>d;
c = mpn_add_1(rp, tp, l, c);
c+= mpn_addmul_1(rp, tp+l, l, CNST_LIMB(1)<<d);
c+= (b!=0)?rp[l-1]>>b:0;
c = mpn_add_1(rp, rp, l, c); rp[l-1]&= GMP_NUMB_MAX>>d;
}

Gian.