Mersenne: Faster GMP exp mod 2**p-1

[Alexander Kruppa]

Hi,

I saw code posted earlier today that did LL testing with the GMP
library. It used the mpz_powm_ui function for exponentiating modulo x.
For x=2**p-1 this function can be sped up considerably, since
a mod (2^p-1) = a % 2^p + a / 2^p (integer operations) and these can
be done by bit shifting/masking.

I use this function below for my P-1 factorer (not released (yet?)).
It gave me a factor 3 speedup. It�s a pretty straigtforward
exponentiation algo 
exp(b,e) = if (e & 1) b*exp(b**2,e >>1) else exp(b**2,e >>1)
with the modulo by bit shifting.
It�s still worlds slower than the FFT code of George, of course. One day
I�ll have to deal with this FFT business.

To Will Edington: I think I saw in the P-1 factorer of the mers package
that it used mpz_powm_ui, too, perhaps you may want to replace it.

Ciao,
  Alex.


void mpz_powmmers_ui(mpz_t r,mpz_t b,unsigned long int e,unsigned long
int n) 
// thats (result, base, exponent, exponent of mersenne to modulo by
(2^n-1))
{
  mpz_t t,u;

  mpz_init(u);
  mpz_init_set(t,b);
  while (mpz_sizeinbase(t,2) > n)
  {
    mpz_tdiv_q_2exp(u,t,n);
    mpz_fdiv_r_2exp(t,t,n);
    mpz_add(t,t,u);
  }  


  if (e & 1)
    { mpz_set(r,t); }
  else
    { mpz_set_ui(r,1); }
  e = e >> 1;
  
  while (e)
  {
    mpz_mul(t,t,t);
    while (mpz_sizeinbase(t,2) > n)
    {
      mpz_tdiv_q_2exp(u,t,n);
      mpz_fdiv_r_2exp(t,t,n);
      mpz_add(t,t,u);
    }
    
    if (e & 1)
    {
      mpz_mul(r,r,t);
      while (mpz_sizeinbase(r,2) > n)
      {
        mpz_tdiv_q_2exp(u,r,n);
        mpz_fdiv_r_2exp(r,r,n);
        mpz_add(r,r,u);
      }
    }    
    e = e >> 1;
  }
}