1 /* Copyright (C) 2012-2020 IBM Corp.
2 * This program is Licensed under the Apache License, Version 2.0
3 * (the "License"); you may not use this file except in compliance
4 * with the License. You may obtain a copy of the License at
5 * http://www.apache.org/licenses/LICENSE-2.0
6 * Unless required by applicable law or agreed to in writing, software
7 * distributed under the License is distributed on an "AS IS" BASIS,
8 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
9 * See the License for the specific language governing permissions and
10 * limitations under the License. See accompanying LICENSE file.
11 */
12 #include <helib/EvalMap.h>
13 #include <helib/apiAttributes.h>
14
15 // needed to get NTL's TraceMap functions...needed for ThinEvalMap
16 #include <NTL/lzz_pXFactoring.h>
17 #include <NTL/GF2XFactoring.h>
18
19 namespace helib {
20
21 // Forward declarations
22 static BlockMatMul1D* buildStep1Matrix(const EncryptedArray& ea,
23 std::shared_ptr<CubeSignature> sig,
24 const NTL::Vec<long>& reps,
25 long dim,
26 long cofactor,
27 bool invert,
28 bool normal_basis);
29 static MatMul1D* buildStep2Matrix(const EncryptedArray& ea,
30 std::shared_ptr<CubeSignature> sig,
31 const NTL::Vec<long>& reps,
32 long dim,
33 long cofactor,
34 bool invert);
35 static void init_representatives(NTL::Vec<long>& representatives,
36 long dim,
37 const NTL::Vec<long>& mvec,
38 const PAlgebra& zMStar);
39
40 // Constructor: initializing tables for the evaluation-map transformations
41
EvalMap(const EncryptedArray & _ea,bool minimal,const NTL::Vec<long> & mvec,bool _invert,bool build_cache,bool normal_basis)42 EvalMap::EvalMap(const EncryptedArray& _ea,
43 bool minimal,
44 const NTL::Vec<long>& mvec,
45 bool _invert,
46 bool build_cache,
47 bool normal_basis) :
48 ea(_ea), invert(_invert)
49 {
50 const PAlgebra& zMStar = ea.getPAlgebra();
51
52 long p = zMStar.getP();
53 long d = zMStar.getOrdP();
54
55 // FIXME: we should check that ea was initialized with
56 // G == factors[0], but this is a slight pain to check
57 // currently
58
59 // NOTE: this code is derived from a more general setting, and
60 // could certainly be greatly simplified
61
62 nfactors = mvec.length();
63 assertTrue(nfactors > 0, "Invalid argument: mvec must not be empty");
64
65 for (long i = 0; i < nfactors; i++) {
66 for (long j = i + 1; j < nfactors; j++) {
67 assertEq(NTL::GCD(mvec[i], mvec[j]),
68 1l,
69 "Invalid argument: mvec elements must be pairwise co-prime");
70 }
71 }
72
73 long m = computeProd(mvec);
74 assertEq(m,
75 (long)zMStar.getM(),
76 "Invalid argument: Product of mvec elements does not match "
77 "ea.zMStar.getM()");
78
79 NTL::Vec<long> phivec(NTL::INIT_SIZE, nfactors);
80 for (long i = 0; i < nfactors; i++)
81 phivec[i] = phi_N(mvec[i]);
82 long phim = computeProd(phivec);
83
84 NTL::Vec<long> dprodvec(NTL::INIT_SIZE, nfactors + 1);
85 dprodvec[nfactors] = 1;
86
87 for (long i = nfactors - 1; i >= 0; i--)
88 dprodvec[i] =
89 dprodvec[i + 1] *
90 multOrd(NTL::PowerMod(p % mvec[i], dprodvec[i + 1], mvec[i]), mvec[i]);
91
92 NTL::Vec<long> dvec(NTL::INIT_SIZE, nfactors);
93 for (long i = 0; i < nfactors; i++)
94 dvec[i] = dprodvec[i] / dprodvec[i + 1];
95
96 long nslots = phim / d;
97 assertEq(d, dprodvec[0], "dprodvec must start with d");
98 assertEq(nslots,
99 (long)zMStar.getNSlots(),
100 "Slot count mismatch between ea and phi(m)/d");
101
102 long inertPrefix = 0;
103 for (long i = 0; i < nfactors && dvec[i] == 1; i++) {
104 inertPrefix++;
105 }
106
107 if (inertPrefix != nfactors - 1)
108 throw LogicError("EvalMap: case not handled: bad inertPrefix");
109
110 NTL::Vec<NTL::Vec<long>> local_reps(NTL::INIT_SIZE, nfactors);
111 for (long i = 0; i < nfactors; i++)
112 init_representatives(local_reps[i], i, mvec, zMStar);
113
114 NTL::Vec<long> crtvec(NTL::INIT_SIZE, nfactors);
115 for (long i = 0; i < nfactors; i++)
116 crtvec[i] = (m / mvec[i]) * NTL::InvMod((m / mvec[i]) % mvec[i], mvec[i]);
117
118 NTL::Vec<long> redphivec(NTL::INIT_SIZE, nfactors);
119 for (long i = 0; i < nfactors; i++)
120 redphivec[i] = phivec[i] / dvec[i];
121
122 CubeSignature redphisig(redphivec);
123
124 NTL::Vec<std::shared_ptr<CubeSignature>> sig_sequence;
125 sig_sequence.SetLength(nfactors + 1);
126 sig_sequence[nfactors] = std::make_shared<CubeSignature>(phivec);
127
128 NTL::Vec<long> reduced_phivec = phivec;
129
130 for (long dim = nfactors - 1; dim >= 0; dim--) {
131 reduced_phivec[dim] /= dvec[dim];
132 sig_sequence[dim] = std::make_shared<CubeSignature>(reduced_phivec);
133 }
134
135 long dim = nfactors - 1;
136 std::unique_ptr<BlockMatMul1D> mat1_data;
137 mat1_data.reset(buildStep1Matrix(ea,
138 sig_sequence[dim],
139 local_reps[dim],
140 dim,
141 m / mvec[dim],
142 invert,
143 normal_basis));
144 mat1.reset(new BlockMatMul1DExec(*mat1_data, minimal));
145
146 matvec.SetLength(nfactors - 1);
147 for (dim = nfactors - 2; dim >= 0; --dim) {
148 std::unique_ptr<MatMul1D> mat_data;
149
150 mat_data.reset(buildStep2Matrix(ea,
151 sig_sequence[dim],
152 local_reps[dim],
153 dim,
154 m / mvec[dim],
155 invert));
156 matvec[dim].reset(new MatMul1DExec(*mat_data, minimal));
157 }
158
159 if (build_cache)
160 upgrade();
161 }
162
upgrade()163 void EvalMap::upgrade()
164 {
165 mat1->upgrade();
166 for (long i = 0; i < matvec.length(); i++)
167 matvec[i]->upgrade();
168 }
169
170 // Applying the evaluation (or its inverse) map to a ciphertext
apply(Ctxt & ctxt) const171 void EvalMap::apply(Ctxt& ctxt) const
172 {
173 if (!invert) { // forward direction
174 mat1->mul(ctxt);
175
176 for (long i = matvec.length() - 1; i >= 0; i--)
177 matvec[i]->mul(ctxt);
178 } else { // inverse transformation
179 for (long i = 0; i < matvec.length(); i++)
180 matvec[i]->mul(ctxt);
181
182 mat1->mul(ctxt);
183 }
184 }
185
init_representatives(NTL::Vec<long> & representatives,long dim,const NTL::Vec<long> & mvec,const PAlgebra & zMStar)186 static void init_representatives(NTL::Vec<long>& representatives,
187 long dim,
188 const NTL::Vec<long>& mvec,
189 const PAlgebra& zMStar)
190 {
191 assertInRange(dim,
192 0l,
193 mvec.length(),
194 "Invalid argument: dim must be between 0 and mvec.length()");
195
196 // special case
197 if (dim >= LONG(zMStar.numOfGens())) {
198 representatives.SetLength(1);
199 representatives[0] = 1;
200 return;
201 }
202
203 long m = mvec[dim];
204 long D = zMStar.OrderOf(dim);
205 long g = NTL::InvMod(zMStar.ZmStarGen(dim) % m, m);
206
207 representatives.SetLength(D);
208 for (long i = 0; i < D; i++)
209 representatives[i] = NTL::PowerMod(g, i, m);
210 }
211
212 // The callback interface for the matrix-multiplication routines.
213
214 //! \cond FALSE (make doxygen ignore these classes)
215 template <typename type>
216 class Step2Matrix : public MatMul1D_derived<type>
217 {
218 PA_INJECT(type)
219
220 const EncryptedArray& base_ea;
221 std::shared_ptr<CubeSignature> sig;
222 long dim;
223 NTL::Mat<RX> A;
224
225 public:
226 // constructor
Step2Matrix(const EncryptedArray & _ea,std::shared_ptr<CubeSignature> _sig,const NTL::Vec<long> & reps,long _dim,long cofactor,bool invert=false)227 Step2Matrix(const EncryptedArray& _ea,
228 std::shared_ptr<CubeSignature> _sig,
229 const NTL::Vec<long>& reps,
230 long _dim,
231 long cofactor,
232 bool invert = false) :
233 base_ea(_ea), sig(_sig), dim(_dim)
234 {
235 long sz = sig->getDim(dim);
236 assertEq(sz,
237 reps.length(),
238 "Invalid argument: sig->getDim(dim) must equal reps.length()");
239
240 const EncryptedArrayDerived<type>& ea = _ea.getDerived(type());
241 RBak bak;
242 bak.save();
243 _ea.getAlMod().restoreContext();
244 const RX& G = ea.getG();
245
246 NTL::Vec<RX> points(NTL::INIT_SIZE, sz);
247 for (long j = 0; j < sz; j++)
248 points[j] = RX(reps[j] * cofactor, 1) % G;
249
250 A.SetDims(sz, sz);
251 for (long j = 0; j < sz; j++)
252 A[0][j] = 1;
253
254 for (long i = 1; i < sz; i++)
255 for (long j = 0; j < sz; j++)
256 A[i][j] = (A[i - 1][j] * points[j]) % G;
257
258 if (invert) {
259 REBak ebak;
260 ebak.save();
261 ea.restoreContextForG();
262
263 mat_RE A1, A2;
264 conv(A1, A);
265
266 long p = _ea.getAlMod().getZMStar().getP();
267 long r = _ea.getAlMod().getR();
268
269 ppInvert(A2, A1, p, r);
270 conv(A, A2);
271 }
272 }
273
get(RX & out,long i,long j,UNUSED long k) const274 bool get(RX& out, long i, long j, UNUSED long k) const override
275 {
276 out = A[i][j];
277 return false;
278 }
279
getEA() const280 const EncryptedArray& getEA() const override { return base_ea; }
multipleTransforms() const281 bool multipleTransforms() const override { return false; }
getDim() const282 long getDim() const override { return dim; }
283 };
284
buildStep2Matrix(const EncryptedArray & ea,std::shared_ptr<CubeSignature> sig,const NTL::Vec<long> & reps,long dim,long cofactor,bool invert)285 static MatMul1D* buildStep2Matrix(const EncryptedArray& ea,
286 std::shared_ptr<CubeSignature> sig,
287 const NTL::Vec<long>& reps,
288 long dim,
289 long cofactor,
290 bool invert)
291 {
292 switch (ea.getTag()) {
293 case PA_GF2_tag:
294 return new Step2Matrix<PA_GF2>(ea, sig, reps, dim, cofactor, invert);
295
296 case PA_zz_p_tag:
297 return new Step2Matrix<PA_zz_p>(ea, sig, reps, dim, cofactor, invert);
298
299 default:
300 return 0;
301 }
302 }
303
304 template <typename type>
305 class Step1Matrix : public BlockMatMul1D_derived<type>
306 {
307 PA_INJECT(type)
308
309 const EncryptedArray& base_ea;
310 std::shared_ptr<CubeSignature> sig;
311 long dim;
312 NTL::Mat<mat_R> A;
313
314 public:
315 // constructor
Step1Matrix(const EncryptedArray & _ea,std::shared_ptr<CubeSignature> _sig,const NTL::Vec<long> & reps,long _dim,long cofactor,bool invert,bool normal_basis)316 Step1Matrix(const EncryptedArray& _ea,
317 std::shared_ptr<CubeSignature> _sig,
318 const NTL::Vec<long>& reps,
319 long _dim,
320 long cofactor,
321 bool invert,
322 bool normal_basis) :
323 base_ea(_ea), sig(_sig), dim(_dim)
324 {
325 const EncryptedArrayDerived<type>& ea = _ea.getDerived(type());
326 RBak bak;
327 bak.save();
328 _ea.getAlMod().restoreContext();
329 const RX& G = ea.getG();
330 long d = deg(G);
331
332 long sz = sig->getDim(dim);
333 assertEq(sz,
334 reps.length(),
335 "Invalid argument: sig->getDim(dim) must equal reps.length()");
336 assertEq(dim,
337 sig->getNumDims() - 1,
338 "Invalid argument: dim must be one less than sig->getNumDims()");
339 assertEq(sig->getSize(), ea.size(), "sig and ea do not have matching size");
340
341 // so sz == phi(m_last)/d, where d = deg(G) = order of p mod m
342
343 NTL::Vec<RX> points(NTL::INIT_SIZE, sz);
344 for (long j = 0; j < sz; j++)
345 points[j] = RX(reps[j] * cofactor, 1) % G;
346
347 NTL::Mat<RX> AA(NTL::INIT_SIZE, sz * d, sz);
348 for (long j = 0; j < sz; j++)
349 AA[0][j] = 1;
350
351 for (long i = 1; i < sz * d; i++)
352 for (long j = 0; j < sz; j++)
353 AA[i][j] = (AA[i - 1][j] * points[j]) % G;
354
355 A.SetDims(sz, sz);
356 for (long i = 0; i < sz; i++)
357 for (long j = 0; j < sz; j++) {
358 A[i][j].SetDims(d, d);
359 for (long k = 0; k < d; k++)
360 VectorCopy(A[i][j][k], AA[i * d + k][j], d);
361 }
362
363 if (invert) {
364 mat_R A1, A2;
365 A1.SetDims(sz * d, sz * d);
366 for (long i = 0; i < sz * d; i++)
367 for (long j = 0; j < sz * d; j++)
368 A1[i][j] = A[i / d][j / d][i % d][j % d];
369
370 long p = _ea.getAlMod().getZMStar().getP();
371 long r = _ea.getAlMod().getR();
372
373 ppInvert(A2, A1, p, r);
374
375 for (long i = 0; i < sz * d; i++)
376 for (long j = 0; j < sz * d; j++)
377 A[i / d][j / d][i % d][j % d] = A2[i][j];
378
379 if (normal_basis) {
380 const NTL::Mat<R>& CB = ea.getNormalBasisMatrix();
381
382 // multiply each entry of A on the right by CB
383 for (long i = 0; i < sz; i++)
384 for (long j = 0; j < sz; j++)
385 A[i][j] = A[i][j] * CB;
386 } // if (normal_basis)
387 } // if (invert)
388 } // constructor
389
get(mat_R & out,long i,long j,UNUSED long k) const390 bool get(mat_R& out, long i, long j, UNUSED long k) const override
391 {
392 out = A[i][j];
393 return false;
394 }
395
getEA() const396 const EncryptedArray& getEA() const override { return base_ea; }
multipleTransforms() const397 bool multipleTransforms() const override { return false; }
getDim() const398 long getDim() const override { return dim; }
399 };
400
buildStep1Matrix(const EncryptedArray & ea,std::shared_ptr<CubeSignature> sig,const NTL::Vec<long> & reps,long dim,long cofactor,bool invert,bool normal_basis)401 static BlockMatMul1D* buildStep1Matrix(const EncryptedArray& ea,
402 std::shared_ptr<CubeSignature> sig,
403 const NTL::Vec<long>& reps,
404 long dim,
405 long cofactor,
406 bool invert,
407 bool normal_basis)
408 {
409 switch (ea.getTag()) {
410 case PA_GF2_tag:
411 return new Step1Matrix<PA_GF2>(ea,
412 sig,
413 reps,
414 dim,
415 cofactor,
416 invert,
417 normal_basis);
418
419 case PA_zz_p_tag:
420 return new Step1Matrix<PA_zz_p>(ea,
421 sig,
422 reps,
423 dim,
424 cofactor,
425 invert,
426 normal_basis);
427
428 default:
429 return 0;
430 }
431 }
432 //! \endcond
433
434 //=============== ThinEvalMap stuff
435
436 // needed to make generic programming work
437
RelaxedInv(NTL::Mat<NTL::zz_p> & x,const NTL::Mat<NTL::zz_p> & a)438 void RelaxedInv(NTL::Mat<NTL::zz_p>& x, const NTL::Mat<NTL::zz_p>& a)
439 {
440 relaxed_inv(x, a);
441 }
442
RelaxedInv(NTL::Mat<NTL::GF2> & x,const NTL::Mat<NTL::GF2> & a)443 void RelaxedInv(NTL::Mat<NTL::GF2>& x, const NTL::Mat<NTL::GF2>& a)
444 {
445 inv(x, a);
446 }
447
TraceMap(NTL::GF2X & w,const NTL::GF2X & a,long d,const NTL::GF2XModulus & F,const NTL::GF2X & b)448 void TraceMap(NTL::GF2X& w,
449 const NTL::GF2X& a,
450 long d,
451 const NTL::GF2XModulus& F,
452 const NTL::GF2X& b)
453
454 {
455 if (d < 0)
456 throw InvalidArgument("TraceMap: d is negative");
457
458 NTL::GF2X y, z, t;
459
460 z = b;
461 y = a;
462 clear(w);
463
464 while (d) {
465 if (d == 1) {
466 if (IsZero(w))
467 w = y;
468 else {
469 CompMod(w, w, z, F);
470 add(w, w, y);
471 }
472 } else if ((d & 1) == 0) {
473 Comp2Mod(z, t, z, y, z, F);
474 add(y, t, y);
475 } else if (IsZero(w)) {
476 w = y;
477 Comp2Mod(z, t, z, y, z, F);
478 add(y, t, y);
479 } else {
480 Comp3Mod(z, t, w, z, y, w, z, F);
481 add(w, w, y);
482 add(y, t, y);
483 }
484
485 d = d >> 1;
486 }
487 }
488
489 // Forward declarations
490 static MatMul1D* buildThinStep1Matrix(const EncryptedArray& ea,
491 std::shared_ptr<CubeSignature> sig,
492 const NTL::Vec<long>& reps,
493 long dim,
494 long cofactor);
495 static MatMul1D* buildThinStep2Matrix(const EncryptedArray& ea,
496 std::shared_ptr<CubeSignature> sig,
497 const NTL::Vec<long>& reps,
498 long dim,
499 long cofactor,
500 bool invert,
501 bool inflate = false);
502 static void init_representatives(NTL::Vec<long>& representatives,
503 long dim,
504 const NTL::Vec<long>& mvec,
505 const PAlgebra& zMStar);
506
507 // Constructor: initializing tables for the evaluation-map transformations
508
ThinEvalMap(const EncryptedArray & _ea,bool minimal,const NTL::Vec<long> & mvec,bool _invert,bool build_cache)509 ThinEvalMap::ThinEvalMap(const EncryptedArray& _ea,
510 bool minimal,
511 const NTL::Vec<long>& mvec,
512 bool _invert,
513 bool build_cache) :
514 ea(_ea), invert(_invert)
515 {
516 const PAlgebra& zMStar = ea.getPAlgebra();
517
518 long p = zMStar.getP();
519 long d = zMStar.getOrdP();
520 long sz = zMStar.numOfGens();
521
522 // FIXME: we should check that ea was initialized with
523 // G == factors[0], but this is a slight pain to check
524 // currently
525
526 // NOTE: this code is derived from a more general setting, and
527 // could certainly be greatly simplified
528
529 nfactors = mvec.length();
530 assertTrue(nfactors > 0, "Invalid argument: mvec must have positive length");
531
532 for (long i = 0; i < nfactors; i++) {
533 for (long j = i + 1; j < nfactors; j++) {
534 assertEq(NTL::GCD(mvec[i], mvec[j]),
535 1l,
536 "Invalid argument: mvec must have pairwise-disjoint entries");
537 }
538 }
539
540 long m = computeProd(mvec);
541 assertEq(m,
542 (long)zMStar.getM(),
543 "Invalid argument: mvec's product does not match ea's m");
544
545 NTL::Vec<long> phivec(NTL::INIT_SIZE, nfactors);
546 for (long i = 0; i < nfactors; i++)
547 phivec[i] = phi_N(mvec[i]);
548 long phim = computeProd(phivec);
549
550 NTL::Vec<long> dprodvec(NTL::INIT_SIZE, nfactors + 1);
551 dprodvec[nfactors] = 1;
552
553 for (long i = nfactors - 1; i >= 0; i--)
554 dprodvec[i] =
555 dprodvec[i + 1] *
556 multOrd(NTL::PowerMod(p % mvec[i], dprodvec[i + 1], mvec[i]), mvec[i]);
557
558 NTL::Vec<long> dvec(NTL::INIT_SIZE, nfactors);
559 for (long i = 0; i < nfactors; i++)
560 dvec[i] = dprodvec[i] / dprodvec[i + 1];
561
562 long nslots = phim / d;
563 assertEq(d, dprodvec[0], "d must match the first entry of dprodvec");
564 assertEq(nslots,
565 (long)zMStar.getNSlots(),
566 "Invalid argument: mismatch of number of slots");
567
568 long inertPrefix = 0;
569 for (long i = 0; i < nfactors && dvec[i] == 1; i++) {
570 inertPrefix++;
571 }
572
573 if (inertPrefix != nfactors - 1)
574 throw LogicError("ThinEvalMap: case not handled: bad inertPrefix");
575
576 NTL::Vec<NTL::Vec<long>> local_reps(NTL::INIT_SIZE, nfactors);
577 for (long i = 0; i < nfactors; i++)
578 init_representatives(local_reps[i], i, mvec, zMStar);
579
580 NTL::Vec<long> crtvec(NTL::INIT_SIZE, nfactors);
581 for (long i = 0; i < nfactors; i++)
582 crtvec[i] = (m / mvec[i]) * NTL::InvMod((m / mvec[i]) % mvec[i], mvec[i]);
583
584 NTL::Vec<long> redphivec(NTL::INIT_SIZE, nfactors);
585 for (long i = 0; i < nfactors; i++)
586 redphivec[i] = phivec[i] / dvec[i];
587
588 CubeSignature redphisig(redphivec);
589
590 NTL::Vec<std::shared_ptr<CubeSignature>> sig_sequence;
591 sig_sequence.SetLength(nfactors + 1);
592 sig_sequence[nfactors] = std::make_shared<CubeSignature>(phivec);
593
594 NTL::Vec<long> reduced_phivec = phivec;
595
596 for (long dim = nfactors - 1; dim >= 0; dim--) {
597 reduced_phivec[dim] /= dvec[dim];
598 sig_sequence[dim] = std::make_shared<CubeSignature>(reduced_phivec);
599 }
600
601 matvec.SetLength(nfactors);
602
603 if (invert) {
604 long dim = nfactors - 1;
605 std::unique_ptr<MatMul1D> mat1_data;
606 mat1_data.reset(buildThinStep1Matrix(ea,
607 sig_sequence[dim],
608 local_reps[dim],
609 dim,
610 m / mvec[dim]));
611 matvec[dim].reset(new MatMul1DExec(*mat1_data, minimal));
612 } else if (sz == nfactors) {
613 long dim = nfactors - 1;
614 std::unique_ptr<MatMul1D> mat1_data;
615 mat1_data.reset(buildThinStep2Matrix(ea,
616 sig_sequence[dim],
617 local_reps[dim],
618 dim,
619 m / mvec[dim],
620 invert,
621 /*inflate=*/true));
622 matvec[dim].reset(new MatMul1DExec(*mat1_data, minimal));
623 }
624
625 for (long dim = nfactors - 2; dim >= 0; --dim) {
626 std::unique_ptr<MatMul1D> mat_data;
627
628 mat_data.reset(buildThinStep2Matrix(ea,
629 sig_sequence[dim],
630 local_reps[dim],
631 dim,
632 m / mvec[dim],
633 invert));
634 matvec[dim].reset(new MatMul1DExec(*mat_data, minimal));
635 }
636
637 if (build_cache)
638 upgrade();
639 }
640
upgrade()641 void ThinEvalMap::upgrade()
642 {
643 for (long i = 0; i < matvec.length(); i++)
644 if (matvec[i])
645 matvec[i]->upgrade();
646 }
647
648 // Applying the evaluation (or its inverse) map to a ciphertext
apply(Ctxt & ctxt) const649 void ThinEvalMap::apply(Ctxt& ctxt) const
650 {
651 if (!invert) { // forward direction
652 for (long i = matvec.length() - 1; i >= 0; i--)
653 if (matvec[i])
654 matvec[i]->mul(ctxt);
655 } else { // inverse transformation
656 for (long i = 0; i < matvec.length(); i++)
657 matvec[i]->mul(ctxt);
658 traceMap(ctxt);
659 }
660 }
661
662 // The callback interface for the matrix-multiplication routines.
663
664 //! \cond FALSE (make doxygen ignore these classes)
665 template <typename type>
666 class ThinStep2Matrix : public MatMul1D_derived<type>
667 {
668 PA_INJECT(type)
669
670 const EncryptedArray& base_ea;
671 std::shared_ptr<CubeSignature> sig;
672 long dim;
673 NTL::Mat<RX> A;
674
675 public:
676 // constructor
ThinStep2Matrix(const EncryptedArray & _ea,std::shared_ptr<CubeSignature> _sig,const NTL::Vec<long> & reps,long _dim,long cofactor,bool invert,bool inflate)677 ThinStep2Matrix(const EncryptedArray& _ea,
678 std::shared_ptr<CubeSignature> _sig,
679 const NTL::Vec<long>& reps,
680 long _dim,
681 long cofactor,
682 bool invert,
683 bool inflate) :
684 base_ea(_ea), sig(_sig), dim(_dim)
685 {
686 long sz = sig->getDim(dim);
687 assertEq(sz,
688 reps.length(),
689 "Invalid argument: sig and reps have inconsistent "
690 "dimension");
691
692 const EncryptedArrayDerived<type>& ea = _ea.getDerived(type());
693 RBak bak;
694 bak.save();
695 _ea.getAlMod().restoreContext();
696 const RX& G = ea.getG();
697 long d = deg(G);
698
699 NTL::Vec<RX> points(NTL::INIT_SIZE, sz);
700 for (long j = 0; j < sz; j++) {
701 points[j] = RX(reps[j] * cofactor, 1) % G;
702 if (inflate)
703 points[j] = NTL::PowerMod(points[j], d, G);
704 }
705
706 A.SetDims(sz, sz);
707 for (long j = 0; j < sz; j++)
708 A[0][j] = 1;
709
710 for (long i = 1; i < sz; i++)
711 for (long j = 0; j < sz; j++)
712 A[i][j] = (A[i - 1][j] * points[j]) % G;
713
714 if (invert) {
715 REBak ebak;
716 ebak.save();
717 ea.restoreContextForG();
718
719 mat_RE A1, A2;
720 conv(A1, A);
721
722 long p = _ea.getAlMod().getZMStar().getP();
723 long r = _ea.getAlMod().getR();
724
725 ppInvert(A2, A1, p, r);
726 conv(A, A2);
727 }
728 }
729
get(RX & out,long i,long j,UNUSED long k) const730 bool get(RX& out, long i, long j, UNUSED long k) const override
731 {
732 out = A[i][j];
733 return false;
734 }
735
getEA() const736 const EncryptedArray& getEA() const override { return base_ea; }
multipleTransforms() const737 bool multipleTransforms() const override { return false; }
getDim() const738 long getDim() const override { return dim; }
739 };
740
buildThinStep2Matrix(const EncryptedArray & ea,std::shared_ptr<CubeSignature> sig,const NTL::Vec<long> & reps,long dim,long cofactor,bool invert,bool inflate)741 static MatMul1D* buildThinStep2Matrix(const EncryptedArray& ea,
742 std::shared_ptr<CubeSignature> sig,
743 const NTL::Vec<long>& reps,
744 long dim,
745 long cofactor,
746 bool invert,
747 bool inflate)
748 {
749 switch (ea.getTag()) {
750 case PA_GF2_tag:
751 return new ThinStep2Matrix<PA_GF2>(ea,
752 sig,
753 reps,
754 dim,
755 cofactor,
756 invert,
757 inflate);
758
759 case PA_zz_p_tag:
760 return new ThinStep2Matrix<PA_zz_p>(ea,
761 sig,
762 reps,
763 dim,
764 cofactor,
765 invert,
766 inflate);
767
768 default:
769 return 0;
770 }
771 }
772
773 template <typename type>
774 class ThinStep1Matrix : public MatMul1D_derived<type>
775 {
776 PA_INJECT(type)
777
778 const EncryptedArray& base_ea;
779 std::shared_ptr<CubeSignature> sig;
780 long dim;
781 NTL::Mat<RX> A_deflated;
782
783 public:
784 // constructor
ThinStep1Matrix(const EncryptedArray & _ea,std::shared_ptr<CubeSignature> _sig,const NTL::Vec<long> & reps,long _dim,long cofactor)785 ThinStep1Matrix(const EncryptedArray& _ea,
786 std::shared_ptr<CubeSignature> _sig,
787 const NTL::Vec<long>& reps,
788 long _dim,
789 long cofactor) :
790 base_ea(_ea), sig(_sig), dim(_dim)
791 {
792 const EncryptedArrayDerived<type>& ea = _ea.getDerived(type());
793 RBak bak;
794 bak.save();
795 _ea.getAlMod().restoreContext();
796 const RXModulus G(ea.getG());
797 long d = deg(G);
798
799 long p = _ea.getAlMod().getZMStar().getP();
800 long r = _ea.getAlMod().getR();
801
802 long sz = sig->getDim(dim);
803 assertEq(sz,
804 reps.length(),
805 "Invalid argument: sig and reps have inconsistent "
806 "dimension");
807 assertEq(dim,
808 sig->getNumDims() - 1,
809 "Invalid argument: dim must be one less than "
810 "sig->getNumDims()");
811 assertEq(sig->getSize(), ea.size(), "sig and ea do not have matching size");
812
813 // so sz == phi(m_last)/d, where d = deg(G) = order of p mod m
814
815 NTL::Vec<RX> points(NTL::INIT_SIZE, sz);
816 for (long j = 0; j < sz; j++)
817 points[j] = RX(reps[j] * cofactor, 1) % G;
818
819 NTL::Mat<RX> AA(NTL::INIT_SIZE, sz * d, sz);
820 for (long j = 0; j < sz; j++)
821 AA[0][j] = 1;
822
823 for (long i = 1; i < sz * d; i++)
824 for (long j = 0; j < sz; j++)
825 AA[i][j] = (AA[i - 1][j] * points[j]) % G;
826
827 NTL::Mat<mat_R> A;
828 A.SetDims(sz, sz);
829 for (long i = 0; i < sz; i++)
830 for (long j = 0; j < sz; j++) {
831 A[i][j].SetDims(d, d);
832 for (long k = 0; k < d; k++)
833 VectorCopy(A[i][j][k], AA[i * d + k][j], d);
834 }
835
836 // if (invert) {
837 // NOTE: this version is only used for the inverse matrix
838 mat_R A1, A2;
839 A1.SetDims(sz * d, sz * d);
840 for (long i = 0; i < sz * d; i++)
841 for (long j = 0; j < sz * d; j++)
842 A1[i][j] = A[i / d][j / d][i % d][j % d];
843
844 ppInvert(A2, A1, p, r);
845
846 for (long i = 0; i < sz * d; i++)
847 for (long j = 0; j < sz * d; j++)
848 A[i / d][j / d][i % d][j % d] = A2[i][j];
849 // }
850
851 A_deflated.SetDims(sz, sz);
852 vec_R v, w;
853 v.SetLength(d);
854 w.SetLength(d);
855
856 RX h; // set h = X^p mod G
857 PowerXMod(h, p, G);
858
859 NTL::Vec<R> trace_vec;
860 trace_vec.SetLength(2 * d - 1);
861 // set trace_vec[i] = trace(X^i mod G)
862 for (long i = 0; i < 2 * d - 1; i++) {
863 RX trace_val;
864 TraceMap(trace_val, (RX(i, 1) % G), d, G, h);
865 assertTrue(deg(trace_val) <= 0, "trace_val is positive");
866 trace_vec[i] = ConstTerm(trace_val);
867 }
868
869 NTL::Mat<R> trace_mat;
870 trace_mat.SetDims(d, d);
871 // set trace_mat[i][j] = trace(X^{i+j} mod G)
872 for (long i = 0; i < d; i++)
873 for (long j = 0; j < d; j++)
874 trace_mat[i][j] = trace_vec[i + j];
875
876 NTL::Mat<R> trace_mat_inv;
877 RelaxedInv(trace_mat_inv, trace_mat);
878
879 for (long i = 0; i < sz; i++)
880 for (long j = 0; j < sz; j++) {
881 for (long i1 = 0; i1 < d; i1++)
882 v[i1] = A[i][j][i1][0];
883 mul(w, v, trace_mat_inv);
884 conv(A_deflated[i][j], w);
885 }
886 } // constructor
887
get(RX & out,long i,long j,UNUSED long k) const888 bool get(RX& out, long i, long j, UNUSED long k) const override
889 {
890 out = A_deflated[i][j];
891 return false;
892 }
893
getEA() const894 const EncryptedArray& getEA() const override { return base_ea; }
multipleTransforms() const895 bool multipleTransforms() const override { return false; }
getDim() const896 long getDim() const override { return dim; }
897 };
898
buildThinStep1Matrix(const EncryptedArray & ea,std::shared_ptr<CubeSignature> sig,const NTL::Vec<long> & reps,long dim,long cofactor)899 static MatMul1D* buildThinStep1Matrix(const EncryptedArray& ea,
900 std::shared_ptr<CubeSignature> sig,
901 const NTL::Vec<long>& reps,
902 long dim,
903 long cofactor)
904 {
905 switch (ea.getTag()) {
906 case PA_GF2_tag:
907 return new ThinStep1Matrix<PA_GF2>(ea, sig, reps, dim, cofactor);
908
909 case PA_zz_p_tag:
910 return new ThinStep1Matrix<PA_zz_p>(ea, sig, reps, dim, cofactor);
911
912 default:
913 return 0;
914 }
915 }
916 //! \endcond
917
918 } // namespace helib
919