src/.svn/text-base/add_explicit.C.svn-base

   1 /* Genus 2 divisor group law via explicit formulas.
   2
   3 Divisor addition: Algorithm 14.19 of "Handbook of EAHCC"
   4
   5 Divisor doubling: Algorithm 14.21 of "Handbook of EAHCC"
   6
   7 Special cases are handled by calling add_cantor()
   8
   9 */
  10
  11 #if HAVE_CONFIG_H
  12 #include <config.h>
  13 #endif
  14
  15 #include <assert.h>
  16 #include <g2hec_Genus2_ops.h>
  17
  18 NS_G2_START_IMPL
  19
  20 static bool_t add_diff(divisor& x, const divisor& a, const divisor& b)
  21      // x = a + b where a != b via explicit formulas
  22      // c.f. Algorithm 14.9 of HOEHCC
  23      // Assumes upolys of a and b are both of degree 2 and have no common
  24      // factor.
  25 {
  26   /* Variables
  27 We choose to use local variables for now to make the function thread safe.
  28 Note that this may affect the performance of the program as opposed to
  29 the method of passing along the temporary variables.
  30   */
  31
  32   bool_t OK = TRUE;
  33
  34   const field_t f4 = coeff(a.get_curve().get_f(), 4),
  35                 h1 = coeff(a.get_curve().get_h(), 1),
  36                 h2 = coeff(a.get_curve().get_h(), 2),
  37                 h0 = coeff(a.get_curve().get_h(), 0),
  38     u11 = coeff(a.get_upoly(), 1), u10 = coeff(a.get_upoly(), 0),
  39     u21 = coeff(b.get_upoly(), 1), u20 = coeff(b.get_upoly(), 0),
  40     v11 = coeff(a.get_vpoly(), 1), v10 = coeff(a.get_vpoly(), 0),
  41     v21 = coeff(b.get_vpoly(), 1), v20 = coeff(b.get_vpoly(), 0);
  42
  43   field_t z1, z2, z3, z4, r, inv1, inv0, w0, w1, w2, w3, s1p, s0p,
  44     w4, w5, s0pp, l2p, l1p, l0p, u0p, u1p, v1p, v0p;
  45
  46   poly_t up, vp;
  47
  48   /* Explicit formulas */
  49
  50   /* Step 1:
  51      Compute r = Res(u1, u2).
  52    z1 = u11 - u21, z2 = u20 - u10, z3 = u11*z1 + z2, z4 = z1*u10,
  53    r = z2*z3 + z4*z1
  54   */
  55   z1 = u11 - u21;
  56   z2 = u20 - u10;
  57   z3 = u11*z1 + z2;
  58   z4 = z1*u10;
  59   r = z2*z3 + z4*z1;
  60
  61   assert(!IsZero(r));  // Res(u1, u2)<>0, i.e., GCD(u1, u2) = 1
  62
  63   /* Step 2:
  64      inv1 = z1, inv0 = z3
  65   */
  66   inv1 = z1;
  67   inv0 = z3;
  68
  69   /* Step 3:
  70      w0 = v10 - v20, w1 = v11 - v21, w2 = inv0*w0, w3 = inv1*w1,
  71      s1p = (inv0 + inv1)*(w0 + w1) - w2 - w3*(1 + u11), s0p = w2 - z4*w1
  72      If s1p = 0, handled by Cantor's algorithm.
  73   */
  74   w0 = v10 - v20;
  75   w1 = v11 - v21;
  76   w2 = inv0*w0;
  77   w3 = inv1*w1;
  78   s1p = (inv0 + inv1)*(w0 + w1) - w2 - w3*(1 + u11);
  79
  80   // Special case
  81   if (IsZero(s1p))
  82     return add_cantor(x, a, b);
  83
  84   s0p = w2 - z4*w1;
  85
  86
  87   /* Step 4:
  88      w1 = 1/(r*s1p), w2 = r*w1, w3 = s1p^2*w1, w4 = r*w2, w5 = w4^2,
  89      s0pp = s0p*w2
  90   */
  91   w1 = 1/(r*s1p);
  92   w2 = r*w1;
  93   w3 = sqr(s1p)*w1;
  94   w4 = r*w2;
  95   w5 = sqr(w4);
  96   s0pp = s0p*w2;
  97
  98   /* Step 5:
  99      l2p = u21 + s0pp, l1p = u21*s0pp + u20, l0p = u20*s0pp
 100   */
 101   l2p = u21 + s0pp;
 102   l1p = u21*s0pp + u20;
 103   l0p = u20*s0pp;
 104
 105   /* Step 6:
 106      u0p = (s0pp - u11)*(s0pp - z1 + h2*w4) - u10,
 107      u0p = u0p + l1p + (h1 + 2*v21)*w4 + (2*u21 + z1 - f4)*w5,
 108      u1p = 2*s0pp - z1 + h2*w4 - w5
 109   */
 110   u0p = (s0pp - u11)*(s0pp - z1 + h2*w4) - u10;
 111   u0p = u0p + l1p + (h1 + 2*v21)*w4 + (2*u21 + z1 - f4)*w5;
 112   u1p = 2*s0pp - z1 + h2*w4 - w5;
 113
 114   /* Step 7:
 115      w1 = l2p - u1p, w2 = u1p*w1 + u0p - l1p, v1p = w2*w3 - v21 - h1 + h2*u1p,
 116      w2 = u0p*w1 - l0p, v0p = w2*w3 - v20 - h0 + h2*u0p
 117   */
 118   w1 = l2p - u1p;
 119   w2 = u1p*w1 + u0p - l1p;
 120   v1p = w2*w3 - v21 - h1 + h2*u1p;
 121   w2 = u0p*w1 - l0p;
 122   v0p = w2*w3 - v20 - h0 + h2*u0p;
 123
 124   /* Step 8:
 125      return [up, vp] = [x^2 + u1p*x + u0p, v1p*x + v0p]
 126   */
 127   SetCoeff(up, 2, 1);
 128   SetCoeff(up, 1, u1p);
 129   SetCoeff(up, 0, u0p);
 130
 131   SetCoeff(vp, 1, v1p);
 132   SetCoeff(vp, 0, v0p);
 133
 134   x.set_upoly(up);
 135   x.set_vpoly(vp);
 136
 137   x.update();
 138
 139 #if DEBUG_LEVEL > 1
 140   cout << "Entered add_diff()" << endl;
 141   cout << x << endl;
 142 #endif
 143
 144   assert(OK = ( OK && x.is_valid_divisor()) );
 145
 146   return OK;
 147 }
 148
 149 static bool_t doubling(divisor& x, const divisor& a)
 150      // x = a + a via explicit formulas
 151
 152 {  /* Variables
 153 We choose to use local variables for now to make the function thread safe.
 154 Note that this may affect the performance of the program as opposed to
 155 the method of passing along the temporary variables.
 156   */
 157
 158   bool_t OK = TRUE;
 159
 160
 161   const field_t f4 = coeff(a.get_curve().get_f(), 4),
 162                 f3 = coeff(a.get_curve().get_f(), 3),
 163                 f2 = coeff(a.get_curve().get_f(), 2),
 164                 h1 = coeff(a.get_curve().get_h(), 1),
 165                 h2 = coeff(a.get_curve().get_h(), 2),
 166                 h0 = coeff(a.get_curve().get_h(), 0),
 167     u1 = coeff(a.get_upoly(), 1), u0 = coeff(a.get_upoly(), 0),
 168     v1 = coeff(a.get_vpoly(), 1), v0 = coeff(a.get_vpoly(), 0);
 169
 170   field_t v1t, v0t, w0, w1, w2, w3, r, inv1p, inv0p, w4, t1p, t0p, s1p,
 171     s0p, w5, s0pp, l2p, l1p, l0p, u0p, u1p, v1p, v0p;
 172
 173   poly_t up, vp;
 174
 175   /* Explicit formulas */
 176
 177   /* Step 1:
 178      v1t = h1 + 2*v1 - h2*u1, v0t = h0 + 2*v0 - h2*u0
 179   */
 180   v1t = h1 + 2*v1 - h2*u1;
 181   v0t = h0 + 2*v0 - h2*u0;
 182
 183   /* Step 2:
 184      w0 = v1^2, w1 = u1^2, w2 = v1t^2, w3 = u1*v1t,
 185      r = u0*w2 + v0t*(v0t - w3)
 186   */
 187   w0 = sqr(v1);
 188   w1 = sqr(u1);
 189   w2 = sqr(v1t);
 190   w3 = u1*v1t;
 191   r = u0*w2 + v0t*(v0t - w3);
 192   assert(!IsZero(r));
 193
 194   /* Step 3:
 195      inv1p = -v1t, inv0p = v0t - w3
 196   */
 197    inv1p = -v1t;
 198    inv0p = v0t - w3;
 199
 200   /* Step 4:
 201      w3 = f3 + w1, w4 = 2*u0, t1p = 2*(w1 - f4*u1) + w3 - w4 - h2*v1,
 202      t0p = u1*(2*w4 - w3 + f4*u1 + h2*v1) + f2 - w0 - 2*f4*u0 - h1*v1 - h2*v0
 203   */
 204    w3 = f3 + w1;
 205    w4 = 2*u0;
 206    t1p = 2*(w1 - f4*u1) + w3 - w4 - h2*v1;
 207    t0p = u1*(2*w4 - w3 + f4*u1 + h2*v1) + f2 - w0 - 2*f4*u0 - h1*v1 - h2*v0;
 208
 209   /* Step 5:
 210      w0 = t0p*inv0p, w1 = t1p*inv1p,
 211      s1p = (inv0p + inv1p)*(t0p + t1p) - w0 - w1*(1+u1),
 212      s0p = w0 - u0*w1.
 213
 214      If s1p = 0, call add_cantor().
 215   */
 216    w0 = t0p*inv0p;
 217    w1 = t1p*inv1p;
 218    s1p = (inv0p + inv1p)*(t0p + t1p) - w0 - w1*(1+u1);
 219
 220 #if DEBUG_LEVEL > 1
 221    cout << "r = " << r << endl;
 222    cout << "s1p = " << s1p << endl;
 223 #endif
 224
 225   // Special case
 226   if (IsZero(s1p))
 227     return add_cantor(x, a, a);
 228
 229    s0p = w0 - u0*w1;
 230
 231   /* Step 6:
 232      w1 = 1/(r*s1p), w2 = r*w1, w3 = s1p^2*w1, w4 = r*w2, w5 = w4^2,
 233      s0pp = s0p*w2
 234    */
 235   w1 = 1/(r*s1p);
 236   w2 = r*w1;
 237   w3 = sqr(s1p)*w1;
 238   w4 = r*w2;
 239   w5 = sqr(w4);
 240   s0pp = s0p*w2;
 241
 242   /* Step 7:
 243      l2p = u1 + s0pp, l1p = u1*s0pp + u0, l0p = u0*s0pp
 244   */
 245    l2p = u1 + s0pp;
 246    l1p = u1*s0pp + u0;
 247    l0p = u0*s0pp;
 248
 249   /* Step 8:
 250      u0p = s0pp^2 + w4*(h2*(s0pp - u1) + 2*v1 + h1) + w5*(2*u1 - f4),
 251      u1p = 2*s0pp + h2*w4 - w5
 252   */
 253    u0p = sqr(s0pp) + w4*(h2*(s0pp - u1) + 2*v1 + h1) + w5*(2*u1 - f4);
 254    u1p = 2*s0pp + h2*w4 - w5;
 255
 256   /* Step 9:
 257      w1 = l2p - u1p, w2 = u1p*w1 + u0p - l1p, v1p = w2*w3 - v1 - h1 + h2*u1p,
 258      w2 = u0p*w1 - l0p, v0p = w2*w3 - v0 - h0 + h2*u0p
 259   */
 260    w1 = l2p - u1p;
 261    w2 = u1p*w1 + u0p - l1p;
 262    v1p = w2*w3 - v1 - h1 + h2*u1p;
 263    w2 = u0p*w1 - l0p;
 264    v0p = w2*w3 - v0 - h0 + h2*u0p;
 265
 266   /* Step 10:
 267      return [up, vp]
 268   */
 269   SetCoeff(up, 2, 1);
 270   SetCoeff(up, 1, u1p);
 271   SetCoeff(up, 0, u0p);
 272
 273   SetCoeff(vp, 1, v1p);
 274   SetCoeff(vp, 0, v0p);
 275
 276   x.set_upoly(up);
 277   x.set_vpoly(vp);
 278
 279   x.update();
 280
 281 #if DEBUG_LEVEL > 1
 282   cout << "Entered doubling(). " << endl;
 283   //  cout << x << endl;
 284 #endif
 285
 286   assert(OK = OK && x.is_valid_divisor());
 287   return OK;
 288 }
 289
 290 bool_t add(divisor& x, const divisor& a, const divisor& b)
 291      // This subroutine wraps other functions that does the actual
 292      // divisor arithmetic. It checks the validity of input divisors
 293      // so that other subroutines it calls do not need to do so.
 294 {
 295   bool_t OK = TRUE;
 296
 297
 298   /* Reduce overhead of checking with NDEBUG flag */
 299   assert(OK = OK && a.is_valid_divisor() && b.is_valid_divisor());
 300
 301   if (deg(a.get_upoly()) == genus && deg(b.get_upoly()) == genus) {
 302
 303     if (a == - b) {
 304       x.set_unit();
 305       OK = TRUE;
 306       return OK;
 307     }
 308
 309     if (a != b && IsOne(GCD(a.get_upoly(), b.get_upoly()))) {
 310       // Addition
 311       OK = OK && add_diff(x, a, b);
 312       return OK;
 313
 314     } else if (a == b && IsOne(GCD(a.get_curve().get_h() +
 315                             2*a.get_vpoly(), a.get_upoly())) ) {
 316       // Doubling
 317       // Exclude the case when one point of the divisor is equal to
 318       // its opposite
 319
 320       OK = OK && doubling(x, a);
 321       return OK;
 322     }
 323
 324
 325   }
 326
 327    // Call add_cantor() to handle other cases
 328   OK = OK && add_cantor(x, a, b);
 329
 330   return OK;
 331 }
 332
 333 NS_G2_END_IMPL