1 //==============================================================================================
2 //
3 // This file is part of LiDIA --- a library for computational number theory
4 //
5 // Copyright (c) 1994--2001 the LiDIA Group. All rights reserved.
6 //
7 // See http://www.informatik.tu-darmstadt.de/TI/LiDIA/
8 //
9 //----------------------------------------------------------------------------------------------
10 //
11 // $Id$
12 //
13 // Author : Nigel Smart (NS)
14 // Adaption of John Cremona's code
15 // Changes : See CVS log
16 //
17 //==============================================================================================
18
19
20 #ifdef HAVE_CONFIG_H
21 # include "config.h"
22 #endif
23 #include "LiDIA/rational_factorization.h"
24 #include "LiDIA/elliptic_curves/ec_arith.h"
25
26
27
28 #ifdef LIDIA_NAMESPACE
29 namespace LiDIA {
30 #endif
31
32
33
convert(bigint & a,const bigrational & b)34 void convert(bigint& a, const bigrational& b)
35 {
36 if (!b.denominator().is_one()) {
37 lidia_error_handler("convert", "Cannot convert bigrational to bigint");
38 }
39 a.assign(b.numerator());
40 }
41
42
43
44 base_vector< bigint >
int_roots_cubic(const bigint & a,const bigint & b,const bigint & c)45 int_roots_cubic(const bigint& a,
46 const bigint& b,
47 const bigint& c)
48 {
49 bigcomplex za(a), zb(b), zc(c);
50 bigcomplex* croots = solve_cubic(za, zb, zc);
51 base_vector< bigint > iroots(3, 0);
52 int niroots = 0;
53 bigint x, cx;
54 for (int i = 0; i < 3; i++) {
55 croots[i].real().bigintify(x);
56 if (x.is_zero()) {
57 if (c.is_zero()) {
58 int newone = 1;
59 for (int j = 0; (j < niroots) && newone; j++)
60 newone = (x != iroots[j]);
61 if (newone)
62 iroots[niroots++] = x;
63 }
64 }
65 else {
66 cx = c/x;
67 if (x*cx == c)
68 if (((x+a)*x+b+cx).is_zero()) {
69 //x is a root, so add to list if not there already:
70 int newone = 1;
71 for (int j = 0; (j < niroots) && newone; j++)
72 newone = (x != iroots[j]);
73 if (newone)
74 iroots[niroots++] = x;
75 }
76 }
77 }
78 delete[] croots;
79 return iroots;
80 }
81
82
83
84 bigcomplex
discriminant(const bigcomplex & b,const bigcomplex & c,const bigcomplex & d)85 discriminant(const bigcomplex& b,
86 const bigcomplex& c,
87 const bigcomplex& d)
88 {
89 bigcomplex bb = b*b, cc = c*c, bc = b*c;
90 return bigfloat(27.0)*d*d-bc*bc+bigfloat(4.0)*bb*b*d-bigfloat(18.0)*bc*d+bigfloat(4.0)*c*cc;
91 }
92
93
94
95 bigcomplex*
solve_cubic(const bigcomplex & c1,const bigcomplex & c2,const bigcomplex & c3)96 solve_cubic(const bigcomplex& c1,
97 const bigcomplex& c2,
98 const bigcomplex& c3)
99 {
100 bigfloat three(3.0), two(2.0), one(1.0);
101 bigfloat third = one/three;
102 bigcomplex w = bigfloat(0.5)*bigcomplex(-1, sqrt(three));
103 bigcomplex disc = discriminant(c1, c2, c3);
104 bigcomplex p3 = three*c2 - c1*c1;
105 bigcomplex mc1over3 = -c1*third;
106 bigcomplex *roots = new bigcomplex[3];
107
108 if (abs(disc).is_approx_zero()) {
109 if (abs(p3).is_approx_zero()) {
110 // triple root
111 roots[0] = roots[1] = roots[2] = mc1over3;
112 }
113 else {
114 // double root
115 roots[0] = roots[1] = (c1*c2 - bigfloat(9.0)*c3)/(p3+p3);
116 roots[2] = -(roots[0] + roots[0]+c1);
117 }
118 }
119 else {
120 // distinct roots
121 bigcomplex q = (((mc1over3+c1)*mc1over3 +c2)*mc1over3 +c3);
122 // = F(mc1over3);
123 if (abs(p3).is_approx_zero()) {
124 // pure cubic
125 roots[0] = power(-q, third);
126 roots[1] = w*roots[0];
127 roots[2] = w*roots[1];
128 roots[0] += mc1over3;
129 roots[1] += mc1over3;
130 roots[2] += mc1over3;
131 }
132 else {
133 bigcomplex d = bigfloat(729.0)*q*q+ bigfloat(4.0)*p3*p3*p3;
134 bigcomplex t1cubed = bigfloat(0.5)*(sqrt(d)- bigfloat(27.0)*q);
135 bigcomplex t1 = power(t1cubed, third);
136 bigcomplex t2 = t1*w;
137 bigcomplex t3 = t2*w;
138 roots[0] = (-c1+t1-p3/t1)* third;
139 roots[1] = (-c1+t2-p3/t2)* third;
140 roots[2] = (-c1+t3-p3/t3)* third;
141 }
142 }
143 long i, iter, niter = 2; // number of iterations in Newton refinement of roots
144 for (i = 0; i < 3; i++) {
145 bigcomplex z = roots[i], fz, fdashz;
146 for (iter = 0; iter < niter; iter++) {
147 fz = ((z+c1)*z+c2)*z+c3;
148 fdashz = (3*z+2*c1)*z+c2;
149 if (!fdashz.is_approx_zero()) {
150 z -= fz/fdashz;
151 }
152 }
153 roots[i] = z;
154 }
155 return roots;
156 }
157
158
159
160 bigcomplex*
solve_quartic(const bigcomplex & a,const bigcomplex & b,const bigcomplex & c,const bigcomplex & d)161 solve_quartic(const bigcomplex& a,
162 const bigcomplex& b,
163 const bigcomplex& c,
164 const bigcomplex& d)
165 {
166 bigcomplex p, q, r, aa, e, f1, f2;
167 bigcomplex a4 = a/4;
168 bigcomplex* roots = new bigcomplex[4];
169 if (d.is_approx_zero()) {
170 roots[0] = 0;
171 bigcomplex* cuberoots = solve_cubic(a, b, c);
172 roots[1] = cuberoots[0];
173 roots[2] = cuberoots[1];
174 roots[3] = cuberoots[2];
175 delete[] cuberoots;
176 }
177 else {
178 p = b - 3*a*a / 8;
179 q = ((a / 2) * (a*a4 - b)) + c;
180 r = ((256*d) - (64*a*c) + (16*a*a*b) - (3*a*a*a*a)) / 256;
181 if (q.is_approx_zero()) {
182 bigcomplex s = sqrt(p*p - 4*r);
183 roots[0] = sqrt((-p + s) / 2) - a4;
184 roots[1] = -sqrt((-p + s) / 2) - a4;
185 roots[2] = sqrt((-p - s) / 2) - a4;
186 roots[3] = -sqrt((-p - s) / 2) - a4;
187 }
188 else {
189 bigcomplex* aaroots = solve_cubic(-p/2, -r, ((p*r)/2-(q*q)/8));
190 aa = aaroots[0];
191 if (aa.is_approx_zero()) aa = 0;
192 e = sqrt(-p + 2*aa);
193 f1 = (aa + q / (2*e));
194 f2 = (aa - q / (2*e));
195 bigcomplex s1 = sqrt(e*e - 4*f1);
196 bigcomplex s2 = sqrt(e*e - 4*f2);
197 roots[0] = ((e + s1) / 2) - a4;
198 roots[1] = ((e - s1) / 2) - a4;
199 roots[2] = ((-e + s2) / 2) - a4;
200 roots[3] = ((-e - s2) / 2) - a4;
201 delete[] aaroots;
202 }
203 }
204 return roots;
205 }
206
207
208
209 // puts in decreasing order
210 void
orderreal(bigfloat & e1,bigfloat & e2,bigfloat & e3)211 orderreal(bigfloat& e1, bigfloat& e2, bigfloat& e3)
212 {
213 if (e1 < e3)
214 swap(e1, e3);
215 if (e1 < e2)
216 swap(e1, e2);
217 else if (e2 < e3)
218 swap(e2, e3);
219 }
220
221
222
223 bigcomplex*
solve_real_quartic(const bigfloat & a,const bigfloat & b,const bigfloat & c,const bigfloat & d,const bigfloat & e)224 solve_real_quartic(const bigfloat& a,
225 const bigfloat& b,
226 const bigfloat& c,
227 const bigfloat& d,
228 const bigfloat& e)
229 {
230 bigfloat ii = 12*a*e - 3*b*d + c*c;
231 bigfloat jj = (72*a*e + 9*b*d - 2*c*c) * c - 27*(a*d*d + b*b*e);
232 bigfloat disc = 4*ii*ii*ii-jj*jj;
233 bigfloat H = 8*a*c - 3*b*b;
234 bigfloat R = b*b*b + 8*d*a*a - 4*a*b*c;
235 bigfloat Q = H*H-16*a*a*ii; // = 3*Q in JC's standard notation
236 int type, nrr;
237
238 if (disc < 0) {
239 type = 3;
240 nrr = 2;
241 } // 2 real roots
242 else {
243 if ((H< 0) && (Q > 0)) {
244 type = 2;
245 nrr = 4;
246 } // 4 real roots
247 else {
248 type = 1;
249 nrr = 0;
250 } // 0 real roots
251 }
252 bigcomplex c1(0), c2(-3*ii), c3(jj);
253 bigcomplex* cphi = solve_cubic(c1, c2, c3);
254 bigcomplex* roots = new bigcomplex[4];
255 bigfloat a4 = 4*a;
256 bigfloat one(1);
257 bigfloat oneover4a = one/a4;
258 bigfloat athird = one/bigfloat(3.0);
259
260 if (type < 3) {
261 // all the phi are real; order them so that a*phi[i] decreases
262 bigfloat phi1 = real(cphi[0]);
263 bigfloat phi2 = real(cphi[1]);
264 bigfloat phi3 = real(cphi[2]);
265
266 if (a > 0)
267 orderreal(phi1, phi2, phi3);
268 else
269 orderreal(phi3, phi2, phi1);
270
271 if (type == 2) {
272 // all roots are real
273 bigfloat r1 = sqrt((a4*phi1-H)*athird);
274 bigfloat r2 = sqrt((a4*phi2-H)*athird);
275 bigfloat r3 = sqrt((a4*phi3-H)*athird);
276 if (R < 0) r3 = -r3;
277 roots[0] = bigcomplex((r1 + r2 - r3 - b) * oneover4a);
278 roots[1] = bigcomplex((r1 - r2 + r3 - b) * oneover4a);
279 roots[2] = bigcomplex((-r1 + r2 + r3 - b) * oneover4a);
280 roots[3] = bigcomplex((-r1 - r2 - r3 - b) * oneover4a);
281 // Those are all real and in descending order of size
282 }
283 else {
284 // no roots are real
285 bigfloat r1 = sqrt((a4*phi1-H)/3);
286 bigfloat ir2 = sqrt(-((a4*phi2-H)/3));
287 bigfloat ir3 = sqrt(-((a4*phi3-H)/3));
288 if (R > 0) r1 = -r1;
289 roots[0] = bigcomplex(r1-b, ir2 - ir3) * oneover4a;
290 roots[1] = conj(roots[0]);
291 roots[2] = bigcomplex(-r1-b, ir2 + ir3) * oneover4a;
292 roots[3] = conj(roots[2]);
293 }
294 }
295 else {
296 // disc < 0
297 bigfloat realphi; // will hold the real root, which will be cphi[2]
298 if (cphi[1].is_approx_real()) {
299 realphi = real(cphi[1]);
300 cphi[1] = cphi[2];
301 cphi[2] = realphi;
302 }
303 else
304 if (cphi[2].is_approx_real()) {
305 realphi = real(cphi[2]);
306 }
307 else {
308 realphi = real(cphi[0]);
309 cphi[0] = cphi[2];
310 cphi[2] = realphi;
311 }
312 bigcomplex r1 = sqrt((a4*cphi[0]-H)/3);
313 bigfloat r3 = sqrt((a4*realphi-H)/3);
314 if (R < 0) r3 = -r3;
315 roots[0] = bigcomplex(r3 - b, 2*imag(r1)) * oneover4a;
316 roots[1] = conj(roots[0]);
317 roots[2] = bigcomplex((2*real(r1) - r3 - b)) * oneover4a;
318 roots[3] = bigcomplex((-2*real(r1) - r3 - b)) * oneover4a;
319 // roots[2] and roots[3] are real
320 }
321 delete[] cphi;
322 return roots;
323 }
324
325
326
327 base_vector< bigint >
square_divs(const bigint & N)328 square_divs(const bigint& N)
329 {
330 rational_factorization plist(N);
331 plist.factor();
332
333 int np = plist.no_of_comp();
334 int i, e = 0, nu = 1, nd = nu;
335
336 int* elist = new int[np];
337 bigint p;
338 for (i = 0; i < np; i++) {
339 p = plist.base(i);
340 elist[i] = e = valuation(p, N)/2;
341 nd *= (1+e);
342 }
343 base_vector< bigint > dlist(nd, 0);
344 dlist[0] = 1;
345 nd = nu;
346 for (i = 0; i < np; i++) {
347 p = plist.base(i);
348 e = elist[i];
349 for (int j = 0; j < e; j++)
350 for (int k = 0; k < nd; k++)
351 dlist[nd*(j+1)+k] = p*dlist[nd*j+k];
352 nd *= (e+1);
353 }
354 delete[] elist;
355 return dlist;
356 }
357
358
359
360 bigcomplex
cagm(const bigcomplex & a,const bigcomplex & b)361 cagm(const bigcomplex& a,
362 const bigcomplex& b)
363 {
364 bigcomplex x = a, y = b, z;
365 bigfloat theta, two(2), pi = Pi();
366 while (!(abs(x-y)).is_approx_zero()) {
367 z = (x+y)/two;
368 y = sqrt(x*y);
369 theta = two*arg(y/z);
370 if (theta > pi || theta <= -pi)
371 y = -y;
372 x = z;
373 }
374 return x;
375 }
376
377
378
379 bigcomplex
normalize(bigcomplex & w1,bigcomplex & w2)380 normalize(bigcomplex& w1, bigcomplex& w2)
381 {
382 bigcomplex tau = w1/w2, w3;
383 if (tau.imag() < 0) {
384 w1 = -w1;
385 tau = -tau;
386 }
387 w1 = w1-w2*round(tau.real());
388 tau = w1/w2;
389
390 for (int i = 1; i < 50 && (abs(tau) < 1); i++) {
391 // {Just to stop infinite loop due to rounding}
392 w3 = -w1;
393 w1 = w2;
394 w2 = w3;
395 tau = w1/w2;
396 w1 = w1-w2*round(tau.real());
397 tau = w1/w2;
398 }
399 return tau;
400 }
401
402
403
404 void
getc4c6(const bigcomplex & w1,const bigcomplex & w2,bigcomplex & c4,bigcomplex & c6)405 getc4c6(const bigcomplex& w1,
406 const bigcomplex& w2,
407 bigcomplex& c4,
408 bigcomplex &c6)
409 {
410 bigcomplex tau = w1/w2;
411 bigfloat zero(0), one(1), two(2);
412 bigfloat pi(Pi());
413 bigfloat x = two*pi*tau.real();
414 bigfloat y = two*pi*tau.imag();
415 bigcomplex q = exp(-y) * bigcomplex(cos(x), sin(x));
416 bigcomplex f = two*pi/w2;
417 bigcomplex f2 = f*f;
418 bigcomplex f4 = f2*f2;
419
420 bigcomplex term = bigcomplex(one);
421 bigcomplex qpower = bigcomplex(one);
422 bigcomplex sum4 = bigcomplex(zero);
423 bigcomplex sum6 = bigcomplex(zero);
424
425 bigfloat n, n2;
426
427 for (n = 1; !term.is_approx_zero(); n += 1) {
428 n2 = n*n;
429 qpower *= q;
430 term = n*n2*qpower/(one-qpower);
431 sum4 += term;
432 term *= n2;
433 sum6 += term;
434 }
435 c4 = (one + bigfloat(240.0)*sum4)*f4;
436 c6 = (one - bigfloat(504.0)*sum6)*f4*f2;
437 }
438
439
440
441 #ifdef LIDIA_NAMESPACE
442 } // end of namespace LiDIA
443 #endif
444