Line data Source code
1 : /* Implements algorithm polyeval and remainder tree using middle product.
2 :
3 : Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2009 Laurent Fousse,
4 : Alexander Kruppa, Paul Zimmermann.
5 :
6 : This file is part of the ECM Library.
7 :
8 : The ECM Library is free software; you can redistribute it and/or modify
9 : it under the terms of the GNU Lesser General Public License as published by
10 : the Free Software Foundation; either version 3 of the License, or (at your
11 : option) any later version.
12 :
13 : The ECM Library is distributed in the hope that it will be useful, but
14 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 : License for more details.
17 :
18 : You should have received a copy of the GNU Lesser General Public License
19 : along with the ECM Library; see the file COPYING.LIB. If not, see
20 : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 :
23 : #include <stdlib.h>
24 : #include <string.h> /* for strlen */
25 : #include "ecm-impl.h"
26 :
27 : #ifdef HAVE_UNISTD_H
28 : # include <unistd.h> /* for unlink */
29 : #endif
30 :
31 :
32 : #ifndef MAX
33 : #define MAX(a,b) (((a) > (b)) ? (a) : (b))
34 : #endif
35 :
36 : /* #define DEBUG_TREEDATA */
37 :
38 : extern unsigned int Fermat;
39 :
40 : #if defined(DEBUG) || defined(DEBUG_TREEDATA)
41 : void
42 : print_vect (listz_t t, unsigned int l)
43 : {
44 : unsigned int i;
45 :
46 : fprintf (ECM_STDOUT, "[");
47 : for (i = 0; i < l; i++)
48 : {
49 : mpz_out_str (ECM_STDOUT, 10, t[i]);
50 : if (i != l - 1)
51 : fprintf (ECM_STDOUT, ", ");
52 : else
53 : fprintf (ECM_STDOUT, "]");
54 : }
55 : }
56 : #endif
57 :
58 : /* Computes TUpTree as described in ref[1]. k is the degree of the
59 : * polynomial at the root of the tree. sh is the shift we need to
60 : * apply to find the actual coefficients of the polynomial at the root
61 : * of the tree.
62 : */
63 :
64 : void
65 2890170 : TUpTree (listz_t b, listz_t *Tree, unsigned int k, listz_t tmp, int dolvl,
66 : unsigned int sh, mpz_t n, FILE *TreeFile)
67 : {
68 : unsigned int m, l;
69 :
70 2890170 : m = k / 2;
71 2890170 : l = k - m;
72 :
73 2890170 : if (k == 1)
74 295965 : return;
75 :
76 : #ifdef DEBUG
77 : fprintf (ECM_STDOUT, "In TupTree, k = %d.\n", k);
78 :
79 : fprintf (ECM_STDOUT, "b = ");
80 : print_vect (b, k);
81 : fprintf (ECM_STDOUT, "\nThe polynomials at that level are: ");
82 : print_vect (Tree[0] + sh, k);
83 : fprintf (ECM_STDOUT, "\n");
84 : #endif
85 :
86 2594205 : if (dolvl == 0 || dolvl == -1)
87 : {
88 1447765 : if (TreeFile != NULL)
89 : {
90 56713 : list_inp_raw (tmp + k, TreeFile, l);
91 : #ifdef DEBUG_TREEDATA
92 : printf ("Read from file: ");
93 : print_vect (tmp + k, l);
94 : #endif
95 56713 : TMulGen (tmp + l, m - 1, tmp + k, l - 1, b, k - 1, tmp + k + l, n);
96 56713 : list_inp_raw (tmp + k, TreeFile, m);
97 : #ifdef DEBUG_TREEDATA
98 : print_vect (tmp + k, m);
99 : printf ("\n");
100 : #endif
101 56713 : TMulGen (tmp, l - 1, tmp + k, m - 1, b, k - 1, tmp + k + m, n);
102 : }
103 : else
104 : {
105 : #ifdef DEBUG_TREEDATA
106 : printf ("Got from Tree: ");
107 : print_vect (Tree[0] + sh, l);
108 : print_vect (Tree[0] + sh + l, m);
109 : printf ("\n");
110 : #endif
111 1391052 : TMulGen (tmp + l, m - 1, Tree[0] + sh, l - 1, b, k - 1, tmp + k, n);
112 1391052 : TMulGen (tmp, l - 1, Tree[0] + sh + l, m - 1, b, k - 1, tmp + k, n);
113 : }
114 :
115 : #if defined(DEBUG) || defined (DEBUG_TREEDATA)
116 : fprintf (ECM_STDOUT, "And the result at that level (before correction) is:");
117 : print_vect (tmp, k);
118 : fprintf (ECM_STDOUT, "\n");
119 : #endif
120 :
121 : /* GMP-ECM specific: leading coefficients in the product tree
122 : * are implicit ones, so we need some extra work here.
123 : */
124 :
125 1447765 : list_add (tmp, tmp, b + m, l);
126 1447765 : list_add (tmp + l, tmp + l, b + l, m);
127 :
128 1447765 : list_mod (b, tmp, k, n); /* reduce both parts simultaneously */
129 :
130 : #ifdef DEBUG
131 : fprintf (ECM_STDOUT, "And the result at this level is:");
132 : print_vect (b, k);
133 : fprintf (ECM_STDOUT, "\n");
134 : #endif
135 : }
136 :
137 2594205 : if (dolvl > 0 || dolvl == -1)
138 : {
139 1442102 : if (dolvl > 0)
140 1146440 : dolvl--;
141 1442102 : TUpTree (b, Tree + 1, l, tmp, dolvl, sh, n, TreeFile);
142 1442102 : TUpTree (b + l, Tree + 1, m, tmp, dolvl, sh + l, n, TreeFile);
143 : }
144 : }
145 :
146 : static unsigned int
147 2434 : TUpTree_space (unsigned int k)
148 : {
149 :
150 : unsigned int m, l;
151 : unsigned int r1, r2;
152 :
153 2434 : m = k / 2;
154 2434 : l = k - m;
155 :
156 2434 : if (k == 1)
157 635 : return 0;
158 :
159 1799 : r1 = TMulGen_space (l - 1, m - 1, k - 1) + l;
160 1799 : if (m != l)
161 : {
162 514 : r2 = TMulGen_space (m - 1, l - 1, k - 1) + k;
163 514 : r1 = MAX (r1, r2);
164 : }
165 :
166 1799 : r2 = TUpTree_space (l);
167 1799 : r1 = MAX (r1, r2);
168 :
169 1799 : if (m != l)
170 : {
171 514 : r2 = TUpTree_space (m);
172 514 : r1 = MAX (r1, r2);
173 : }
174 :
175 1799 : return r1;
176 : }
177 :
178 : /* This is the documentation of the (now removed) polyeval() function.
179 : Algorithm polyeval from section 3.7 of Peter Montgomery's dissertation.
180 : Input:
181 : G - an array of k elements of R, G[i], 0 <= i < k
182 : representing the coefficients of a polynomial G(x) of degree < k
183 : Tree - the product tree produced by PolyFromRoots
184 : Tree[0][0..k-1] (degree k/2)
185 : Tree[1][0..k-1] (degree k/4), ...,
186 : Tree[lgk-1][0..k-1] (degree 1)
187 : Output: the sequence of values of G(a[i]) are stored in G[i] for 0 <= i < k
188 : Remark: we need an auxiliary (k+1)-th cell G[k] in G.
189 : The memory used is M(k) = max(3*floor(k/2)+list_mul_mem(floor(k/2)),
190 : k+list_mul_mem(ceil(k/2)),
191 : floor(k/2) + M(ceil(k/2))).
192 : Since list_mul_mem(k) >= 2*k, the maximum is the 1st.
193 : */
194 :
195 : /* Same as polyeval. Needs invF as extra argument.
196 : Return non-zero iff an error occurred.
197 : */
198 : int
199 121 : polyeval_tellegen (listz_t b, unsigned int k, listz_t *Tree, listz_t tmp,
200 : unsigned int sizeT, listz_t invF, mpz_t n,
201 : char *TreeFilename)
202 : {
203 : unsigned int tupspace;
204 : unsigned int tkspace;
205 121 : int allocated = 0,
206 121 : r = 0; /* return value, 0 = no error */
207 : listz_t T;
208 :
209 : ASSERT(Tree != NULL || TreeFilename != NULL);
210 :
211 121 : tupspace = TUpTree_space (k) + k;
212 121 : tkspace = 2 * k - 1 + list_mul_mem (k);
213 :
214 121 : tupspace = MAX (tupspace, tkspace);
215 :
216 121 : if (TreeFilename != NULL)
217 2 : tupspace += (k + 1) / 2;
218 :
219 121 : if (sizeT >= tupspace)
220 116 : T = tmp;
221 : else
222 : {
223 5 : outputf (OUTPUT_DEVVERBOSE, "polyeval_tellegen: allocating extra temp"
224 : " space, want %d but T has only %d\n", tupspace, sizeT);
225 5 : T = init_list (tupspace);
226 5 : if (T == NULL)
227 0 : return ECM_ERROR;
228 5 : allocated = 1;
229 : }
230 :
231 : #ifdef TELLEGEN_DEBUG
232 : fprintf (ECM_STDOUT, "In polyeval_tellegen, k = %d.\n", k);
233 : fprintf (ECM_STDOUT, "Required memory: %d.\n",
234 : TMulGen_space (k - 1, k - 1, k - 1));
235 : #endif
236 :
237 121 : if (Fermat)
238 : {
239 : /* Schoenhage-Strassen can't do a half product faster than a full */
240 4 : F_mul (T, invF, b, k, DEFAULT, Fermat, T + 2 * k);
241 4 : list_mod (T, T + k - 1, k, n);
242 : }
243 : else
244 : {
245 : /* need space 2k-1+list_mul_mem(k) in T */
246 117 : list_mul_high (T, invF, b, k);
247 117 : list_mod (T, T + k - 1, k, n);
248 : }
249 121 : list_revert (T, k);
250 121 : if (TreeFilename != NULL)
251 : {
252 : unsigned int lgk, i;
253 : FILE *TreeFile;
254 2 : char *fullname = (char *) malloc (strlen (TreeFilename) + 1 + 2 + 1);
255 2 : ASSERT_ALWAYS(fullname != NULL);
256 :
257 2 : lgk = ceil_log2 (k);
258 18 : for (i = 0; i < lgk; i++)
259 : {
260 16 : sprintf (fullname, "%s.%d", TreeFilename, i);
261 :
262 16 : TreeFile = fopen (fullname, "rb");
263 16 : if (TreeFile == NULL)
264 : {
265 0 : outputf (OUTPUT_ERROR,
266 : "Error opening file %s for product tree of F\n",
267 : fullname);
268 0 : r = ECM_ERROR;
269 0 : goto clear_T;
270 : }
271 16 : TUpTree (T, NULL, k, T + k, i, 0, n, TreeFile);
272 16 : fclose (TreeFile);
273 16 : unlink (fullname);
274 : }
275 2 : free (fullname);
276 : }
277 : else
278 119 : TUpTree (T, Tree, k, T + k, -1, 0, n, NULL);
279 121 : list_swap (b, T, k); /* more efficient than list_set, since T is not
280 : needed anymore */
281 :
282 121 : clear_T:
283 121 : if (allocated)
284 5 : clear_list (T, tupspace);
285 :
286 121 : return r;
287 : }
|
