modulo.H

Go to the documentation of this file.
00001 // Modulo-Operations for unsigned int
00002 // written by Thorsten Reinecke (1999-03-12)
00003 // last change: 2007-06-21
00004 
00005 // some of the following functions are necessary because a*b%m will
00006 // return false results if a*b overflows unsigned-int-range.
00007 // these routines are okay for unsigned-int-range/2
00008 
00023 #ifndef MODULO_H_INCLUDED
00024 #define MODULO_H_INCLUDED
00025 
00026 #include <cstdlib>
00027 #include <iostream>
00028 
00030 namespace numtheory
00031 {
00032 
00033 
00042 inline unsigned int normalized_signed_mod(const signed int x, const int m)
00043 {
00044 #if defined(ASM_CMOV)
00045   // *** inline x86 assembler using conditional mov operations ***
00046 #ifdef DEBUG
00047  #warning "using cmov-operation for normalized_signed_mod(const signed int,int)"
00048 #endif
00049   // damned sign!!
00050   // x%m will not be normalized by default, but often you depend on it!
00051   // -1%2 -> -1 , 1%2 -> 1,
00052   // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
00053   register unsigned int result;
00054   register unsigned int clobbered;
00055   asm ( \
00056    "cdq # prepare signed division \n\t"
00057    "idivl %[mod] # remainder -> edx \n\t" \
00058    "mov %%edx,%%eax \n\t" \
00059    "addl %[mod],%%edx # this is a positive value \n\t" \
00060    "cmovnc %%eax,%%edx # and now it is normalized (minimal) \n\t"
00061    : "=a" (clobbered), "=&d" (result) : "a" (x), [mod] "rm" (m) : "cc");
00062   return result;
00063 #elif defined(ASM_386)
00064   // *** inline i386 assembler ***
00065   // damned sign!!
00066   // x%m will not be normalized by default, but often you depend on it!
00067   // -1%2 -> -1 , 1%2 -> 1,
00068   // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
00069   register unsigned int result;
00070   register unsigned int clobbered;
00071   asm ( \
00072    "cdq # prepare signed division \n\t"
00073    "idivl %[mod] # remainder -> edx \n\t" \
00074    "#bt $31,%%edx # signed? \n\t" \
00075    "#sbb %%eax,%%eax # compute instead of branch! \n\t" \
00076    "#and %[mod],%%eax \n\t" \
00077    "#add %%eax,%%edx # and normalize to positive value \n\t"
00078    "testl %%edx,%%edx \n\t" \
00079    "jns 1f # signed? \n\t"
00080    "addl %[mod],%%edx # normalize to positive value \n\t"
00081    "1:"
00082    : "=a" (clobbered), "=&d" (result) : "a" (x), [mod] "rm" (m) : "cc");
00083   return result;
00084 #else
00085   // *** generic code ***
00086   // damned sign!!
00087   // x%m will not be normalized by default, but often you depend on it!
00088   // -1%2 -> -1 , 1%2 -> 1,
00089   // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
00090   return (x%m<0) ? (x%m)+m : x%m;
00091 #endif
00092 }
00093 
00094 
00099 inline unsigned int reciprocal(register const unsigned int m)
00100 {
00101 #if 1 && defined(ASM_386)
00102  #ifdef DEBUG
00103   #warning "activated 32 bit reciprocal"
00104  #endif
00105   register unsigned int result;
00106   register unsigned int clobbered;
00107   asm ( \
00108    "mov $1,%%edx \n\t" \
00109    "mov %[mod],%%eax \n\t" \
00110    "divl %[mod] \n\t" \
00111    : [h] "=&d" (clobbered), [res] "=&a" (result)
00112    : [mod] "r" (m)
00113    : "cc");
00114   return result;
00115 #else
00116   return static_cast<unsigned int>(0x100000000LL/m)+1;
00117 #endif
00118 }
00119 
00120 
00133 inline unsigned int normalized_signed_mod(register const signed int x, register const int m, const int &recip_m)
00134 {
00135 #if 1 && defined(ASM_X86_64)
00136  #ifdef DEBUG
00137   #warning "activated normalized signed mod using reciprocals, X86_64"
00138  #endif
00139   register unsigned int result;
00140   register unsigned int clobbered;
00141   asm ( \
00142    "imull %[a] \n\t" \
00143    "movl %[a],%[res] \n\t" \
00144    "imull %[mod],%[h] \n\t" \
00145    "subl %[h],%[res] \n\t" \
00146    "mov $0,%[h] \n\t" \
00147    "cmovs %[mod],%[h] \n\t" \
00148    "addl %[h],%[res] # normalize step1 \n\t"
00149    "movl %[res],%[h] \n\t" \
00150    "subl %[mod],%[res] \n\t"
00151    "cmovb %[h],%[res] # step2 \n\t" \
00152    : [h] "=&d" (clobbered), [res] "=a" (result)
00153    : [a] "r" (x), [mod] "r" (m), [recip] "a" (recip_m)
00154    : "cc");
00155 #if 0 || defined(DEBUG)
00156   if (result!=normalized_signed_mod(x,m))
00157    {
00158      std::cerr << "error in normalized_signed_mod (using reciprocal)..." << std::endl;
00159      std::cerr << "x=" << x << "  m=" << m << std::endl;
00160      std::cerr << "res=" << static_cast<signed int>(result) << " instead of " << normalized_signed_mod(x,m) << std::endl;
00161      exit(1);
00162    }
00163 #endif
00164   return result;
00165 #elif 1 && defined(ASM_CMOV)
00166  #ifdef DEBUG
00167   #warning "activated normalized signed mod using reciprocals, cmov"
00168  #endif
00169   register unsigned int result;
00170   register unsigned int clobbered;
00171   asm ( \
00172    "imull %[a] \n\t" \
00173    "movl %[a],%[res] \n\t" \
00174    "imull %[mod],%[h] \n\t" \
00175    "subl %[h],%[res] \n\t" \
00176    "mov $0,%[h] \n\t" \
00177    "cmovs %[mod],%[h] \n\t" \
00178    "addl %[h],%[res] # normalize step1 \n\t"
00179    "movl %[res],%[h] \n\t" \
00180    "subl %[mod],%[res] \n\t"
00181    "cmovb %[h],%[res] # step2 \n\t" \
00182    : [h] "=&d" (clobbered), [res] "=a" (result)
00183    : [a] "r" (x), [mod] "r" (m), [recip] "a" (recip_m)
00184    : "cc");
00185 #if 0 || defined(DEBUG)
00186   if (result!=normalized_signed_mod(x,m))
00187    {
00188      std::cerr << "error in normalized_signed_mod (using reciprocal)..." << std::endl;
00189      std::cerr << "x=" << x << "  m=" << m << std::endl;
00190      std::cerr << "res=" << static_cast<signed int>(result) << " instead of " << normalized_signed_mod(x,m) << std::endl;
00191      exit(1);
00192    }
00193 #endif
00194   return result;
00195 #elif 1 && defined(ASM_386)
00196  #ifdef DEBUG
00197   #warning "activated normalized signed mod using reciprocals, i386"
00198  #endif
00199   register unsigned int result;
00200   register unsigned int clobbered;
00201   asm ( \
00202    "imull %[a] \n\t" \
00203    "movl %[a],%[res] \n\t" \
00204    "imull %[mod],%[h] \n\t" \
00205    "subl %[h],%[res] \n\t" \
00206    "jns 0f \n\t" \
00207    "addl %[mod],%[res] # normalize step1 \n\t"
00208    "0: cmp %[mod],%[res] \n\t"
00209    "jb 0f \n\t" \
00210    "subl %[mod],%[res] # step2 \n\t" \
00211    "0: \n\t" \
00212    : [h] "=&d" (clobbered), [res] "=a" (result)
00213    : [a] "r" (x), [mod] "r" (m), [recip] "a" (recip_m)
00214    : "cc");
00215 #if 0 || defined(DEBUG)
00216   if (result!=normalized_signed_mod(x,m))
00217    {
00218      std::cerr << "error in normalized_signed_mod (using reciprocal)..." << std::endl;
00219      std::cerr << "x=" << x << "  m=" << m << std::endl;
00220      std::cerr << "res=" << static_cast<signed int>(result) << " instead of " << normalized_signed_mod(x,m) << std::endl;
00221      exit(1);
00222    }
00223 #endif
00224   return result;
00225 #else
00226  #warning "normalized signed mod fallback code (ignore reciprocals)"
00227   return normalized_signed_mod(x,m);
00228 #endif
00229 }
00230 
00231 
00232 
00233 #if defined(ASM_386) || defined(ASM_CMOV) || defined(ASM_X86_64)
00234 // *** i386 specific assembler code ***
00235 // runs much faster than the generic code and has less constraints :)
00236 
00237 #ifdef DEBUG
00238  #warning "i386-assembler code optimizations for mulmod and squaremod enabled"
00239 #endif
00240 
00252 inline unsigned int mulmod(const unsigned int a, const unsigned int b, const unsigned int m)
00253 {
00254   // returns (a*b) mod m
00255   unsigned int x;
00256   register unsigned int clobbered;
00257   asm ( \
00258    "mull %[factor2] # multiply the factors \n\t" \
00259    "divl %[mod] # remainder -> edx \n\t"
00260    : "=a" (clobbered), "=&d" (x) : "a" (a), [factor2] "rm" (b), [mod] "rm" (m) : "cc");
00261 
00262   // IMPORTANT REMARK
00263   // If you request a register for an input operand and this register gets modified,
00264   // you have to tell GCC that this register is invalid for this operand after leaving the
00265   // assembler code...
00266   // You cannot put the register on the clobberlist (because GCC forbids operands to use
00267   // clobbered registers), 
00268   // instead you have to declare a dummy output operand that uses the same register...
00269   // Be careful: if all output operands are dummies, you have to declare the assembler code
00270   // as volatile!
00271 
00272   return x;
00273 }
00274 
00286 inline unsigned int squaremod(const unsigned int a, const unsigned int m)
00287 {
00288   // returns (a^2) mod m
00289   unsigned int x;
00290   register unsigned int clobbered;
00291   asm ( \
00292    "mull %%eax # square the factor \n\t" \
00293    "divl %[mod] # remainder -> edx \n\t"
00294    : "=a" (clobbered), "=&d" (x) : "a" (a), [mod] "g" (m) : "cc");
00295   return x;
00296 }
00297 
00298 #else
00299 // *** generic variants for mulmod and squaremod ***
00300 
00309 inline unsigned int mulmod(register unsigned int a, register unsigned int b, const unsigned int m)
00310 {
00311   // assumes a<m, b<m !!
00312   // hint: multiplication runs faster if b<=a
00313   // returns (a*b) mod m
00314   //a%=m; b%=m; // (only necessary if a>=m or b>=m)
00315   register unsigned int x = (b&1) ? a : 0;
00316   while (b>>=1)
00317    {
00318     a<<=1; if (a>=m) a-=m;
00319     if (b&1) { x+=a; if (x>=m) x-=m; }
00320    }
00321   return x;
00322 }
00323 
00332 inline unsigned int squaremod(register unsigned int a, const unsigned int m)
00333 {
00334   // assumes a<m !!
00335   // returns (a^2) mod m
00336   //a%=m; // only necessary if a>=m
00337   register unsigned int b=a;
00338   register unsigned int x = (b&1) ? a : 0;
00339   while (b>>=1)
00340    {
00341     a<<=1; if (a>=m) a-=m;
00342     if (b&1) { x+=a; if (x>=m) x-=m; }
00343    }
00344   return x;
00345 }
00346 
00347 #endif
00348 
00349 
00358 inline unsigned int normalized_signed_mulmod(signed int x1, signed int x2, const int m)
00359 {
00360 #if defined(ASM_CMOV) || defined(ASM_X86_64)
00361   // *** inline x86 assembler using conditional mov operations ***
00362 #ifdef DEBUG
00363  #warning "using cmov-operation for normalized_signed_mulmod(const signed int,const signed int,int)"
00364 #endif
00365   // damned sign!!
00366   // x%m will not be normalized by default, but often you depend on it!
00367   // -1%2 -> -1 , 1%2 -> 1,
00368   // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
00369   register unsigned int result;
00370   register unsigned int clobbered;
00371   asm ( \
00372    "imull %[x2] # signed multiplication \n\t"
00373    "idivl %[mod] # remainder -> edx \n\t" \
00374    "mov %%edx,%%eax \n\t" \
00375    "addl %[mod],%%edx # this is a positive value \n\t" \
00376    "cmovnc %%eax,%%edx # and now it is normalized (minimal) \n\t"
00377    : "=a" (clobbered), "=&d" (result) : "a" (x1), [x2] "g" (x2), [mod] "g" (m) : "cc");
00378   return result;
00379 #elif defined(ASM_386)
00380   // *** inline i386 assembler ***
00381   // damned sign!!
00382   // x%m will not be normalized by default, but often you depend on it!
00383   // -1%2 -> -1 , 1%2 -> 1,
00384   // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
00385   register unsigned int result;
00386   register unsigned int clobbered;
00387   asm ( \
00388    "imull %[x2] # prepare signed division \n\t"
00389    "idivl %[mod] # remainder -> edx \n\t" \
00390    "#bt $31,%%edx # signed? \n\t" \
00391    "#sbb %%eax,%%eax # compute instead of branch! \n\t" \
00392    "#and %[mod],%%eax \n\t" \
00393    "#add %%eax,%%edx # and normalize to positive value \n\t" 
00394    "testl %%edx,%%edx \n\t" \
00395    "jns 1f # signed? \n\t"
00396    "addl %[mod],%%edx # normalize to positive value \n\t"
00397    "1:"
00398    : "=a" (clobbered), "=&d" (result) : "a" (x1), [x2] "g" (x2), [mod] "g" (m) : "cc");
00399   return result;
00400 #else
00401   // *** generic code ***
00402   // damned sign!!
00403   // x%m will not be normalized by default, but often you depend on it!
00404   // -1%2 -> -1 , 1%2 -> 1,
00405   // therefore (-1%2 != 1%2) but one surely expects -1=1 (mod 2) !!!
00406   bool sign = (x1<0)^(x2<0);
00407   if (x1<0) x1=-x1;
00408   if (x2<0) x2=-x2;
00409   register unsigned int r = mulmod(x1,x2,m);
00410   return sign ? m-r : r;
00411 #endif
00412 }
00413 
00414 
00415 #if defined(ASM_CMOV) || defined(ASM_X86_64)
00416 
00417 #ifdef DEBUG
00418  #warning "using cmov-operation for powmod(unsigned int, unsigned int, const unsigned int)"
00419 #endif
00420 
00428 inline unsigned int powmod(register unsigned int a, register unsigned int pow, const unsigned int m)
00429 {
00430   // returns (a^pow) mod m
00431   // assumes a<m
00432   register unsigned int x;
00433   unsigned int temp;
00434   asm ( \
00435    "movl $1,%[X]       \n\t" \
00436    "shrl $1,%[pow]     \n\t" \
00437    "cmovc %%eax,%[X]   \n\t" \
00438    "jz 1f              \n\t" \
00439    "0: mull %%eax      \n\t" \
00440    "divl %[M]          \n\t" \
00441    "shrl $1,%[pow]     \n\t" \
00442    "movl %%edx,%%eax   \n\t" \
00443    "ja 0b  # CF=0 and ZF=0 \n\t" \
00444    "movl %%edx,%[tmp]  \n\t" \
00445    "mull %[X]          \n\t" \
00446    "divl %[M]          \n\t" \
00447    "cmpl $0,%[pow]     \n\t" \
00448    "movl %[tmp],%%eax  \n\t" \
00449    "movl %%edx,%[X]    \n\t" \
00450    "jne 0b             \n\t" \
00451    "1:  " \
00452 //   : [X] "=&r" (x), "+a" (a), [pow] "+rm" (pow), [tmp] "=&rm" (temp) : [M] "rm" (m) : "cc", "edx");
00453 // FIXME: "can't find a register in class 'AD_REGS' while reloading 'asm'"
00454 // if the above line triggers compilation errors, because gcc refuses to choose an alternative which works,
00455 // then you need to use suboptimal constraints (use always memory operand instead of possible register):
00456    : [X] "=&r" (x), "+a" (a), [pow] "+rm" (pow), [tmp] "=m" (temp) : [M] "rm" (m) : "cc", "edx");
00457 
00458   return x;
00459 }
00460 
00461 #else
00462 // *** generic variant of powmod ***
00463 
00471 inline unsigned int powmod(register unsigned int a, register unsigned int pow, const unsigned int m)
00472 {
00473   // assumes a<m
00474   // returns (a^pow) mod m
00475   unsigned int x = (pow&1) ? a : 1;
00476   while (pow>>=1)
00477    {
00478      a=squaremod(a,m);
00479      if (pow&1) x=mulmod(x,a,m);
00480    }
00481   return x;
00482 }
00483 #endif
00484 
00485 
00491 bool strong_probab_prime(const unsigned int p, const unsigned int base);
00492 
00493 
00502 bool probab_prime(const unsigned int p);
00503 
00504 
00512 bool is_prime(register const int p);
00513 
00514 
00518 inline unsigned int gcd(register unsigned int a, register unsigned int b)
00519 {
00520   // returns greatest common divisor
00521   while (a&&b) if (a>b) a%=b; else b%=a;
00522   return a ? a : b;
00523 }
00524 
00525 inline unsigned short int gcd(register unsigned short int a, register unsigned short int b)
00526 {
00527   // returns greatest common divisor
00528   while (a&&b) if (a>b) a%=b; else b%=a;
00529   return a ? a : b;
00530 }
00531 
00542 inline unsigned int oddgcd(register unsigned int a, register unsigned int b)
00543 {
00544 #if defined(ASM_CMOV) || defined(ASM_X86_64)
00545  // avoid branches since cmovs are faster on random data and enable parallel execution
00546 
00547 #ifdef DEBUG
00548  #warning "using cmov-operation for oddgcd(unsigned int, unsigned int)"
00549 #endif
00550 
00551   asm ( \
00552    "cmpl %[a],%[b] \n\t" \
00553    "je 2f \n\t" \
00554    "1: movl %[a],%%ecx \n\t" \
00555    "shrl $1,%[a] \n\t" \
00556    "cmovcl %%ecx,%[a] \n\t" \
00557    "movl %[b],%%edx \n\t" \
00558    "shrl $1,%[b] \n\t" \
00559    "cmovcl %%edx,%[b] \n\t" \
00560    "movl %[a],%%ecx \n\t" \
00561    "subl %[b],%[a] \n\t" \
00562    "cmovbe %%ecx,%[a] \n\t" \
00563    "movl %[b],%%edx \n\t" \
00564    "subl %[a],%[b] \n\t" \
00565    "cmovbe %%edx,%[b] \n\t" \
00566    "jnz 1b \n\t" \
00567    "2: \n" \
00568    : [a] "+r" (a), [b] "+r" (b) : : "ecx", "edx", "cc");
00569   return a;
00570 #else
00571   // *** generic version ***
00572   // returns greatest (odd) common divisor
00573   if (a==0) return b;
00574   if (b==0) return a;
00575   do
00576    {
00577      //cout << a << "," << b << endl;
00578      while (1&a^1) a>>=1;
00579      while (1&b^1) b>>=1;
00580      if (a>b) a-=b;
00581      if (b>a) b-=a;
00582    } while (a!=b);
00583   //cout << a << endl;
00584   return a;
00585 #endif
00586 }
00587 
00588 
00594 inline bool coprime(register unsigned int a, register unsigned int b)
00595 {
00596   return ((a|b)&1) ? oddgcd(a,b)==1 : false;
00597 }
00598 
00599 
00605 signed int jacobi (int a, unsigned int b);
00606 
00607 
00620 signed int legendre (signed int a, const unsigned int p);
00621 
00622 
00631 unsigned int invmod(const unsigned int x, const unsigned int m) __attribute__ ((const));
00632 
00633 
00643 unsigned int bininvmod(register unsigned int x, register const unsigned int m) __attribute__ ((const));
00644 
00645 
00655 unsigned int montgominvmod(register unsigned int x, unsigned int m) __attribute__ ((const));
00656 
00657 
00670 inline unsigned int fastinvmod(const unsigned int x, const unsigned int m)
00671 {
00672 #if defined(DEBUG)
00673   if (x==0 || x>=m) std::cerr << "fastinvmod: constraint 0<y<m not met: x=" << x << ", m=" << m << std::endl;
00674   if (!coprime(x,m)) std::cerr << "fastinvmod: constraint coprime(x,m) not met: x=" << x << ", m=" << m << std::endl;
00675   if ((m&1)==0) std::cerr << "fastinvmod: constraint m=1 (mod 2)  [=odd m] not met: x=" << x << ", m=" << m << std::endl;
00676 #endif 
00677 
00678 #if defined (ASM_386) || defined(ASM_X86_64)
00679   //return invmod(x,m);
00680   //return bininvmod(x,m);
00681   return montgominvmod(x,m);
00682 #else
00683   //return invmod(x,m);
00684   return bininvmod(x,m);
00685   //return montgominvmod(x,m);
00686 #endif
00687 }
00688 
00689 
00690 #if 0
00691  // use this, if
00692  // a) there is no crafted invmod support for your cpu
00693  // and b) your cpu supports very fast floating point arithmetic
00694  #define USE_JASONS_INVMOD
00695  namespace japinvmod
00696  {
00697    unsigned int fast_invmod(unsigned int a, unsigned int p);
00698  }
00699 #endif
00700 
00701 
00702 inline unsigned int fastinvmod_23bit(const unsigned int x, const unsigned int m)
00703 {
00704 #if defined(DEBUG)
00705   if (x==0 || x>=m) std::cerr << "fastinvmod_23bit: constraint 0<y<m not met: x=" << x << ", m=" << m << std::endl;
00706   if (!coprime(x,m)) std::cerr << "fastinvmod_23bit: constraint coprime(x,m) not met: x=" << x << ", m=" << m << std::endl;
00707   if ((m&1)==0) std::cerr << "fastinvmod_23bit: constraint m=1 (mod 2)  [=odd m] not met: x=" << x << ", m=" << m << std::endl;
00708   if (m>8388608) std::cerr << "fastinvmod_23bit: m exceeds 23 bit: x=" << x << ", m=" << m << std::endl;
00709 #endif 
00710 
00711 #if defined(USE_JASONS_INVMOD)
00712   return japinvmod::fast_invmod(x,m);
00713 #else
00714   //return invmod(x,m);
00715   //return bininvmod(x,m);
00716   return montgominvmod(x,m);
00717 #endif
00718 }
00719 
00720 
00721 
00722 
00736 unsigned int sqrtmod(const unsigned int Radikant, const unsigned int Primzahl);
00737 
00738 
00739 } // namespace numtheory
00740 
00741 #endif /* MODULO_H_INCLUDED */

Generated on Wed Nov 7 23:29:25 2007 for Qsieve by  doxygen 1.5.4