1 /*
2 Copyright (C) 2014 William Hart
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 #define ulong ulongxx /* interferes with system includes */
13 #include <stdlib.h>
14 #include <stdio.h>
15 #include <math.h>
16 #undef ulong
17 #define ulong mp_limb_t
18 #include <gmp.h>
19 #include "flint.h"
20 #include "ulong_extras.h"
21 #include "fmpz.h"
22 #include "fmpz_vec.h"
23
fmpz_lucas_chain(fmpz_t Vm,fmpz_t Vm1,const fmpz_t A,const fmpz_t m,const fmpz_t n)24 void fmpz_lucas_chain(fmpz_t Vm, fmpz_t Vm1, const fmpz_t A,
25 const fmpz_t m, const fmpz_t n)
26 {
27 fmpz_t t;
28 slong i, B = fmpz_sizeinbase(m, 2);
29
30 fmpz_init(t);
31 fmpz_set_ui(Vm, 2);
32 fmpz_set(Vm1, A);
33
34 for (i = B - 1; i >= 0; i--)
35 {
36 if (fmpz_tstbit(m, i)) /* 1 in binary repn */
37 {
38 fmpz_mul(t, Vm, Vm1);
39 fmpz_sub(t, t, A);
40 fmpz_mod(Vm, t, n);
41
42 fmpz_mul(t, Vm1, Vm1);
43 fmpz_sub_ui(t, t, 2);
44 fmpz_mod(Vm1, t, n);
45 } else /* 0 in binary repn */
46 {
47 fmpz_mul(t, Vm, Vm1);
48 fmpz_sub(t, t, A);
49 fmpz_mod(Vm1, t, n);
50
51 fmpz_mul(t, Vm, Vm);
52 fmpz_sub_ui(t, t, 2);
53 fmpz_mod(Vm, t, n);
54 }
55 }
56
57 fmpz_clear(t);
58 }
59
fmpz_lucas_chain_full(fmpz_t Vm,fmpz_t Vm1,const fmpz_t A,const fmpz_t B,const fmpz_t m,const fmpz_t n)60 void fmpz_lucas_chain_full(fmpz_t Vm, fmpz_t Vm1, const fmpz_t A, const fmpz_t B,
61 const fmpz_t m, const fmpz_t n)
62 {
63 fmpz_t t, Q;
64 slong i, bits = fmpz_sizeinbase(m, 2);
65
66 fmpz_init(t);
67 fmpz_init(Q);
68 fmpz_set_ui(Q, 1);
69 fmpz_set_ui(Vm, 2);
70 fmpz_set(Vm1, A);
71
72 for (i = bits - 1; i >= 0; i--)
73 {
74 if (fmpz_tstbit(m, i)) /* 1 in binary repn */
75 {
76 fmpz_mul(t, Vm1, Vm);
77 fmpz_submul(t, Q, A);
78 fmpz_mod(Vm, t, n);
79
80 fmpz_mul(Vm1, Vm1, Vm1);
81 fmpz_mul_ui(t, Q, 2);
82 fmpz_mul(t, t, B);
83 fmpz_sub(Vm1, Vm1, t);
84 fmpz_mod(Vm1, Vm1, n);
85
86 fmpz_mul(Q, Q, Q);
87 fmpz_mul(Q, Q, B);
88 fmpz_mod(Q, Q, n);
89 } else /* 0 in binary repn */
90 {
91 fmpz_mul(t, Vm, Vm1);
92 fmpz_submul(t, Q, A);
93 fmpz_mod(Vm1, t, n);
94
95 fmpz_mul(t, Vm, Vm);
96 fmpz_submul_ui(t, Q, 2);
97 fmpz_mod(Vm, t, n);
98
99 fmpz_mul(Q, Q, Q);
100 fmpz_mod(Q, Q, n);
101 }
102 }
103
104 fmpz_clear(Q);
105 fmpz_clear(t);
106 }
107
108 /* Compute U_{2m}, U_{2m + 1} given U_m, U_{m + 1} */
fmpz_lucas_chain_double(fmpz_t U2m,fmpz_t U2m1,const fmpz_t Um,const fmpz_t Um1,const fmpz_t A,const fmpz_t B,const fmpz_t n)109 void fmpz_lucas_chain_double(fmpz_t U2m, fmpz_t U2m1, const fmpz_t Um,
110 const fmpz_t Um1, const fmpz_t A, const fmpz_t B,
111 const fmpz_t n)
112 {
113 fmpz_t t, t2;
114
115 fmpz_init(t);
116 fmpz_init(t2);
117
118 fmpz_mul_2exp(t, Um1, 1); /* U_m(2U_{m+1) - AU_m) */
119 fmpz_submul(t, Um, A);
120 fmpz_mul(t, t, Um);
121
122 fmpz_mul(U2m1, Um1, Um1); /* U_{m+1}^2 - BU_m^2 */
123 fmpz_mul(t2, Um, Um);
124 fmpz_submul(U2m1, t2, B);
125 fmpz_mod(U2m1, U2m1, n);
126
127 fmpz_mod(U2m, t, n);
128
129 fmpz_clear(t);
130 fmpz_clear(t2);
131 }
132
133 /*
134 Compute U_{m + n}, U_{m + n + 1} given U_m, U_{m + 1} and
135 U_n, U_{n + 1}
136 */
fmpz_lucas_chain_add(fmpz_t Umn,fmpz_t Umn1,const fmpz_t Um,const fmpz_t Um1,const fmpz_t Un,const fmpz_t Un1,const fmpz_t A,const fmpz_t B,const fmpz_t n)137 void fmpz_lucas_chain_add(fmpz_t Umn, fmpz_t Umn1, const fmpz_t Um,
138 const fmpz_t Um1, const fmpz_t Un,
139 const fmpz_t Un1, const fmpz_t A, const fmpz_t B,
140 const fmpz_t n)
141 {
142 fmpz_t t, t2;
143
144 fmpz_init(t);
145 fmpz_init(t2);
146
147 fmpz_mul(t, Un, A); /* U_nU_{m + 1} - BU_m(AU_n - U_{n + 1})/B */
148 fmpz_sub(t, Un1, t);
149 fmpz_mul(t, t, Um);
150 fmpz_addmul(t, Un, Um1);
151
152 fmpz_mul(Umn1, Un1, Um1); /* U_{n + 1}U_{m + 1} - BU_mU_n */
153 fmpz_mul(t2, Um, Un);
154 fmpz_submul(Umn1, t2, B);
155 fmpz_mod(Umn1, Umn1, n);
156
157 fmpz_mod(Umn, t, n);
158
159 fmpz_clear(t);
160 fmpz_clear(t2);
161 }
162
163 /* Compute U_{km}, U_{km + 1} from U_m, U_{m + 1}, k > 0 */
fmpz_lucas_chain_mul(fmpz_t Ukm,fmpz_t Ukm1,const fmpz_t Um,const fmpz_t Um1,const fmpz_t A,const fmpz_t B,const fmpz_t k,const fmpz_t n)164 void fmpz_lucas_chain_mul(fmpz_t Ukm, fmpz_t Ukm1,
165 const fmpz_t Um, const fmpz_t Um1,
166 const fmpz_t A, const fmpz_t B, const fmpz_t k,
167 const fmpz_t n)
168 {
169 slong i = 0, b = fmpz_sizeinbase(k, 2);
170 fmpz_t t, t1;
171
172 fmpz_init(t);
173 fmpz_init(t1);
174
175 fmpz_set(Ukm, Um);
176 fmpz_set(Ukm1, Um1);
177
178 while (!fmpz_tstbit(k, i))
179 {
180 fmpz_lucas_chain_double(Ukm, Ukm1, Ukm, Ukm1, A, B, n);
181 i++;
182 }
183
184 i++;
185
186 if (i < b)
187 {
188 fmpz_set(t, Ukm);
189 fmpz_set(t1, Ukm1);
190 }
191
192 while (i < b)
193 {
194 fmpz_lucas_chain_double(t, t1, t, t1, A, B, n);
195 if (fmpz_tstbit(k, i))
196 fmpz_lucas_chain_add(Ukm, Ukm1, Ukm, Ukm1, t, t1, A, B, n);
197 i++;
198 }
199
200 fmpz_clear(t);
201 fmpz_clear(t1);
202 }
203
204 /* Compute U_m, U_{m + 1} from V_m, V_{m + 1} */
fmpz_lucas_chain_VtoU(fmpz_t Um,fmpz_t Um1,const fmpz_t Vm,const fmpz_t Vm1,const fmpz_t A,const fmpz_t B,const fmpz_t Dinv,const fmpz_t n)205 void fmpz_lucas_chain_VtoU(fmpz_t Um, fmpz_t Um1,
206 const fmpz_t Vm, const fmpz_t Vm1,
207 const fmpz_t A, const fmpz_t B, const fmpz_t Dinv,
208 const fmpz_t n)
209 {
210 fmpz_t t;
211
212 fmpz_init(t);
213
214 fmpz_mul_2exp(t, Vm1, 1); /* (2V_{m + 1} - AV_m) / D */
215 fmpz_submul(t, Vm, A);
216 fmpz_mul(t, t, Dinv);
217
218 fmpz_set(Um1, Vm);
219 fmpz_mod(Um, t, n);
220
221 fmpz_addmul(Um1, Um, A); /* (V_m + AU_m) / 2 */
222 if (!fmpz_is_even(Um1))
223 fmpz_add(Um1, Um1, n);
224 fmpz_tdiv_q_2exp(Um1, Um1, 1);
225
226 fmpz_mod(Um1, Um1, n);
227
228 fmpz_clear(t);
229 }
230