src/modinteger/cl_MI_pow2m1.h

   1 // m > 0, m = 2^m1 - 1 (m1 > 1)
   2
   3 class cl_heap_modint_ring_pow2m1 : public cl_heap_modint_ring {
   4         SUBCLASS_cl_heap_modint_ring()
   5 public:
   6         // Constructor.
   7         cl_heap_modint_ring_pow2m1 (const cl_I& m, uintL m1); // m = 2^m1 - 1
   8         // Virtual destructor.
   9         ~cl_heap_modint_ring_pow2m1 () {}
  10         // Additional information.
  11         uintL m1;
  12 };
  13
  14 static
  15 #if !(defined(__GNUC__) && (__GNUC__ == 2) && (__GNUC_MINOR__ <= 90)) // workaround g++-2.7.2 and egcs-1.0.2-prerelease bug
  16 inline
  17 #endif
  18 const cl_I pow2m1_reduce_modulo (cl_heap_modint_ring* _R, const cl_I& x)
  19 {
  20         var cl_heap_modint_ring_pow2m1* R = (cl_heap_modint_ring_pow2m1*)_R;
  21         // Method:
  22         // If x>=0, split x into pieces of m1 bits and sum them up.
  23         //   x = x0 + 2^m1*x1 + 2^(2*m1)*x2 + ... ==>
  24         //   mod(x,m) = mod(x0+x1+x2+...,m).
  25         // If x<0, apply this to -1-x, and use mod(x,m) = m-1-mod(-1-x,m).
  26  {      Mutable(cl_I,x);
  27         var cl_boolean sign = minusp(x);
  28         if (sign) { x = lognot(x); }
  29         var const uintL m1 = R->m1;
  30         if (x >= R->modulus) {
  31                 x = plus1(x); // avoid staying at x = m
  32                 do {
  33                         var uintL xlen = integer_length(x);
  34                         var cl_I y = ldb(x,cl_byte(m1,0));
  35                         for (var uintL i = m1; i < xlen; i += m1)
  36                                 y = y + ldb(x,cl_byte(m1,i));
  37                         x = y;
  38                 } while (x > R->modulus);
  39                 x = minus1(x);
  40         }
  41         // Now 0 <= x < m.
  42         if (sign) { x = R->modulus - 1 - x; }
  43         return x;
  44 }}
  45
  46 static const _cl_MI pow2m1_canonhom (cl_heap_modint_ring* R, const cl_I& x)
  47 {
  48         return _cl_MI(R, pow2m1_reduce_modulo(R,x));
  49 }
  50
  51 static const _cl_MI pow2m1_mul (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
  52 {
  53         var cl_heap_modint_ring_pow2m1* R = (cl_heap_modint_ring_pow2m1*)_R;
  54         var const uintL m1 = R->m1;
  55         var cl_I zr = x.rep * y.rep;
  56         zr = ldb(zr,cl_byte(m1,m1)) + ldb(zr,cl_byte(m1,0));
  57         return _cl_MI(R, zr >= R->modulus ? zr - R->modulus : zr);
  58 }
  59
  60 static const _cl_MI pow2m1_square (cl_heap_modint_ring* _R, const _cl_MI& x)
  61 {
  62         var cl_heap_modint_ring_pow2m1* R = (cl_heap_modint_ring_pow2m1*)_R;
  63         var const uintL m1 = R->m1;
  64         var cl_I zr = square(x.rep);
  65         zr = ldb(zr,cl_byte(m1,m1)) + ldb(zr,cl_byte(m1,0));
  66         return _cl_MI(R, zr >= R->modulus ? zr - R->modulus : zr);
  67 }
  68
  69 #define pow2m1_addops std_addops
  70 static cl_modint_mulops pow2m1_mulops = {
  71         std_one,
  72         pow2m1_canonhom,
  73         pow2m1_mul,
  74         pow2m1_square,
  75         std_expt_pos,
  76         std_recip,
  77         std_div,
  78         std_expt,
  79         pow2m1_reduce_modulo,
  80         std_retract
  81 };
  82
  83 // Constructor.
  84 inline cl_heap_modint_ring_pow2m1::cl_heap_modint_ring_pow2m1 (const cl_I& m, uintL _m1)
  85         : cl_heap_modint_ring (m, &std_setops, &pow2m1_addops, &pow2m1_mulops), m1 (_m1) {}