modulo.H

Go to the documentation of this file.

// Modulo-Operations for unsigned int
// written by Thorsten Reinecke (1999-03-12)
// last change: 2007-06-21

// some of the following functions are necessary because a*b%m will
// return false results if a*b overflows unsigned-int-range.
// these routines are okay for unsigned-int-range/2

#ifndef MODULO_H_INCLUDED
#define MODULO_H_INCLUDED

#include <cstdlib>
#include <iostream>

namespace numtheory
{


inline unsigned int normalized_signed_mod(const signed int x, const int m)
{
#if defined(ASM_CMOV)
  // *** inline x86 assembler using conditional mov operations ***
#ifdef DEBUG
 #warning "using cmov-operation for normalized_signed_mod(const signed int,int)"
#endif
  // damned sign!!
  // x%m will not be normalized by default, but often you depend on it!
  // -1%2 -> -1 , 1%2 -> 1,
  // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
  register unsigned int result;
  register unsigned int clobbered;
  asm ( \
   "cdq # prepare signed division \n\t"
   "idivl %[mod] # remainder -> edx \n\t" \
   "mov %%edx,%%eax \n\t" \
   "addl %[mod],%%edx # this is a positive value \n\t" \
   "cmovnc %%eax,%%edx # and now it is normalized (minimal) \n\t"
   : "=a" (clobbered), "=&d" (result) : "a" (x), [mod] "rm" (m) : "cc");
  return result;
#elif defined(ASM_386)
  // *** inline i386 assembler ***
  // damned sign!!
  // x%m will not be normalized by default, but often you depend on it!
  // -1%2 -> -1 , 1%2 -> 1,
  // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
  register unsigned int result;
  register unsigned int clobbered;
  asm ( \
   "cdq # prepare signed division \n\t"
   "idivl %[mod] # remainder -> edx \n\t" \
   "#bt $31,%%edx # signed? \n\t" \
   "#sbb %%eax,%%eax # compute instead of branch! \n\t" \
   "#and %[mod],%%eax \n\t" \
   "#add %%eax,%%edx # and normalize to positive value \n\t"
   "testl %%edx,%%edx \n\t" \
   "jns 1f # signed? \n\t"
   "addl %[mod],%%edx # normalize to positive value \n\t"
   "1:"
   : "=a" (clobbered), "=&d" (result) : "a" (x), [mod] "rm" (m) : "cc");
  return result;
#else
  // *** generic code ***
  // damned sign!!
  // x%m will not be normalized by default, but often you depend on it!
  // -1%2 -> -1 , 1%2 -> 1,
  // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
  return (x%m<0) ? (x%m)+m : x%m;
#endif
}


inline unsigned int reciprocal(register const unsigned int m)
{
#if 1 && defined(ASM_386)
 #ifdef DEBUG
  #warning "activated 32 bit reciprocal"
 #endif
  register unsigned int result;
  register unsigned int clobbered;
  asm ( \
   "mov $1,%%edx \n\t" \
   "mov %[mod],%%eax \n\t" \
   "divl %[mod] \n\t" \
   : [h] "=&d" (clobbered), [res] "=&a" (result)
   : [mod] "r" (m)
   : "cc");
  return result;
#else
  return static_cast<unsigned int>(0x100000000LL/m)+1;
#endif
}


inline unsigned int normalized_signed_mod(register const signed int x, register const int m, const int &recip_m)
{
#if 1 && defined(ASM_X86_64)
 #ifdef DEBUG
  #warning "activated normalized signed mod using reciprocals, X86_64"
 #endif
  register unsigned int result;
  register unsigned int clobbered;
  asm ( \
   "imull %[a] \n\t" \
   "movl %[a],%[res] \n\t" \
   "imull %[mod],%[h] \n\t" \
   "subl %[h],%[res] \n\t" \
   "mov $0,%[h] \n\t" \
   "cmovs %[mod],%[h] \n\t" \
   "addl %[h],%[res] # normalize step1 \n\t"
   "movl %[res],%[h] \n\t" \
   "subl %[mod],%[res] \n\t"
   "cmovb %[h],%[res] # step2 \n\t" \
   : [h] "=&d" (clobbered), [res] "=a" (result)
   : [a] "r" (x), [mod] "r" (m), [recip] "a" (recip_m)
   : "cc");
#if 0 || defined(DEBUG)
  if (result!=normalized_signed_mod(x,m))
   {
     std::cerr << "error in normalized_signed_mod (using reciprocal)..." << std::endl;
     std::cerr << "x=" << x << "  m=" << m << std::endl;
     std::cerr << "res=" << static_cast<signed int>(result) << " instead of " << normalized_signed_mod(x,m) << std::endl;
     exit(1);
   }
#endif
  return result;
#elif 1 && defined(ASM_CMOV)
 #ifdef DEBUG
  #warning "activated normalized signed mod using reciprocals, cmov"
 #endif
  register unsigned int result;
  register unsigned int clobbered;
  asm ( \
   "imull %[a] \n\t" \
   "movl %[a],%[res] \n\t" \
   "imull %[mod],%[h] \n\t" \
   "subl %[h],%[res] \n\t" \
   "mov $0,%[h] \n\t" \
   "cmovs %[mod],%[h] \n\t" \
   "addl %[h],%[res] # normalize step1 \n\t"
   "movl %[res],%[h] \n\t" \
   "subl %[mod],%[res] \n\t"
   "cmovb %[h],%[res] # step2 \n\t" \
   : [h] "=&d" (clobbered), [res] "=a" (result)
   : [a] "r" (x), [mod] "r" (m), [recip] "a" (recip_m)
   : "cc");
#if 0 || defined(DEBUG)
  if (result!=normalized_signed_mod(x,m))
   {
     std::cerr << "error in normalized_signed_mod (using reciprocal)..." << std::endl;
     std::cerr << "x=" << x << "  m=" << m << std::endl;
     std::cerr << "res=" << static_cast<signed int>(result) << " instead of " << normalized_signed_mod(x,m) << std::endl;
     exit(1);
   }
#endif
  return result;
#elif 1 && defined(ASM_386)
 #ifdef DEBUG
  #warning "activated normalized signed mod using reciprocals, i386"
 #endif
  register unsigned int result;
  register unsigned int clobbered;
  asm ( \
   "imull %[a] \n\t" \
   "movl %[a],%[res] \n\t" \
   "imull %[mod],%[h] \n\t" \
   "subl %[h],%[res] \n\t" \
   "jns 0f \n\t" \
   "addl %[mod],%[res] # normalize step1 \n\t"
   "0: cmp %[mod],%[res] \n\t"
   "jb 0f \n\t" \
   "subl %[mod],%[res] # step2 \n\t" \
   "0: \n\t" \
   : [h] "=&d" (clobbered), [res] "=a" (result)
   : [a] "r" (x), [mod] "r" (m), [recip] "a" (recip_m)
   : "cc");
#if 0 || defined(DEBUG)
  if (result!=normalized_signed_mod(x,m))
   {
     std::cerr << "error in normalized_signed_mod (using reciprocal)..." << std::endl;
     std::cerr << "x=" << x << "  m=" << m << std::endl;
     std::cerr << "res=" << static_cast<signed int>(result) << " instead of " << normalized_signed_mod(x,m) << std::endl;
     exit(1);
   }
#endif
  return result;
#else
 #warning "normalized signed mod fallback code (ignore reciprocals)"
  return normalized_signed_mod(x,m);
#endif
}



#if defined(ASM_386) || defined(ASM_CMOV) || defined(ASM_X86_64)
// *** i386 specific assembler code ***
// runs much faster than the generic code and has less constraints :)

#ifdef DEBUG
 #warning "i386-assembler code optimizations for mulmod and squaremod enabled"
#endif

inline unsigned int mulmod(const unsigned int a, const unsigned int b, const unsigned int m)
{
  // returns (a*b) mod m
  unsigned int x;
  register unsigned int clobbered;
  asm ( \
   "mull %[factor2] # multiply the factors \n\t" \
   "divl %[mod] # remainder -> edx \n\t"
   : "=a" (clobbered), "=&d" (x) : "a" (a), [factor2] "rm" (b), [mod] "rm" (m) : "cc");

  // IMPORTANT REMARK
  // If you request a register for an input operand and this register gets modified,
  // you have to tell GCC that this register is invalid for this operand after leaving the
  // assembler code...
  // You cannot put the register on the clobberlist (because GCC forbids operands to use
  // clobbered registers), 
  // instead you have to declare a dummy output operand that uses the same register...
  // Be careful: if all output operands are dummies, you have to declare the assembler code
  // as volatile!

  return x;
}

inline unsigned int squaremod(const unsigned int a, const unsigned int m)
{
  // returns (a^2) mod m
  unsigned int x;
  register unsigned int clobbered;
  asm ( \
   "mull %%eax # square the factor \n\t" \
   "divl %[mod] # remainder -> edx \n\t"
   : "=a" (clobbered), "=&d" (x) : "a" (a), [mod] "g" (m) : "cc");
  return x;
}

#else
// *** generic variants for mulmod and squaremod ***

inline unsigned int mulmod(register unsigned int a, register unsigned int b, const unsigned int m)
{
  // assumes a<m, b<m !!
  // hint: multiplication runs faster if b<=a
  // returns (a*b) mod m
  //a%=m; b%=m; // (only necessary if a>=m or b>=m)
  register unsigned int x = (b&1) ? a : 0;
  while (b>>=1)
   {
    a<<=1; if (a>=m) a-=m;
    if (b&1) { x+=a; if (x>=m) x-=m; }
   }
  return x;
}

inline unsigned int squaremod(register unsigned int a, const unsigned int m)
{
  // assumes a<m !!
  // returns (a^2) mod m
  //a%=m; // only necessary if a>=m
  register unsigned int b=a;
  register unsigned int x = (b&1) ? a : 0;
  while (b>>=1)
   {
    a<<=1; if (a>=m) a-=m;
    if (b&1) { x+=a; if (x>=m) x-=m; }
   }
  return x;
}

#endif


inline unsigned int normalized_signed_mulmod(signed int x1, signed int x2, const int m)
{
#if defined(ASM_CMOV) || defined(ASM_X86_64)
  // *** inline x86 assembler using conditional mov operations ***
#ifdef DEBUG
 #warning "using cmov-operation for normalized_signed_mulmod(const signed int,const signed int,int)"
#endif
  // damned sign!!
  // x%m will not be normalized by default, but often you depend on it!
  // -1%2 -> -1 , 1%2 -> 1,
  // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
  register unsigned int result;
  register unsigned int clobbered;
  asm ( \
   "imull %[x2] # signed multiplication \n\t"
   "idivl %[mod] # remainder -> edx \n\t" \
   "mov %%edx,%%eax \n\t" \
   "addl %[mod],%%edx # this is a positive value \n\t" \
   "cmovnc %%eax,%%edx # and now it is normalized (minimal) \n\t"
   : "=a" (clobbered), "=&d" (result) : "a" (x1), [x2] "g" (x2), [mod] "g" (m) : "cc");
  return result;
#elif defined(ASM_386)
  // *** inline i386 assembler ***
  // damned sign!!
  // x%m will not be normalized by default, but often you depend on it!
  // -1%2 -> -1 , 1%2 -> 1,
  // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
  register unsigned int result;
  register unsigned int clobbered;
  asm ( \
   "imull %[x2] # prepare signed division \n\t"
   "idivl %[mod] # remainder -> edx \n\t" \
   "#bt $31,%%edx # signed? \n\t" \
   "#sbb %%eax,%%eax # compute instead of branch! \n\t" \
   "#and %[mod],%%eax \n\t" \
   "#add %%eax,%%edx # and normalize to positive value \n\t" 
   "testl %%edx,%%edx \n\t" \
   "jns 1f # signed? \n\t"
   "addl %[mod],%%edx # normalize to positive value \n\t"
   "1:"
   : "=a" (clobbered), "=&d" (result) : "a" (x1), [x2] "g" (x2), [mod] "g" (m) : "cc");
  return result;
#else
  // *** generic code ***
  // damned sign!!
  // x%m will not be normalized by default, but often you depend on it!
  // -1%2 -> -1 , 1%2 -> 1,
  // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
  bool sign = (x1<0)^(x2<0);
  if (x1<0) x1=-x1;
  if (x2<0) x2=-x2;
  register unsigned int r = mulmod(x1,x2,m);
  return sign ? m-r : r;
#endif
}


#if defined(ASM_CMOV) || defined(ASM_X86_64)

#ifdef DEBUG
 #warning "using cmov-operation for powmod(unsigned int, unsigned int, const unsigned int)"
#endif

inline unsigned int powmod(register unsigned int a, register unsigned int pow, const unsigned int m)
{
  // returns (a^pow) mod m
  // assumes a<m
  register unsigned int x;
  unsigned int temp;
  asm ( \
   "movl $1,%[X]       \n\t" \
   "shrl $1,%[pow]     \n\t" \
   "cmovc %%eax,%[X]   \n\t" \
   "jz 1f              \n\t" \
   "0: mull %%eax      \n\t" \
   "divl %[M]          \n\t" \
   "shrl $1,%[pow]     \n\t" \
   "movl %%edx,%%eax   \n\t" \
   "ja 0b  # CF=0 and ZF=0 \n\t" \
   "movl %%edx,%[tmp]  \n\t" \
   "mull %[X]          \n\t" \
   "divl %[M]          \n\t" \
   "cmpl $0,%[pow]     \n\t" \
   "movl %[tmp],%%eax  \n\t" \
   "movl %%edx,%[X]    \n\t" \
   "jne 0b             \n\t" \
   "1:  " \
//   : [X] "=&r" (x), "+a" (a), [pow] "+rm" (pow), [tmp] "=&rm" (temp) : [M] "rm" (m) : "cc", "edx");
// FIXME: "can't find a register in class 'AD_REGS' while reloading 'asm'"
// if the above line triggers compilation errors, because gcc refuses to choose an alternative which works,
// then you need to use suboptimal constraints (use always memory operand instead of possible register):
   : [X] "=&r" (x), "+a" (a), [pow] "+rm" (pow), [tmp] "=m" (temp) : [M] "rm" (m) : "cc", "edx");

  return x;
}

#else
// *** generic variant of powmod ***

inline unsigned int powmod(register unsigned int a, register unsigned int pow, const unsigned int m)
{
  // assumes a<m
  // returns (a^pow) mod m
  unsigned int x = (pow&1) ? a : 1;
  while (pow>>=1)
   {
     a=squaremod(a,m);
     if (pow&1) x=mulmod(x,a,m);
   }
  return x;
}
#endif


bool strong_probab_prime(const unsigned int p, const unsigned int base);


bool probab_prime(const unsigned int p);


bool is_prime(register const int p);


inline unsigned int gcd(register unsigned int a, register unsigned int b)
{
  // returns greatest common divisor
  while (a&&b) if (a>b) a%=b; else b%=a;
  return a ? a : b;
}

inline unsigned short int gcd(register unsigned short int a, register unsigned short int b)
{
  // returns greatest common divisor
  while (a&&b) if (a>b) a%=b; else b%=a;
  return a ? a : b;
}

inline unsigned int oddgcd(register unsigned int a, register unsigned int b)
{
#if defined(ASM_CMOV) || defined(ASM_X86_64)
 // avoid branches since cmovs are faster on random data and enable parallel execution

#ifdef DEBUG
 #warning "using cmov-operation for oddgcd(unsigned int, unsigned int)"
#endif

  asm ( \
   "cmpl %[a],%[b] \n\t" \
   "je 2f \n\t" \
   "1: movl %[a],%%ecx \n\t" \
   "shrl $1,%[a] \n\t" \
   "cmovcl %%ecx,%[a] \n\t" \
   "movl %[b],%%edx \n\t" \
   "shrl $1,%[b] \n\t" \
   "cmovcl %%edx,%[b] \n\t" \
   "movl %[a],%%ecx \n\t" \
   "subl %[b],%[a] \n\t" \
   "cmovbe %%ecx,%[a] \n\t" \
   "movl %[b],%%edx \n\t" \
   "subl %[a],%[b] \n\t" \
   "cmovbe %%edx,%[b] \n\t" \
   "jnz 1b \n\t" \
   "2: \n" \
   : [a] "+r" (a), [b] "+r" (b) : : "ecx", "edx", "cc");
  return a;
#else
  // *** generic version ***
  // returns greatest (odd) common divisor
  if (a==0) return b;
  if (b==0) return a;
  do
   {
     //cout << a << "," << b << endl;
     while (1&a^1) a>>=1;
     while (1&b^1) b>>=1;
     if (a>b) a-=b;
     if (b>a) b-=a;
   } while (a!=b);
  //cout << a << endl;
  return a;
#endif
}


inline bool coprime(register unsigned int a, register unsigned int b)
{
  return ((a|b)&1) ? oddgcd(a,b)==1 : false;
}


signed int jacobi (int a, unsigned int b);


signed int legendre (signed int a, const unsigned int p);


unsigned int invmod(const unsigned int x, const unsigned int m) __attribute__ ((const));


unsigned int bininvmod(register unsigned int x, register const unsigned int m) __attribute__ ((const));


unsigned int montgominvmod(register unsigned int x, unsigned int m) __attribute__ ((const));


inline unsigned int fastinvmod(const unsigned int x, const unsigned int m)
{
#if defined(DEBUG)
  if (x==0 || x>=m) std::cerr << "fastinvmod: constraint 0<y<m not met: x=" << x << ", m=" << m << std::endl;
  if (!coprime(x,m)) std::cerr << "fastinvmod: constraint coprime(x,m) not met: x=" << x << ", m=" << m << std::endl;
  if ((m&1)==0) std::cerr << "fastinvmod: constraint m=1 (mod 2)  [=odd m] not met: x=" << x << ", m=" << m << std::endl;
#endif 

#if defined (ASM_386) || defined(ASM_X86_64)
  //return invmod(x,m);
  //return bininvmod(x,m);
  return montgominvmod(x,m);
#else
  //return invmod(x,m);
  return bininvmod(x,m);
  //return montgominvmod(x,m);
#endif
}


#if 0
 // use this, if
 // a) there is no crafted invmod support for your cpu
 // and b) your cpu supports very fast floating point arithmetic
 #define USE_JASONS_INVMOD
 namespace japinvmod
 {
   unsigned int fast_invmod(unsigned int a, unsigned int p);
 }
#endif


inline unsigned int fastinvmod_23bit(const unsigned int x, const unsigned int m)
{
#if defined(DEBUG)
  if (x==0 || x>=m) std::cerr << "fastinvmod_23bit: constraint 0<y<m not met: x=" << x << ", m=" << m << std::endl;
  if (!coprime(x,m)) std::cerr << "fastinvmod_23bit: constraint coprime(x,m) not met: x=" << x << ", m=" << m << std::endl;
  if ((m&1)==0) std::cerr << "fastinvmod_23bit: constraint m=1 (mod 2)  [=odd m] not met: x=" << x << ", m=" << m << std::endl;
  if (m>8388608) std::cerr << "fastinvmod_23bit: m exceeds 23 bit: x=" << x << ", m=" << m << std::endl;
#endif 

#if defined(USE_JASONS_INVMOD)
  return japinvmod::fast_invmod(x,m);
#else
  //return invmod(x,m);
  //return bininvmod(x,m);
  return montgominvmod(x,m);
#endif
}




unsigned int sqrtmod(const unsigned int Radikant, const unsigned int Primzahl);


} // namespace numtheory

#endif /* MODULO_H_INCLUDED */

---