1 /*
2
3 nt.c - number theoretic helper functions
4
5 Copyright (C) 2009, 2010, 2015 Andreas Enge
6
7 This file is part of CM.
8
9 CM is free software; you can redistribute it and/or modify it under
10 the terms of the GNU General Public License as published by the
11 Free Software Foundation; either version 3 of the license, or (at your
12 option) any later version.
13
14 CM is distributed in the hope that it will be useful, but WITHOUT ANY
15 WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 for more details.
18
19 You should have received a copy of the GNU General Public License along
20 with CM; see the file COPYING. If not, write to the Free Software
21 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
22 */
23
24 #include "cm_common-impl.h"
25
26 static void nt_elliptic_curve_double (mpz_t P_x, mpz_t P_y, bool *P_infty,
27 mpz_t a, mpz_t p);
28 static void nt_elliptic_curve_add (mpz_t P_x, mpz_t P_y, bool *P_infty,
29 mpz_t Q_x, mpz_t Q_y, const bool Q_infty, mpz_t a, mpz_t p);
30
31 /*****************************************************************************/
32
cm_nt_mod(const long int a,const long int p)33 long int cm_nt_mod (const long int a, const long int p)
34 /* returns a representative of a % p in [0, p-1[ */
35
36 {
37 long int res = a % p;
38 return (res < 0 ? res + p : res);
39 }
40
41 /*****************************************************************************/
42
cm_nt_sqrt(const unsigned long int n)43 long int cm_nt_sqrt (const unsigned long int n)
44 /* returns the positive root of n if it is a perfect square, -1 otherwise */
45
46 {
47 double root;
48 unsigned long int s;
49
50 root = sqrt ((double) n);
51 s = (unsigned long int) (root + 0.5);
52 if (s * s == n)
53 return (long int) s;
54 else
55 return -1;
56 }
57
58 /*****************************************************************************/
59
cm_nt_gcd(long int a,long int b)60 long int cm_nt_gcd (long int a, long int b)
61 /* returns the positive gcd of a and b */
62
63 {
64 long int r, a_local, b_local;
65
66 if (a == 0)
67 return b;
68 else if (b == 0)
69 return a;
70 else
71 {
72 if (a < 0)
73 a_local = -a;
74 else
75 a_local = a;
76 if (b < 0)
77 b_local = -b;
78 else
79 b_local = b;
80 r = a_local % b_local;
81 while (r > 0)
82 {
83 a_local = b_local;
84 b_local = r;
85 r = a_local % b_local;
86 }
87 return b_local;
88 }
89 }
90
91 /*****************************************************************************/
92
cm_nt_gcdext(long int * u,long int * v,long int a,long int b)93 long int cm_nt_gcdext (long int *u, long int *v, long int a, long int b)
94 /* returns the positive gcd d of a and b and u, v such that u a + v b = d */
95
96 {
97 long int r0, r1, r2, u0, u1, u2, v0, v1, v2, q;
98 int sgn_a, sgn_b;
99
100 if (a < 0) {
101 sgn_a = -1;
102 r0 = -a;
103 }
104 else {
105 sgn_a = 1;
106 r0 = a;
107 }
108 if (b < 0) {
109 sgn_b = -1;
110 r1 = -b;
111 }
112 else {
113 sgn_b = 1;
114 r1 = b;
115 }
116 u0 = 1;
117 u1 = 0;
118 v0 = 0;
119 v1 = 1;
120
121 while (r1 != 0) {
122 q = r0 / r1;
123 r2 = r0 % r1;
124 r0 = r1;
125 r1 = r2;
126 u2 = u0 - q * u1;
127 u0 = u1;
128 u1 = u2;
129 v2 = v0 - q * v1;
130 v0 = v1;
131 v1 = v2;
132 }
133
134 *u = sgn_a * u0;
135 *v = sgn_b * v0;
136
137 return r0;
138 }
139
140 /*****************************************************************************/
141
cm_nt_invert(long int a,const long int p)142 long int cm_nt_invert (long int a, const long int p)
143 /* returns a^(-1) mod p */
144
145 {
146 long int r, a_local, b_local, u0, u1, u2;
147
148 if (a == 0)
149 {
150 printf ("*** nt_invert called with 0\n");
151 exit (1);
152 }
153 else if (a < 0)
154 /* we assume that a is reduced mod p by the C function */
155 a += p;
156
157 a_local = a;
158 b_local = p;
159 u0 = 1;
160 u1 = 0;
161 /* loop invariant: u0 * a + something * p = a_local */
162 /* loop invariant: u1 * a + something * p = b_local */
163 while (b_local > 0)
164 {
165 u2 = u0 - (a_local / b_local) * u1;
166 r = a_local % b_local;
167 u0 = u1;
168 u1 = u2;
169 a_local = b_local;
170 b_local = r;
171 }
172 return u0;
173 }
174
175 /*****************************************************************************/
176
cm_nt_kronecker(long int a,long int b)177 int cm_nt_kronecker (long int a, long int b)
178 /* computes the Kronecker symbol (a / b) following Algorithm 1.4.12 */
179 /* in Cohen (by the binary algorithm) */
180 {
181 const int tab [8] = {0, 1, 0, -1, 0, -1, 0, 1};
182 int k, r;
183
184 /* step 1 */
185 if (b == 0)
186 {
187 if (a == 1 || a == -1)
188 return 1;
189 else
190 return 0;
191 }
192
193 /* step 2*/
194 if (a % 2 == 0 && b % 2 == 0)
195 return 0;
196
197 while (b % 4 == 0)
198 b >>= 2;
199 if (b % 2 == 0)
200 {
201 b >>= 1;
202 k = tab [a & 7];
203 }
204 else
205 k = 1;
206
207 if (b < 0)
208 {
209 b = -b;
210 if (a < 0)
211 k = -k;
212 }
213
214 /* step 3 and added test; here, b is already odd */
215 if (a < 0)
216 {
217 a = -a;
218 if (b & 2)
219 k = -k;
220 }
221 a %= b;
222
223 /* steps 4 to 6 */
224 while (a != 0)
225 {
226 while (a % 4 == 0)
227 a >>= 2;
228 if (a % 2 == 0)
229 {
230 a >>= 1;
231 k *= tab [b & 7];
232 }
233 if (b > a)
234 {
235 r = b - a;
236 if (a & b & 2)
237 k = -k;
238 b = a;
239 a = r;
240 }
241 else
242 a -= b;
243 }
244
245 if (b > 1)
246 return 0;
247 else
248 return k;
249 }
250
251 /*****************************************************************************/
252
cm_nt_is_prime(mpz_t a)253 int cm_nt_is_prime (mpz_t a)
254
255 {
256 return (mpz_probab_prime_p (a, 10) > 0);
257 }
258
259 /*****************************************************************************/
260
cm_nt_is_prime_l(const unsigned long int prime)261 int cm_nt_is_prime_l (const unsigned long int prime)
262
263 {
264 int res;
265 mpz_t a;
266
267 mpz_init_set_ui (a, prime);
268 res = cm_nt_is_prime (a);
269
270 mpz_clear (a);
271
272 return res;
273 }
274
275 /*****************************************************************************/
276
cm_nt_next_prime(const unsigned long int n)277 unsigned long int cm_nt_next_prime (const unsigned long int n)
278 /* returns the prime following n */
279
280 {
281 const unsigned long int P [] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
282 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103,
283 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179,
284 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257,
285 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347,
286 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431,
287 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509,
288 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607,
289 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691,
290 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
291 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883,
292 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991,
293 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063,
294 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153,
295 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237,
296 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319,
297 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433,
298 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499,
299 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597,
300 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669,
301 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777,
302 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873,
303 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979,
304 1987, 1993, 1997, 1999 };
305 const unsigned int size = 303; /* number of entries in P */
306
307 if (n < P [size - 1]) {
308 /* search in prime list; loop invariants:
309 left <= right
310 P [0], ..., P [left - 1] <= n
311 P [right], ..., P [size - 1] > n */
312 unsigned int left = 0, right = size - 1, middle;
313 while (left != right) {
314 middle = (left + right) / 2; /* < right */
315 if (P [middle] > n)
316 right = middle; /* becomes smaller */
317 else
318 left = middle + 1; /* becomes bigger, remains <= right */
319 }
320 return P [left];
321 }
322 else {
323 unsigned long int res;
324 mpz_t a;
325 if (n & 1)
326 mpz_init_set_ui (a, n+2);
327 else
328 mpz_init_set_ui (a, n+1);
329 while (!cm_nt_is_prime (a))
330 mpz_add_ui (a, a, 2ul);
331 res = mpz_get_ui (a);
332 mpz_clear (a);
333 return res;
334 }
335 }
336
337 /*****************************************************************************/
338
cm_nt_factor(long int d,unsigned long int * factors,unsigned int * exponents)339 void cm_nt_factor (long int d, unsigned long int *factors,
340 unsigned int *exponents)
341 /* factors the absolute value of d by trial division. The prime factors */
342 /* are stored in "factors", their multiplicities in "exponents", which */
343 /* must provide sufficiently much space. The list of prime factors is */
344 /* terminated by 0, so that 12 entries suffice for a number of 32 bits, */
345 /* and 17 entries for a number of 64 bits. */
346
347 {
348 unsigned long int no, trial, trial2;
349 int i, j;
350
351 if (d < 0)
352 no = -d;
353 else
354 no = d;
355
356 i = 0;
357 j = 0;
358 trial = 0;
359 trial2 = 0;
360 while (trial2 <= no)
361 {
362 if (trial == 0)
363 {
364 trial = 2;
365 trial2 = 4;
366 }
367 else if (trial == 2)
368 {
369 trial = 3;
370 trial2 = 9;
371 }
372 else
373 {
374 trial += 2;
375 trial2 += 4 * (trial - 1);
376 }
377 if (no % trial == 0)
378 {
379 factors [j] = trial;
380 exponents [j] = 1;
381 no /= trial;
382 while (no % trial == 0)
383 {
384 no /= trial;
385 exponents [j]++;
386 }
387 j++;
388 }
389 i++;
390 }
391 if (no != 1)
392 {
393 factors [j] = no;
394 exponents [j] = 1;
395 j++;
396 }
397 factors [j] = 0;
398 }
399
400 /*****************************************************************************/
401
cm_nt_mpz_tonelli_z(mpz_t root,mpz_t a,mpz_t p)402 void cm_nt_mpz_tonelli_z (mpz_t root, mpz_t a, mpz_t p)
403 /* computes a square root of a modulo p by the Tonelli-Shanks algorithm, */
404 /* see Cohen93, Algorithm 1.5. */
405
406 {
407 mpz_t pm1, z;
408 mpz_t a_local, y, x, b, tmp;
409 unsigned int e;
410 unsigned long int r, m;
411
412 mpz_init (a_local);
413 mpz_tdiv_r (a_local, a, p);
414 if (mpz_cmp_ui (a_local, 0ul) == 0) {
415 mpz_set_ui (root, 0ul);
416 mpz_clear (a_local);
417 return;
418 }
419
420 mpz_init (pm1);
421 mpz_init (z);
422 mpz_init (y);
423 mpz_init (x);
424 mpz_init (b);
425 mpz_init (tmp);
426
427 mpz_sub_ui (pm1, p, 1ul);
428 e = 0;
429 while (mpz_divisible_2exp_p (pm1, 1ul) != 0) {
430 mpz_tdiv_q_2exp (pm1, pm1, 1ul);
431 e++;
432 }
433 if (e > 1) {
434 /* find generator of the 2-Sylow group */
435 for (mpz_set_ui (z, 2ul); mpz_legendre (z, p) != -1;
436 mpz_add_ui (z, z, 1ul));
437 mpz_powm (z, z, pm1, p);
438 }
439
440 if (e == 1) /* p=3 (mod 8) */ {
441 mpz_add_ui (tmp, p, 1ul);
442 mpz_tdiv_q_2exp (tmp, tmp, 2ul);
443 mpz_powm (x, a_local, tmp, p);
444 }
445 else {
446 /* initialisation */
447 mpz_set (y, z);
448 r = e;
449 mpz_sub_ui (tmp, pm1, 1ul);
450 mpz_tdiv_q_2exp (tmp, tmp, 1ul);
451 mpz_powm (x, a_local, tmp, p);
452 mpz_powm_ui (b, x, 2ul, p);
453 mpz_mul (b, b, a_local);
454 mpz_mod (b, b, p);
455 mpz_mul (x, x, a_local);
456 mpz_mod (x, x, p);
457 while (mpz_cmp_ui (b, 1ul) != 0) {
458 /* find exponent */
459 m = 1;
460 mpz_powm_ui (tmp, b, 2ul, p);
461 while (mpz_cmp_ui (tmp, 1ul) != 0) {
462 m++;
463 mpz_powm_ui (tmp, tmp, 2ul, p);
464 }
465 if (m == r) {
466 printf ("*** mpz_tonelli called with a = ");
467 mpz_out_str (stdout, 10, a);
468 printf (" and p = ");
469 mpz_out_str (stdout, 10, p);
470 printf (",\nbut a is not a square modulo p!\n");
471 exit (1);
472 }
473 r = 1 << (r - m - 1);
474 mpz_powm_ui (tmp, y, r, p);
475 mpz_powm_ui (y, tmp, 2ul, p);
476 r = m;
477 mpz_mul (x, x, tmp);
478 mpz_mod (x, x, p);
479 mpz_mul (b, b, y);
480 mpz_mod (b, b, p);
481 }
482 }
483
484 mpz_set (root, x);
485
486 mpz_clear (a_local);
487 mpz_clear (pm1);
488 mpz_clear (z);
489 mpz_clear (y);
490 mpz_clear (x);
491 mpz_clear (b);
492 mpz_clear (tmp);
493 }
494
495 /*****************************************************************************/
496
cm_nt_mpz_tonelli(mpz_t root,const long int a,mpz_t p)497 void cm_nt_mpz_tonelli (mpz_t root, const long int a, mpz_t p)
498 /* computes a square root of a modulo p */
499
500 {
501 mpz_t tmp_a;
502
503 mpz_init_set_si (tmp_a, a);
504 cm_nt_mpz_tonelli_z (root, tmp_a, p);
505 mpz_clear (tmp_a);
506 }
507
508 /*****************************************************************************/
509 /*****************************************************************************/
510
nt_elliptic_curve_double(mpz_t P_x,mpz_t P_y,bool * P_infty,mpz_t a,mpz_t p)511 static void nt_elliptic_curve_double (mpz_t P_x, mpz_t P_y, bool *P_infty,
512 mpz_t a, mpz_t p)
513 /* replaces P by 2P on the elliptic curve given by a and the prime field */
514 /* of characteristic p; b is implicit since it is assumed that the point */
515 /* lies on the curve. */
516 /* P_x and P_y are the X- and Y-coordinates of the point, respectively, */
517 /* P_infty is true if the point is the point at infinity. In this case, */
518 /* the values of P_x and P_y are not defined. */
519 /* P_x, P_y and P_infty are modified. */
520
521 {
522 mpz_t factor, x_local, tmp, tmp2;
523
524 mpz_init (factor);
525 mpz_init (x_local);
526 mpz_init (tmp);
527 mpz_init (tmp2);
528
529 if (!(*P_infty))
530 /* there is nothing to do when P is already infinity */
531 {
532 if (mpz_cmp_ui (P_y, 0ul) == 0)
533 *P_infty = true;
534 else
535 {
536 /* factor = (3*P_x^2 + a) / 2*P_y */
537 mpz_pow_ui (tmp, P_x, 2ul);
538 mpz_mod (tmp, tmp, p);
539 mpz_mul_ui (factor, tmp, 3ul);
540 mpz_add (factor, factor, a);
541 mpz_mul_2exp (tmp, P_y, 1ul);
542 mpz_invert (tmp2, tmp, p);
543 mpz_mul (tmp, factor, tmp2);
544 mpz_mod (factor, tmp, p);
545 /* P_x = factor^2 - 2 P_x; save P_x */
546 mpz_set (x_local, P_x);
547 mpz_pow_ui (P_x, factor, 2ul);
548 mpz_submul_ui (P_x, x_local, 2ul);
549 mpz_mod (P_x, P_x, p);
550 /* P_y = factor * (old x - new x) - P_y */
551 mpz_sub (tmp, x_local, P_x);
552 mpz_mul (tmp2, factor, tmp);
553 mpz_sub (P_y, tmp2, P_y);
554 mpz_mod (P_y, P_y, p);
555 }
556 }
557
558 mpz_clear (factor);
559 mpz_clear (x_local);
560 mpz_clear (tmp);
561 mpz_clear (tmp2);
562 }
563
564 /*****************************************************************************/
565
nt_elliptic_curve_add(mpz_t P_x,mpz_t P_y,bool * P_infty,mpz_t Q_x,mpz_t Q_y,const bool Q_infty,mpz_t a,mpz_t p)566 static void nt_elliptic_curve_add (mpz_t P_x, mpz_t P_y, bool *P_infty,
567 mpz_t Q_x, mpz_t Q_y, const bool Q_infty, mpz_t a, mpz_t p)
568 /* replaces P by P+Q on the elliptic curve given by a and the implicit b */
569 /* over the prime field of characteristic p */
570 /* P_x, P_y and P_infty are modified. */
571
572 {
573 mpz_t factor, x_local, tmp, tmp2;
574
575 mpz_init (factor);
576 mpz_init (x_local);
577 mpz_init (tmp);
578 mpz_init (tmp2);
579
580 if (!Q_infty)
581 {
582 if (*P_infty)
583 {
584 mpz_set (P_x, Q_x);
585 mpz_set (P_y, Q_y);
586 *P_infty = false;
587 }
588 else if (mpz_cmp (P_x, Q_x) == 0)
589 {
590 if (mpz_cmp (P_y, Q_y) == 0)
591 nt_elliptic_curve_double (P_x, P_y, P_infty, a, p);
592 else
593 *P_infty = true;
594 }
595 else
596 {
597 /* factor = Delta y / Delta x */
598 mpz_sub (tmp, Q_x, P_x);
599 mpz_invert (tmp2, tmp, p);
600 mpz_sub (tmp, Q_y, P_y);
601 mpz_mul (factor, tmp, tmp2);
602 mpz_mod (factor, factor, p);
603 /* P_x = factor^2 - (P_x + Q_x), keep a copy of P_x */
604 mpz_set (x_local, P_x);
605 mpz_pow_ui (tmp, factor, 2ul);
606 mpz_add (tmp2, P_x, Q_x);
607 mpz_sub (P_x, tmp, tmp2);
608 mpz_mod (P_x, P_x, p);
609 /* P_y = factor * (old P_x - new x) - P_y */
610 mpz_sub (tmp, x_local, P_x);
611 mpz_mul (tmp2, factor, tmp);
612 mpz_sub (P_y, tmp2, P_y);
613 mpz_mod (P_y, P_y, p);
614 }
615 }
616
617 mpz_clear (factor);
618 mpz_clear (x_local);
619 mpz_clear (tmp);
620 mpz_clear (tmp2);
621 }
622
623 /*****************************************************************************/
624
cm_nt_elliptic_curve_multiply(mpz_t P_x,mpz_t P_y,bool * P_infty,mpz_t m,mpz_t a,mpz_t p)625 void cm_nt_elliptic_curve_multiply (mpz_t P_x, mpz_t P_y, bool *P_infty,
626 mpz_t m, mpz_t a, mpz_t p)
627 /* replaces P by mP on the elliptic curve given by a and the implicit b */
628 /* over the prime field of characteristic p */
629 /* P_x, P_y and P_infty are modified. */
630
631 {
632 mpz_t Q_x, Q_y, m_local, m_new;
633 bool Q_infty = true;
634 /* m P is stored in Q; we use a right to left binary exponentiation */
635 /* scheme with the invariant Q + m P = const. */
636
637 mpz_init (Q_x);
638 mpz_init (Q_y);
639 mpz_init (m_local);
640 mpz_init (m_new);
641 *P_infty = false;
642
643 mpz_set (m_local, m);
644 if (mpz_cmp_ui (m_local, 0ul) == 0)
645 *P_infty = true;
646 else
647 {
648 if (mpz_cmp_ui (m_local, 0ul) < 0)
649 {
650 mpz_neg (m_local, m_local);
651 mpz_neg (P_y, P_y);
652 }
653 while (mpz_cmp_ui (m_local, 0ul) > 0)
654 {
655 while (mpz_tdiv_q_ui (m_new, m_local, 2ul) == 0)
656 {
657 mpz_set (m_local, m_new);
658 nt_elliptic_curve_double (P_x, P_y, P_infty, a, p);
659 }
660 mpz_sub_ui (m_local, m_local, 1ul);
661 nt_elliptic_curve_add (Q_x, Q_y, &Q_infty, P_x, P_y, *P_infty, a, p);
662 }
663 /* copy the result */
664 mpz_set (P_x, Q_x);
665 mpz_set (P_y, Q_y);
666 *P_infty = Q_infty;
667 }
668
669 mpz_clear (Q_x);
670 mpz_clear (Q_y);
671 mpz_clear (m_local);
672 mpz_clear (m_new);
673 }
674
675 /*****************************************************************************/
676
cm_nt_elliptic_curve_random(mpz_t P_x,mpz_t P_y,mpz_t cofactor,mpz_t a,mpz_t b,mpz_t p)677 void cm_nt_elliptic_curve_random (mpz_t P_x, mpz_t P_y,
678 mpz_t cofactor, mpz_t a, mpz_t b, mpz_t p)
679 /* creates a point on the elliptic curve given by a and b over the prime */
680 /* field of characteristic p and multiplies it by the cofactor until */
681 /* the result is different from infinity. If the curve order is cofactor */
682 /* times a prime, this results in a point of order this prime. */
683 /* The point is not really random, since successive X-coordinates from */
684 /* 1 on are tested. */
685 /* P_x and P_y are modified. */
686
687 {
688 mpz_t tmp;
689 long unsigned int P_x_long = 0;
690 bool P_infty = true;
691
692 mpz_init (tmp);
693 while (P_infty)
694 {
695 P_x_long++;
696 /* P_y = P_x^3 + a P_x + b */
697 mpz_mul_ui (P_y, a, P_x_long);
698 mpz_add (P_y, P_y, b);
699 mpz_add_ui (P_y, P_y, P_x_long * P_x_long * P_x_long);
700 mpz_mod (P_y, P_y, p);
701 /* try to compute the square root of P_y */
702 if (mpz_jacobi (P_y, p) != -1)
703 {
704 mpz_set_ui (P_x, P_x_long);
705 cm_nt_mpz_tonelli_z (P_y, P_y, p);
706 /* get rid of the cofactor */
707 P_infty = false;
708 cm_nt_elliptic_curve_multiply (P_x, P_y, &P_infty, cofactor, a, p);
709 }
710 }
711 mpz_clear (tmp);
712 }
713
714 /*****************************************************************************/
715 /*****************************************************************************/
716
cm_nt_fget_z(mpz_t out,ftype in)717 bool cm_nt_fget_z (mpz_t out, ftype in)
718 /* tries to round the real to an integer value */
719 /* The return value reflects the success of the operation. */
720
721 {
722 ftype rounded, diff;
723 bool ok;
724 mp_exp_t expo;
725
726 finit (rounded, fget_prec (in));
727 finit (diff, fget_prec (in));
728
729 fround (rounded, in);
730 fsub (diff, in, rounded);
731 if (fsgn (diff) == 0 || (- fget_exp (diff)) >= 10) {
732 expo = fget_z_exp (out, rounded);
733 if (expo > 0)
734 ok = false;
735 else if (expo < 0) {
736 mpz_tdiv_q_2exp (out, out, (unsigned long int) (-expo));
737 ok = true;
738 }
739 else
740 ok = true;
741 }
742 else
743 ok = false;
744
745 fclear (rounded);
746 fclear (diff);
747
748 return ok;
749 }
750
751 /*****************************************************************************/
752 /*****************************************************************************/
753