dcrtpoly.cpp 62.3 KB
Newer Older
Gerard Ryan's avatar
Gerard Ryan committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
/*
 * @file dcrtpoly.cpp - implementation of the integer lattice using double-CRT representations.
 * @author  TPOC: palisade@njit.edu
 *
 * @copyright Copyright (c) 2017, New Jersey Institute of Technology (NJIT)
 * All rights reserved.
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 * list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice, this
 * list of conditions and the following disclaimer in the documentation and/or other
 * materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
 * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN
 * IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */

#include "dcrtpoly.h"
#include <fstream>
#include <memory>
using std::shared_ptr;
using std::string;
#include "../utils/serializablehelper.h"
#include "../utils/debug.h"

namespace lbcrypto
{

/*CONSTRUCTORS*/
Gerard Ryan's avatar
Gerard Ryan committed
39 40
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl()
Gerard Ryan's avatar
Gerard Ryan committed
41 42
{
	m_format = EVALUATION;
Gerard Ryan's avatar
Gerard Ryan committed
43
	m_params.reset( new DCRTPolyImpl::Params(0,1) );
Gerard Ryan's avatar
Gerard Ryan committed
44 45
}

Gerard Ryan's avatar
Gerard Ryan committed
46 47
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl(const shared_ptr<DCRTPolyImpl::Params> dcrtParams, Format format, bool initializeElementToZero)
Gerard Ryan's avatar
Gerard Ryan committed
48 49 50 51 52 53 54 55 56 57 58 59
{
	m_format = format;
	m_params = dcrtParams;

	size_t vecSize = dcrtParams->GetParams().size();
	m_vectors.reserve(vecSize);

	for (usint i = 0; i < vecSize; i++) {
		m_vectors.push_back(std::move(PolyType(dcrtParams->GetParams()[i],format,initializeElementToZero)));
	}
}

Gerard Ryan's avatar
Gerard Ryan committed
60 61
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl(const DCRTPolyImpl &element)
Gerard Ryan's avatar
Gerard Ryan committed
62 63 64 65 66 67
{
	m_format = element.m_format;
	m_vectors = element.m_vectors;
	m_params = element.m_params;
}

Gerard Ryan's avatar
Gerard Ryan committed
68 69 70
template<typename VecType>
const DCRTPolyImpl<VecType>&
DCRTPolyImpl<VecType>::operator=(const PolyLargeType &element)
Gerard Ryan's avatar
Gerard Ryan committed
71 72
{

Gerard Ryan's avatar
Gerard Ryan committed
73
	if( element.GetModulus() > m_params->GetModulus() ) {
Gerard Ryan's avatar
Gerard Ryan committed
74 75 76
		throw std::logic_error("Modulus of element passed to constructor is bigger that DCRT big modulus");
	}

Gerard Ryan's avatar
Gerard Ryan committed
77 78
	size_t vecCount = m_params->GetParams().size();
	m_vectors.clear();
Gerard Ryan's avatar
Gerard Ryan committed
79 80 81 82
	m_vectors.reserve(vecCount);

	// fill up with vectors with the proper moduli
	for(usint i = 0; i < vecCount; i++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
83
		PolyType newvec(m_params->GetParams()[i], m_format, true);
Gerard Ryan's avatar
Gerard Ryan committed
84 85 86 87
		m_vectors.push_back( std::move(newvec) );
	}

	// need big ints out of the little ints for the modulo operations, below
Gerard Ryan's avatar
Gerard Ryan committed
88
	std::vector<Integer> bigmods;
Gerard Ryan's avatar
Gerard Ryan committed
89 90
	bigmods.reserve(vecCount);
	for( usint i = 0; i < vecCount; i++ )
Gerard Ryan's avatar
Gerard Ryan committed
91
		bigmods.push_back( Integer(m_params->GetParams()[i]->GetModulus().ConvertToInt()) );
Gerard Ryan's avatar
Gerard Ryan committed
92 93 94 95 96

	// copy each coefficient mod the new modulus
	for(usint p = 0; p < element.GetLength(); p++ ) {
		for( usint v = 0; v < vecCount; v++ ) {

Gerard Ryan's avatar
Gerard Ryan committed
97
			Integer tmp = element.at(p) % bigmods[v];
Gerard Ryan's avatar
Gerard Ryan committed
98
			m_vectors[v].at(p)= tmp.ConvertToInt();
Gerard Ryan's avatar
Gerard Ryan committed
99 100
		}
	}
Gerard Ryan's avatar
Gerard Ryan committed
101 102

	return *this;
Gerard Ryan's avatar
Gerard Ryan committed
103 104 105
}

/* Construct from a single Poly. The format is derived from the passed in Poly.*/
Gerard Ryan's avatar
Gerard Ryan committed
106 107
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl(const PolyLargeType &element, const shared_ptr<DCRTPolyImpl::Params> params)
Gerard Ryan's avatar
Gerard Ryan committed
108 109 110 111 112
{
	Format format;
	try {
		format = element.GetFormat();
	} catch (const std::exception& e) {
Gerard Ryan's avatar
Gerard Ryan committed
113
		throw std::logic_error("There is an issue with the format of the Poly passed to the constructor of DCRTPolyImpl");
Gerard Ryan's avatar
Gerard Ryan committed
114 115 116 117 118 119 120 121
	}

	if( element.GetCyclotomicOrder() != params->GetCyclotomicOrder() )
		throw std::logic_error("Cyclotomic order mismatch on input vector and parameters");

	m_format = format;
	m_params = params;

Gerard Ryan's avatar
Gerard Ryan committed
122
	*this = element;
Gerard Ryan's avatar
Gerard Ryan committed
123 124 125 126 127
}

/* Construct using a tower of vectors.
 * The params and format for the DCRTPolyImpl will be derived from the towers
 */
Gerard Ryan's avatar
Gerard Ryan committed
128 129
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl(const std::vector<PolyType> &towers)
Gerard Ryan's avatar
Gerard Ryan committed
130 131
{
	usint cyclotomicOrder = towers.at(0).GetCyclotomicOrder();
Gerard Ryan's avatar
Gerard Ryan committed
132
	std::vector<std::shared_ptr<ILNativeParams>> parms;
Gerard Ryan's avatar
Gerard Ryan committed
133 134 135 136 137 138 139
	for (usint i = 0; i < towers.size(); i++) {
		if ( towers[i].GetCyclotomicOrder() != cyclotomicOrder ) {
			throw std::logic_error("Polys provided to constructor must have the same ring dimension");
		}
		parms.push_back( towers[i].GetParams() );
	}

Gerard Ryan's avatar
Gerard Ryan committed
140
	shared_ptr<DCRTPolyImpl::Params> p( new DCRTPolyImpl::Params(cyclotomicOrder, parms) );
Gerard Ryan's avatar
Gerard Ryan committed
141 142 143 144 145 146 147

	m_params = p;
	m_vectors = towers;
	m_format = m_vectors[0].GetFormat();
}

/*The dgg will be the seed to populate the towers of the DCRTPolyImpl with random numbers. The algorithm to populate the towers can be seen below.*/
Gerard Ryan's avatar
Gerard Ryan committed
148 149
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl(const DggType& dgg, const shared_ptr<DCRTPolyImpl::Params> dcrtParams, Format format)
Gerard Ryan's avatar
Gerard Ryan committed
150 151 152 153 154 155 156 157
{
	m_format = format;
	m_params = dcrtParams;

	size_t vecSize = dcrtParams->GetParams().size();
	m_vectors.reserve(vecSize);

	//dgg generating random values
Gerard Ryan's avatar
Gerard Ryan committed
158
	std::shared_ptr<int32_t> dggValues = dgg.GenerateIntVector(dcrtParams->GetRingDimension());
Gerard Ryan's avatar
Gerard Ryan committed
159 160 161

	for(usint i = 0; i < vecSize; i++) {

Gerard Ryan's avatar
Gerard Ryan committed
162
		NativeVector ilDggValues(dcrtParams->GetRingDimension(), dcrtParams->GetParams()[i]->GetModulus());
Gerard Ryan's avatar
Gerard Ryan committed
163 164 165 166 167 168 169 170 171 172 173 174 175

		for(usint j = 0; j < dcrtParams->GetRingDimension(); j++) {
			uint64_t	entry;
			// if the random generated value is less than zero, then multiply it by (-1) and subtract the modulus of the current tower to set the coefficient
			int64_t k = (dggValues.get())[j];
			if(k < 0) {
				k *= (-1);
				entry = (uint64_t)dcrtParams->GetParams()[i]->GetModulus().ConvertToInt() - (uint64_t)k;
			}
			//if greater than or equal to zero, set it the value generated
			else {
				entry = k;
			}
Gerard Ryan's avatar
Gerard Ryan committed
176
			ilDggValues.at(j)=entry;
Gerard Ryan's avatar
Gerard Ryan committed
177 178 179 180 181 182 183 184 185 186 187
		}

		PolyType ilvector(dcrtParams->GetParams()[i]);
		ilvector.SetValues(ilDggValues, Format::COEFFICIENT); // the random values are set in coefficient format
		if(m_format == Format::EVALUATION) {  // if the input format is evaluation, then once random values are set in coefficient format, switch the format to achieve what the caller asked for.
			ilvector.SwitchFormat();
		}
		m_vectors.push_back(ilvector);
	}
}

Gerard Ryan's avatar
Gerard Ryan committed
188 189
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl(DugType& dug, const shared_ptr<DCRTPolyImpl::Params> dcrtParams, Format format)
Gerard Ryan's avatar
Gerard Ryan committed
190 191 192 193 194 195 196 197 198 199 200
{

	m_format = format;
	m_params = dcrtParams;

	size_t numberOfTowers = dcrtParams->GetParams().size();
	m_vectors.reserve(numberOfTowers);

	for (usint i = 0; i < numberOfTowers; i++) {

		dug.SetModulus(dcrtParams->GetParams()[i]->GetModulus());
Gerard Ryan's avatar
Gerard Ryan committed
201
		NativeVector vals(dug.GenerateVector(dcrtParams->GetRingDimension()));
Gerard Ryan's avatar
Gerard Ryan committed
202 203 204 205 206 207 208 209 210 211
		PolyType ilvector(dcrtParams->GetParams()[i]);

		ilvector.SetValues(vals, Format::COEFFICIENT); // the random values are set in coefficient format
		if (m_format == Format::EVALUATION) {  // if the input format is evaluation, then once random values are set in coefficient format, switch the format to achieve what the caller asked for.
			ilvector.SwitchFormat();
		}
		m_vectors.push_back(ilvector);
	}
}

Gerard Ryan's avatar
Gerard Ryan committed
212 213
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl(const TugType& tug, const shared_ptr<DCRTPolyImpl::Params> dcrtParams, Format format)
Gerard Ryan's avatar
Gerard Ryan committed
214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253
{

	m_format = format;
	m_params = dcrtParams;

	size_t numberOfTowers = dcrtParams->GetParams().size();
	m_vectors.reserve(numberOfTowers);

	//tug generating random values
	std::shared_ptr<int32_t> tugValues = tug.GenerateIntVector(dcrtParams->GetRingDimension());

	for (usint i = 0; i < numberOfTowers; i++) {

		NativeVector ilTugValues(dcrtParams->GetRingDimension(), dcrtParams->GetParams()[i]->GetModulus());

		for(usint j = 0; j < dcrtParams->GetRingDimension(); j++) {
			uint64_t	entry;
			// if the random generated value is less than zero, then multiply it by (-1) and subtract the modulus of the current tower to set the coefficient
			int64_t k = (tugValues.get())[j];
			if(k < 0) {
				k *= (-1);
				entry = (uint64_t)dcrtParams->GetParams()[i]->GetModulus().ConvertToInt() - (uint64_t)k;
			}
			//if greater than or equal to zero, set it the value generated
			else {
				entry = k;
			}
			ilTugValues.at(j)=entry;
		}

		PolyType ilvector(dcrtParams->GetParams()[i]);
		ilvector.SetValues(ilTugValues, Format::COEFFICIENT); // the random values are set in coefficient format
		if(m_format == Format::EVALUATION) {  // if the input format is evaluation, then once random values are set in coefficient format, switch the format to achieve what the caller asked for.
			ilvector.SwitchFormat();
		}
		m_vectors.push_back(ilvector);
	}

}

Gerard Ryan's avatar
Gerard Ryan committed
254
/*Move constructor*/
Gerard Ryan's avatar
Gerard Ryan committed
255 256
template<typename VecType>
DCRTPolyImpl<VecType>::DCRTPolyImpl(const DCRTPolyImpl &&element)
Gerard Ryan's avatar
Gerard Ryan committed
257 258 259 260 261 262
{
	m_format = element.m_format;
	m_vectors = std::move(element.m_vectors);
	m_params = std::move(element.m_params);
}

Gerard Ryan's avatar
Gerard Ryan committed
263 264
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::CloneParametersOnly() const
Gerard Ryan's avatar
Gerard Ryan committed
265 266 267 268 269 270
{

	DCRTPolyImpl res(this->m_params, this->m_format);
	return std::move(res);
}

Gerard Ryan's avatar
Gerard Ryan committed
271 272
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::CloneWithNoise(const DiscreteGaussianGeneratorImpl<VecType> &dgg, Format format) const
Gerard Ryan's avatar
Gerard Ryan committed
273 274 275 276
{

	DCRTPolyImpl res = CloneParametersOnly();

Gerard Ryan's avatar
Gerard Ryan committed
277
	VecType randVec = dgg.GenerateVector(m_params->GetCyclotomicOrder() / 2, m_params->GetModulus());
Gerard Ryan's avatar
Gerard Ryan committed
278 279 280

	// create an Element to pull from
	// create a dummy parm to use in the Poly world
Gerard Ryan's avatar
Gerard Ryan committed
281
	shared_ptr<ILParamsImpl<Integer>> parm( new ILParamsImpl<Integer>(m_params->GetCyclotomicOrder(), m_params->GetModulus(), 1) );
Gerard Ryan's avatar
Gerard Ryan committed
282
	PolyLargeType element( parm );
Gerard Ryan's avatar
Gerard Ryan committed
283 284
	element.SetValues( randVec, m_format );

Gerard Ryan's avatar
Gerard Ryan committed
285
	res = element;
Gerard Ryan's avatar
Gerard Ryan committed
286 287 288 289 290 291

	return std::move(res);
}

// DESTRUCTORS

Gerard Ryan's avatar
Gerard Ryan committed
292 293
template<typename VecType>
DCRTPolyImpl<VecType>::~DCRTPolyImpl() {}
Gerard Ryan's avatar
Gerard Ryan committed
294 295

// GET ACCESSORS
Gerard Ryan's avatar
Gerard Ryan committed
296 297
template<typename VecType>
const typename DCRTPolyImpl<VecType>::PolyType& DCRTPolyImpl<VecType>::GetElementAtIndex (usint i) const
Gerard Ryan's avatar
Gerard Ryan committed
298 299 300 301 302 303 304 305
{
	if(m_vectors.empty())
		throw std::logic_error("DCRTPolyImpl's towers are not initialized.");
	if(i > m_vectors.size()-1)
		throw std::logic_error("Index: " + std::to_string(i) + " is out of range.");
	return m_vectors[i];
}

Gerard Ryan's avatar
Gerard Ryan committed
306 307
template<typename VecType>
usint DCRTPolyImpl<VecType>::GetNumOfElements() const
Gerard Ryan's avatar
Gerard Ryan committed
308 309 310 311
{
	return m_vectors.size();
}

Gerard Ryan's avatar
Gerard Ryan committed
312 313
template<typename VecType>
const std::vector<typename DCRTPolyImpl<VecType>::PolyType>& DCRTPolyImpl<VecType>::GetAllElements() const
Gerard Ryan's avatar
Gerard Ryan committed
314 315 316 317
{
	return m_vectors;
}

Gerard Ryan's avatar
Gerard Ryan committed
318 319
template<typename VecType>
Format DCRTPolyImpl<VecType>::GetFormat() const
Gerard Ryan's avatar
Gerard Ryan committed
320 321 322 323
{
	return m_format;
}

Gerard Ryan's avatar
Gerard Ryan committed
324 325
template<typename VecType>
std::vector<DCRTPolyImpl<VecType>> DCRTPolyImpl<VecType>::BaseDecompose(usint baseBits, bool evalModeAnswer) const
Gerard Ryan's avatar
Gerard Ryan committed
326
{
Gerard Ryan's avatar
Gerard Ryan committed
327 328 329
	bool dbg_flag = false;
	DEBUG("...::BaseDecompose" );
	DEBUG("baseBits=" << baseBits );
Gerard Ryan's avatar
Gerard Ryan committed
330

Gerard Ryan's avatar
Gerard Ryan committed
331
	PolyLargeType v( CRTInterpolate() );
Gerard Ryan's avatar
Gerard Ryan committed
332

Gerard Ryan's avatar
Gerard Ryan committed
333 334
	DEBUG("<v>" << std::endl << v << "</v>" );

Gerard Ryan's avatar
Gerard Ryan committed
335
	std::vector<PolyLargeType> bdV = v.BaseDecompose(baseBits, false);
Gerard Ryan's avatar
Gerard Ryan committed
336

Gerard Ryan's avatar
Gerard Ryan committed
337 338 339 340 341
	DEBUG("<bdV>" );
	for( auto i : bdV )
		DEBUG(i );
	DEBUG("</bdV>" );

Gerard Ryan's avatar
Gerard Ryan committed
342
	std::vector<DCRTPolyImpl<VecType>> result;
Gerard Ryan's avatar
Gerard Ryan committed
343 344 345

	// populate the result by converting each of the big vectors into a VectorArray
	for( usint i=0; i<bdV.size(); i++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
346
		DCRTPolyImpl<VecType> dv(bdV[i], this->GetParams());
Gerard Ryan's avatar
Gerard Ryan committed
347 348 349 350 351
		if( evalModeAnswer )
			dv.SwitchFormat();
		result.push_back( std::move(dv) );
	}

Gerard Ryan's avatar
Gerard Ryan committed
352 353 354 355 356 357 358 359
	DEBUG("<BaseDecompose.result>" );
	for( auto i : result )
		DEBUG(i );
	DEBUG("</BaseDecompose.result>" );

	return std::move(result);
}

Gerard Ryan's avatar
Gerard Ryan committed
360 361
template<typename VecType>
std::vector<DCRTPolyImpl<VecType>> DCRTPolyImpl<VecType>::CRTDecompose(uint32_t baseBits) const
Gerard Ryan's avatar
Gerard Ryan committed
362 363
{

Gerard Ryan's avatar
Gerard Ryan committed
364 365 366 367 368 369 370 371 372 373 374
	uint32_t nWindows = 1;

	if (baseBits > 0) {
		uint32_t nBits = m_vectors[0].GetModulus().GetLengthForBase(2);

		nWindows = nBits / baseBits;
		if (nBits % baseBits > 0)
			nWindows++;
	}

	std::vector<DCRTPolyType> result(m_vectors.size()*nWindows);
Gerard Ryan's avatar
Gerard Ryan committed
375 376 377 378 379 380

	DCRTPolyType input = this->Clone();

	if (input.GetFormat() == EVALUATION)
		input.SwitchFormat();

Gerard Ryan's avatar
Gerard Ryan committed
381
#pragma omp parallel for
Gerard Ryan's avatar
Gerard Ryan committed
382 383
	for( usint i=0; i<m_vectors.size(); i++ ) {

Gerard Ryan's avatar
Gerard Ryan committed
384 385 386 387 388 389 390 391 392 393 394 395
		if (baseBits == 0)
		{
			DCRTPolyType currentDCRTPoly = input.Clone();

			for ( usint k=0; k<m_vectors.size(); k++ ){
				PolyType temp(input.m_vectors[i]);
				if (i!=k)
					temp.SwitchModulus(input.m_vectors[k].GetModulus(),input.m_vectors[k].GetRootOfUnity());
				currentDCRTPoly.m_vectors[k] = temp;
			}

			currentDCRTPoly.SwitchFormat();
Gerard Ryan's avatar
Gerard Ryan committed
396

Gerard Ryan's avatar
Gerard Ryan committed
397
			result[i] = std::move(currentDCRTPoly);
Gerard Ryan's avatar
Gerard Ryan committed
398
		}
Gerard Ryan's avatar
Gerard Ryan committed
399 400 401 402 403 404 405 406 407 408 409 410 411 412 413
		else
		{

			vector<PolyType> decomposed = input.m_vectors[i].BaseDecompose(baseBits,false);

			for (size_t j = 0; j < decomposed.size();  j++) {

				DCRTPolyType currentDCRTPoly = input.Clone();

				for ( usint k=0; k<m_vectors.size(); k++ ){
					PolyType temp(decomposed[j]);
					if (i!=k)
						temp.SwitchModulus(input.m_vectors[k].GetModulus(),input.m_vectors[k].GetRootOfUnity());
					currentDCRTPoly.m_vectors[k] = temp;
				}
Gerard Ryan's avatar
Gerard Ryan committed
414

Gerard Ryan's avatar
Gerard Ryan committed
415
				currentDCRTPoly.SwitchFormat();
Gerard Ryan's avatar
Gerard Ryan committed
416

Gerard Ryan's avatar
Gerard Ryan committed
417 418 419 420 421
				result[j + i*nWindows] = std::move(currentDCRTPoly);

			}

		}
Gerard Ryan's avatar
Gerard Ryan committed
422 423
	}

Gerard Ryan's avatar
Gerard Ryan committed
424 425 426
	return std::move(result);
}

Gerard Ryan's avatar
Gerard Ryan committed
427 428
template<typename VecType>
std::vector<DCRTPolyImpl<VecType>> DCRTPolyImpl<VecType>::PowersOfBase(usint baseBits) const
Gerard Ryan's avatar
Gerard Ryan committed
429
{
Gerard Ryan's avatar
Gerard Ryan committed
430
	bool dbg_flag = false;
Gerard Ryan's avatar
Gerard Ryan committed
431

Gerard Ryan's avatar
Gerard Ryan committed
432
	std::vector<DCRTPolyImpl<VecType>> result;
Gerard Ryan's avatar
Gerard Ryan committed
433 434 435 436 437 438 439 440 441 442

	usint nBits = m_params->GetModulus().GetLengthForBase(2);

	usint nWindows = nBits / baseBits;
	if (nBits % baseBits > 0)
		nWindows++;

	result.reserve(nWindows);

	// prepare for the calculations by gathering a big integer version of each of the little moduli
Gerard Ryan's avatar
Gerard Ryan committed
443
	std::vector<Integer> mods(m_params->GetParams().size());
Gerard Ryan's avatar
Gerard Ryan committed
444
	for( usint i = 0; i < m_params->GetParams().size(); i++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
445
		mods[i] = Integer(m_params->GetParams()[i]->GetModulus().ConvertToInt());
Gerard Ryan's avatar
Gerard Ryan committed
446 447 448
		DEBUG("DCRTPolyImpl::PowersOfBase.mods[" << i << "] = " << mods[i] );
	}

Gerard Ryan's avatar
Gerard Ryan committed
449 450 451 452

	for( usint i = 0; i < nWindows; i++ ) {
		DCRTPolyType x( m_params, m_format );

Gerard Ryan's avatar
Gerard Ryan committed
453 454
		// Shouldn't this be Integer twoPow ( Integer::ONE << (i*baseBits)  ??
		Integer twoPow( Integer(2).Exp( i*baseBits ) );
Gerard Ryan's avatar
Gerard Ryan committed
455
		DEBUG("DCRTPolyImpl::PowersOfBase.twoPow (" << i << ") = " << twoPow );
Gerard Ryan's avatar
Gerard Ryan committed
456
		for( usint t = 0; t < m_params->GetParams().size(); t++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
457 458
			DEBUG("@(" << i << ", " << t << ")" );
			DEBUG("twoPow= " << twoPow << ", mods[" << t << "]" << mods[t] );
Gerard Ryan's avatar
Gerard Ryan committed
459
			Integer pI (twoPow % mods[t]);
Gerard Ryan's avatar
Gerard Ryan committed
460 461 462
			DEBUG("twoPow= " << twoPow << ", mods[" << t << "]" << mods[t] << ";   pI.ConvertToInt=" << pI.ConvertToInt() << ";   pI=" << pI );
			DEBUG("m_vectors= " << m_vectors[t] );

Gerard Ryan's avatar
Gerard Ryan committed
463
			x.m_vectors[t] = m_vectors[t] * pI.ConvertToInt();
Gerard Ryan's avatar
Gerard Ryan committed
464
			DEBUG("DCRTPolyImpl::PowersOfBase.x.m_vectors[" << t << ", " << i << "]" << x.m_vectors[t] );
Gerard Ryan's avatar
Gerard Ryan committed
465 466 467 468 469 470 471 472 473
		}
		result.push_back( x );
	}

	return std::move(result);
}

/*VECTOR OPERATIONS*/

Gerard Ryan's avatar
Gerard Ryan committed
474 475
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::MultiplicativeInverse() const
Gerard Ryan's avatar
Gerard Ryan committed
476
{
Gerard Ryan's avatar
Gerard Ryan committed
477
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
478 479 480 481 482 483 484

	for (usint i = 0; i < m_vectors.size(); i++) {
		tmp.m_vectors[i] = m_vectors[i].MultiplicativeInverse();
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
485 486
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::ModByTwo() const
Gerard Ryan's avatar
Gerard Ryan committed
487
{
Gerard Ryan's avatar
Gerard Ryan committed
488
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
489 490 491 492 493 494 495

	for (usint i = 0; i < m_vectors.size(); i++) {
		tmp.m_vectors[i] = m_vectors[i].ModByTwo();
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
496 497
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::Plus(const DCRTPolyImpl &element) const
Gerard Ryan's avatar
Gerard Ryan committed
498 499 500 501
{
	if( m_vectors.size() != element.m_vectors.size() ) {
		throw std::logic_error("tower size mismatch; cannot add");
	}
Gerard Ryan's avatar
Gerard Ryan committed
502
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
503 504 505 506 507 508 509

	for (usint i = 0; i < tmp.m_vectors.size(); i++) {
		tmp.m_vectors[i] += element.GetElementAtIndex (i);
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
510 511
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::Negate() const
Gerard Ryan's avatar
Gerard Ryan committed
512
{
Gerard Ryan's avatar
Gerard Ryan committed
513
	DCRTPolyImpl<VecType> tmp(this->CloneParametersOnly());
Gerard Ryan's avatar
Gerard Ryan committed
514 515 516 517 518 519 520 521 522
	tmp.m_vectors.clear();

	for (usint i = 0; i < this->m_vectors.size(); i++) {
		tmp.m_vectors.push_back(std::move(this->m_vectors.at(i).Negate()));
	}

	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
523 524
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::Minus(const DCRTPolyImpl &element) const
Gerard Ryan's avatar
Gerard Ryan committed
525 526 527 528
{
	if( m_vectors.size() != element.m_vectors.size() ) {
		throw std::logic_error("tower size mismatch; cannot subtract");
	}
Gerard Ryan's avatar
Gerard Ryan committed
529
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
530 531 532 533 534 535 536

	for (usint i = 0; i < tmp.m_vectors.size(); i++) {
		tmp.m_vectors[i] -= element.GetElementAtIndex (i);
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
537 538
template<typename VecType>
const DCRTPolyImpl<VecType>& DCRTPolyImpl<VecType>::operator+=(const DCRTPolyImpl &rhs)
Gerard Ryan's avatar
Gerard Ryan committed
539 540
{
	for (usint i = 0; i < this->GetNumOfElements(); i++) {
Gerard Ryan's avatar
Gerard Ryan committed
541
		this->m_vectors[i] += rhs.m_vectors[i];
Gerard Ryan's avatar
Gerard Ryan committed
542 543 544 545 546
	}
	return *this;

}

Gerard Ryan's avatar
Gerard Ryan committed
547 548
template<typename VecType>
const DCRTPolyImpl<VecType>& DCRTPolyImpl<VecType>::operator-=(const DCRTPolyImpl &rhs)
Gerard Ryan's avatar
Gerard Ryan committed
549 550
{
	for (usint i = 0; i < this->GetNumOfElements(); i++) {
Gerard Ryan's avatar
Gerard Ryan committed
551
		this->m_vectors.at(i) -= rhs.m_vectors[i];
Gerard Ryan's avatar
Gerard Ryan committed
552 553 554 555 556
	}
	return *this;

}

Gerard Ryan's avatar
Gerard Ryan committed
557 558
template<typename VecType>
const DCRTPolyImpl<VecType>& DCRTPolyImpl<VecType>::operator*=(const DCRTPolyImpl &element)
Gerard Ryan's avatar
Gerard Ryan committed
559 560 561 562 563 564 565 566 567
{
	for (usint i = 0; i < this->m_vectors.size(); i++) {
		this->m_vectors.at(i) *= element.m_vectors.at(i);
	}

	return *this;

}

Gerard Ryan's avatar
Gerard Ryan committed
568 569
template<typename VecType>
bool DCRTPolyImpl<VecType>::operator==(const DCRTPolyImpl &rhs) const
Gerard Ryan's avatar
Gerard Ryan committed
570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589
{

	if( GetCyclotomicOrder() != rhs.GetCyclotomicOrder() )
		return false;

	if( GetModulus() != rhs.GetModulus() )
		return false;

	if (m_format != rhs.m_format) {
		return false;
	}

	if (m_vectors.size() != rhs.m_vectors.size()) {
		return false;
	}

	//check if the towers are the same
	else return (m_vectors == rhs.GetAllElements());
}

Gerard Ryan's avatar
Gerard Ryan committed
590 591
template<typename VecType>
const DCRTPolyImpl<VecType> & DCRTPolyImpl<VecType>::operator=(const DCRTPolyImpl & rhs)
Gerard Ryan's avatar
Gerard Ryan committed
592 593 594 595 596 597 598 599 600
{
	if (this != &rhs) {
		m_vectors = rhs.m_vectors;
		m_format = rhs.m_format;
		m_params = rhs.m_params;
	}
	return *this;
}

Gerard Ryan's avatar
Gerard Ryan committed
601 602
template<typename VecType>
const DCRTPolyImpl<VecType> & DCRTPolyImpl<VecType>::operator=(DCRTPolyImpl&& rhs)
Gerard Ryan's avatar
Gerard Ryan committed
603 604 605 606 607 608 609 610 611
{
	if (this != &rhs) {
		m_vectors = std::move(rhs.m_vectors);
		m_format = std::move(rhs.m_format);
		m_params = std::move(rhs.m_params);
	}
	return *this;
}

Gerard Ryan's avatar
Gerard Ryan committed
612 613
template<typename VecType>
DCRTPolyImpl<VecType>& DCRTPolyImpl<VecType>::operator=(std::initializer_list<uint64_t> rhs)
Gerard Ryan's avatar
Gerard Ryan committed
614 615 616 617 618 619 620 621
{
	usint len = rhs.size();
	static PolyType::Integer ZERO(0);
	if(!IsEmpty()) {
		usint vectorLength = this->m_vectors[0].GetLength();
		for(usint i = 0; i < m_vectors.size(); ++i) { // this loops over each tower
			for(usint j = 0; j < vectorLength; ++j) { // loops within a tower
				if(j<len) {
Gerard Ryan's avatar
Gerard Ryan committed
622
				  this->m_vectors[i].at(j)= *(rhs.begin()+j);
Gerard Ryan's avatar
Gerard Ryan committed
623
				} else {
Gerard Ryan's avatar
Gerard Ryan committed
624
				  this->m_vectors[i].at(j)= ZERO;
Gerard Ryan's avatar
Gerard Ryan committed
625 626 627 628
				}
			}
		}
	} else {
Gerard Ryan's avatar
Gerard Ryan committed
629 630
		for(size_t i=0; i<m_vectors.size(); i++) {
			NativeVector temp(m_params->GetRingDimension());
Gerard Ryan's avatar
Gerard Ryan committed
631 632 633 634 635 636 637 638 639
			temp.SetModulus(m_vectors.at(i).GetModulus());
			temp = rhs;
			m_vectors.at(i).SetValues(std::move(temp),m_format);
		}

	}
	return *this;
}

Gerard Ryan's avatar
Gerard Ryan committed
640
// Used only inside a Matrix object; so an allocator already initializes the values
Gerard Ryan's avatar
Gerard Ryan committed
641 642
template<typename VecType>
DCRTPolyImpl<VecType>& DCRTPolyImpl<VecType>::operator=(uint64_t val)
Gerard Ryan's avatar
Gerard Ryan committed
643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661
{
	if (!IsEmpty()) {
		for (usint i = 0; i < m_vectors.size(); i++) {
			m_vectors[i] = val;
		}
	}
	else {
		for (usint i = 0; i<m_vectors.size(); i++) {
			NativeVector temp(m_params->GetRingDimension());
			temp.SetModulus(m_vectors.at(i).GetModulus());
			temp = val;
			m_vectors.at(i).SetValues(std::move(temp), m_format);
		}
	}

	return *this;
}

// Used only inside a Matrix object; so an allocator already initializes the values
Gerard Ryan's avatar
Gerard Ryan committed
662 663
template<typename VecType>
DCRTPolyImpl<VecType>& DCRTPolyImpl<VecType>::operator=(std::vector<int64_t> val)
Gerard Ryan's avatar
Gerard Ryan committed
664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684
{
	if (!IsEmpty()) {
		for (usint i = 0; i < m_vectors.size(); i++) {
			m_vectors[i] = val;
		}
	}
	else {
		for (usint i = 0; i<m_vectors.size(); i++) {
			NativeVector temp(m_params->GetRingDimension());
			temp.SetModulus(m_vectors.at(i).GetModulus());
			m_vectors.at(i).SetValues(std::move(temp), m_format);
			m_vectors[i] = val;
		}
	}

	m_format = COEFFICIENT;

	return *this;
}

// Used only inside a Matrix object; so an allocator already initializes the values
Gerard Ryan's avatar
Gerard Ryan committed
685 686
template<typename VecType>
DCRTPolyImpl<VecType>& DCRTPolyImpl<VecType>::operator=(std::vector<int32_t> val)
Gerard Ryan's avatar
Gerard Ryan committed
687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707
{
	if (!IsEmpty()) {
		for (usint i = 0; i < m_vectors.size(); i++) {
			m_vectors[i] = val;
		}
	}
	else {
		for (usint i = 0; i<m_vectors.size(); i++) {
			NativeVector temp(m_params->GetRingDimension());
			temp.SetModulus(m_vectors.at(i).GetModulus());
			m_vectors.at(i).SetValues(std::move(temp), m_format);
			m_vectors[i] = val;
		}
	}

	m_format = COEFFICIENT;

	return *this;
}


Gerard Ryan's avatar
Gerard Ryan committed
708 709
/*SCALAR OPERATIONS*/

Gerard Ryan's avatar
Gerard Ryan committed
710 711
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::Plus(const Integer &element) const
Gerard Ryan's avatar
Gerard Ryan committed
712
{
Gerard Ryan's avatar
Gerard Ryan committed
713
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
714 715 716 717 718 719 720

	for (usint i = 0; i < tmp.m_vectors.size(); i++) {
		tmp.m_vectors[i] += element.ConvertToInt();
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
721 722
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::Minus(const Integer &element) const
Gerard Ryan's avatar
Gerard Ryan committed
723
{
Gerard Ryan's avatar
Gerard Ryan committed
724
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
725 726 727 728 729 730 731

	for (usint i = 0; i < tmp.m_vectors.size(); i++) {
		tmp.m_vectors[i] -= element.ConvertToInt();
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
732 733
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::Times(const DCRTPolyImpl & element) const
Gerard Ryan's avatar
Gerard Ryan committed
734 735 736 737
{
	if( m_vectors.size() != element.m_vectors.size() ) {
		throw std::logic_error("tower size mismatch; cannot multiply");
	}
Gerard Ryan's avatar
Gerard Ryan committed
738
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
739

Gerard Ryan's avatar
Gerard Ryan committed
740
#pragma omp parallel for
Gerard Ryan's avatar
Gerard Ryan committed
741 742
	for (usint i = 0; i < m_vectors.size(); i++) {
		//ModMul multiplies and performs a mod operation on the results. The mod is the modulus of each tower.
Gerard Ryan's avatar
Gerard Ryan committed
743
		tmp.m_vectors[i] *= element.m_vectors[i];
Gerard Ryan's avatar
Gerard Ryan committed
744 745 746 747
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
748 749
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::Times(const Integer &element) const
Gerard Ryan's avatar
Gerard Ryan committed
750
{
Gerard Ryan's avatar
Gerard Ryan committed
751
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
752

Gerard Ryan's avatar
Gerard Ryan committed
753
#pragma omp parallel for
Gerard Ryan's avatar
Gerard Ryan committed
754
	for (usint i = 0; i < m_vectors.size(); i++) {
Gerard Ryan's avatar
Gerard Ryan committed
755
		tmp.m_vectors[i] = tmp.m_vectors[i] * element.ConvertToInt(); // (element % Integer((*m_params)[i]->GetModulus().ConvertToInt())).ConvertToInt();
Gerard Ryan's avatar
Gerard Ryan committed
756 757 758 759
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
760 761
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::Times(
Gerard Ryan's avatar
Gerard Ryan committed
762 763
		const std::vector<NativeInteger> &element) const
{
Gerard Ryan's avatar
Gerard Ryan committed
764
	DCRTPolyImpl<VecType> tmp(*this);
Gerard Ryan's avatar
Gerard Ryan committed
765

Gerard Ryan's avatar
Gerard Ryan committed
766
#pragma omp parallel for
Gerard Ryan's avatar
Gerard Ryan committed
767
	for (usint i = 0; i < m_vectors.size(); i++) {
Gerard Ryan's avatar
Gerard Ryan committed
768
		tmp.m_vectors[i] *= element[i]; // (element % Integer((*m_params)[i]->GetModulus().ConvertToInt())).ConvertToInt();
Gerard Ryan's avatar
Gerard Ryan committed
769 770 771 772
	}
	return std::move(tmp);
}

Gerard Ryan's avatar
Gerard Ryan committed
773 774
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::MultiplyAndRound(const Integer &p, const Integer &q) const
Gerard Ryan's avatar
Gerard Ryan committed
775 776 777 778 779 780
{
	std::string errMsg = "Operation not implemented yet";
	throw std::runtime_error(errMsg);
	return *this;
}

Gerard Ryan's avatar
Gerard Ryan committed
781 782
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::DivideAndRound(const Integer &q) const
Gerard Ryan's avatar
Gerard Ryan committed
783 784 785 786 787 788
{
	std::string errMsg = "Operation not implemented yet";
	throw std::runtime_error(errMsg);
	return *this;
}

Gerard Ryan's avatar
Gerard Ryan committed
789 790
template<typename VecType>
const DCRTPolyImpl<VecType>& DCRTPolyImpl<VecType>::operator*=(const Integer &element)
Gerard Ryan's avatar
Gerard Ryan committed
791 792
{
	for (usint i = 0; i < this->m_vectors.size(); i++) {
Gerard Ryan's avatar
Gerard Ryan committed
793
		this->m_vectors.at(i) *= element.ConvertToInt(); //this->m_vectors.at(i) * (element % Integer((*m_params)[i]->GetModulus().ConvertToInt())).ConvertToInt();
Gerard Ryan's avatar
Gerard Ryan committed
794 795 796 797 798 799
	}

	return *this;
}


Gerard Ryan's avatar
Gerard Ryan committed
800 801
template<typename VecType>
  void DCRTPolyImpl<VecType>::SetValuesToZero()
Gerard Ryan's avatar
Gerard Ryan committed
802 803 804 805 806 807 808
  {
  	for(usint i = 0; i < m_vectors.size(); i++) {
  		m_vectors[i].SetValuesToZero();
  	}
  }
/*OTHER FUNCTIONS*/

Gerard Ryan's avatar
Gerard Ryan committed
809

Gerard Ryan's avatar
Gerard Ryan committed
810
  
Gerard Ryan's avatar
Gerard Ryan committed
811 812
template<typename VecType>
void DCRTPolyImpl<VecType>::AddILElementOne()
Gerard Ryan's avatar
Gerard Ryan committed
813 814
{
	if(m_format != Format::EVALUATION)
Gerard Ryan's avatar
Gerard Ryan committed
815
		throw std::runtime_error("DCRTPolyImpl<VecType>::AddILElementOne cannot be called on a DCRTPolyImpl in COEFFICIENT format.");
Gerard Ryan's avatar
Gerard Ryan committed
816 817 818 819 820
	for(usint i = 0; i < m_vectors.size(); i++) {
		m_vectors[i].AddILElementOne();
	}
}

Gerard Ryan's avatar
Gerard Ryan committed
821 822
template<typename VecType>
void DCRTPolyImpl<VecType>::MakeSparse(const uint32_t &wFactor)
Gerard Ryan's avatar
Gerard Ryan committed
823 824 825 826 827 828 829 830
{
	for(usint i = 0; i < m_vectors.size(); i++) {
		m_vectors[i].MakeSparse(wFactor);
	}
}

// This function modifies PolyArrayImpl to keep all the even indices in the tower.
// It reduces the ring dimension of the tower by half.
Gerard Ryan's avatar
Gerard Ryan committed
831 832
template<typename VecType>
void DCRTPolyImpl<VecType>::Decompose()
Gerard Ryan's avatar
Gerard Ryan committed
833 834 835 836 837 838 839 840 841 842 843 844
{

	if(m_format != Format::COEFFICIENT) {
		std::string errMsg = "DCRTPolyImpl not in COEFFICIENT format to perform Decompose.";
		throw std::runtime_error(errMsg);
	}

	for( size_t i = 0; i < m_vectors.size(); i++) {
		m_vectors[i].Decompose();
	}

	// the individual vectors parms have changed, so change the DCRT parms
Gerard Ryan's avatar
Gerard Ryan committed
845
	std::vector<std::shared_ptr<ILNativeParams>> vparms(m_vectors.size());
Gerard Ryan's avatar
Gerard Ryan committed
846 847
	for( size_t i = 0; i < m_vectors.size(); i++)
		vparms[i] = m_vectors[i].GetParams();
Gerard Ryan's avatar
Gerard Ryan committed
848
	m_params.reset( new DCRTPolyImpl::Params(vparms[0]->GetCyclotomicOrder(), vparms) );
Gerard Ryan's avatar
Gerard Ryan committed
849 850
}

Gerard Ryan's avatar
Gerard Ryan committed
851 852
template<typename VecType>
bool DCRTPolyImpl<VecType>::IsEmpty() const
Gerard Ryan's avatar
Gerard Ryan committed
853 854 855 856 857 858 859 860
{
	for(size_t i=0; i<m_vectors.size(); i++) {
		if(!m_vectors.at(i).IsEmpty())
			return false;
	}
	return true;
}

Gerard Ryan's avatar
Gerard Ryan committed
861 862
template<typename VecType>
void DCRTPolyImpl<VecType>::DropLastElement()
Gerard Ryan's avatar
Gerard Ryan committed
863 864 865 866 867 868
{
	if(m_vectors.size() == 0) {
		throw std::out_of_range("Last element being removed from empty list");
	}

	m_vectors.resize(m_vectors.size() - 1);
Gerard Ryan's avatar
Gerard Ryan committed
869
	DCRTPolyImpl::Params *newP = new DCRTPolyImpl::Params( *m_params );
Gerard Ryan's avatar
Gerard Ryan committed
870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886
	newP->PopLastParam();
	m_params.reset(newP);
}

/**
* This function performs ModReduce on ciphertext element and private key element. The algorithm can be found from this paper:
* D.Cousins, K. Rohloff, A Scalabale Implementation of Fully Homomorphic Encyrption Built on NTRU, October 2014, Financial Cryptography and Data Security
* http://link.springer.com/chapter/10.1007/978-3-662-44774-1_18
*
* Modulus reduction reduces a ciphertext from modulus q to a smaller modulus q/qi. The qi is generally the largest. In the code below,
* ModReduce is written for DCRTPolyImpl and it drops the last tower while updating the necessary parameters.
* The steps taken here are as follows:
* 1. compute a short d in R such that d = c mod q
* 2. compute a short delta in R such that delta = (vq′−1)·d mod (pq′). E.g., all of delta’s integer coefficients can be in the range [−pq′/2, pq′/2).
* 3. let d′ = c + delta mod q. By construction, d′ is divisible by q′.
* 4. output (d′/q′) in R(q/q′).
*/
Gerard Ryan's avatar
Gerard Ryan committed
887 888
template<typename VecType>
void DCRTPolyImpl<VecType>::ModReduce(const Integer &plaintextModulus)
Gerard Ryan's avatar
Gerard Ryan committed
889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953
{
	bool dbg_flag = false;
	if(m_format != Format::EVALUATION) {
		throw std::logic_error("Mod Reduce function expects EVAL Formatted DCRTPolyImpl. It was passed COEFF Formatted DCRTPolyImpl.");
	}
	this->SwitchFormat();

	usint lastTowerIndex = m_vectors.size() - 1;

	DEBUG("ModReduce(" << plaintextModulus << ") on tower size " << m_vectors.size()<< " m=" << GetCyclotomicOrder());

	PolyType towerT(m_vectors[lastTowerIndex]); //last tower that will be dropped
	PolyType d(towerT);

	//precomputations
	typename PolyType::Integer ptm(plaintextModulus.ConvertToInt());
	typename PolyType::Integer qt(m_vectors[lastTowerIndex].GetModulus());
	DEBUG("qt: "<< qt);
	DEBUG("plaintextModulus: "<< ptm);
	typename PolyType::Integer v(qt.ModInverse(ptm));
	DEBUG("v: "<< v);
	typename PolyType::Integer a((v * qt).ModSub(1, ptm*qt));
	DEBUG("a:	"<<a);

	// Since only positive values are being used for Discrete gaussian generator, a call to switch modulus needs to be done
	d.SwitchModulus( ptm*qt, d.GetRootOfUnity() );
	// FIXME NOT CHANGING ROOT OF UNITY-TODO: What to do with SwitchModulus and is it necessary to pass rootOfUnity

	// Calculating delta, step 2
	PolyType delta(d.Times(a));

	// Calculating d' = c + delta mod q (step 3)
	// no point in going to size() since the last tower's being dropped
	for(usint i=0; i<m_vectors.size(); i++) {
		PolyType temp(delta);
		temp.SwitchModulus(m_vectors[i].GetModulus(), m_vectors[i].GetRootOfUnity());
		m_vectors[i] += temp;
	}

	//step 4
	DropLastElement();

	std::vector<PolyType::Integer> qtInverseModQi(m_vectors.size());
	for(usint i=0; i<m_vectors.size(); i++) {
		const PolyType::Integer& mod = m_vectors[i].GetModulus();
		qtInverseModQi[i] = qt.ModInverse(mod);
		m_vectors[i] = qtInverseModQi[i].ConvertToInt() * m_vectors[i];
	}

	SwitchFormat();
}

/*
 * This method applies the Chinese Remainder Interpolation on an ILVectoArray2n and produces an Poly
* How the Algorithm works:
* Consider the DCRTPolyImpl as a 2-dimensional matrix M, with dimension ringDimension * Number of Towers.
* For brevity , lets say this is r * t
* Let qt denote the bigModulus (all the towers' moduli multiplied together) and qi denote the modulus of a particular tower.
* Let V be a BigVector of size tower (tower size). Each coefficient of V is calculated as follows:
* for every r
*   calculate: V[j]= {Sigma(i = 0 --> t-1) ValueOf M(r,i) * qt/qi *[ (qt/qi)^(-1) mod qi ]}mod qt
*
* Once we have the V values, we construct an Poly from V, use qt as it's modulus, and calculate a root of unity
* for parameter selection of the Poly.
*/
Gerard Ryan's avatar
Gerard Ryan committed
954 955
template<typename VecType>
typename DCRTPolyImpl<VecType>::PolyLargeType DCRTPolyImpl<VecType>::CRTInterpolate() const
Gerard Ryan's avatar
Gerard Ryan committed
956 957 958 959 960 961
{
	bool dbg_flag = false;

	usint ringDimension = GetRingDimension();
	usint nTowers = m_vectors.size();

Gerard Ryan's avatar
Gerard Ryan committed
962
	DEBUG("in Interpolate ring " << ringDimension << " towers " << nTowers);
Gerard Ryan's avatar
Gerard Ryan committed
963 964 965 966

	for( usint vi = 0; vi < nTowers; vi++ )
		DEBUG("tower " << vi << " is " << m_vectors[vi]);

Gerard Ryan's avatar
Gerard Ryan committed
967
	Integer bigModulus(GetModulus()); // qT
Gerard Ryan's avatar
Gerard Ryan committed
968 969 970 971

	DEBUG("bigModulus " << bigModulus);

	// this is the resulting vector of coefficients
Gerard Ryan's avatar
Gerard Ryan committed
972
	VecType coefficients(ringDimension, bigModulus);
Gerard Ryan's avatar
Gerard Ryan committed
973 974 975 976

	// this will finally be  V[j]= {Sigma(i = 0 --> t-1) ValueOf M(r,i) * qt/qj *[ (qt/qj)^(-1) mod qj ]}modqt

	// first, precompute qt/qj factors
Gerard Ryan's avatar
Gerard Ryan committed
977
	vector<Integer> multiplier(nTowers);
Gerard Ryan's avatar
Gerard Ryan committed
978
	for( usint vi = 0 ; vi < nTowers; vi++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
979 980 981
		Integer qj(m_vectors[vi].GetModulus().ConvertToInt());
		Integer divBy = bigModulus / qj;
		Integer modInv = divBy.ModInverse(qj).Mod(qj);
Gerard Ryan's avatar
Gerard Ryan committed
982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002
		multiplier[vi] = divBy * modInv;

		DEBUG("multiplier " << vi << " " << qj << " " << multiplier[vi]);
	}

	// if the vectors are not in COEFFICIENT form, they need to be, so we will need to make a copy
	// of them and switchformat on them... otherwise we can just use what we have
	const std::vector<PolyType> *vecs = &m_vectors;
	std::vector<PolyType> coeffVecs;
	if( m_format == EVALUATION ) {
		for( usint i=0; i<m_vectors.size(); i++ ) {
			PolyType vecCopy(m_vectors[i]);
			vecCopy.SetFormat(COEFFICIENT);
			coeffVecs.push_back( std::move(vecCopy) );
		}
		vecs = &coeffVecs;
	}

	for( usint vi = 0; vi < nTowers; vi++ )
		DEBUG("tower " << vi << " is " << (*vecs)[vi]);

Gerard Ryan's avatar
Gerard Ryan committed
1003
	//Precompute the Barrett mu parameter
Gerard Ryan's avatar
Gerard Ryan committed
1004
	Integer mu = ComputeMu<Integer>(bigModulus);
Gerard Ryan's avatar
Gerard Ryan committed
1005

Gerard Ryan's avatar
Gerard Ryan committed
1006
	// now, compute the values for the vector
Gerard Ryan's avatar
Gerard Ryan committed
1007
#pragma omp parallel for
Gerard Ryan's avatar
Gerard Ryan committed
1008 1009 1010
	for( usint ri = 0; ri < ringDimension; ri++ ) {
		coefficients[ri] = 0;
		for( usint vi = 0; vi < nTowers; vi++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
1011
			coefficients[ri] += (Integer((*vecs)[vi].GetValues()[ri].ConvertToInt()) * multiplier[vi]);
Gerard Ryan's avatar
Gerard Ryan committed
1012 1013
		}
		DEBUG( (*vecs)[0].GetValues()[ri] << " * " << multiplier[0] << " == " << coefficients[ri] );
Gerard Ryan's avatar
Gerard Ryan committed
1014
		coefficients[ri].ModBarrettInPlace(bigModulus,mu);
Gerard Ryan's avatar
Gerard Ryan committed
1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026
	}

	DEBUG("passed loops");
	DEBUG(coefficients);

	// Create an Poly for this BigVector

	DEBUG("elementing after vectoring");
	DEBUG("m_cyclotomicOrder " << GetCyclotomicOrder());
	DEBUG("modulus "<< bigModulus);

	// Setting the root of unity to ONE as the calculation is expensive and not required.
Gerard Ryan's avatar
Gerard Ryan committed
1027
	typename DCRTPolyImpl<VecType>::PolyLargeType polynomialReconstructed( shared_ptr<ILParamsImpl<Integer>>( new ILParamsImpl<Integer>(GetCyclotomicOrder(), bigModulus, 1) ) );
Gerard Ryan's avatar
Gerard Ryan committed
1028 1029 1030 1031 1032 1033 1034
	polynomialReconstructed.SetValues(coefficients,COEFFICIENT);

	DEBUG("answer: " << polynomialReconstructed);

	return std::move( polynomialReconstructed );
}

Gerard Ryan's avatar
Gerard Ryan committed
1035
// todo can we be smarter with this method?
Gerard Ryan's avatar
Gerard Ryan committed
1036 1037
template<typename VecType>
NativePoly DCRTPolyImpl<VecType>::DecryptionCRTInterpolate(PlaintextModulus ptm) const {
Gerard Ryan's avatar
Gerard Ryan committed
1038 1039 1040
	return this->CRTInterpolate().DecryptionCRTInterpolate(ptm);
}

Gerard Ryan's avatar
Gerard Ryan committed
1041
//Source: Halevi S., Polyakov Y., and Shoup V. An Improved RNS Variant of the BFV Homomorphic Encryption Scheme. Cryptology ePrint Archive, Report 2018/117. (https://eprint.iacr.org/2018/117)
Gerard Ryan's avatar
Gerard Ryan committed
1042 1043 1044 1045 1046 1047 1048
//
//Computes Round(p/q*x) mod p as [\sum_i x_i*alpha_i + Round(\sum_i x_i*beta_i)] mod p for fast rounding in RNS
// vectors alpha and beta are precomputed as
// alpha_i = Floor[(p*[(q/qi)^{-1}]_qi)/qi]_p
// beta_i = ((p*[(q/qi)^{-1}]_qi)%qi)/qi in (0,1)
// used in decryption of BFVrns

Gerard Ryan's avatar
Gerard Ryan committed
1049 1050 1051 1052 1053
template<typename VecType>
PolyImpl<NativeVector>
DCRTPolyImpl<VecType>::ScaleAndRound(const NativeInteger &p,
		const std::vector<NativeInteger> &alpha, const std::vector<double> &beta,
		const std::vector<NativeInteger> &alphaPrecon, const std::vector<QuadFloat> &quadBeta,
Gerard Ryan's avatar
Gerard Ryan committed
1054
		const std::vector<long double> &extBeta) const {
Gerard Ryan's avatar
Gerard Ryan committed
1055 1056 1057 1058

	usint ringDimension = GetRingDimension();
	usint nTowers = m_vectors.size();

Gerard Ryan's avatar
Gerard Ryan committed
1059
	typename PolyType::Vector coefficients(ringDimension, p.ConvertToInt());
Gerard Ryan's avatar
Gerard Ryan committed
1060

Gerard Ryan's avatar
Gerard Ryan committed
1061 1062 1063 1064 1065
	if(m_vectors[0].GetModulus().GetMSB() < 45)
	{
#pragma omp parallel for
		for( usint ri = 0; ri < ringDimension; ri++ ) {
			double curFloatSum = 0.0;
Gerard Ryan's avatar
Gerard Ryan committed
1066
			NativeInteger curIntSum = 0;
Gerard Ryan's avatar
Gerard Ryan committed
1067
			for( usint vi = 0; vi < nTowers; vi++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
1068
				const NativeInteger &xi = m_vectors[vi].GetValues()[ri];
Gerard Ryan's avatar
Gerard Ryan committed
1069 1070 1071 1072 1073 1074 1075 1076 1077

				// We assume that that the value of p is smaller than 64 bits (like 58)
				// Thus we do not make additional curIntSum.Mod(p) calls for each value of vi
				//curIntSum += xi.ModMul(alpha[vi],p);
				curIntSum += xi.ModMulPreconOptimized(alpha[vi],p,alphaPrecon[vi]);

				curFloatSum += (double)(xi.ConvertToInt())*beta[vi];
			}

Gerard Ryan's avatar
Gerard Ryan committed
1078
			coefficients[ri] = (curIntSum + NativeInteger(std::llround(curFloatSum))).Mod(p);
Gerard Ryan's avatar
Gerard Ryan committed
1079 1080 1081 1082 1083 1084 1085
		}
	}
	else if (m_vectors[0].GetModulus().GetMSB() < 58)
	{
#pragma omp parallel for
		for( usint ri = 0; ri < ringDimension; ri++ ) {
			long double curFloatSum = 0.0;
Gerard Ryan's avatar
Gerard Ryan committed
1086
			NativeInteger curIntSum = 0;
Gerard Ryan's avatar
Gerard Ryan committed
1087
			for( usint vi = 0; vi < nTowers; vi++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
1088
				const NativeInteger &xi = m_vectors[vi].GetValues()[ri];
Gerard Ryan's avatar
Gerard Ryan committed
1089 1090 1091 1092 1093 1094 1095 1096 1097

				// We assume that that the value of p is smaller than 64 bits (like 58)
				// Thus we do not make additional curIntSum.Mod(p) calls for each value of vi
				//curIntSum += xi.ModMul(alpha[vi],p);
				curIntSum += xi.ModMulPreconOptimized(alpha[vi],p,alphaPrecon[vi]);

				curFloatSum += (long double)(xi.ConvertToInt())*extBeta[vi];
			}

Gerard Ryan's avatar
Gerard Ryan committed
1098
			coefficients[ri] = (curIntSum + NativeInteger(std::llround(curFloatSum))).Mod(p);
Gerard Ryan's avatar
Gerard Ryan committed
1099 1100 1101 1102
		}
	}
	else
	{
Gerard Ryan's avatar
Gerard Ryan committed
1103

Gerard Ryan's avatar
Gerard Ryan committed
1104 1105 1106
		if (nTowers > 16) // handles the case when curFloatSum exceeds 2^63 (causing an an overflow in int)
			{
			QuadFloat pFloat = quadFloatFromInt64(p.ConvertToInt());
Gerard Ryan's avatar
Gerard Ryan committed
1107

Gerard Ryan's avatar
Gerard Ryan committed
1108 1109 1110
#pragma omp parallel for
			for( usint ri = 0; ri < ringDimension; ri++ ) {
				QuadFloat curFloatSum = QuadFloat(0);
Gerard Ryan's avatar
Gerard Ryan committed
1111
				NativeInteger curIntSum = 0;
Gerard Ryan's avatar
Gerard Ryan committed
1112
				for( usint vi = 0; vi < nTowers; vi++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
1113
					const NativeInteger &xi = m_vectors[vi].GetValues()[ri];
Gerard Ryan's avatar
Gerard Ryan committed
1114 1115 1116 1117 1118 1119 1120 1121 1122

					// We assume that that the value of p is smaller than 64 bits (like 58)
					// Thus we do not make additional curIntSum.Mod(p) calls for each value of vi
					//curIntSum += xi.ModMul(alpha[vi],p);
					curIntSum += xi.ModMulPreconOptimized(alpha[vi],p,alphaPrecon[vi]);

					curFloatSum += quadFloatFromInt64(xi.ConvertToInt())*quadBeta[vi];
				}

Gerard Ryan's avatar
Gerard Ryan committed
1123
				coefficients[ri] = (curIntSum + NativeInteger(quadFloatRound(curFloatSum - pFloat*floor(curFloatSum/pFloat)))).Mod(p);
Gerard Ryan's avatar
Gerard Ryan committed
1124
			}
Gerard Ryan's avatar
Gerard Ryan committed
1125
		}
Gerard Ryan's avatar
Gerard Ryan committed
1126 1127 1128 1129 1130
		else
		{
#pragma omp parallel for
			for( usint ri = 0; ri < ringDimension; ri++ ) {
				QuadFloat curFloatSum = QuadFloat(0);
Gerard Ryan's avatar
Gerard Ryan committed
1131
				NativeInteger curIntSum = 0;
Gerard Ryan's avatar
Gerard Ryan committed
1132
				for( usint vi = 0; vi < nTowers; vi++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
1133
					const NativeInteger &xi = m_vectors[vi].GetValues()[ri];
Gerard Ryan's avatar
Gerard Ryan committed
1134 1135 1136 1137 1138 1139 1140 1141

					// We assume that that the value of p is smaller than 64 bits (like 58)
					// Thus we do not make additional curIntSum.Mod(p) calls for each value of vi
					//curIntSum += xi.ModMul(alpha[vi],p);
					curIntSum += xi.ModMulPreconOptimized(alpha[vi],p,alphaPrecon[vi]);

					curFloatSum += quadFloatFromInt64(xi.ConvertToInt())*quadBeta[vi];
				}
Gerard Ryan's avatar
Gerard Ryan committed
1142

Gerard Ryan's avatar
Gerard Ryan committed
1143
				coefficients[ri] = (curIntSum + NativeInteger(quadFloatRound(curFloatSum))).Mod(p);
Gerard Ryan's avatar
Gerard Ryan committed
1144 1145
			}
		}
Gerard Ryan's avatar
Gerard Ryan committed
1146 1147 1148 1149
	}

	// Setting the root of unity to ONE as the calculation is expensive
	// It is assumed that no polynomial multiplications in evaluation representation are performed after this
Gerard Ryan's avatar
Gerard Ryan committed
1150
	PolyType result( shared_ptr<typename PolyType::Params>( new typename PolyType::Params(GetCyclotomicOrder(), p.ConvertToInt(), 1) ) );
Gerard Ryan's avatar
Gerard Ryan committed
1151 1152 1153 1154 1155 1156 1157
	result.SetValues(coefficients,COEFFICIENT);

	return std::move(result);

}

/*
Gerard Ryan's avatar
Gerard Ryan committed
1158
 * Source: Halevi S., Polyakov Y., and Shoup V. An Improved RNS Variant of the BFV Homomorphic Encryption Scheme. Cryptology ePrint Archive, Report 2018/117. (https://eprint.iacr.org/2018/117)
Gerard Ryan's avatar
Gerard Ryan committed
1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182
 *
 * The goal is to switch the basis of x from Q to S
 *
 * Let us write x as
 * x = \sum_i [xi (q/qi)^{-1}]_qi \times q/qi - alpha*q,
 * where alpha is a number between 0 and k-1 (assuming we iterate over i \in [0,k-1]).
 *
 * Now let us take mod s_i (to go to the S basis).
 * mod s_i = \sum_i [xi (q/qi)^{-1}]_qi \times q/qi mod s_i - alpha*q mod s_i
 *
 * The main problem is that we need to find alpha.
 * If we know alpha, we can compute x mod s_i (assuming that q mod s_i is precomputed).
 *
 * We compute x mod s_i in two steps:
 * 	(1) find x' mod s_i = \sum_k [xi (q/qi)^{-1}]_qi \times q/qi mod s_i and find alpha when computing this sum;
 * 	(2) subtract alpha*q mod s_i from x' mod s_i.
 *
 * We compute lyam_i =  [xi (q/qi)^{-1}]_qi/qi, which is a floating-point number between 0 and 1, during the summation in step 1.
 * Then we compute alpha as Round(\sum_i lyam_i).
 *
 * Finally, we evaluate (x' - alpha q) mod s_i to get the CRT basis of x with respect to S.
 *
 */

Gerard Ryan's avatar
Gerard Ryan committed
1183 1184 1185 1186 1187
template<typename VecType>
DCRTPolyImpl<VecType> DCRTPolyImpl<VecType>::SwitchCRTBasis(
		const shared_ptr<DCRTPolyImpl::Params> params, const std::vector<NativeInteger> &qInvModqi,
		const std::vector<std::vector<NativeInteger>> &qDivqiModsi, const std::vector<NativeInteger> &qModsi,
		const std::vector<DoubleNativeInteger> &siModulimu, const std::vector<NativeInteger> &qInvModqiPrecon) const{
Gerard Ryan's avatar
Gerard Ryan committed
1188 1189 1190 1191 1192 1193 1194

	DCRTPolyType ans(params,m_format,true);

	usint ringDimension = GetRingDimension();
	usint nTowers = m_vectors.size();
	usint nTowersNew = ans.m_vectors.size();

Gerard Ryan's avatar
Gerard Ryan committed
1195
#pragma omp parallel for
Gerard Ryan's avatar
Gerard Ryan committed
1196 1197
	for( usint rIndex = 0; rIndex < ringDimension; rIndex++ ) {

Gerard Ryan's avatar
Gerard Ryan committed
1198
		std::vector<NativeInteger> xInvVector(nTowers);
Gerard Ryan's avatar
Gerard Ryan committed
1199
		double nu = 0.0;
Gerard Ryan's avatar
Gerard Ryan committed
1200 1201 1202

		// Compute alpha and vector of x_i terms
		for( usint vIndex = 0; vIndex < nTowers; vIndex++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
1203 1204
			const NativeInteger &xi = m_vectors[vIndex].GetValues()[rIndex];
			const NativeInteger &qi = m_vectors[vIndex].GetModulus();
Gerard Ryan's avatar
Gerard Ryan committed
1205 1206

			//computes [xi (q/qi)^{-1}]_qi
Gerard Ryan's avatar
Gerard Ryan committed
1207
			xInvVector[vIndex] = xi.ModMulPreconOptimized(qInvModqi[vIndex],qi,qInvModqiPrecon[vIndex]);
Gerard Ryan's avatar
Gerard Ryan committed
1208 1209

			//computes [xi (q/qi)^{-1}]_qi / qi to keep track of the number of q-overflows
Gerard Ryan's avatar
Gerard Ryan committed
1210
			nu += (double)xInvVector[vIndex].ConvertToInt()/(double)qi.ConvertToInt();
Gerard Ryan's avatar
Gerard Ryan committed
1211 1212 1213
		}

		// alpha corresponds to the number of overflows
Gerard Ryan's avatar
Gerard Ryan committed
1214
		NativeInteger alpha = std::llround(nu);
Gerard Ryan's avatar
Gerard Ryan committed
1215

Gerard Ryan's avatar
Gerard Ryan committed
1216 1217
		for (usint newvIndex = 0; newvIndex < nTowersNew; newvIndex ++ ) {

Gerard Ryan's avatar
Gerard Ryan committed
1218
			DoubleNativeInteger curValue = 0;
Gerard Ryan's avatar
Gerard Ryan committed
1219

Gerard Ryan's avatar
Gerard Ryan committed
1220
			const NativeInteger &si = ans.m_vectors[newvIndex].GetModulus();
Gerard Ryan's avatar
Gerard Ryan committed
1221 1222 1223

			//first round - compute "fast conversion"
			for( usint vIndex = 0; vIndex < nTowers; vIndex++ ) {
Gerard Ryan's avatar
Gerard Ryan committed
1224
				curValue += Mul128(xInvVector[vIndex].ConvertToInt(),qDivqiModsi[newvIndex][vIndex].ConvertToInt());
Gerard Ryan's avatar
Gerard Ryan committed
1225 1226
			}

Gerard Ryan's avatar
Gerard Ryan committed
1227
			const NativeInteger &curNativeValue = NativeInteger(BarrettUint128ModUint64( curValue, si.ConvertToInt(), siModulimu[newvIndex]));
Gerard Ryan's avatar
Gerard Ryan committed
1228 1229

			//second round - remove q-overflows
Gerard Ryan's avatar
Gerard Ryan committed
1230
			ans.m_vectors[newvIndex].at(rIndex) = curNativeValue.ModSubFast(alpha.ModMulFastOptimized(qModsi[newvIndex],si),si);
Gerard Ryan's avatar
Gerard Ryan committed
1231 1232 1233 1234 1235 1236 1237 1238 1239

		}

	}

	return std::move(ans);

}

Gerard Ryan's avatar
Gerard Ryan committed
1240
// Source: Halevi S., Polyakov Y., and Shoup V. An Improved RNS Variant of the BFV Homomorphic Encryption Scheme. Cryptology ePrint Archive, Report 2018/117. (https://eprint.iacr.org/2018/117)
Gerard Ryan's avatar
Gerard Ryan committed
1241 1242 1243 1244
//
// @brief Expands polynomial in CRT basis Q = q1*q2*...*qn to a larger CRT basis Q*S, where S = s1*s2*...*sn;
// uses SwichCRTBasis as a subroutine; Outputs the resulting polynomial in EVALUATION representation

Gerard Ryan's avatar
Gerard Ryan committed
1245 1246 1247 1248 1249
template<typename VecType>
void DCRTPolyImpl<VecType>::ExpandCRTBasis(const shared_ptr<DCRTPolyImpl::Params> paramsExpanded,
		const shared_ptr<DCRTPolyImpl::Params> params, const std::vector<NativeInteger> &qInvModqi,
		const std::vector<std::vector<NativeInteger>> &qDivqiModsi, const std::vector<NativeInteger> &qModsi,
		const std::vector<DoubleNativeInteger> &siModulimu, const std::vector<NativeInteger> &qInvModqiPrecon) {
Gerard Ryan's avatar
Gerard Ryan committed
1250 1251 1252 1253 1254 1255 1256 1257 1258

	std::vector<PolyType> polyInNTT;

	// if the input polynomial is in evaluation representation, store it for later use to reduce the number of NTTs
	if (this->GetFormat() == EVALUATION) {
		polyInNTT = m_vectors;
		this->SwitchFormat();
	}

Gerard Ryan's avatar
Gerard Ryan committed
1259
	DCRTPolyType polyWithSwitchedCRTBasis = SwitchCRTBasis(params,qInvModqi,qDivqiModsi,qModsi, siModulimu, qInvModqiPrecon);
Gerard Ryan's avatar
Gerard Ryan committed
1260 1261 1262 1263 1264 1265

	size_t size = m_vectors.size();
	size_t newSize = polyWithSwitchedCRTBasis.m_vectors.size() + size;

	m_vectors.resize(newSize);

Gerard Ryan's avatar
Gerard Ryan committed
1266
#pragma omp parallel for
Gerard Ryan's avatar
Gerard Ryan committed
1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279
	// populate the towers corresponding to CRT basis S and convert them to evaluation representation
	for (size_t i = 0; i < polyWithSwitchedCRTBasis.m_vectors.size(); i++ ) {
		m_vectors[size + i] = polyWithSwitchedCRTBasis.GetElementAtIndex(i);
		m_vectors[size + i].SwitchFormat();
	}

	if (polyInNTT.size() > 0) // if the input polynomial was in evaluation representation, use the towers for Q from it
	{
		for (size_t i = 0; i < size; i++ )
			m_vectors[i] = polyInNTT[i];
	}
	else
	{ // else call NTT for the towers for Q
Gerard Ryan's avatar
Gerard Ryan committed
1280
#pragma omp parallel for
Gerard Ryan's avatar
Gerard Ryan committed
1281 1282 1283 1284 1285 1286 1287 1288 1289 1290
		for (size_t i = 0; i <size; i++ )
			m_vectors[i].SwitchFormat();
	}

	m_format = EVALUATION;

	m_params = paramsExpanded;

}

Gerard Ryan's avatar
Gerard Ryan committed
1291 1292 1293 1294 1295 1296 1297 1298
//Source: Jean-Claude Bajard and Julien Eynard and Anwar Hasan and Vincent Zucca. A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes. Cryptology ePrint Archive: Report 2016/510. (https://eprint.iacr.org/2016/510)
//
// Computes Round(t/q*x) mod t for fast rounding in RNS
// vector qDivqiModqiTable are precomputed as (q/qi)^-1 mod qi
// matrix qDivqiModtgammaTable are precomputed as (q/qi) mod {t U gamma}, we assume t is stored first in the vector
// GCD(t, gamma) = 1
// used in decryption of BFVrnsB

Gerard Ryan's avatar
Gerard Ryan committed
1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312
template<typename VecType>
PolyImpl<NativeVector>
DCRTPolyImpl<VecType>::ScaleAndRound(
		const std::vector<NativeInteger> &qModuliTable,
		const NativeInteger &gamma,
		const NativeInteger &t,
		const NativeInteger &gammaInvModt,
		const NativeInteger &gammaInvModtPrecon,
		const std::vector<NativeInteger> &negqInvModtgammaTable,
		const std::vector<NativeInteger> &negqInvModtgammaPreconTable,
		const std::vector<NativeInteger> &tgammaqDivqiModqiTable,
		const std::vector<NativeInteger> &tgammaqDivqiModqiPreconTable,
		const std::vector<std::vector<NativeInteger>> &qDivqiModtgammaTable,
		const std::vector<std::vector<NativeInteger>> &qDivqiModtgammaPreconTable) const {
Gerard Ryan's avatar
Gerard Ryan committed
1313 1314 1315 1316

	usint n = GetRingDimension();
	usint numq = m_vectors.size();

Gerard Ryan's avatar
Gerard Ryan committed
1317
	typename PolyType::Vector coefficients(n, t.ConvertToInt());
Gerard Ryan's avatar
Gerard Ryan committed
1318 1319 1320 1321

#pragma omp parallel for
	for (usint k = 0; k < n; k++)
	{
Gerard Ryan's avatar
Gerard Ryan committed
1322
		NativeInteger sgamma = 0, st = 0, tmp, tmpt, tmpgamma;
Gerard Ryan's avatar
Gerard Ryan committed
1323 1324
		for (usint i = 0; i < numq; i++)
		{
Gerard Ryan's avatar
Gerard Ryan committed
1325 1326
			const NativeInteger &qi = qModuliTable[i];
			const NativeInteger &xi = m_vectors[i][k];
Gerard Ryan's avatar
Gerard Ryan committed
1327 1328 1329 1330 1331 1332
			tmp = xi;
			tmp.ModMulPreconOptimizedEq( tgammaqDivqiModqiTable[i], qi, tgammaqDivqiModqiPreconTable[i] ); // xi*t*gamma*(q/qi)^-1 mod qi

			tmpt = tmp.ModMulPreconOptimized( qDivqiModtgammaTable[i][0], t, qDivqiModtgammaPreconTable[i][0] ); // mod t
			tmpgamma = tmp.ModMulPreconOptimized( qDivqiModtgammaTable[i][1], gamma, qDivqiModtgammaPreconTable[i][1] ); // mod gamma

Gerard Ryan's avatar
Gerard Ryan committed
1333
			st.ModAddFastOptimizedEq( tmpt, t.ConvertToInt() );
Gerard Ryan's avatar
Gerard Ryan committed
1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345
			sgamma.ModAddFastOptimizedEq( tmpgamma, gamma );
		}

		// mul by -q^-1
		st.ModMulPreconOptimizedEq(negqInvModtgammaTable[0], t, negqInvModtgammaPreconTable[0]);
		sgamma.ModMulPreconOptimizedEq( negqInvModtgammaTable[1], gamma, negqInvModtgammaPreconTable[1] );

		if ( sgamma > (gamma >> 1) )
			sgamma = sgamma.ModSubFast( gamma, t );

		tmp = st.ModSub( sgamma, t );

Gerard Ryan's avatar
Gerard Ryan committed
1346
		coefficients[k] = tmp.ModMulPreconOptimized( gammaInvModt, t, gammaInvModtPrecon ).ConvertToInt();
Gerard Ryan's avatar
Gerard Ryan committed
1347 1348 1349 1350
	}

	// Setting the root of unity to ONE as the calculation is expensive
	// It is assumed that no polynomial multiplications in evaluation representation are performed after this
Gerard Ryan's avatar
Gerard Ryan committed
1351
	PolyType result( shared_ptr<typename PolyType::Params>( new typename PolyType::Params(GetCyclotomicOrder(), t.ConvertToInt(), 1) ) );
Gerard Ryan's avatar
Gerard Ryan committed
1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362
	result.SetValues(coefficients,COEFFICIENT);

	return std::move(result);
}

//Source: Jean-Claude Bajard and Julien Eynard and Anwar Hasan and Vincent Zucca. A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes. Cryptology ePrint Archive: Report 2016/510. (https://eprint.iacr.org/2016/510)
// Almost equivalent to "ExpandCRTBasis"
// @brief Expands polynomial in CRT basis q to a larger CRT basis {Bsk U mtilde}, mtilde is a redundant modulus used to remove q overflows generated from fast conversion.
// Outputs the resulting polynomial in CRT/RNS representation in basis {q U Bsk}
// used in EvalMult of BFVrnsB

Gerard Ryan's avatar
Gerard Ryan committed
1363 1364 1365 1366 1367
template<typename VecType>
void DCRTPolyImpl<VecType>::FastBaseConvqToBskMontgomery(
		const shared_ptr<DCRTPolyImpl::Params> paramsBsk,
		const std::vector<NativeInteger> &qModuli,
		const std::vector<NativeInteger> &BskmtildeModuli,
Gerard Ryan's avatar
Gerard Ryan committed
1368
		const std::vector<DoubleNativeInteger> &BskmtildeModulimu,
Gerard Ryan's avatar
Gerard Ryan committed
1369 1370 1371 1372 1373 1374 1375 1376 1377
		const std::vector<NativeInteger> &mtildeqDivqiModqi,
		const std::vector<NativeInteger> &mtildeqDivqiModqiPrecon,
		const std::vector<std::vector<NativeInteger>> &qDivqiModBj,
		const std::vector<NativeInteger> &qModBski,
		const std::vector<NativeInteger> &qModBskiPrecon,
		const NativeInteger &negqInvModmtilde,
		const NativeInteger &negqInvModmtildePrecon,
		const std::vector<NativeInteger> &mtildeInvModBskiTable,
		const std::vector<NativeInteger> &mtildeInvModBskiPreconTable) {
Gerard Ryan's avatar
Gerard Ryan committed
1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404

	// Input: poly in basis q
	// Output: poly in base Bsk = {B U msk}

	//computing steps 0 and 1 in Algorithm 3 in source paper.

	std::vector<PolyType> polyInNTT;

	// if the input polynomial is in evaluation representation, store it for later use to reduce the number of NTTs
	if (this->GetFormat() == EVALUATION) {
		polyInNTT = m_vectors;
		this->SwitchFormat();
	}

	size_t numq = qModuli.size();
	size_t numBsk = BskmtildeModuli.size() - 1;
	size_t newSize = numq + BskmtildeModuli.size();

	m_vectors.resize(newSize);

	uint32_t n = GetLength();

	m_params = paramsBsk;

	// ----------------------- step 0 -----------------------

	// first we twist xi by mtilde*(q/qi)^-1 mod qi
Gerard Ryan's avatar
Gerard Ryan committed
1405
	NativeInteger *ximtildeqiDivqModqi = new NativeInteger[n*numq];
Gerard Ryan's avatar
Gerard Ryan committed
1406 1407
    for (uint32_t i = 0; i < numq; i++)
    {
Gerard Ryan's avatar
Gerard Ryan committed
1408 1409
    	const NativeInteger &currentmtildeqDivqiModqi = mtildeqDivqiModqi[i];
    	const NativeInteger &currentmtildeqDivqiModqiPrecon = mtildeqDivqiModqiPrecon[i];
Gerard Ryan's avatar
Gerard Ryan committed
1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438

#pragma omp parallel for
        for (uint32_t k = 0; k < n; k++)
        {
        	ximtildeqiDivqModqi[i*n + k] = m_vectors[i][k].ModMulPreconOptimized( currentmtildeqDivqiModqi, qModuli[i], currentmtildeqDivqiModqiPrecon);
        }
    }

    for (uint32_t j = 0; j < numBsk + 1; j++)
    {
    	if( j < numBsk)
    	{
    		// TODO check this
			PolyType newvec(m_params->GetParams()[j], m_format, true);
			m_vectors[numq+j] = std::move(newvec);
    	}
    	else
    	{
    		// the mtilde vector (params not important)
    		PolyType newvec(m_params->GetParams()[0], m_format, true);
			m_vectors[numq+j] = std::move(newvec);
    	}

#pragma omp parallel for
    	for ( uint32_t k = 0; k < n; k++ )
    	{
    		DoubleNativeInteger result = 0;
    		for (uint32_t i = 0; i < numq; i++)
    		{
Gerard Ryan's avatar
Gerard Ryan committed
1439
    			const NativeInteger &qDivqiModBjValue = qDivqiModBj[i][j];
Gerard Ryan's avatar
Gerard Ryan committed
1440 1441 1442 1443 1444 1445 1446 1447 1448
    			result += Mul128( ximtildeqiDivqModqi[i*n+k].ConvertToInt(), qDivqiModBjValue.ConvertToInt() );
    		}
    		m_vectors[numq+j][k] = BarrettUint128ModUint64( result, BskmtildeModuli[j].ConvertToInt(), BskmtildeModulimu[j] );
    	}
    }

    // now we have input in Basis (q U Bsk U mtilde)
    // next we perform Small Motgomery Reduction mod q
    // ----------------------- step 1 -----------------------
Gerard Ryan's avatar
Gerard Ryan committed
1449
    const NativeInteger &mtilde = BskmtildeModuli[numBsk];
Gerard Ryan's avatar
Gerard Ryan committed
1450

Gerard Ryan's avatar
Gerard Ryan committed
1451
    NativeInteger *r_m_tildes = new NativeInteger[n];
Gerard Ryan's avatar
Gerard Ryan committed
1452 1453 1454 1455 1456 1457 1458 1459 1460 1461

#pragma omp parallel for
    for ( uint32_t k = 0; k < n; k++ )
	{
    	r_m_tildes[k] = m_vectors[numq+numBsk][k]; // c``_mtilde
		r_m_tildes[k].ModMulPreconOptimizedEq(negqInvModmtilde, mtilde, negqInvModmtildePrecon); // c``_mtilde*-1/q mod mtilde
	}

    for (uint32_t i = 0; i < numBsk; i++)
    {
Gerard Ryan's avatar
Gerard Ryan committed
1462 1463
    	const NativeInteger &currentqModBski = qModBski[i];
    	const NativeInteger &currentqModBskiPrecon = qModBskiPrecon[i];
Gerard Ryan's avatar
Gerard Ryan committed
1464 1465 1466 1467 1468

#pragma omp parallel for
    	for ( uint32_t k = 0; k < n; k++ )
		{
    		// collapsing
Gerard Ryan's avatar
Gerard Ryan committed
1469
    		NativeInteger r_m_tilde = r_m_tildes[k]; // m_tilde < than all Bsk_i
Gerard Ryan's avatar
Gerard Ryan committed
1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509
    		r_m_tilde.ModMulPreconOptimizedEq( currentqModBski, BskmtildeModuli[i], currentqModBskiPrecon ); // (r_mtilde) * q mod Bski
    		r_m_tilde.ModAddFastOptimizedEq( m_vectors[numq+i][k], BskmtildeModuli[i] ); // (c``_m + (r_mtilde* q)) mod Bski
    		m_vectors[numq+i][k] = r_m_tilde.ModMulPreconOptimized( mtildeInvModBskiTable[i], BskmtildeModuli[i], mtildeInvModBskiPreconTable[i] ); // (c``_m + (r_mtilde* q)) * mtilde mod Bski
		}
    }

    // remove mtilde residue
    m_vectors.erase(m_vectors.begin()+numq+numBsk);

    if (polyInNTT.size() > 0) // if the input polynomial was in evaluation representation, use the towers for Q from it
	{
		for (size_t i = 0; i < numq; i++ )
			m_vectors[i] = polyInNTT[i];
	}
	else
	{ // else call NTT for the towers for q
#pragma omp parallel for
		for (size_t i = 0; i <numq; i++ )
			m_vectors[i].SwitchFormat();
	}

#pragma omp parallel for
    for (uint32_t i = 0; i < numBsk; i++)
		m_vectors[numq+i].SwitchFormat();


    m_format = EVALUATION;

    delete[] r_m_tildes;
    delete[] ximtildeqiDivqModqi;
    r_m_tildes = nullptr;
    ximtildeqiDivqModqi = nullptr;
}

//Source: Jean-Claude Bajard and Julien Eynard and Anwar Hasan and Vincent Zucca. A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes. Cryptology ePrint Archive: Report 2016/510. (https://eprint.iacr.org/2016/510)
// Almost equivalent to "ScaleAndRound"
// @brief Scales polynomial in CRT basis {q U Bsk} by scalar t/q.
// Outputs the resulting polynomial in CRT/RNS representation in basis {q U Bsk}. Note that the actual result is basically in basis {Bsk}.
// used in EvalMult of BFVrnsB

Gerard Ryan's avatar
Gerard Ryan committed
1510 1511 1512 1513 1514
template<typename VecType>
void DCRTPolyImpl<VecType>::FastRNSFloorq(
		const NativeInteger &t,
		const std::vector<NativeInteger> &qModuli,
		const std::vector<NativeInteger> &BskModuli,
Gerard Ryan's avatar
Gerard Ryan committed
1515
		const std::vector<DoubleNativeInteger> &BskModulimu,
Gerard Ryan's avatar
Gerard Ryan committed
1516 1517 1518 1519 1520
		const std::vector<NativeInteger> &tqDivqiModqi,
		const std::vector<NativeInteger> &tqDivqiModqiPrecon,
		const std::vector<std::vector<NativeInteger>> &qDivqiModBj,
		const std::vector<NativeInteger> &qInvModBi,
		const std::vector<NativeInteger> &qInvModBiPrecon) {
Gerard Ryan's avatar
Gerard Ryan committed
1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533

	// Input: poly in basis {q U Bsk}
	// Output: approximateFloor(t/q*poly) in basis Bsk

	// --------------------- step 3 ---------------------
	// approximate rounding

	size_t numq = qModuli.size();
	size_t numBsk = BskModuli.size();

	uint32_t n = GetLength();

	// Twist xi by t*(q/qi)^-1 mod qi
Gerard Ryan's avatar
Gerard Ryan committed
1534
	NativeInteger *txiqiDivqModqi = new NativeInteger[n*numBsk];
Gerard Ryan's avatar
Gerard Ryan committed
1535 1536 1537

	for (uint32_t i = 0; i < numq; i++)
	{
Gerard Ryan's avatar
Gerard Ryan committed
1538 1539
		const NativeInteger &currenttqDivqiModqi = tqDivqiModqi[i];
		const NativeInteger &currenttqDivqiModqiPrecon = tqDivqiModqiPrecon[i];
Gerard Ryan's avatar
Gerard Ryan committed
1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556

#pragma omp parallel for
		for (uint32_t k = 0; k < n; k++)
		{
			// multiply by t*(q/qi)^-1 mod qi
			m_vectors[i][k].ModMulPreconOptimizedEq(currenttqDivqiModqi, qModuli[i], currenttqDivqiModqiPrecon);
		}
	}

	for (uint32_t j = 0; j < numBsk; j++)
	{
#pragma omp parallel for
		for ( uint32_t k = 0; k < n; k++ )
		{
			DoubleNativeInteger aq = 0;
			for (uint32_t i = 0; i < numq; i++)
			{
Gerard Ryan's avatar
Gerard Ryan committed
1557 1558
				const NativeInteger &qDivqiModBjValue = qDivqiModBj[i][j];
				NativeInteger &xi = m_vectors[i][k];
Gerard Ryan's avatar
Gerard Ryan committed
1559 1560 1561 1562 1563 1564 1565 1566 1567 1568
				aq += Mul128( xi.ConvertToInt(), qDivqiModBjValue.ConvertToInt() );
			}
			txiqiDivqModqi[j*n + k] = BarrettUint128ModUint64( aq, BskModuli[j].ConvertToInt(), BskModulimu[j] );
		}
	}

	// now we have FastBaseConv( |t*ct|q, q, Bsk ) in txiqiDivqModqi

    for (uint32_t i = 0; i < numBsk; i++)
    {
Gerard Ryan's avatar
Gerard Ryan committed
1569 1570
        const NativeInteger &currentqInvModBski = qInvModBi[i];
        const NativeInteger &currentqInvModBskiPrecon = qInvModBiPrecon[i];
Gerard Ryan's avatar
Gerard Ryan committed
1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589
#pragma omp parallel for
        for (uint32_t k = 0; k < n; k++)
        {
        	// Not worthy to use lazy reduction here
        	m_vectors[i+numq][k].ModMulFastEq(t, BskModuli[i]);
        	m_vectors[i+numq][k].ModSubEq( txiqiDivqModqi[i*n+k], BskModuli[i] );
        	m_vectors[i+numq][k].ModMulPreconOptimizedEq( currentqInvModBski, BskModuli[i], currentqInvModBskiPrecon );
        }
    }
	delete[] txiqiDivqModqi;
	txiqiDivqModqi = nullptr;
}

//Source: Jean-Claude Bajard and Julien Eynard and Anwar Hasan and Vincent Zucca. A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes. Cryptology ePrint Archive: Report 2016/510. (https://eprint.iacr.org/2016/510)
// // Almost qeuivalent to "SwitchCRTBasis"
// @brief Converts fast polynomial in CRT basis {q U Bsk} to basis {q} using Shenoy Kumaresan method.
// Outputs the resulting polynomial in CRT/RNS representation in basis q. Note that the actual result is basically in basis {Bsk}.
// used in EvalMult of BFVrnsB

Gerard Ryan's avatar
Gerard Ryan committed
1590 1591 1592
template<typename VecType>
void DCRTPolyImpl<VecType>::FastBaseConvSK(
			const std::vector<NativeInteger> &qModuli,
Gerard Ryan's avatar
Gerard Ryan committed
1593
			const std::vector<DoubleNativeInteger> &qModulimu,
Gerard Ryan's avatar
Gerard Ryan committed
1594
			const std::vector<NativeInteger> &BskModuli,
Gerard Ryan's avatar
Gerard Ryan committed
1595
			const std::vector<DoubleNativeInteger> &BskModulimu,
Gerard Ryan's avatar
Gerard Ryan committed
1596 1597 1598 1599 1600 1601 1602 1603
			const std::vector<NativeInteger> &BDivBiModBi,
			const std::vector<NativeInteger> &BDivBiModBiPrecon,
			const std::vector<NativeInteger> &BDivBiModmsk,
			const NativeInteger &BInvModmsk,
			const NativeInteger &BInvModmskPrecon,
			const std::vector<std::vector<NativeInteger>> &BDivBiModqj,
			const std::vector<NativeInteger> &BModqi,
			const std::vector<NativeInteger> &BModqiPrecon
Gerard Ryan's avatar
Gerard Ryan committed
1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616
			)
{
	// Input: poly in basis Bsk
	// Output: poly in basis q

	// FastBaseconv(x, B, q)
	size_t numq = qModuli.size();
	size_t numBsk = BskModuli.size();

	uint32_t n = GetLength();

    for (uint32_t i = 0; i < numBsk-1; i++) // exclude msk residue
    {
Gerard Ryan's avatar
Gerard Ryan committed
1617 1618
        const NativeInteger &currentBDivBiModBi = BDivBiModBi[i];
        const NativeInteger &currentBDivBiModBiPrecon = BDivBiModBiPrecon[i];
Gerard Ryan's avatar
Gerard Ryan committed
1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634
#pragma omp parallel for
        for (uint32_t k = 0; k < n; k++)
        {
            m_vectors[numq+i][k].ModMulPreconOptimizedEq( currentBDivBiModBi, BskModuli[i], currentBDivBiModBiPrecon);
        }
    }

    for (uint32_t j = 0; j < numq; j++)
	{
#pragma omp parallel for
		for (uint32_t k = 0; k < n; k++)
		{

			DoubleNativeInteger result = 0;
			for (uint32_t i = 0; i < numBsk-1; i++) // exclude msk residue
			{
Gerard Ryan's avatar
Gerard Ryan committed
1635 1636
				const NativeInteger &currentBDivBiModqj = BDivBiModqj[i][j];
				const NativeInteger &xi = m_vectors[numq+i][k];