1 /*
2 ref_mul.c: reference implementations for polynomial multiplication,
3 middle product, scalar multiplication, integer middle product
4
5 Copyright (C) 2007, 2008, David Harvey
6
7 This file is part of the zn_poly library (version 0.9).
8
9 This program is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 2 of the License, or
12 (at your option) version 3 of the License.
13
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with this program. If not, see <http://www.gnu.org/licenses/>.
21
22 */
23
24 #include "support.h"
25 #include "zn_poly_internal.h"
26 #include <string.h>
27
28
29 /*
30 Sets x = op[0] + op[1]*B + ... + op[n-1]*B^(n-1), where B = 2^b.
31
32 Running time is soft-linear in output length.
33 */
34 void
pack(mpz_t x,const ulong * op,size_t n,unsigned b)35 pack (mpz_t x, const ulong* op, size_t n, unsigned b)
36 {
37 ZNP_ASSERT (n >= 1);
38
39 if (n == 1)
40 {
41 // base case
42 mpz_set_ui (x, op[0]);
43 }
44 else
45 {
46 // recursively split into top and bottom halves
47 mpz_t y;
48 mpz_init (y);
49 pack (x, op, n / 2, b);
50 pack (y, op + n / 2, n - n / 2, b);
51 mpz_mul_2exp (y, y, (n / 2) * b);
52 mpz_add (x, x, y);
53 mpz_clear (y);
54 }
55 }
56
57
58
59 /*
60 Inverse operation of pack(), with output coefficients reduced mod m.
61
62 Running time is soft-linear in output length.
63 */
64 void
unpack(ulong * res,const mpz_t op,size_t n,unsigned b,ulong m)65 unpack (ulong* res, const mpz_t op, size_t n, unsigned b, ulong m)
66 {
67 ZNP_ASSERT (n >= 1);
68 ZNP_ASSERT (mpz_sizeinbase (op, 2) <= n * b);
69
70 mpz_t y;
71 mpz_init(y);
72
73 if (n == 1)
74 {
75 // base case
76 mpz_set (y, op);
77 mpz_fdiv_r_ui (y, y, m);
78 *res = mpz_get_ui (y);
79 }
80 else
81 {
82 // recursively split into top and bottom halves
83 mpz_tdiv_q_2exp (y, op, (n / 2) * b);
84 unpack (res + n / 2, y, n - n / 2, b, m);
85 mpz_tdiv_r_2exp (y, op, (n / 2) * b);
86 unpack (res, y, n / 2, b, m);
87 }
88
89 mpz_clear (y);
90 }
91
92
93
94 /*
95 Reference implementation of zn_array_mul().
96 Very simple Kronecker substitution, uses GMP for multiplication.
97 */
98 void
ref_zn_array_mul(ulong * res,const ulong * op1,size_t n1,const ulong * op2,size_t n2,const zn_mod_t mod)99 ref_zn_array_mul (ulong* res,
100 const ulong* op1, size_t n1,
101 const ulong* op2, size_t n2,
102 const zn_mod_t mod)
103 {
104 ZNP_ASSERT (n2 >= 1);
105 ZNP_ASSERT (n1 >= n2);
106
107 mpz_t x, y;
108 mpz_init (x);
109 mpz_init (y);
110
111 unsigned b = 2 * mod->bits + ceil_lg (n2);
112 unsigned words = CEIL_DIV (b, ULONG_BITS);
113
114 pack (x, op1, n1, b);
115 pack (y, op2, n2, b);
116 mpz_mul (x, x, y);
117 unpack (res, x, n1 + n2 - 1, b, mod->m);
118
119 mpz_clear (y);
120 mpz_clear (x);
121 }
122
123
124
125 /*
126 Reference implementation of zn_array_mulmid().
127 Just calls ref_zn_array_mul() and extracts relevant part of output.
128 */
129 void
ref_zn_array_mulmid(ulong * res,const ulong * op1,size_t n1,const ulong * op2,size_t n2,const zn_mod_t mod)130 ref_zn_array_mulmid (ulong* res,
131 const ulong* op1, size_t n1,
132 const ulong* op2, size_t n2,
133 const zn_mod_t mod)
134 {
135 ZNP_ASSERT (n2 >= 1);
136 ZNP_ASSERT (n1 >= n2);
137
138 ulong* temp = (ulong*) malloc ((n1 + n2 - 1) * sizeof (ulong));
139 ref_zn_array_mul (temp, op1, n1, op2, n2, mod);
140 ulong i;
141 for (i = 0; i < n1 - n2 + 1; i++)
142 res[i] = temp[i + n2 - 1];
143 free (temp);
144 }
145
146
147
148 /*
149 Reference implementation of negacyclic multiplication.
150
151 Multiplies op1[0, n) by op2[0, n) negacyclically, puts result into
152 res[0, n).
153 */
154 void
ref_zn_array_negamul(ulong * res,const ulong * op1,const ulong * op2,size_t n,const zn_mod_t mod)155 ref_zn_array_negamul (ulong* res, const ulong* op1, const ulong* op2,
156 size_t n, const zn_mod_t mod)
157 {
158 ulong* temp = (ulong*) malloc (sizeof (ulong) * 2 * n);
159
160 ref_zn_array_mul (temp, op1, n, op2, n, mod);
161 temp[2 * n - 1] = 0;
162
163 mpz_t x;
164 mpz_init (x);
165
166 size_t i;
167 for (i = 0; i < n; i++)
168 {
169 mpz_set_ui (x, temp[i]);
170 mpz_sub_ui (x, x, temp[i + n]);
171 mpz_mod_ui (x, x, mod->m);
172 res[i] = mpz_get_ui (x);
173 }
174
175 mpz_clear (x);
176 free (temp);
177 }
178
179
180
181 /*
182 Reference implementation of scalar multiplication.
183
184 Multiplies op[0, n) by x, puts result in res[0, n).
185 */
186 void
ref_zn_array_scalar_mul(ulong * res,const ulong * op,size_t n,ulong x,const zn_mod_t mod)187 ref_zn_array_scalar_mul (ulong* res, const ulong* op, size_t n,
188 ulong x, const zn_mod_t mod)
189 {
190 mpz_t y;
191 mpz_init (y);
192
193 size_t i;
194 for (i = 0; i < n; i++)
195 {
196 mpz_set_ui (y, op[i]);
197 mpz_mul_ui (y, y, x);
198 mpz_mod_ui (y, y, mod->m);
199 res[i] = mpz_get_ui (y);
200 }
201
202 mpz_clear (y);
203 }
204
205
206
207 /*
208 Reference implementation of mpn_smp.
209
210 Computes SMP(op1[0, n1), op2[0, n2)), stores result at res[0, n1 - n2 + 3).
211 */
212 void
ref_mpn_smp(mp_limb_t * res,const mp_limb_t * op1,size_t n1,const mp_limb_t * op2,size_t n2)213 ref_mpn_smp (mp_limb_t* res,
214 const mp_limb_t* op1, size_t n1,
215 const mp_limb_t* op2, size_t n2)
216 {
217 ZNP_ASSERT (n2 >= 1);
218 ZNP_ASSERT (n1 >= n2);
219
220 mp_limb_t* prod = (mp_limb_t*) malloc (sizeof (mp_limb_t) * (n1 + n2));
221
222 // first compute the ordinary product
223 mpn_mul (prod, op1, n1, op2, n2);
224
225 // now we want to remove the cross-terms that could possibly interfere
226 // with the result we want, i.e. in the following diagram, we want only
227 // contributions from O's, but mpn_mul has given us all of O, A and X,
228 // and we will remove the A's.
229
230 // OOOOAAXX
231 // AOOOOAAX
232 // AAOOOOAA
233 // XAAOOOOA
234 // XXAAOOOO
235
236 int which; // 0 == bottom-left corner, 1 == top-right corner
237 size_t diag; // 0 == closest to diagonal, 1 == next diagonal
238 size_t i, x, y, off;
239 mp_limb_t lo, hi;
240
241 for (which = 0; which <= 1; which++)
242 for (diag = 0; diag < ZNP_MIN (n2 - 1, 2); diag++)
243 for (i = 0; i < n2 - 1 - diag; i++)
244 {
245 x = n2 - 2 - i - diag;
246 y = i;
247 if (which)
248 {
249 x = n1 - 1 - x;
250 y = n2 - 1 - y;
251 }
252 off = x + y;
253 hi = mpn_mul_1 (&lo, op1 + x, 1, op2[y]);
254 mpn_sub_1 (prod + off, prod + off, n1 + n2 - off, lo);
255 mpn_sub_1 (prod + off + 1, prod + off + 1, n1 + n2 - off - 1, hi);
256 }
257
258 // copy the result to the output array
259 memcpy (res, prod + n2 - 1, sizeof (mp_limb_t) * (n1 - n2 + 2));
260 res[n1 - n2 + 2] = (n2 > 1) ? prod[n1 + 1] : 0;
261 free (prod);
262 }
263
264
265 /*
266 Reference implementation of mpn_mulmid.
267
268 Let P = product op1 * op2. Computes P[n2 + 1, n1), stores result at
269 res[2, n1 - n2 + 1).
270 */
271 void
ref_mpn_mulmid(mp_limb_t * res,const mp_limb_t * op1,size_t n1,const mp_limb_t * op2,size_t n2)272 ref_mpn_mulmid (mp_limb_t* res, const mp_limb_t* op1, size_t n1,
273 const mp_limb_t* op2, size_t n2)
274 {
275 ZNP_ASSERT (n2 >= 1);
276 ZNP_ASSERT (n1 >= n2);
277
278 mp_limb_t* prod = (mp_limb_t*) malloc (sizeof (mp_limb_t) * (n1 + n2));
279
280 // compute the ordinary product
281 mpn_mul (prod, op1, n1, op2, n2);
282
283 // copy relevant segment to output
284 if (n1 > n2)
285 memcpy (res + 2, prod + n2 + 1, sizeof (mp_limb_t) * (n1 - n2 - 1));
286
287 free (prod);
288 }
289
290
291 // end of file ****************************************************************
292