1 // curvedata.cc -- implementation of Curvedata class
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
25 #include <eclib/curve.h>
26 #include <eclib/cubic.h>
27
Curvedata(const bigint & aa1,const bigint & aa2,const bigint & aa3,const bigint & aa4,const bigint & aa6,int min_on_init)28 Curvedata::Curvedata(const bigint& aa1, const bigint& aa2,
29 const bigint& aa3, const bigint& aa4,
30 const bigint& aa6, int min_on_init)
31 : Curve(aa1,aa2,aa3,aa4,aa6), minimal_flag(0), ntorsion(0)
32 {
33 b2 = a1*a1 + 4*a2; b4 = 2*a4 + a1*a3;
34 b6 = a3*a3 + 4*a6; b8 = (b2*b6 - b4*b4) / 4;
35 c4 = b2*b2 - 24*b4; c6 = -b2*b2*b2 + 36*b2*b4 - 216*b6;
36 discr = (c4*c4*c4 - c6*c6) / 1728;
37 discr_factored=0;
38 if(sign(discr)==0) // singular curve, replace by null
39 {
40 a1=0;a2=0;a3=0;a4=0;a6=0;
41 b2=0;b4=0;b6=0;b8=0;
42 c4=0;c6=0;
43 conncomp=0;
44 }
45 else
46 {
47 conncomp = sign(discr)>0 ? 2 : 1;
48 if (min_on_init) minimalize(); // which sets discr
49 }
50 }
51
52 //#define DEBUG_Q_INPUT
Curvedata(const vector<bigrational> & qai,bigint & scale)53 Curvedata::Curvedata(const vector<bigrational>& qai, bigint& scale)
54 : minimal_flag(0), ntorsion(0)
55 {
56 bigrational qa1(qai[0]), qa2(qai[1]), qa3(qai[2]), qa4(qai[3]), qa6(qai[4]);
57 scale=BIGINT(1);
58 a1=num(qa1); a2=num(qa2); a3=num(qa3); a4=num(qa4); a6=num(qa6);
59 #ifdef DEBUG_Q_INPUT
60 cout<<"In Curvedata constructor with ["<<qa1<<","<<qa2<<","<<qa3<<","<<qa4<<","<<qa6<<"]"<<endl;
61 #endif
62 vector<bigint> plist=pdivs(den(qa1));
63 plist=vector_union(plist,pdivs(den(qa2)));
64 plist=vector_union(plist,pdivs(den(qa3)));
65 plist=vector_union(plist,pdivs(den(qa4)));
66 plist=vector_union(plist,pdivs(den(qa6)));
67 #ifdef DEBUG_Q_INPUT
68 cout<<"Denominator primes: "<<plist<<endl;
69 #endif
70 for(vector<bigint>::const_iterator pp=plist.begin(); pp!=plist.end(); pp++)
71 {
72 bigint p=*pp;
73 #ifdef DEBUG_Q_INPUT
74 cout<<"p = "<<p<<endl;
75 #endif
76 long e=val(p,den(qa1));
77 e=max(e,ceil(rational(val(p,den(qa2)),2)));
78 e=max(e,ceil(rational(val(p,den(qa3)),3)));
79 e=max(e,ceil(rational(val(p,den(qa4)),4)));
80 e=max(e,ceil(rational(val(p,den(qa6)),6)));
81 #ifdef DEBUG_Q_INPUT
82 cout<<"e = "<<e<<endl;
83 #endif
84 if(e>0)
85 {
86 bigint pe=pow(p,e);
87 scale *= pe;
88 bigint pei=pe;
89 a1*=pei; pei*=pe;
90 a2*=pei; pei*=pe;
91 a3*=pei; pei*=pe;
92 a4*=pei; pei*=pe; pei*=pe;
93 a6*=pei;
94 }
95 }
96 a1/=den(qa1); a2/=den(qa2); a3/=den(qa3); a4/=den(qa4); a6/=den(qa6);
97 #ifdef DEBUG_Q_INPUT
98 cout<<"After scaling, coeffs are ["<<a1<<","<<a2<<","<<a3<<","<<a4<<","<<a6<<"]"<<endl;
99 #endif
100 b2 = a1*a1 + 4*a2; b4 = 2*a4 + a1*a3;
101 b6 = a3*a3 + 4*a6; b8 = (b2*b6 - b4*b4) / 4;
102 c4 = b2*b2 - 24*b4; c6 = -b2*b2*b2 + 36*b2*b4 - 216*b6;
103 discr = (c4*c4*c4 - c6*c6) / 1728;
104 discr_factored=0;
105 if(sign(discr)==0) // singular curve, replace by null
106 {
107 a1=0;a2=0;a3=0;a4=0;a6=0;
108 b2=0;b4=0;b6=0;b8=0;
109 c4=0;c6=0;
110 conncomp=0;
111 }
112 else
113 {
114 conncomp = sign(discr)>0 ? 2 : 1;
115 }
116 }
117
Curvedata(const Curve & c,int min_on_init)118 Curvedata::Curvedata(const Curve& c, int min_on_init)
119 : Curve(c), minimal_flag(0), ntorsion(0)
120 {
121 b2 = a1*a1 + 4*a2; b4 = 2*a4 + a1*a3;
122 b6 = a3*a3 + 4*a6; b8 = (b2*b6 - b4*b4) / 4;
123 c4 = b2*b2 - 24*b4; c6 = -b2*b2*b2 + 36*b2*b4 - 216*b6;
124 discr = (c4*c4*c4 - c6*c6) / 1728;
125 discr_factored=0;
126 if(sign(discr)==0) // singular curve, replace by null
127 {
128 a1=0;a2=0;a3=0;a4=0;a6=0;
129 b2=0;b4=0;b6=0;b8=0;
130 c4=0;c6=0;
131 conncomp=0;
132 }
133 else
134 {
135 conncomp = sign(discr)>0 ? 2 : 1;
136 if (min_on_init) minimalize(); // which sets discr
137 }
138 }
139
Curvedata(const Curvedata & c)140 Curvedata::Curvedata(const Curvedata& c)
141 : Curve(c), b2(c.b2), b4(c.b4), b6(c.b6), b8(c.b8), c4(c.c4),
142 c6(c.c6), discr(c.discr), minimal_flag(c.minimal_flag),
143 discr_factored(c.discr_factored), conncomp(c.conncomp),
144 ntorsion(c.ntorsion)
145 {
146 if(discr_factored) the_bad_primes=c.the_bad_primes;
147 }
148
Curvedata(const Curvedata & c,int min_on_init)149 Curvedata::Curvedata(const Curvedata& c, int min_on_init)
150 : Curve(c), b2(c.b2), b4(c.b4), b6(c.b6), b8(c.b8), c4(c.c4),
151 c6(c.c6), discr(c.discr), minimal_flag(c.minimal_flag),
152 discr_factored(c.discr_factored), conncomp(c.conncomp),
153 ntorsion(c.ntorsion)
154 {
155 if(discr_factored) the_bad_primes=c.the_bad_primes;
156 if (min_on_init) minimalize();
157 }
158
Curvedata(const bigint & cc4,const bigint & cc6,int min_on_init)159 Curvedata::Curvedata(const bigint& cc4, const bigint& cc6, int min_on_init)
160 :minimal_flag(0), discr_factored(0), ntorsion(0)
161 {
162 if (valid_invariants(cc4, cc6))
163 {
164 c4=cc4; c6=cc6;
165 c4c6_to_ai(cc4, cc6, a1, a2, a3, a4, a6, b2, b4, b6, b8);
166 /*
167 cout<<"a1="<<a1<<"\t";
168 cout<<"a2="<<a2<<"\t";
169 cout<<"a3="<<a3<<"\t";
170 cout<<"a4="<<a4<<"\t";
171 cout<<"a6="<<a6<<"\n";
172 */
173 if (min_on_init) minimalize();
174 else
175 {
176 discr = (c4*c4*c4 - c6*c6) / 1728;
177 }
178 conncomp = sign(discr)>0 ? 2 : 1;
179 }
180 else
181 {
182 cout << " ## attempt to call Curve constructor\n"
183 << " with invalid invariants c4 = "<<cc4<<", cc6 = "<<c6
184 << ": reading as null curve\n";
185 a1=0; a2=0; a3=0; a4=0; a6=0;
186 b2=0; b4=0; b6=0; b8=0;
187 c4=0; c6=0; discr=0;
188 }
189 }
190
operator =(const Curvedata & c)191 void Curvedata::operator=(const Curvedata& c)
192 {
193 a1 = c.a1, a2 = c.a2, a3 = c.a3, a4 = c.a4, a6 = c.a6;
194 b2 = c.b2, b4 = c.b4, b6 = c.b6, b8 = c.b8;
195 c4 = c.c4, c6 = c.c6, discr = c.discr;
196 minimal_flag = c.minimal_flag;
197 discr_factored = c.discr_factored;
198 if(discr_factored) the_bad_primes = c.the_bad_primes;
199 conncomp = c.conncomp;
200 ntorsion = c.ntorsion;
201 }
202
203 // Minimalizing, reducing, and standardizing a general curve
204 // based on Laska--Kraus--Connell algorithm
205
minimalize()206 void Curvedata::minimalize()
207 {
208 if (minimal_flag) return;
209 if ( isnull() ) {minimal_flag = 1; return; }
210
211 // else we are ready for Laska-Kraus-Connell reduction
212
213 bigint newc4, newc6, newdiscr, u;
214 //cout<<"minimising c4, c6 = "<<c4<<", "<<c6<<"\n";
215 minimise_c4c6(c4,c6,discr,newc4,newc6,newdiscr,u);
216 //cout<<"minimal c4, c6 = "<<newc4<<", "<<newc6<<"\t"<<"( u = "<<u<<")\n";
217
218 // now compute minimal equation
219 if ( u > 1) { c4 = newc4; c6 = newc6; }
220 discr = newdiscr;
221
222 if(discr_factored)
223 {
224 if(u>1) // trim list of bad primes
225 {
226 bigint p;
227 vector<bigint> new_bad_primes;
228 vector<bigint>::iterator pr = the_bad_primes.begin();
229 while(pr!=the_bad_primes.end())
230 {
231 bigint p=*pr++;
232 if(div(p,discr)) new_bad_primes.push_back(p);
233 }
234 the_bad_primes=new_bad_primes;
235 }
236 }
237 else
238 the_bad_primes = pdivs(discr);
239
240 // cout<<"After Curvedata::minimalize(): discr = "<<discr<<", ";
241 // cout<<"with bad primes "<<the_bad_primes<<endl;
242 c4c6_to_ai(c4,c6,a1,a2,a3,a4,a6,b2,b4,b6,b8);
243 minimal_flag = 1;
244 }
245
minimalize(bigint & u,bigint & r,bigint & s,bigint & t) const246 Curvedata Curvedata::minimalize(bigint& u, bigint& r, bigint& s, bigint& t) const
247 {
248 if (minimal_flag)
249 {
250 Curvedata newc(*this);
251 r=0; s=0; t=0; u=1;
252 return newc;
253 }
254 // else we use Laska-Kraus-Connell reduction
255
256 bigint newc4, newc6, newdiscr, u2;
257 //cout<<"minimising c4, c6 = "<<c4<<", "<<c6<<"\n";
258 minimise_c4c6(c4,c6,discr,newc4,newc6,newdiscr,u);
259 //cout<<"minimal c4, c6 = "<<newc4<<", "<<newc6<<"\t"<<"( u = "<<u<<")\n";
260
261 Curvedata newc(newc4, newc6, 0); // no need to re-minimalise
262 //cout<<"minimal curve = "<<newc<<endl;
263
264 s = (u*newc.a1 - a1)/2; u2=u*u;
265 r = (u2*newc.a2 - a2 + s*a1 + s*s)/3;
266 t = (u2*u*newc.a3 - a3 - r*a1)/2;
267 // cout<<"r,s,t="<<r<<","<<s<<","<<t<<endl;
268 return newc;
269 }
270
transform(const bigint & r,const bigint & s,const bigint & t)271 void Curvedata::transform(const bigint& r, const bigint& s, const bigint& t) //NB u=1;
272 {
273 a6 += r*(a4 + r*(a2 + r)) - t*(a3 + r*a1 + t);
274 a4 += -s*a3 + 2*r*a2 - (t + r*s)*a1 + 3*r*r - 2*s*t;
275 a3 += r*a1 +t+t;
276 a2 += -s*a1 + 3*r - s*s;
277 a1 += s+s;
278 b2 = a1*a1 + 4*a2;
279 b4 = a4+a4 + a1*a3;
280 b6 = a3*a3 + 4*a6;
281 b8 = (b2*b6 - b4*b4) / 4;
282 }
283
input(istream & is)284 void Curvedata::input(istream& is)
285 {
286 Curve::input(is);
287 b2=a1*a1 + 4*a2; b4=2*a4 + a1*a3;
288 b6=a3*a3 + 4*a6; b8= (b2*b6 - b4*b4) / 4;
289 c4=b2*b2 - 24*b4; c6=-b2*b2*b2 + 36*b2*b4 - 216*b6;
290 discr= (c4*c4*c4 - c6*c6) / 1728;
291 minimal_flag=0;
292 discr_factored=0;
293 conncomp= sign(discr)>0 ? 2 : 1;
294 ntorsion=0;
295 }
296
output(ostream & os) const297 void Curvedata::output(ostream& os) const
298 {
299 Curve::output(os);
300 if (isnull()) {os << " --singular\n"; return; }
301 if (minimal_flag) os << " (reduced minimal model)";
302 os << endl;
303 os << "b2 = " << b2 << "\t "
304 << "b4 = " << b4 << "\t "
305 << "b6 = " << b6 << "\t "
306 << "b8 = " << b8 << endl;
307 os << "c4 = " << c4 << "\t\t"
308 << "c6 = " << c6 << endl;
309 os << "disc = " << discr << "\t(";
310 if (minimal_flag&&discr_factored)
311 os << "bad primes: " << the_bad_primes << ";\t";
312 os << "# real components = " << conncomp << ")" << endl;
313 if (ntorsion) os<<"#torsion = "<<ntorsion<<endl;
314 else os<<"#torsion not yet computed"<<endl;
315 }
316
317
opt_x_shift(const Curvedata & C,bigint & k)318 Curvedata opt_x_shift(const Curvedata& C, bigint& k)
319 {
320 bigint b2,b4,b6,b8,four,zero; four=BIGINT(4); zero=BIGINT(0);
321 C.getbi(b2,b4,b6,b8);
322 cubic b_cubic(four,b2,2*b4,b6);
323 k = b_cubic.shift_reduce();
324 Curvedata CC(C);
325 CC.transform(k,zero,zero);
326 return CC;
327 }
328
329 #if(0)
330
331 // next the invariants package
332
agm(bigfloat a,bigfloat b)333 bigfloat agm( bigfloat a, bigfloat b)
334 {
335 bigfloat a0, b0 ;
336 do{
337 a0 = (a + b) / 2.0 ;
338 b0 = sqrt( a * b ) ;
339 a = a0 ; b = b0 ;
340 } while ( !is_approx_zero(a-b) ) ;
341 return ( a ) ;
342 }
343
cuberoot(bigfloat x)344 bigfloat cuberoot(bigfloat x ) /* computes x^(1.0/3.0) for all x */
345 {
346 bigfloat third(1); third/=3;
347 if (x >= 0.0)
348 return (pow(x, third) ) ;
349 else
350 return( - pow(-x, third) );
351 }
352
sort(bigfloat roots[3])353 void sort(bigfloat roots[3] )
354 {
355 bigfloat x0, x1, x2, temp ;
356 x0 = roots[0] ; x1 = roots[1] ; x2 = roots[2] ;
357 if (x0 > x1) { temp = x0; x0 = x1 ; x1 = temp ; } /* now x0 < x1 */
358 if (x1 > x2) { temp = x1; x1 = x2 ; x2 = temp ; } /* now x1 < x2, x0 < x2 */
359 if (x0 > x1) { temp = x0; x0 = x1 ; x1 = temp ; } /* now x0 < x1 < x2 */
360 roots[0] = x0 ; roots[1] = x1; roots[2] = x2 ;
361 }
362
CurvedataExtra(const Curvedata & c)363 CurvedataExtra::CurvedataExtra(const Curvedata& c)
364 : Curvedata(c)
365 {
366 if ( isnull() ) {roots[0]=roots[1]=roots[2]=period=0; return; }
367 bigfloat fb2 = I2bigfloat(b2);
368 bigfloat fc6 = I2bigfloat(c6);
369 bigfloat fdiscr = I2bigfloat(discr);
370 if (conncomp == 1) {
371 bigfloat sqdisc = 24.0 * sqrt(-3.0 * fdiscr) ;
372 roots[2] = ( -fb2 + cuberoot(fc6 - sqdisc) +
373 cuberoot(fc6 + sqdisc) ) / 12.0 ;
374 } else {
375 bigfloat theta = atan2( 24.0 * sqrt(3.0*fdiscr), fc6) ;
376 bigfloat fc4 = I2bigfloat(c4);
377 bigfloat sqc4 = sqrt( fc4 ) ;
378 roots[0] = (-fb2 + 2.0 * sqc4 * cos( theta/3.0 )) / 12.0 ;
379 roots[1] = (-fb2 + 2.0 * sqc4 * cos( (theta+2.0*Pi())/3.0)) / 12.0;
380 roots[2] = (-fb2 + 2.0 * sqc4 * cos( (theta+4.0*Pi())/3.0)) / 12.0;
381 sort(roots) ;
382 }
383 // note that roots[2] (the THIRD one!) is always the largest root
384 // now compute the period
385 bigfloat a, b ; // parameters to feed to agm()
386 if (conncomp == 2){
387 a = sqrt(roots[2] - roots[0]) ;
388 b = sqrt(roots[2] - roots[1]) ;
389 } else {
390 bigfloat fb4 = I2bigfloat(b4); // bigfloat version of b4
391 bigfloat fp = sqrt( fb4/2.0 + roots[2] * (fb2/2.0 + 3.0 * roots[2]) );
392 // f', where f is the cubic
393 a = 2.0 * sqrt( fp ) ;
394 b = sqrt( 3.0 * roots[2] + fb2/4.0 + 2.0 * fp ) ;
395 }
396 period = (2.0 * Pi() / agm(a, b) );
397 }
398
output(ostream & os) const399 void CurvedataExtra::output(ostream& os) const
400 {
401 Curvedata::output(os);
402 if(conncomp == 1)
403 os << "The real two-division point is " << roots[2] << endl ;
404 else
405 if(conncomp == 2)
406 os << "The real two-division points are "
407 << roots[0] << ", " << roots[1] << ", "<< roots[2] << endl ;
408 os << "The real period is " << period << endl;
409 }
410 #endif
411
412