src/modinteger/cl_MI_pow2.h

   1 // m > 0, m = 2^m1
   2
   3 class cl_heap_modint_ring_pow2 : public cl_heap_modint_ring {
   4         SUBCLASS_cl_heap_modint_ring()
   5 public:
   6         // Constructor.
   7         cl_heap_modint_ring_pow2 (const cl_I& m, uintL m1); // m = 2^m1
   8         // Virtual destructor.
   9         ~cl_heap_modint_ring_pow2 () {}
  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 pow2_reduce_modulo (cl_heap_modint_ring* _R, const cl_I& x)
  19 {
  20         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
  21         return ldb(x,cl_byte(R->m1,0));
  22 }
  23
  24 static const _cl_MI pow2_canonhom (cl_heap_modint_ring* R, const cl_I& x)
  25 {
  26         return _cl_MI(R, pow2_reduce_modulo(R,x));
  27 }
  28
  29 static const _cl_MI pow2_plus (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
  30 {
  31         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
  32         var cl_I zr = x.rep + y.rep;
  33         return _cl_MI(R, ldb(zr,cl_byte(R->m1,0)));
  34 }
  35
  36 static const _cl_MI pow2_minus (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
  37 {
  38         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
  39         var cl_I zr = x.rep - y.rep;
  40         return _cl_MI(R, ldb(zr,cl_byte(R->m1,0)));
  41 }
  42
  43 static const _cl_MI pow2_uminus (cl_heap_modint_ring* _R, const _cl_MI& x)
  44 {
  45         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
  46         var cl_I zr = - x.rep;
  47         return _cl_MI(R, ldb(zr,cl_byte(R->m1,0)));
  48 }
  49
  50 static const _cl_MI pow2_one (cl_heap_modint_ring* _R)
  51 {
  52         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
  53         return _cl_MI(R, R->m1==0 ? 0 : 1);
  54 }
  55
  56 static const _cl_MI pow2_mul (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
  57 {
  58         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
  59         var cl_I zr = x.rep * y.rep;
  60         return _cl_MI(R, ldb(zr,cl_byte(R->m1,0)));
  61 }
  62
  63 static const _cl_MI pow2_square (cl_heap_modint_ring* _R, const _cl_MI& x)
  64 {
  65         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
  66         var cl_I zr = square(x.rep);
  67         return _cl_MI(R, ldb(zr,cl_byte(R->m1,0)));
  68 }
  69
  70 // Timing comparison with std_recip, on a i486 33 MHz running Linux:
  71 // timeMIpow2recip N inverts an (N*32)-bit number.
  72 //   N  std_recip pow2_recip
  73 //   10   0.0030  0.00017
  74 //   25   0.011   0.00068
  75 //   50   0.035   0.0024
  76 //  100   0.124   0.0090
  77 //  250   0.71    0.055
  78 //  500   2.76    0.193
  79 // 1000  11.0     0.61
  80 // 2500  68.7     2.2
  81 // 5000 283       5.0
  82 static const cl_MI_x pow2_recip (cl_heap_modint_ring* _R, const _cl_MI& x)
  83 {
  84         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
  85         var const cl_I& xr = x.rep;
  86         if (!oddp(xr)) {
  87                 if (R->m1 == 0)
  88                         return cl_MI(R, 0);
  89                 if (zerop(xr))
  90                         cl_error_division_by_0();
  91                 else
  92                         return cl_notify_composite(R,xr);
  93         } else
  94                 return cl_MI(R, cl_recip2adic(R->m1,xr));
  95 }
  96
  97 // Timing comparison with std_div, on a i486 33 MHz running Linux:
  98 // timeMIpow2div N divides two (N*32)-bit numbers.
  99 //   N   std_div  pow2_div
 100 //   10   0.0035  0.00017
 101 //   25   0.0136  0.00068
 102 //   50   0.043   0.0024
 103 //  100   0.151   0.0090
 104 //  250   0.85    0.054
 105 //  500   3.3     0.21
 106 // 1000  12.3     0.86
 107 static const cl_MI_x pow2_div (cl_heap_modint_ring* _R, const _cl_MI& x, const _cl_MI& y)
 108 {
 109         var cl_heap_modint_ring_pow2* R = (cl_heap_modint_ring_pow2*)_R;
 110         var const cl_I& yr = y.rep;
 111         if (!oddp(yr)) {
 112                 if (R->m1 == 0)
 113                         return cl_MI(R, 0);
 114                 if (zerop(yr))
 115                         cl_error_division_by_0();
 116                 else
 117                         return cl_notify_composite(R,yr);
 118         } else
 119                 return cl_MI(R, cl_div2adic(R->m1,x.rep,yr));
 120 }
 121
 122 static cl_modint_addops pow2_addops = {
 123         std_zero,
 124         std_zerop,
 125         pow2_plus,
 126         pow2_minus,
 127         pow2_uminus
 128 };
 129 static cl_modint_mulops pow2_mulops = {
 130         pow2_one,
 131         pow2_canonhom,
 132         pow2_mul,
 133         pow2_square,
 134         std_expt_pos,
 135         pow2_recip,
 136         pow2_div,
 137         std_expt,
 138         pow2_reduce_modulo,
 139         std_retract
 140 };
 141
 142 // Constructor.
 143 inline cl_heap_modint_ring_pow2::cl_heap_modint_ring_pow2 (const cl_I& m, uintL _m1)
 144         : cl_heap_modint_ring (m, &std_setops, &pow2_addops, &pow2_mulops), m1 (_m1) {}