sqrt_modulo.cc

Go to the documentation of this file.

// The following functions compute square roots modulo prime numbers.
// please note:
//  - primality is not checked! 
//    (if necessary, you have to do this beforehand, eg. using is_probab_prime)
//  - no check is done, whether the square root exists!
//    (if necessary, check it using legendre(Radikant,Primenumber)==1) 

// written by Thorsten Reinecke 1998,1999
// last change: 2005-02-28

// now it should be threadsafe


#include "modulo.H"

namespace numtheory
{

class Clucas_capsule
{
 private:
  static const int lucas_cache_size = 128;
  unsigned int lucas_p, lucas_q, lucas_p_inv;
  int lucas_cache_index;
  unsigned int lucas_cache [lucas_cache_size][2];
  unsigned int lucasv(const unsigned int Primenumber, const unsigned int m);
 public:
  inline Clucas_capsule(void) : lucas_cache_index(0) { }
  unsigned int lucas(const unsigned int Radikant, const unsigned int Primenumber);  
};


unsigned int Clucas_capsule::lucasv(const unsigned int Primenumber, const unsigned int m)
{
  using std::cerr;
  using std::endl;
  if (m==0) return 2;
  if (m==1) return lucas_p;
  
  // value already known?
  for (int i=lucas_cache_index-1; i>=0; --i)
    {
      if (lucas_cache[i][0]==m) return lucas_cache[i][1];
    }
  
  unsigned int wert; // return value for lucas_v
  // printf("lucasv(%i)\n",m);
  
  if (m&1) // odd?
    {
      const unsigned int h1=lucasv(Primenumber,m+1);
      const unsigned int h2=lucasv(Primenumber,m-1);
      register unsigned int xx;
      xx=mulmod(h2,lucas_q,Primenumber)+h1;
      wert=mulmod(xx%Primenumber,lucas_p_inv,Primenumber);
    }
  else // even
    {
      const unsigned int h1=lucasv(Primenumber,m>>1);
      unsigned int xx;
      xx=(powmod(lucas_q,m>>1,Primenumber)<<1)%Primenumber;
      wert=(squaremod(h1,Primenumber)+Primenumber-xx)%Primenumber;
  
      // write value to cache...
      if (lucas_cache_index<lucas_cache_size)
        {
          lucas_cache[lucas_cache_index][0]=m;
          lucas_cache[lucas_cache_index++][1]=wert;
        } else cerr << "Lucas-Cache should be increased!" << endl;
    }
  return wert;
}


unsigned int Clucas_capsule::lucas(const unsigned int Radikant, const unsigned int Primenumber)
{
  using std::cerr;
  using std::endl;
  if ((Primenumber&3)!=1)
    {
      cerr << "Error in Lucassequence Primenumber%4<>1!: " << Primenumber%4 << endl;
      exit(1);
    }
  
  // compute lucas_p
  lucas_q=Radikant;
  lucas_p=1;
 
  {
    int h1,h2;
    do
     {
       lucas_p++;
       h1=squaremod(lucas_p,Primenumber);
       h2=mulmod(4,lucas_q,Primenumber);
       if (h1>=h2) h1-=h2; else h1=Primenumber-(h2-h1);
     } while (legendre(h1,Primenumber)!=-1);
  }

  // compute inverse of lucas_p mod Primenumber
  lucas_p_inv=fastinvmod(lucas_p,Primenumber);
  
  //cout << "Lucas-Cache was " << lucas_cache_index << endl;
  lucas_cache_index=0; // clear cache!
  
  // now compute the squareroot of n
  unsigned int v=lucasv(Primenumber,(Primenumber+1)>>1);
  // divide v by 2 (but do it modulo Primenumber!!)
  if (v&1) v+=Primenumber;
  return v>>1;
}


unsigned int sqrtmod(const unsigned int Radikant, const unsigned int Primenumber)
{
  // we could perform a primality check here, but we omit it for the sake of efficiency.  

  if ((Primenumber&3)==3)
    {
      /* we can can compute our result easily... */
      //cout << "sqrtmod: p=3 (mod 4)" << endl;
      return powmod(Radikant,(Primenumber+1)>>2,Primenumber);
    }
  else
    if ((Primenumber&7)==5)
      {
        /* result can be computed rather easily, too... */
        //cout << "sqrtmod: p=5 (mod 8)" << endl;
        unsigned int y = powmod(Radikant,(Primenumber+3)>>3,Primenumber);
        if (squaremod(y,Primenumber)==Radikant) return y;
        else 
          {
            unsigned int w=powmod(2,Primenumber>>2,Primenumber);
            return mulmod(y,w,Primenumber);
          }
      }
    else
     { 
       if (Radikant<=1) return Radikant; // special cases: 0^2=0, 1^2=1
       /* efficient solution by using Lucas-Sequences for Primenumber=1 (mod 4) */
       //cout << "sqrtmod: p=1 (mod 4) -> Lucassequences" << endl;
       Clucas_capsule lucas_capsule;
       return lucas_capsule.lucas(Radikant,Primenumber);
     }
}

} // namespace numtheory

---