1 /*
2     Copyright (C) 2015 Kushagra Singh
3 
4     This file is part of FLINT.
5 
6     FLINT is free software: you can redistribute it and/or modify it under
7     the terms of the GNU Lesser General Public License (LGPL) as published
8     by the Free Software Foundation; either version 2.1 of the License, or
9     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
10 */
11 
12 #include <gmp.h>
13 #include "flint.h"
14 #include "fmpz.h"
15 #include "mpn_extras.h"
16 
17 /* P (x0 : z0) <- kP using Montgomery ladder algorithm */
18 
19 /* tstbit uses 0 based indexing */
20 
21 void
fmpz_factor_ecm_mul_montgomery_ladder(mp_ptr x,mp_ptr z,mp_ptr x0,mp_ptr z0,mp_limb_t k,mp_ptr n,ecm_t ecm_inf)22 fmpz_factor_ecm_mul_montgomery_ladder(mp_ptr x, mp_ptr z, mp_ptr x0, mp_ptr z0,
23                                       mp_limb_t k, mp_ptr n, ecm_t ecm_inf)
24 {
25     mp_ptr x1, z1, x2, z2;      /* Q (x1 : z1), P (x2 : z2) */
26     mp_limb_t len;
27 
28     TMP_INIT;
29 
30     if (k == 0)
31     {
32         mpn_zero(x, ecm_inf->n_size);
33         mpn_zero(z, ecm_inf->n_size);
34         return;
35     }
36 
37     if (k == 1)
38     {
39         flint_mpn_copyi(x, x0, ecm_inf->n_size);
40         flint_mpn_copyi(z, z0, ecm_inf->n_size);
41         return;
42     }
43 
44     TMP_START;
45     x1 = TMP_ALLOC(ecm_inf->n_size * sizeof(mp_limb_t));
46     z1 = TMP_ALLOC(ecm_inf->n_size * sizeof(mp_limb_t));
47     x2 = TMP_ALLOC(ecm_inf->n_size * sizeof(mp_limb_t));
48     z2 = TMP_ALLOC(ecm_inf->n_size * sizeof(mp_limb_t));
49 
50 
51     flint_mpn_copyi(x1, x0, ecm_inf->n_size);    /* Q <- P0 */
52     flint_mpn_copyi(z1, z0, ecm_inf->n_size);
53     mpn_zero(x2, ecm_inf->n_size);
54     mpn_zero(z2, ecm_inf->n_size);
55 
56     fmpz_factor_ecm_double(x2, z2, x0, z0, n, ecm_inf);   /* P <- 2P0 */
57 
58     len = n_sizeinbase(k, 2) - 2;
59 
60     while (1)
61     {
62         if (((UWORD(1) << len) & k) != 0)       /* ith bit is 1 */
63         {
64             /* Q <- P + Q */
65             fmpz_factor_ecm_add(x1, z1, x1, z1, x2, z2, x0, z0, n, ecm_inf);
66             /* P <- 2 * P */
67             fmpz_factor_ecm_double(x2, z2, x2, z2, n, ecm_inf);
68         }
69         else
70         {
71             /* P <- P + Q */
72             fmpz_factor_ecm_add(x2, z2, x1, z1, x2, z2, x0, z0, n, ecm_inf);
73             /* Q <- 2 * Q */
74             fmpz_factor_ecm_double(x1, z1, x1, z1, n, ecm_inf);
75         }
76 
77         if (len == 0)
78             break;
79         else
80             len -= 1;
81     }
82 
83     flint_mpn_copyi(x, x1, ecm_inf->n_size);
84     flint_mpn_copyi(z, z1, ecm_inf->n_size);
85 
86     TMP_END;
87 }
88