guru/sing.c

   1 /*
   2  * Example of a singular curve, similar to 19.c
   3  * but the Tate pairing degenerates
   4  *
   5  * Consider the curve E: y^2 = x^3 + x^2 over F_19:
   6  * E_ns(F_19) is a cyclic group of order 18.
   7  */
   8
   9 #include "pbc.h"
  10 #include "pbc_singular.h"
  11 #include "pbc_fp.h"
  12
  13 static void miller(element_t res, element_t P, element_t Q, element_t R, int n)
  14 {
  15     //collate divisions
  16     mp_bitcnt_t m;
  17     element_t v, vd;
  18     element_t Z;
  19     element_t a, b, c;
  20     element_t e0, e1;
  21     mpz_t q;
  22     element_ptr Zx, Zy;
  23     const element_ptr Px = curve_x_coord(P);
  24     const element_ptr Py = curve_y_coord(P);
  25     const element_ptr numx = curve_x_coord(Q);
  26     const element_ptr numy = curve_y_coord(Q);
  27     const element_ptr denomx = curve_x_coord(R);
  28     const element_ptr denomy = curve_y_coord(R);
  29
  30     #define do_vertical(e, edenom) {         \
  31         element_sub(e0, numx, Zx);           \
  32         element_mul((e), (e), e0);           \
  33                                              \
  34         element_sub(e0, denomx, Zx);         \
  35         element_mul((edenom), (edenom), e0); \
  36     }
  37
  38     #define do_tangent(e, edenom) {                  \
  39         /*a = -slope_tangent(A.x, A.y);              \
  40           b = 1;                                     \
  41           c = -(A.y + a * A.x);                      \
  42           but we multiply by 2*A.y to avoid division \
  43                                                      \
  44           a = -Ax * (Ax + Ax + Ax + twicea_2) - a_4; \
  45           This curve is special:                     \
  46           a = -(3 Ax^2 + 2Ax)                        \
  47           b = 2 * Ay                                 \
  48           c = -(2 Ay^2 + a Ax);                   */ \
  49                                                      \
  50         if (element_is0(Zy)) {                       \
  51             do_vertical((e), (edenom));              \
  52         } else {                                     \
  53           element_square(a, Zx);                     \
  54           element_mul_si(a, a, 3);                   \
  55           element_add(a, a, Zx);                     \
  56           element_add(a, a, Zx);                     \
  57           element_neg(a, a);                         \
  58                                                      \
  59           element_add(b, Zy, Zy);                    \
  60                                                      \
  61           element_mul(e0, b, Zy);                    \
  62           element_mul(c, a, Zx);                     \
  63           element_add(c, c, e0);                     \
  64           element_neg(c, c);                         \
  65                                                      \
  66           element_mul(e0, a, numx);                  \
  67           element_mul(e1, b, numy);                  \
  68           element_add(e0, e0, e1);                   \
  69           element_add(e0, e0, c);                    \
  70           element_mul((e), (e), e0);                 \
  71                                                      \
  72           element_mul(e0, a, denomx);                \
  73           element_mul(e1, b, denomy);                \
  74           element_add(e0, e0, e1);                   \
  75           element_add(e0, e0, c);                    \
  76           element_mul((edenom), (edenom), e0);       \
  77         }                                            \
  78     }
  79
  80     #define do_line(e, edenom) {               \
  81         if (!element_cmp(Zx, Px)) {            \
  82             if (!element_cmp(Zy, Py)) {        \
  83                 do_tangent((e), (edenom));     \
  84             } else {                           \
  85                 do_vertical((e), (edenom));    \
  86             }                                  \
  87         } else {                               \
  88           element_sub(b, Px, Zx);              \
  89           element_sub(a, Zy, Py);              \
  90           element_mul(c, Zx, Py);              \
  91           element_mul(e0, Zy, Px);             \
  92           element_sub(c, c, e0);               \
  93                                                \
  94           element_mul(e0, a, numx);            \
  95           element_mul(e1, b, numy);            \
  96           element_add(e0, e0, e1);             \
  97           element_add(e0, e0, c);              \
  98           element_mul((e), (e), e0);           \
  99                                                \
 100           element_mul(e0, a, denomx);          \
 101           element_mul(e1, b, denomy);          \
 102           element_add(e0, e0, e1);             \
 103           element_add(e0, e0, c);              \
 104           element_mul((edenom), (edenom), e0); \
 105        }                                       \
 106     }
 107
 108     element_init(a, res->field);
 109     element_init(b, res->field);
 110     element_init(c, res->field);
 111     element_init(e0, res->field);
 112     element_init(e1, res->field);
 113
 114     element_init(v, res->field);
 115     element_init(vd, res->field);
 116     element_init(Z, P->field);
 117
 118     element_set(Z, P);
 119     Zx = curve_x_coord(Z);
 120     Zy = curve_y_coord(Z);
 121
 122     element_set1(v);
 123     element_set1(vd);
 124
 125     mpz_init(q);
 126     mpz_set_ui(q, n);
 127     m = (mp_bitcnt_t)mpz_sizeinbase(q, 2);
 128     m = (m > 2 ? m - 2 : 0);
 129
 130     for (;;) {
 131         element_square(v, v);
 132         element_square(vd, vd);
 133         do_tangent(v, vd);
 134         element_double(Z, Z);
 135         do_vertical(vd, v);
 136
 137         if (mpz_tstbit(q, m)) {
 138             do_line(v, vd);
 139             element_add(Z, Z, P);
 140             if (m) {
 141                 do_vertical(vd, v);
 142             }
 143         }
 144         if (!m) break;
 145         m--;
 146     }
 147
 148     mpz_clear(q);
 149
 150     element_invert(vd, vd);
 151     element_mul(res, v, vd);
 152
 153     element_clear(v);
 154     element_clear(vd);
 155     element_clear(Z);
 156     element_clear(a);
 157     element_clear(b);
 158     element_clear(c);
 159     element_clear(e0);
 160     element_clear(e1);
 161     #undef do_vertical
 162     #undef do_tangent
 163     #undef do_line
 164 }
 165
 166 static void tate_3(element_ptr out, element_ptr P, element_ptr Q, element_ptr R)
 167 {
 168     mpz_t six;
 169
 170     mpz_init(six);
 171     mpz_set_ui(six, 6);
 172     element_t QR;
 173     element_t e0;
 174
 175     element_init(QR, P->field);
 176     element_init(e0, out->field);
 177
 178     element_add(QR, Q, R);
 179
 180     //for subgroup size 3, -2P = P, hence
 181     //the tangent line at P has divisor 3(P) - 3(O)
 182
 183     miller(out, P, QR, R, 3);
 184
 185     element_pow_mpz(out, out, six);
 186     element_clear(QR);
 187     element_clear(e0);
 188     mpz_clear(six);
 189 }
 190
 191 static void tate_9(element_ptr out, element_ptr P, element_ptr Q, element_ptr R)
 192 {
 193     element_t QR;
 194     element_init(QR, P->field);
 195
 196     element_add(QR, Q, R);
 197
 198     miller(out, P, QR, R, 9);
 199
 200     element_square(out, out);
 201
 202     element_clear(QR);
 203 }
 204
 205 int main(void)
 206 {
 207     field_t c;
 208     field_t Z19;
 209     element_t P, Q, R;
 210     mpz_t q, z;
 211     element_t a;
 212     int i;
 213
 214     mpz_init(q);
 215     mpz_init(z);
 216
 217     mpz_set_ui(q, 19);
 218
 219     field_init_fp(Z19, q);
 220     element_init(a, Z19);
 221
 222     field_init_curve_singular_with_node(c, Z19);
 223
 224     element_init(P, c);
 225     element_init(Q, c);
 226     element_init(R, c);
 227
 228     //(3,+/-6) is a generator
 229     //we have an isomorphism from E_ns to F_19^*
 230     // (3,6) --> 3
 231     //(generally (x,y) --> (y+x)/(y-x)
 232
 233     curve_set_si(R, 3, 6);
 234
 235     for (i=1; i<=18; i++) {
 236         mpz_set_si(z, i);
 237         element_mul_mpz(Q, R, z);
 238         element_printf("%dR = %B\n", i, Q);
 239     }
 240
 241     mpz_set_ui(z, 6);
 242     element_mul_mpz(P, R, z);
 243     //P has order 3
 244     element_printf("P = %B\n", P);
 245
 246     for (i=1; i<=3; i++) {
 247         mpz_set_si(z, i);
 248         element_mul_mpz(Q, R, z);
 249         tate_3(a, P, Q, R);
 250         element_printf("e_3(P,%dP) = %B\n", i, a);
 251     }
 252
 253     element_double(P, R);
 254     //P has order 9
 255     element_printf("P = %B\n", P);
 256     for (i=1; i<=9; i++) {
 257         mpz_set_si(z, i);
 258         element_mul_mpz(Q, P, z);
 259         tate_9(a, P, Q, R);
 260         element_printf("e_9(P,%dP) = %B\n", i, a);
 261     }
 262
 263     return 0;
 264 }