1 // points.cc: implementations for Point class for points on elliptic curves
2 //////////////////////////////////////////////////////////////////////////
3 //
4 // Copyright 1990-2012 John Cremona
5 //
6 // This file is part of the eclib package.
7 //
8 // eclib is free software; you can redistribute it and/or modify it
9 // under the terms of the GNU General Public License as published by the
10 // Free Software Foundation; either version 2 of the License, or (at your
11 // option) any later version.
12 //
13 // eclib is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 // for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with eclib; if not, write to the Free Software Foundation,
20 // Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
21 //
22 //////////////////////////////////////////////////////////////////////////
23
24 // originally adapted from Elliptic.cc by Oisin McGuiness
25
26 #include <eclib/points.h> // which includes curve.h
27 #include <eclib/cperiods.h>
28 #include <NTL/RR.h> // for the realify_point() function
29
30 //#define DEBUG_TORSION
31
32 //
33 // Point member functions
34 //
35
36 // Shifts P via T = (u, r, s, t) on to newc.
37 // N.B. Assumes that (1) T(P.E) = newc
38 // (2) T(P) is on newc
39 // "back" means use reverse transform
40
transform(const Point & P,Curvedata * newc,const bigint & u,const bigint & r,const bigint & s,const bigint & t,int back)41 Point transform(const Point& P, Curvedata* newc,
42 const bigint& u,
43 const bigint& r, const bigint& s, const bigint& t,
44 int back)
45 {
46 if(P.is_zero()) return Point(newc);
47 if(!P.isvalid())
48 cout << "Attempting to trabsform the point " << P
49 << "which is not a valid point on its curve " << P.getcurve() << "!\n";
50 Point Q(newc,transform(P,u,r,s,t,back));
51 if(!Q.isvalid())
52 {
53 cout << "Result of transforming the point " << P << " on curve "<<(Curve)*(P.E) << " via [u,r,s,t]=["<<u<<","<<r<<","<<s<<","<<t<<"]";
54 if(back) cout<<" (inverse) ";
55 cout << " is " << Q << " which is not a valid point on its curve "
56 << (Curve)(*newc) << "!\n";
57 }
58 return Q;
59 }
60
61 // Genuine elliptic curve point functions:
62
63 // addition and subtraction
operator +=(const Point & Q)64 void Point::operator+=(const Point& Q) // P += Q;
65 {
66 Point sum = (*this) + Q ;
67 E = sum.E; X = sum.X; Y = sum.Y; Z = sum.Z; reduce();
68 ord = 0; height = -1;
69 }
70
operator -=(const Point & Q)71 void Point::operator-=(const Point& Q) // P -= Q ;
72 {
73 Point sum = -Q; sum+=(*this);
74 E = sum.E; X = sum.X; Y = sum.Y; Z = sum.Z; reduce();
75 ord = 0; height = -1;
76 }
77
operator +(const Point & Q) const78 Point Point::operator+(const Point& Q) const // P + Q
79 {
80 Point ans(E);
81 // make sure that both points are on the same curve
82 if(E != Q.E){
83 cerr << "## Can't add points on different curves!" << endl;
84 return ans;
85 }
86 // NOTE: we don't do any reductions here, but rely on assignment to do it
87 // check special cases first
88 if(Z==0) return Q ;
89 if(Q.is_zero()) return *this ;
90 if(eq(*this , Q)) {return Q.twice();}
91 Point minusQ=-Q;
92 if(eq(*this , minusQ)) {return ans;} // zero
93
94 // we now have genuine work to do
95 // let's set up some local variables to avoid repeated references
96 // coefficients
97 bigint A1,A2,A3,A4,A6; E->getai(A1,A2,A3,A4,A6);
98 // coordinates of P
99 const bigint& X1 = X ;
100 const bigint& Y1 = Y ;
101 const bigint& Z1 = Z ;
102 // coordinates of Q
103 const bigint& X2 = Q.X ;
104 const bigint& Y2 = Q.Y ;
105 const bigint& Z2 = Q.Z ;
106 const bigint& Z12 = Z1 * Z2 ;
107
108 const bigint& L = - Y2 * Z1 + Y1 * Z2 ; /* lambda */
109 const bigint& M = - X2 * Z1 + X1 * Z2 ; /* mu */
110 const bigint& N = - Y1 * X2 + Y2 * X1 ; /* nu */
111
112 const bigint& Mz = M * M * Z12 ;
113
114 const bigint& t = L * L * Z12 + M * ( A1 * L * Z12 - M * (A2 * Z12
115 + X1 * Z2 + X2 * Z1 ) ) ;
116 const bigint& newX = M * t ;
117 const bigint& newY = - ( t * (L + A1 * M) + Mz * (N + A3 * M )) ;
118 const bigint& newZ = M * Mz ;
119 ans.init(E, newX, newY, newZ);
120 return ans;
121 }
122
operator -(const Point & Q) const123 Point Point::operator-(const Point& Q) const // P - Q
124 {
125 Point ans = (-Q) ;
126 ans += (*this);
127 return ans;
128 }
129
twice(void) const130 Point Point::twice(void) const // doubles P
131 {
132 // do trivial cases
133 Point ans(E);
134 if( Z==0) return ans;
135 bigint A1,A2,A3,A4,A6; E->getai(A1,A2,A3,A4,A6);
136 Point minusthis = -(*this);
137 if(eq(*this,minusthis)) return ans; // order 2
138
139 const bigint& Zsq = Z * Z ;
140
141 const bigint& L = 3 * X * X + 2 * A2 * X * Z +
142 A4 * Zsq - A1 * Y * Z ;
143 const bigint& M = 2 * Y + A1 * X + A3 * Z ;
144
145 const bigint& Mz = M * Z ;
146 const bigint& N = - X * X * X -A3 * Y * Zsq + A4 * X *Zsq +
147 2 * A6 * Z *Zsq ;
148
149 const bigint& t = L * L + Mz * ( A1 * L - M * ( A2 * Z + 2 * X ) );
150
151 const bigint& newX = t * Mz ;
152 const bigint& newY = - ( L * t + Mz * ( A1 * t+ M * ( N + A3 * Mz * Z) ) );
153 const bigint& newZ = Mz * Mz * Mz;
154
155 ans.init(E, newX, newY, newZ) ;
156 return ans;
157 }
158
operator -(void) const159 Point Point::operator-(void) const // -Q
160 {
161 Point negation(*this);
162 negation.Y = -Y - (E->a3)*Z - (E->a1)*X;
163 return negation;
164 }
165
166 // calculates nP
operator *(int n,const Point & P)167 Point operator*(int n, const Point& P)
168 {
169 Point ans(*(P.E));
170 if(P.is_zero() || n == 0) return ans;
171 int negative = (n < 0) ;
172 if(negative) n = - n ;
173 if(n == 1) {
174 ans = P;
175 if(negative) ans = -P;
176 return ans ;
177 }
178 // now n >= 2
179 if(n == 2){
180 ans = P.twice() ;
181 if(negative) ans = -ans;
182 return ans ;
183 }
184 // now n >= 3
185 if(n&1) ans = P ; // else ans is ZERO
186 Point Q = P ;
187 while(n > 1){
188 Q = Q.twice() ; // 2^k P
189 n /= 2 ;
190 if(n&1) ans = ans + Q ;
191 }
192 if(negative) ans = -ans ;
193 return ans ;
194 }
195
order(Point & p)196 int order(Point& p)
197 {
198 // ASSUME that point is valid; check before calling if unknown
199 if (p.ord) {return p.ord;}
200 bigint eight, z=p.getZ(); eight=8;
201 if (is_zero(z)) {p.ord = 1; return 1; }
202 if (z>eight) {p.ord =-1; return -1;}
203 Point q = p; long ord=1;
204 while ( (sign(z)!=0) && (z<=eight) ) // (worst denom of a torsion point)
205 {ord++; q+=p; z = q.getZ(); }
206 if (sign(z)!=0) ord = -1;
207 p.ord = ord;
208 return ord;
209 }
210
order(Point & p,vector<Point> & multiples)211 int order(Point& p, vector<Point>& multiples)
212 {
213 #ifdef DEBUG_TORSION
214 cout<<"In order() with p = "<<p<<endl;
215 #endif
216 // ASSUME that point is valid; check before calling if unknown
217 multiples.clear(); multiples.reserve(13);
218 multiples.push_back(Point(*(p.E)));
219 if (p.is_zero()) {p.ord=1; return 1; }
220 multiples.push_back(p);
221 Point q = p;
222 bigint eight; eight=8;
223 while ( (!q.is_zero()) && (q.getZ()<=eight)
224 && (multiples.size()<13) ) // 12 is max poss order
225 {
226 q+=p;
227 if (!q.is_zero()) multiples.push_back(q);
228 }
229 if (q.is_zero())
230 p.ord = multiples.size();
231 else
232 p.ord = -1;
233 #ifdef DEBUG_TORSION
234 cout<<"Returning ord="<<p.ord<<", multiples = "<<multiples<<endl;
235 #endif
236 return p.ord;
237 }
238
isvalid() const239 int Point::isvalid() const // P on its curve ?
240 {
241 if(E == 0){
242 cerr << "## Bad point: null curve pointer!"<<endl;
243 return 0 ;
244 }
245 // Null point, useful for terminating input:
246 if((sign(X)==0)&&(sign(Y)==0)&&(sign(Z)==0)) return 0;
247 // Point at infinity is on a curve
248 if((sign(X)==0)&&(sign(Z)==0)) return 1 ;
249 else{
250 // should calculate
251 // y^2 +a1 x y + a3 y
252 // and
253 // x^3 + a2 x^2 + a4 x + a6
254 // where
255 // x = X/Z, y = Y/Z
256 // and verify equality.
257 //
258 // In homogeneous coordinates:
259 //
260 // Lhs = Y^2 Z + a1 XYZ + a3 YZ^2 = (YZ) *(Y + a1 X + a3 Z)
261 //
262 //
263 // Rhs = X^3 +a2 X^2 Z + a4 X Z^2 + a6 Z^3
264 //
265 bigint A1,A2,A3,A4,A6; E->getai(A1,A2,A3,A4,A6);
266 const bigint& Lhs = Y*Z*(Y + A1*X + A3*Z) ;
267 const bigint& Rhs = A6*pow(Z,3) + X*(A4*Z*Z + X*(A2*Z + X)) ;
268 return Lhs == Rhs ;
269 }
270 }
271
272 // Find all points with a given rational x-coordinate:
points_from_x(Curvedata & E,const bigrational & x)273 vector<Point> points_from_x(Curvedata &E, const bigrational& x)
274 {
275 // cout<<"Trying to construct points with x="<<x<<endl;
276 bigint a1,a2,a3,a4,a6,b2,b4,b6,b8;
277 E.getai(a1,a2,a3,a4,a6);
278 E.getbi(b2,b4,b6,b8);
279 vector<Point> ans;
280 bigint xn = num(x), xd2=den(x), xd, xd4, s, t, yn;
281 // cout<<"xd2 = "<<xd2<<endl;
282 if(isqrt(xd2,xd)) // xd2=xd^2
283 {
284 // cout<<"xd = "<<xd<<endl;
285 xd4 = xd2*xd2;
286 s = ((4*xn+b2*xd2)*xn+(2*b4*xd4))*xn+b6*xd4*xd2;
287 // cout<<"s = "<<s<<endl;
288 if(isqrt(s,t)) // s=t^2
289 {
290 // cout<<"t = "<<t<<endl;
291 yn = t - (a1*xn+a3*xd2)*xd;
292 divide_exact(yn,BIGINT(2),yn);
293 // cout<<"yn = "<<yn<<endl;
294 Point P(E,xn*xd,yn,xd2*xd);
295 // cout<<"point="<<P<<endl;
296 ans.push_back(P);
297 if (!is_zero(t)) ans.push_back(-P);
298 }
299 }
300 return ans;
301 }
302
get_ntorsion()303 long Curvedata::get_ntorsion()
304 {
305 if(ntorsion==0)
306 {
307 #ifdef DEBUG_TORSION
308 cout<<"Calling torsion_points() on "<<(Curve)(*this)<<endl;
309 #endif
310 vector<Point> ans = torsion_points(*this);
311 ntorsion = ans.size();
312 #ifdef DEBUG_TORSION
313 cout<<"torsion_points() returns "<<ans<<" of size "<<ntorsion<<endl;
314 #endif
315 }
316 #ifdef DEBUG_TORSION
317 cout<<"Curvedata::get_ntorsion() returns "<<ntorsion<<endl;
318 #endif
319 return ntorsion;
320 }
321
322 // Cperiods is a class containing a basis for the period lattice L;
323 // it knows how to compute points from z mod L; so this function
324 // effectively does the same as PARI's pointell() (called ellzto point
325 // in the next PARI release)
326 //
make_tor_pt(Curvedata & E,Cperiods & per,const bigfloat & ra1,const bigfloat & ra2,const bigfloat & ra3,const bigcomplex & z)327 Point make_tor_pt(Curvedata& E, Cperiods& per,
328 const bigfloat& ra1, const bigfloat& ra2, const bigfloat& ra3,
329 const bigcomplex& z)
330 {
331 bigcomplex cx,cy;
332 #ifdef DEBUG_TORSION
333 cout<<"In make_tor_pt() with z = " << z << endl;
334 #endif
335 per.XY_coords(cx,cy,z);
336 cx = cx-(ra1*ra1+4*ra2)/to_bigfloat(12);
337 cy = (cy - ra1*cx - ra3)/to_bigfloat(2);
338 #ifdef DEBUG_TORSION
339 cout<<"(x,y) = ("<<(cx)<<","<<(cy)<<")\n";
340 #endif
341 bigint x=Iround(real(cx)), y=Iround(real(cy));
342 Point P(E, x, y);
343 return P;
344 }
345
346
two_torsion(Curvedata & E)347 vector<Point> two_torsion(Curvedata& E)
348 {
349 #ifdef DEBUG_TORSION
350 cout<<"\nIn two_torsion() with curve "<<(Curve)E<<"\n";
351 #endif
352 bigint a1, a2, a3, a4, a6, b2, b4, b6, b8;
353 E.getai(a1,a2,a3,a4,a6);
354 E.getbi(b2,b4,b6,b8);
355 int scaled=0;
356 if (odd(a1) || odd(a3))
357 {
358 b4*=BIGINT(8);
359 b6*=BIGINT(16);
360 scaled=1;
361 }
362 else
363 {
364 b2=a2; b4=a4; b6=a6;
365 }
366 vector<bigint> xlist = Introotscubic(b2,b4,b6);
367 int n, n2tors = xlist.size();
368
369 // If there are 3 points of order 2, we order them for consistency:
370 if(n2tors==3) sort(xlist.begin(),xlist.end());
371
372 vector<Point> two_tors;
373 two_tors.push_back(Point(E)) ; // zero point
374 for(n=0; n<n2tors; n++)
375 {
376 bigint ei = xlist[n];
377 if(scaled)
378 two_tors.push_back(Point(E,2*ei,-a1*ei-4*a3,BIGINT(8)));
379 else
380 two_tors.push_back(Point(E,ei,BIGINT(0),BIGINT(1)));
381 }
382 #ifdef DEBUG_TORSION
383 cout<<"\ntwo_torsion() returning "<<two_tors<<"\n";
384 #endif
385 return two_tors;
386 }
387
388 // Returns vector of x values such that x/3 is the rational x-coord of
389 // a point in E[3], possibly with y not rational. If y is rational
390 // then these x will in fact be multiples of 3 since rational
391 // 3-torsion is integral
three_torsion_x(Curvedata & E)392 vector<bigint> three_torsion_x(Curvedata& E)
393 {
394 #ifdef DEBUG_TORSION
395 cout<<"\nIn three_torsion_x() with curve "<<(Curve)E<<"\n";
396 #endif
397 bigint b2, b4, b6, b8;
398 E.getbi(b2,b4,b6,b8);
399 vector<bigint> xlist = Introotsquartic(b2,9*b4,27*b6,27*b8);
400 // NB The implementation of Introosquartic() in marith.cc does not
401 // fix the order of the roots, which depends on the order of the
402 // factors in NTL's Z[X] factorization routine. HENCE the order of
403 // the 3-torsion point x-coordinates (when there are two) is not
404 // well-defined without the sorting done here.
405 #ifdef DEBUG_TORSION
406 cout<<"\nthree_torsion_x() finds unsorted xlist = "<<xlist<<"\n";
407 #endif
408 if(xlist.size()==2)
409 {
410 sort(xlist.begin(),xlist.end());
411 #ifdef DEBUG_TORSION
412 cout<<"\nthree_torsion_x() returns sorted xlist = "<<xlist<<"\n";
413 #endif
414 }
415 return xlist;
416 }
417
three_torsion(Curvedata & E)418 vector<Point> three_torsion(Curvedata& E)
419 {
420 #ifdef DEBUG_TORSION
421 cout<<"\nIn three_torsion() with curve "<<(Curve)E<<"\n";
422 #endif
423 bigint a1, a2, a3, a4, a6, b2, b4, b6, b8, xi, d, rd;
424 E.getai(a1,a2,a3,a4,a6);
425 E.getbi(b2,b4,b6,b8);
426 vector<bigint> xlist = three_torsion_x(E);
427 vector<Point> three_tors;
428 three_tors.push_back(Point(E)) ; // zero point
429 for(unsigned int n=0; n<xlist.size(); n++)
430 {
431 xi = xlist[n];
432 if(xi%3==0) // 3-torsion must be integral
433 {
434 xi/=3;
435 d = ((4*xi+b2)*xi+(2*b4))*xi+b6;
436 if(isqrt(d,rd))
437 {
438 Point P(E,2*xi,rd-(a1*xi+a3),BIGINT(2));
439 three_tors.push_back(P);
440 three_tors.push_back(-P);
441 }
442 }
443 }
444 #ifdef DEBUG_TORSION
445 cout<<"\nthree_torsion() returning "<<three_tors<<"\n";
446 #endif
447 return three_tors;
448 }
449
450 // New torsion point routine using Mazur's theorem to limit possibilities
451 // and computing possible real torsion from the period lattice. Suggestion
452 // of Darrin Doud.
453 //
torsion_points(Curvedata & E)454 vector<Point> torsion_points(Curvedata& E) // After Darrin Doud, adapted by JC
455 {
456 if ( E.isnull() ) return vector<Point>(0);
457 //
458 // table[i][] contains a list of possible maximal orders for a point,
459 // given that the 2-torsion subgroup has order i
460 //
461 static long table[5][5] = {{},{5,7,9,3},{12,6,8,4,10},{},{8,6,4}};
462 static long nt[5] = {0,4,5,0,3};
463 #ifdef DEBUG_TORSION
464 cout<<"\nIn torsion_points() with curve "<<(Curve)E<<"\n";
465 #endif
466 bigint a1,a2,a3,a4,a6,sa2,sa4,sa6, x, y ;
467 E.getai(a1,a2,a3,a4,a6);
468 bigfloat ra1=I2bigfloat(a1), ra2=I2bigfloat(a2), ra3=I2bigfloat(a3);
469 long i,j,ntp=1,nt2;
470 vector<Point> points;
471 vector<Point> cycle;
472 Cperiods per(E); bigcomplex w1,w2,z,z2;
473 per.getwRI(w1,w2); z2=w2/to_bigfloat(2);
474 int ncc = getconncomp(E);
475 int found;
476 #ifdef DEBUG_TORSION
477 cout<<"Periods: "<<per<<"\nReal Period = "<<w1<<endl;
478 cout<<ncc<<" real component(s)"<<endl;
479 #endif
480 points.push_back(Point(E)) ; // zero point
481 Point p, q;
482
483 // We find the two-torsion algebraically
484
485 vector<Point> two_tors = two_torsion(E);
486 nt2=two_tors.size();
487 #ifdef DEBUG_TORSION
488 cout<<"Size of 2-torsion subgroup = " << nt2 << endl;
489 cout << two_tors << endl;
490 #endif
491
492 for(i=0, found=0; (i<nt[nt2])&&(!found); i++)
493 {
494 long ni=table[nt2][i];
495 #ifdef DEBUG_TORSION
496 cout<<"Looking for a point of order " << ni << "\n";
497 #endif
498 if(ni==3)
499 {
500 p=Point(E);
501 vector<Point> p3 = three_torsion(E);
502 if(p3.size()>1) {p=p3[1];}
503 }
504 else
505 {
506 z=w1/to_bigfloat(ni);
507 p = make_tor_pt(E,per,ra1,ra2,ra3,z);
508 }
509 #ifdef DEBUG_TORSION
510 cout<<"p = " << p <<"?\n";
511 #endif
512 found=(p.isvalid())&&(order(p,cycle)==ni);
513 if(!found&&(ncc==2)&&even(ni))
514 {
515 p = make_tor_pt(E,per,ra1,ra2,ra3,z+z2);
516 #ifdef DEBUG_TORSION
517 cout<<"p = " << p <<"?\n";
518 #endif
519 found=(p.isvalid())&&(order(p,cycle)==ni);
520 if(!found&&((ni%4)==2))
521 {
522 p = make_tor_pt(E,per,ra1,ra2,ra3,z+z+z2);
523 #ifdef DEBUG_TORSION
524 cout<<"p = " << p <<"?\n";
525 #endif
526 found=(p.isvalid())&&(order(p,cycle)==ni);
527 }
528 }
529 if(found)
530 {
531 ntp=ni;
532 points=cycle;
533 #ifdef DEBUG_TORSION
534 cout<<"Found a point " << p << " of order " << ni << endl;
535 cout<<"generating subgroup "<<cycle<<endl;
536 #endif
537 }
538 #ifdef DEBUG_TORSION
539 else
540 {
541 cout<<"none\n";
542 }
543 #endif
544 }
545
546 #ifdef DEBUG_TORSION
547 cout<<"Number of points in cyclic part = " << ntp << endl;
548 cout<<points<<endl;
549 #endif
550
551 if(ntp==1) {return two_tors;} // C1, C2 or C2xC2
552 if(nt2==4)
553 // non-cyclic case, C2xC4, C2xC6 or C2xC8
554 // Find a point of order 2 not in the second factor and add it in
555 {
556 for(i=1; i<4; i++) // 0 is first point in two_tors !
557 {
558 p=two_tors[i];
559 if(find(points.begin(),points.end(),p)==points.end())
560 {
561 #ifdef DEBUG_TORSION
562 cout<<"Using 2-torsion point "<<p<<" as coset rep\n";
563 #endif
564 for(j=0; j<ntp; j++)
565 {
566 points.push_back(points[j]+p);
567 }
568 ntp*=2;
569 break;
570 }
571 #ifdef DEBUG_TORSION
572 cout<<"2-torsion point "<<p<<" is in subgroup\n";
573 #endif
574 }
575 }
576 return points;
577 }
578
579 //#define DEBUG_DIVISION_POINTS 1
division_points(int m)580 vector<Point> Point::division_points(int m) // list of Q s.t. n*Q=this
581 {
582 vector<Point> ans;
583 vector<Point> Qs;
584 Point Q;
585 if (is_torsion())
586 {
587 Qs = torsion_points(*E);
588 for (vector<Point>::iterator Qi = Qs.begin(); Qi!=Qs.end(); Qi++)
589 {
590 Q = *Qi;
591 if (m * Q == *this)
592 ans.push_back(Q);
593 }
594 return ans;
595 }
596 ZPoly pol = division_points_X_pol(E->a1,E->a2,E->a3,E->a4,E->a6, m, X, Z);
597 #ifdef DEBUG_DIVISION_POINTS
598 cout << "division poly = " << pol << endl;
599 #endif
600
601 vector<bigrational> xQs = roots(pol);
602 #ifdef DEBUG_DIVISION_POINTS
603 cout << " with roots " << xQs << endl;
604 #endif
605 for(vector<bigrational>::iterator xQi = xQs.begin(); xQi!=xQs.end(); xQi++)
606 {
607 Qs = points_from_x(*E, *xQi);
608 // will have length 0 or 2 since non-torsion, and we only want
609 // exctly one when there two so, must check which works
610 if (Qs.size()>0)
611 {
612 Q = Qs[0];
613 if (m*Q == *this)
614 {
615 ans.push_back(Q);
616 }
617 else
618 {
619 ans.push_back(-Q);
620 }
621 }
622 }
623 return ans;
624 }
625
626 // end of file: points.cc
627