src/modinteger/cl_MI_pow2.h

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