1 /* ---------------------------------------------------------------------
2 *
3 * -- Integer Matrix Library (IML)
4 * (C) Copyright 2004, 2006 All Rights Reserved
5 *
6 * -- IML routines -- Version 1.0.1 -- November, 2006
7 *
8 * Author : Zhuliang Chen
9 * Contributor(s) : Arne Storjohann
10 * University of Waterloo -- School of Computer Science
11 * Waterloo, Ontario, N2L3G1 Canada
12 *
13 * ---------------------------------------------------------------------
14 *
15 * -- Copyright notice and Licensing terms:
16 *
17 * Redistribution and use in source and binary forms, with or without
18 * modification, are permitted provided that the following conditions
19 * are met:
20 *
21 * 1. Redistributions of source code must retain the above copyright
22 * notice, this list of conditions and the following disclaimer.
23 * 2. Redistributions in binary form must reproduce the above copyright
24 * notice, this list of conditions, and the following disclaimer in
25 * the documentation and/or other materials provided with the distri-
26 * bution.
27 * 3. The name of the University, the IML group, or the names of its
28 * contributors may not be used to endorse or promote products deri-
29 * ved from this software without specific written permission.
30 *
31 * -- Disclaimer:
32 *
33 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
34 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
35 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
36 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
37 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPE-
38 * CIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
39 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
40 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEO-
41 * RY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (IN-
42 * CLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
43 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
44 *
45 */
46
47
48 #include "RNSop.h"
49
50 /*
51 *
52 * Calling Sequence:
53 * basisExt(len, n, p, basis, cmbasis, cumprod, bdcoeff, R, RE)
54 *
55 * Summary:
56 * Given a representation of a matrix/vector in some RNS, extend to compute
57 * the representation in another positive integer
58 *
59 * Description:
60 * Let R be the representation of a matrix/vector M in a residue basis
61 * 'basis', i.e., R[i] = mod(M, basis[i]) (i = 0..len-1). The function
62 * computes the representation of M in another positive integer p,
63 * RE = mod(M, p) using Garner's algorithm. 'mod' represents positive modular
64 * operation.
65 *
66 * Let q be product of all elements in the RNS basis. Every entry m in M
67 * satisfies -(q-1)/2 <= m <= (q-1)/2. That is, M has both positive entries
68 * and negative entries.
69 *
70 * To avoid repeat computations, the function takes some precomputed
71 * informations as input, which are listed below.
72 *
73 * Input:
74 * len: long, dimension of RNS basis
75 * n: long, length of array RE
76 * p: FiniteField, modulus
77 * basis: 1-dim FiniteField array length len, RNS basis
78 * cmbasis: 1-dim FiniteField array length len, computed by function
79 * combBasis, inverses of special combination of RNS basis
80 * cumprod: double, computed by function cumProd,
81 * (-basis[0]*basis[1]*...*basis[len-1]) mod p
82 * bdcoeff: 1-dim FiniteField array length len, computed by function repBound
83 * R: 1-dim Double array length n*len, representation of a len x n
84 * matrix, R[i]=mod(M, basis[i]) (i=0..len-1)
85 *
86 * Output:
87 * RE: 1-dim Double array length n, RE = mod(M, p), the space of RE
88 * should be allocated before calling the function.
89 *
90 * Precondition:
91 * t <= 2^53-1, where t is the maximal intermidiate value arised in this
92 * function,
93 * t = max(2*(basis[len-1]-1)^2, (p-1)^2+basis[len-1]-1)
94 *
95 */
96
97 void
basisExt(const long len,const long n,const FiniteField p,const FiniteField * basis,const FiniteField * cmbasis,const double cumprod,const FiniteField * bdcoeff,Double ** R,Double * RE)98 basisExt (const long len, const long n, const FiniteField p, \
99 const FiniteField *basis, const FiniteField *cmbasis, \
100 const double cumprod, const FiniteField *bdcoeff, Double **R, \
101 Double *RE)
102 {
103 long i, j;
104 double temp;
105 Double **U;
106 const FiniteField *q, *qinv;
107
108 q = basis;
109 qinv = cmbasis;
110
111 /* if p = q[i] then just copy the corresponding column to RE */
112 for (i = 0; i < len; i++)
113 { if (p == q[i]) { cblas_dcopy(n, R[i], 1, RE, 1); return; } }
114 U = XMALLOC(Double *, len);
115 for (i = 0; i < len; i++) { U[i] = XMALLOC(Double, n); }
116
117 /* compute the coefficients of mix radix in positive representation by
118 inplacing modular matrix U */
119 cblas_dcopy(n, R[0], 1, U[0], 1);
120 for (i = 1; i < len; i++)
121 {
122 cblas_dcopy(n, U[i-1], 1, U[i], 1);
123 for (j = i-2; j >= 0; j--)
124 {
125 cblas_dscal(n, (double)(q[j] % q[i]), U[i], 1);
126 cblas_daxpy(n, 1.0, U[j], 1, U[i], 1);
127 Dmod((double)q[i], U[i], 1, n, n);
128 }
129 temp = (double)qinv[i]*(double)(q[i]-1);
130 temp = fmod(temp, (double)q[i]);
131 cblas_dscal(n, temp, U[i], 1);
132 cblas_daxpy(n, (double)qinv[i], R[i], 1, U[i], 1);
133 Dmod((double)q[i], U[i], 1, n, n);
134 }
135
136 /* compute mod(r, p) in positive representation and store into RE */
137 cblas_dcopy(n, U[len-1], 1, RE, 1);
138 Dmod((double)p, RE, 1, n, n);
139 for (i = len-2; i >= 0; i--)
140 {
141 cblas_dscal(n, (double)(q[i] % p), RE, 1);
142 cblas_daxpy(n, 1.0, U[i], 1, RE, 1);
143 Dmod((double)p, RE, 1, n, n);
144 }
145
146 /* convert to symmetric representation */
147 for (i = 0; i < n; i++)
148 for (j = len-1; j >= 0; j--)
149 {
150 if (U[j][i] > bdcoeff[j])
151 {
152 RE[i] = fmod(RE[i]+cumprod, (double)p);
153 break;
154 }
155 else if (U[j][i] < bdcoeff[j]) { break; }
156 }
157 for (i = 0; i < len; i++) { XFREE(U[i]); } { XFREE(U); }
158
159 return;
160 }
161
162
163 /*
164 *
165 * Calling Sequence:
166 * basisExtPos(len, n, p, basis, cmbasis, R, RE)
167 *
168 * Summary:
169 * Given a representation of a non-negative matrix/vector in some RNS,
170 * extend to compute the representation in another positive integer
171 *
172 * Description:
173 * Let R be the representation of a matrix/vector M in a residue basis
174 * 'basis', i.e., R[i] = mod(M, basis[i]) (i = 0..len-1). The function
175 * computes the representation of M in another positive integer p,
176 * RE = mod(M, p) using Garner's algorithm. 'mod' represents positive modular
177 * operation.
178 *
179 * Let q be product of all elements in the RNS basis. Every entry m in M
180 * satisfies 0 <= m <= q-1. That is, M only contains non-negative entries.
181 *
182 * To avoid repeat computations, the function takes some precomputed
183 * informations as input, which are listed below.
184 *
185 * Input:
186 * len: long, dimension of RNS basis
187 * n: long, length of array RE
188 * p: FiniteField, modulus
189 * basis: 1-dim FiniteField array length len, RNS basis
190 * cmbasis: 1-dim FiniteField array length len, computed by function
191 * combBasis, inverses of special combination of RNS basis
192 * R: 1-dim Double array length n*len, representation of a len x n
193 * matrix, R[i]=mod(M, basis[i]) (i=0..len-1)
194 *
195 * Output:
196 * RE: 1-dim Double array length n, RE = mod(M, p), the space of RE
197 * should be allocated before calling the function.
198 *
199 * Precondition:
200 * t <= 2^53-1, where t is the maximal intermidiate value arised in this
201 * function, t = max(2*(basis[len-1]-1)^2, (p-1)^2+basis[len-1]-1)
202 *
203 */
204
205 void
basisExtPos(const long len,const long n,const FiniteField p,const FiniteField * basis,const FiniteField * cmbasis,Double ** R,Double * RE)206 basisExtPos (const long len, const long n, const FiniteField p, \
207 const FiniteField *basis, const FiniteField *cmbasis, \
208 Double **R, Double *RE)
209 {
210 long i, j;
211 double temp;
212 Double **U;
213 const FiniteField *q, *qinv;
214
215 q = basis;
216 qinv = cmbasis;
217
218 /* if p = q[i] then just copy the corresponding column to RE */
219 for (i = 0; i < len; i++)
220 { if (p == q[i]) { cblas_dcopy(n, R[i], 1, RE, 1); return; } }
221 U = XMALLOC(Double *, len);
222 for (i = 0; i < len; i++) { U[i] = XMALLOC(Double, n); }
223
224 /* compute the coefficients of mix radix in positive representation by
225 inplacing modular matrix U */
226 cblas_dcopy(n, R[0], 1, U[0], 1);
227 for (i = 1; i < len; i++)
228 {
229 cblas_dcopy(n, U[i-1], 1, U[i], 1);
230 for (j = i-2; j >= 0; j--)
231 {
232 cblas_dscal(n, (double)(q[j] % q[i]), U[i], 1);
233 cblas_daxpy(n, 1.0, U[j], 1, U[i], 1);
234 Dmod((double)q[i], U[i], 1, n, n);
235 }
236 temp = (double)qinv[i]*(double)(q[i]-1);
237 temp = fmod(temp, (double)q[i]);
238 cblas_dscal(n, temp, U[i], 1);
239 cblas_daxpy(n, (double)qinv[i], R[i], 1, U[i], 1);
240 Dmod((double)q[i], U[i], 1, n, n);
241 }
242
243 /* compute mod(r, p) in positive representation and store into RE */
244 cblas_dcopy(n, U[len-1], 1, RE, 1);
245 Dmod((double)p, RE, 1, n, n);
246 for (i = len-2; i >= 0; i--)
247 {
248 cblas_dscal(n, (double)(q[i] % p), RE, 1);
249 cblas_daxpy(n, 1.0, U[i], 1, RE, 1);
250 Dmod((double)p, RE, 1, n, n);
251 }
252 for (i = 0; i < len; i++) { XFREE(U[i]); } { XFREE(U); }
253
254 return;
255 }
256
257
258
259 /*
260 *
261 * Calling Sequence:
262 * basisProd(len, basis, mp_prod)
263 *
264 * Summary:
265 * Compute the product of elements of a RNS basis
266 *
267 * Description:
268 * Let a RNS basis be 'basis'. The function computes the product of its
269 * elements basis[0]*basis[1]*...*basis[len-1].
270 *
271 * Input:
272 * len: long, dimension of RNS basis
273 * basis: 1-dim FiniteField array length len, RNS basis
274 *
275 * Output:
276 * mp_prod: mpz_t, product of elements in 'basis'
277 *
278 */
279
280 void
basisProd(const long len,const FiniteField * basis,mpz_t mp_prod)281 basisProd (const long len, const FiniteField *basis, mpz_t mp_prod)
282 {
283 long i;
284
285 mpz_set_ui(mp_prod, basis[0]);
286 for (i = 1; i < len; i++) { mpz_mul_ui(mp_prod, mp_prod, basis[i]); }
287 }
288
289
290
291 /*
292 *
293 * Calling Sequence:
294 * ChineseRemainder(len, mp_prod, basis, cmbasis, bdcoeff, Ac, mp_Ac)
295 *
296 * Summary:
297 * Given a representation of an integer in some RNS, use Chinese Remainder
298 * Algorithm to reconstruct the integer
299 *
300 * Description:
301 * Let A be an integer, and Ac contains the representation of A in a RNS
302 * basis 'basis', i.e. Ac[i] = mod(A, basis[i]), (i = 0..len). Here 'mod'
303 * is in positive representation. The function reconstructs the integer A
304 * given the RNS basis 'basis' and Ac.
305 *
306 * To avoid repeat computations, the function takes some precomputed
307 * informations as input, which are listed below.
308 *
309 * Input:
310 * len: long, dimension of RNS basis
311 * mp_prod: mpz_t, computed by function basisProd, product of RNS basis
312 * basis: 1-dim FiniteField array length len, RNS basis
313 * cmbasis: 1-dim FiniteField array length len, computed by function
314 * combBasis, inverses of special combination of RNS basis
315 * bdcoeff: 1-dim FiniteField array length len, computed by function repBound
316 * Ac: 1-dim Double array length n, representation of A in RNS
317 *
318 * Output:
319 * mp_Ac: mpz_t, reconstructed integer A
320 *
321 * Precondition:
322 * Let q be product of all elements in the RNS basis. Then A must satisfy
323 * -(q-1)/2 <= A <= (q-1)/2.
324 *
325 */
326
327 void
ChineseRemainder(const long len,const mpz_t mp_prod,const FiniteField * basis,const FiniteField * cmbasis,const FiniteField * bdcoeff,Double * Ac,mpz_t mp_Ac)328 ChineseRemainder (const long len, const mpz_t mp_prod, \
329 const FiniteField *basis, const FiniteField *cmbasis, \
330 const FiniteField *bdcoeff, Double *Ac, mpz_t mp_Ac)
331 {
332 long i, j;
333 double temp, tempq, tempqinv;
334 Double *U;
335
336 U = XMALLOC(Double, len);
337
338 /* compute the coefficients of mix radix in positive representation by
339 inplacing modular matrix U */
340 U[0] = Ac[0];
341 for (i = 1; i < len; i++)
342 {
343 U[i] = U[i-1];
344 tempq = (double)basis[i];
345 tempqinv = (double)cmbasis[i];
346 for (j = i-2; j >= 0; j--)
347 {
348 U[i] = U[j] + U[i]*fmod((double)basis[j], tempq);
349 U[i] = fmod(U[i], tempq);
350 }
351 temp = fmod(tempqinv*(double)(basis[i]-1), tempq);
352 U[i] = fmod(tempqinv*Ac[i]+temp*U[i], tempq);
353 }
354 /* compute Ac in positive representation */
355 mpz_set_d(mp_Ac, U[len-1]);
356 for (j = len-2; j >= 0; j--)
357 {
358 mpz_mul_ui(mp_Ac, mp_Ac, basis[j]);
359 mpz_add_ui(mp_Ac, mp_Ac, (FiniteField)U[j]);
360 }
361 /* transfer from positive representation to symmetric representation */
362 for (j = len-1; j >= 0; j--)
363 {
364 if (U[j] > bdcoeff[j])
365 {
366 mpz_sub(mp_Ac, mp_Ac, mp_prod);
367 break;
368 }
369 else if (U[j] < bdcoeff[j]) { break; }
370 }
371 XFREE(U);
372
373 return;
374 }
375
376
377
378 /*
379 *
380 * Calling Sequence:
381 * ChineseRemainderPos(len, basis, cmbasis, Ac, mp_Ac)
382 *
383 * Summary:
384 * Given a representation of a non-negative integer in some RNS, use Chinese
385 * Remainder Algorithm to reconstruct the integer
386 *
387 * Description:
388 * Let A be a non-negative integer, and Ac contains the representation of A
389 * in a RNS basis 'basis', i.e. Ac[i] = mod(A, basis[i]), (i = 0..len).
390 * Here 'mod' is in positive representation. The function reconstructs the
391 * integer A given the RNS basis 'basis' and Ac.
392 *
393 * To avoid repeat computations, the function takes some precomputed
394 * informations as input, which are listed below.
395 *
396 * Input:
397 * len: long, dimension of RNS basis
398 * basis: 1-dim FiniteField array length len, RNS basis
399 * cmbasis: 1-dim FiniteField array length len, computed by function
400 * combBasis, inverses of special combination of RNS basis
401 * Ac: 1-dim Double array length n, representation of A in RNS
402 *
403 * Output:
404 * mp_Ac: mpz_t, reconstructed integer A
405 *
406 * Precondition:
407 * Let q be product of all elements in the RNS basis. Then A must satisfy
408 * 0 <= A <= q-1.
409 *
410 */
411
412 void
ChineseRemainderPos(const long len,const FiniteField * basis,const FiniteField * cmbasis,Double * Ac,mpz_t mp_Ac)413 ChineseRemainderPos (const long len, const FiniteField *basis, \
414 const FiniteField *cmbasis, Double *Ac, mpz_t mp_Ac)
415 {
416 long i, j;
417 double temp, tempq, tempqinv;
418 Double *U;
419
420 U = XMALLOC(Double, len);
421
422 /* compute the coefficients of mix radix in positive representation by
423 inplacing modular matrix U */
424 U[0] = Ac[0];
425 for (i = 1; i < len; i++)
426 {
427 U[i] = U[i-1];
428 tempq = (double)basis[i];
429 tempqinv = (double)cmbasis[i];
430 for (j = i-2; j >= 0; j--)
431 {
432 U[i] = U[j] + U[i]*fmod((double)basis[j], tempq);
433 U[i] = fmod(U[i], tempq);
434 }
435 temp = fmod(tempqinv*(double)(basis[i]-1), tempq);
436 U[i] = fmod(tempqinv*Ac[i]+temp*U[i], tempq);
437
438 }
439 /* compute Ac in positive representation */
440 mpz_set_d(mp_Ac, U[len-1]);
441 for (j = len-2; j >= 0; j--)
442 {
443 mpz_mul_ui(mp_Ac, mp_Ac, basis[j]);
444 mpz_add_ui(mp_Ac, mp_Ac, (FiniteField)U[j]);
445 }
446 XFREE(U);
447
448 return;
449 }
450
451
452
453 /*
454 *
455 * Calling Sequence:
456 * cmbasis <-- combBasis(basislen, basis)
457 *
458 * Summary:
459 * Compute the special combination of a RNS basis
460 *
461 * Description:
462 * Let 'basis' be RNS basis. The function computes an array cmbasis
463 * satisfying
464 * cmbasis[0] = 0, cmbasis[i] = mod(1/(basis[0]*...*basis[i-1]), basis[i])
465 * (i = 1..basislen-1)
466 *
467 * Input:
468 * basislen: long, dimension of RNS basis
469 * basis: 1-dim FiniteField array length basislen, RNS basis
470 *
471 * Return:
472 * cmbasis: 1-dim FiniteField array length basislen, shown as above
473 *
474 */
475
476 FiniteField *
combBasis(const long basislen,const FiniteField * basis)477 combBasis (const long basislen, const FiniteField *basis)
478 {
479 long i, j;
480 double dtemp;
481 mpz_t mp_prod, mp_q;
482 FiniteField *cmbasis;
483
484 cmbasis = XMALLOC(FiniteField, basislen);
485 cmbasis[0] = 0;
486 mpz_init(mp_prod);
487 mpz_init(mp_q);
488 for (i = 1; i < basislen; i++)
489 {
490 dtemp = fmod((double)basis[0], (double)basis[i]);
491 for (j = 1; j <= i-1; j++)
492 dtemp = fmod(dtemp*(double)basis[j], (double)basis[i]);
493 mpz_set_ui(mp_q, basis[i]);
494 mpz_set_d(mp_prod, dtemp);
495 mpz_invert(mp_prod, mp_prod, mp_q);
496 cmbasis[i] = mpz_get_ui(mp_prod);
497 }
498 mpz_clear(mp_prod);
499 mpz_clear(mp_q);
500
501 return cmbasis;
502 }
503
504
505
506 /*
507 *
508 * Calling Sequence:
509 * cumprod <-- cumProd(basislen, basis, extbasislen, extbasis)
510 *
511 * Summary:
512 * Compute the representation of the combination of elements of one RNS basis
513 * in another RNS basis
514 *
515 * Description:
516 * Let 'basis' be one RNS basis with dimension basislen, and 'extbasis' be
517 * another RNS basis with dimension extbasislen. The function computes an
518 * array cumprod length extbasislen satisfying
519 * cumprod[i] = modp(-basis[0]*...*basis[basislen-1], extbasis[i]),
520 * i = 0..extbasislen-1
521 *
522 * Input:
523 * basislen: long, dimension of RNS basis 'basis'
524 * basis: 1-dim FiniteField array length basislen, one RNS basis
525 * extbasislen: long, dimension of RNS basis 'extbasis'
526 * extbasis: 1-dim FiniteField array length basislen, another RNS basis
527 *
528 * Return:
529 * cumprod: 1-dim double array length extbasislen, shown above
530 *
531 */
532
533 double *
cumProd(const long basislen,const FiniteField * basis,const long extbasislen,const FiniteField * extbasis)534 cumProd (const long basislen, const FiniteField *basis, \
535 const long extbasislen, const FiniteField *extbasis)
536 {
537 long i, j;
538 double dtemp, dextbasis;
539 double *cumprod;
540
541 cumprod = XMALLOC(double, extbasislen);
542 for (i = 0; i < extbasislen; i++)
543 {
544 dextbasis = (double)extbasis[i];
545 cumprod[i] = fmod((double)basis[0], dextbasis);
546 for (j = 1; j < basislen; j++)
547 {
548 dtemp = fmod((double)basis[j], dextbasis);
549 cumprod[i] = fmod(cumprod[i]*dtemp, dextbasis);
550 }
551 cumprod[i] = dextbasis-cumprod[i];
552 }
553
554 return cumprod;
555 }
556
557
558 /*
559 *
560 * Calling Sequence:
561 * basiscmb <-- findRNS(RNS_bound, mp_maxInter, len)
562 *
563 * Summary:
564 * Find a RNS basis and its special combination
565 *
566 * Description:
567 * Given RNS_bound, the upper bound of the RNS basis, and mp_maxInter, the
568 * function finds a best RNS basis and a combination of that basis.
569 *
570 * The RNS basis 'basis' has the property:
571 * - its elements are all primes
572 * - basis[0] is the largest prime among all the primes at most RNS_bound
573 * - basis[i+1] is the next prime smaller than basis[i] (i = 0..len-2)
574 * - basis[0]*basis[1]*...*basis[len-1] >= mp_maxInter
575 *
576 * After finding 'basis', the functions also computes the combination of
577 * 'basis' as the operations in function combBasis.
578 *
579 * Input:
580 * RNS_bound: FiniteField, the upper bound of the RNS basis
581 * mp_maxInter: mpz_t, the lower bound for the product of elements of basis
582 *
583 * Return:
584 * basiscmb: 2-dim FiniteField array, dimension 2 x len, where
585 * - basiscmb[0] represents the RNS basis
586 * - basiscmb[1] represents the special combination of basis
587 *
588 * Output:
589 * len: pointer to a long int, storing the dimension of the computed
590 * RNS basis
591 *
592 */
593
594 FiniteField **
findRNS(const FiniteField RNS_bound,const mpz_t mp_maxInter,long * length)595 findRNS (const FiniteField RNS_bound, const mpz_t mp_maxInter, long *length)
596 {
597 long i, j, len=0;
598 double prod;
599 mpz_t mp_l, mp_prod, mp_q;
600 FiniteField **qqinv;
601
602 mpz_init_set_ui(mp_prod, 1);
603 mpz_init_set_ui(mp_l, RNS_bound);
604 qqinv = XMALLOC(FiniteField *, 2);
605 qqinv[0] = NULL;
606 while (mpz_cmp(mp_maxInter, mp_prod) > 0)
607 {
608 ++len;
609 qqinv[0] = XREALLOC(FiniteField, qqinv[0], len);
610 while (mpz_probab_prime_p(mp_l, 10) == 0) { mpz_sub_ui(mp_l, mp_l, 1); }
611 qqinv[0][len-1] = mpz_get_ui(mp_l);
612 mpz_sub_ui(mp_l, mp_l, 1);
613 mpz_mul_ui(mp_prod, mp_prod, qqinv[0][len-1]);
614 }
615 mpz_clear(mp_prod);
616 mpz_clear(mp_l);
617 qqinv[1] = XMALLOC(FiniteField, len);
618 qqinv[1][0] = 0;
619 mpz_init(mp_prod);
620 mpz_init(mp_q);
621 for (i = 1; i < len; i++)
622 {
623 prod = (double)(qqinv[0][0] % qqinv[0][i]);
624 for (j = 1; j <= i-1; j++)
625 prod = fmod(prod*(double)qqinv[0][j], (double)qqinv[0][i]);
626 mpz_set_ui(mp_q, qqinv[0][i]);
627 mpz_set_d(mp_prod, prod);
628 mpz_invert(mp_prod, mp_prod, mp_q);
629 qqinv[1][i] = mpz_get_ui(mp_prod);
630 }
631 mpz_clear(mp_prod);
632 mpz_clear(mp_q);
633 *length = len;
634
635 return qqinv;
636 }
637
638
639
640 /*
641 *
642 * Calling Sequence:
643 * maxInter(mp_prod, mp_alpha, n, mp_b)
644 *
645 * Summary:
646 * Compute the maximum interval of positive and negative results of a
647 * matrix-matrix or matrix-vector product
648 *
649 * Description:
650 * Let mp_alpha be the maximum magnitude of a m x n matrix A, mp_prod-1 be
651 * the maximum magnitude of a n x k matrix C. The function computes the
652 * maximum interval of positive and negative entries of A.C. That is, the
653 * function computes mp_b satisfying
654 * (mp_b-1)/2 = n*mp_alpha*(mp_prod-1)
655 *
656 * Input:
657 * mp_prod: mpz_t, mp_prod-1 be the maximum magnitude of matrix C
658 * mp_alpha: mpz_t, maximum magnitude of matrix A
659 * n: long, column dimension of A
660 *
661 * Output:
662 * mp_b: mpz_t, shown above
663 *
664 */
665
666 void
maxInter(const mpz_t mp_prod,const mpz_t mp_alpha,const long n,mpz_t mp_b)667 maxInter (const mpz_t mp_prod, const mpz_t mp_alpha, const long n, mpz_t mp_b)
668 {
669 mpz_t mp_temp;
670
671 mpz_init(mp_temp);
672 mpz_sub_ui(mp_temp, mp_prod, 1);
673 mpz_set(mp_b, mp_alpha);
674 mpz_mul_ui(mp_b, mp_b, n);
675 mpz_mul(mp_b, mp_b, mp_temp);
676 mpz_mul_ui(mp_b, mp_b, 2);
677 mpz_add_ui(mp_b, mp_b, 1);
678 mpz_clear(mp_temp);
679 }
680
681
682
683 /*
684 *
685 * Calling Sequence:
686 * maxExtInter(mp_alpha, n, mp_b)
687 *
688 * Summary:
689 * Compute the maximum interval of positive and negative results for
690 * lifting
691 *
692 * Description:
693 * Let mp_alpha be the maximum magnitude of a m x n matrix A,
694 * The function computes the mp_b satisfying
695 * (mp_b-1)/2 = n*mp_alpha+1
696 *
697 * Input:
698 * mp_alpha: mpz_t, maximum magnitude of matrix A
699 * n: long, column dimension of A
700 *
701 * Output:
702 * mp_b: mpz_t, shown above
703 *
704 */
705
706 void
maxExtInter(const mpz_t mp_alpha,const long n,mpz_t mp_b)707 maxExtInter (const mpz_t mp_alpha, const long n, mpz_t mp_b)
708 {
709 mpz_set_ui(mp_b, 1);
710 mpz_addmul_ui(mp_b, mp_alpha, n);
711 mpz_mul_ui(mp_b, mp_b, 2);
712 mpz_add_ui(mp_b, mp_b, 1);
713 }
714
715
716
717 /*
718 *
719 * Calling Sequence:
720 * bdcoeff <-- repBound(len, basis, cmbasis)
721 *
722 * Summary:
723 * Compute the mix radix coefficients of a special integer in a RNS basis
724 *
725 * Description:
726 * Given a RNS basis, suppose the product of elements in the basis be q,
727 * then this RNS basis is able to represent integers lying in
728 * [-(q-1)/2, (q-1)/2] and [0, q-1] respectively with symmetric
729 * representation and positive representation. To transfer the result from
730 * positive representation to symmetric representation, the function
731 * computes the mix radix coefficients of the boundary value (q-1)/2 in the
732 * positive representation.
733 *
734 * Let RNS basis be P. The function computes coefficient array U, such that
735 * (q-1)/2 = U[0] + U[1]*P[0] + U[2]*P[0]*P[1] +...+ U[len-1]*P[0]*...*P[len-2]
736 *
737 * Input:
738 * len: long, dimension of RNS basis
739 * basis: 1-dim FiniteField array length len, RNS basis
740 * cmbasis: 1-dim FiniteField array length len, computed by function
741 * combBasis, inverses of special combination of RNS basis
742 *
743 * Output:
744 * bdcoeff: 1-dim FiniteField array length len, the coefficient array U above
745 *
746 */
747
748 FiniteField *
repBound(const long len,const FiniteField * basis,const FiniteField * cmbasis)749 repBound (const long len, const FiniteField *basis, const FiniteField *cmbasis)
750 {
751 long i, j;
752 double dtemp;
753 mpz_t mp_bd, mp_prod;
754 FiniteField *bdcoeff;
755 const FiniteField *q, *qinv;
756
757 q = basis;
758 qinv = cmbasis;
759
760 /* set the bound of transformation from positive to negative */
761 mpz_init(mp_prod);
762 basisProd(len, q, mp_prod);
763 mpz_init(mp_bd);
764 mpz_sub_ui(mp_bd, mp_prod, 1);
765 mpz_divexact_ui(mp_bd, mp_bd, 2);
766
767 /* compute the coeffcients of bound of mix radix and store in bdcoeff */
768 bdcoeff = XMALLOC(FiniteField, len);
769 bdcoeff[0] = mpz_fdiv_ui(mp_bd, q[0]);
770 for (i = 1; i < len; i++)
771 {
772 dtemp = (double)bdcoeff[i-1];
773 for (j = i-2; j >= 0; j--)
774 {
775 dtemp = fmod(dtemp*q[j], (double)q[i]);
776 dtemp = fmod(dtemp+bdcoeff[j], (double)q[i]);
777 }
778 bdcoeff[i] = mpz_fdiv_ui(mp_bd, q[i]);
779 dtemp = fmod((double)bdcoeff[i]-dtemp, (double)q[i]);
780 if (dtemp < 0) { dtemp = q[i]+dtemp; }
781 bdcoeff[i] = (FiniteField)fmod(dtemp*qinv[i], (double)q[i]);
782 }
783 mpz_clear(mp_bd);
784 mpz_clear(mp_prod);
785
786 return bdcoeff;
787 }
788
789
790
791 /*
792 * Calling Sequence:
793 * bd <-- RNSbound(n)
794 *
795 * Summary:
796 * Compute the upper bound of a RNS basis
797 *
798 * Description:
799 * Given a m x n mod p matrix A, and a n x k mod p matrix B, the maximum
800 * magnitude of A.B is n*(p-1)^2. To use BLAS, it is needed that
801 * n*(p-1)^2 <= 2^53-1 to make the result of product correct.
802 *
803 * The function computes an integer bd, such that
804 * n*(bd-1)^2 <= 2^53-1 and n*((bd+1)-1)^2 > 2^53-1
805 *
806 * Input:
807 * n: long, column dimension of matrix A
808 *
809 * Output:
810 * bd: FiniteField, shown above
811 *
812 */
813
814 FiniteField
RNSbound(const long n)815 RNSbound (const long n)
816 {
817 FiniteField bd;
818 mpz_t mp_n, mp_d, mp_q;
819
820 mpz_init(mp_n);
821 mpz_init_set_ui(mp_d, n);
822 mpz_init(mp_q);
823 mpz_ui_pow_ui(mp_n, 2, 53);
824 mpz_sub_ui(mp_n, mp_n, 1);
825 mpz_fdiv_q(mp_q, mp_n, mp_d);
826 mpz_sqrt(mp_q, mp_q);
827 bd = mpz_get_ui(mp_q)+1;
828 mpz_clear(mp_n);
829 mpz_clear(mp_d);
830 mpz_clear(mp_q);
831
832 return bd;
833 }
834
835
836
837
838
839
840
841