text: //==================================================================================
  // BSD 2-Clause License
  //
  // Copyright (c) 2014-2022, NJIT, Duality Technologies Inc. and other contributors
  //
  // All rights reserved.
  //
  // Author TPOC: contact@openfhe.org
  //
  // 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.
  //==================================================================================
  
  /*
    This code provides number theory utilities
   */
  
  #define _USE_MATH_DEFINES
  
  #include <time.h>
  #include <cmath>
  
  #include <chrono>
  #include <sstream>
  
  #include "math/distributiongenerator.h"
  #include "math/nbtheory.h"
  
  #include "utils/debug.h"
  
  namespace lbcrypto {
  
  /*
   Generates a random number between 0 and n.
   Input: BigInteger n.
   Output: Randomly generated BigInteger  between 0 and n.
   */
  template <typename IntType>
  static IntType RNG(const IntType& modulus) {
      // static parameters for the 32-bit unsigned integers used for multiprecision
      // random number generation
      static const usint chunk_min   = 0;
      static const usint chunk_width = std::numeric_limits<uint32_t>::digits;
      static const usint chunk_max   = std::numeric_limits<uint32_t>::max();
  
      static std::uniform_int_distribution<uint32_t> distribution(chunk_min, chunk_max);
  
      // Update values that depend on modulus.
      usint modulusWidth = modulus.GetMSB();
      // Get the number of chunks in the modulus
      // 1 is subtracted to make sure the last chunk is fully used by the modulus
      usint chunksPerValue = (modulusWidth - 1) / chunk_width;
  
      // result is initialized to 0
      IntType result;
  
      // temp is used for intermediate multiprecision computations
      IntType temp;
  
      // stores current random number generated by built-in C++ 11 uniform generator
      // (used for 32-bit unsigned integers)
      uint32_t value;
  
      do {
          result = 0;
          // Generate random uint32_t "limbs" of the BigInteger
          for (usint i = 0; i < chunksPerValue; i++) {
              // Generate an unsigned long integer
              value = distribution(PseudoRandomNumberGenerator::GetPRNG());
              // converts value into IntType
              temp = value;
              // Move it to the appropriate chunk of the big integer
              temp <<= i * chunk_width;
              // Add it to the current big integer storing the result
              result += temp;
          }
  
          // work with the remainder - after all 32-bit chunks were processed
          temp = modulus >> chunksPerValue * chunk_width;
  
          // Generate a uniform number for the remainder
          // If not 1, i.e., the modulus is either 1 or a power of 2*CHUNK_WIDTH
          if (temp.GetMSB() != 1) {
              uint32_t bound = temp.ConvertToInt();
  
              // default generator for the most significant chunk of the multiprecision
              // number
              std::uniform_int_distribution<uint32_t> distribution2 =
                  std::uniform_int_distribution<uint32_t>(chunk_min, bound);
  
              value = distribution2(PseudoRandomNumberGenerator::GetPRNG());
              // converts value into IntType
              temp = value;
              // Move it to the appropriate chunk of the big integer
              temp <<= chunksPerValue * chunk_width;
              // Add it to the current big integer storing the result
              result += temp;
          }
      } while (result > modulus);  // deals with the rare scenario when the bits in
                                   // the most significant chunk are the same
      // and the bits in the following chunk of the result are larger than in the
      // modulus
      return result;
  }
  /*
   A witness function used for the Miller-Rabin Primality test.
   Inputs: a is a randomly generated witness between 2 and p-1,
   p is the number to be tested for primality,
   s and d satisfy p-1 = ((2^s) * d), d is odd.
   Output: true if p is composite,
   false if p is likely prime
   */
  template <typename IntType>
  static bool WitnessFunction(const IntType& a, const IntType& d, usint s, const IntType& p) {
      OPENFHE_DEBUG_FLAG(false);
      OPENFHE_DEBUG("calling modexp a " << a << " d " << d << " p " << p);
      IntType mod = a.ModExp(d, p);
      OPENFHE_DEBUG("mod " << mod);
      bool prevMod = false;
      for (usint i = 1; i < s + 1; i++) {
          OPENFHE_DEBUG("wf " << i);
          if (mod != IntType(1) && mod != p - IntType(1)) {
              prevMod = true;
          }
          else {
              prevMod = false;
          }
          mod.ModMulFastEq(mod, p);
          if (mod == IntType(1) && prevMod) {
              return true;
          }
      }
      return (mod != IntType(1));
  }
  
  /*
   A helper function to RootOfUnity function. This finds a generator for a given
   prime q. Input: BigInteger q which is a prime. Output: A generator of prime q
   */
  template <typename IntType>
  static IntType FindGenerator(const IntType& q) {
      OPENFHE_DEBUG_FLAG(false);
      std::set<IntType> primeFactors;
      OPENFHE_DEBUG("FindGenerator(" << q << "),calling PrimeFactorize");
  
      IntType qm1 = q - IntType(1);
      IntType qm2 = q - IntType(2);
      PrimeFactorize<IntType>(qm1, primeFactors);
      OPENFHE_DEBUG("prime factors of " << qm1);
  #if !defined(NDEBUG)
      for (auto& v : primeFactors)
          OPENFHE_DEBUG(v << " ");
  #endif
      bool generatorFound = false;
      IntType gen;
      while (!generatorFound) {
          usint count = 0;
  
          // gen = RNG(qm2).ModAdd(IntType::ONE, q); //modadd note needed
          gen = RNG(qm2) + IntType(1);
          OPENFHE_DEBUG("generator " << gen);
          OPENFHE_DEBUG("cycling thru prime factors");
  
          for (auto it = primeFactors.begin(); it != primeFactors.end(); ++it) {
              OPENFHE_DEBUG(qm1 << " / " << *it << " " << gen.ModExp(qm1 / (*it), q));
  
              if (gen.ModExp(qm1 / (*it), q) == IntType(1))
                  break;
              else
                  count++;
          }
          if (count == primeFactors.size())
              generatorFound = true;
      }
      return gen;
  }
  
  /*
   A helper function for arbitrary cyclotomics. This finds a generator for any
   composite q (cyclic group). Input: BigInteger q (cyclic group). Output: A
   generator of q
   */
  template <typename IntType>
  IntType FindGeneratorCyclic(const IntType& q) {
      OPENFHE_DEBUG_FLAG(false);
      std::set<IntType> primeFactors;
      OPENFHE_DEBUG("calling PrimeFactorize");
  
      IntType phi_q    = IntType(GetTotient(q.ConvertToInt()));
      IntType phi_q_m1 = IntType(GetTotient(q.ConvertToInt()));
  
      PrimeFactorize<IntType>(phi_q, primeFactors);
      OPENFHE_DEBUG("done");
      bool generatorFound = false;
      IntType gen;
      while (!generatorFound) {
          usint count = 0;
          OPENFHE_DEBUG("count " << count);
  
          gen = RNG(phi_q_m1) + IntType(1);  // gen is random in [1, phi(q)]
          if (GreatestCommonDivisor<IntType>(gen, q) != IntType(1)) {
              // Generator must lie in the group!
              continue;
          }
  
          // Order of a generator cannot divide any co-factor
          for (auto it = primeFactors.begin(); it != primeFactors.end(); ++it) {
              OPENFHE_DEBUG("in set");
              OPENFHE_DEBUG("divide " << phi_q << " by " << *it);
  
              if (gen.ModExp(phi_q / (*it), q) == IntType(1))
                  break;
              else
                  count++;
          }
  
          if (count == primeFactors.size())
              generatorFound = true;
      }
      return gen;
  }
  
  /*
   A helper function for arbitrary cyclotomics. Checks if g is a generator of q
   (supports any cyclic group, not just prime-modulus groups) Input: Candidate
   generator g and modulus q Output: returns true if g is a generator for q
   */
  template <typename IntType>
  bool IsGenerator(const IntType& g, const IntType& q) {
      OPENFHE_DEBUG_FLAG(false);
      std::set<IntType> primeFactors;
      OPENFHE_DEBUG("calling PrimeFactorize");
  
      IntType qm1 = IntType(GetTotient(q.ConvertToInt()));
  
      PrimeFactorize<IntType>(qm1, primeFactors);
      OPENFHE_DEBUG("done");
  
      usint count = 0;
  
      for (auto it = primeFactors.begin(); it != primeFactors.end(); ++it) {
          OPENFHE_DEBUG("in set");
          OPENFHE_DEBUG("divide " << qm1 << " by " << *it);
  
          if (g.ModExp(qm1 / (*it), q) == IntType(1))
              break;
          else
              count++;
      }
  
      if (count == primeFactors.size())
          return true;
      else
          return false;
  }
  
  /*
   finds roots of unity for given input.  Assumes the the input is a power of two.
   Mostly likely does not give correct results otherwise. input:  m as number
   which is cyclotomic(in format of int), modulo which is used to find generator
   (in format of BigInteger)
  
   output:  root of unity (in format of BigInteger)
   */
  template <typename IntType>
  IntType RootOfUnity(usint m, const IntType& modulo) {
      OPENFHE_DEBUG_FLAG(false);
      OPENFHE_DEBUG("in Root of unity m :" << m << " modulo " << modulo.ToString());
      IntType M(m);
      if ((modulo - IntType(1)).Mod(M) != IntType(0)) {
          std::string errMsg =
              "Please provide a primeModulus(q) and a cyclotomic number(m) "
              "satisfying the condition: (q-1)/m is an integer. The values of "
              "primeModulus = " +
              modulo.ToString() + " and m = " + std::to_string(m) + " do not satisfy this condition";
          OPENFHE_THROW(math_error, errMsg);
      }
      IntType result;
      OPENFHE_DEBUG("calling FindGenerator");
      IntType gen = FindGenerator(modulo);
      OPENFHE_DEBUG("gen = " << gen.ToString());
  
      OPENFHE_DEBUG("calling gen.ModExp( " << ((modulo - IntType(1)).DividedBy(M)).ToString() << ", modulus "
                                   << modulo.ToString());
      result = gen.ModExp((modulo - IntType(1)).DividedBy(M), modulo);
      OPENFHE_DEBUG("result = " << result.ToString());
      if (result == IntType(1)) {
          OPENFHE_DEBUG("LOOP?");
          result = RootOfUnity(m, modulo);
      }
  
      /*
     * At this point, result contains a primitive root of unity. However,
     * we want to return the minimum root of unity, to avoid different
     * crypto contexts having different roots of unity for the same
     * cyclotomic order and moduli. Therefore, we are going to cycle over
     * all primitive roots of unity and select the smallest one (minRU).
     *
     * To cycle over all primitive roots of unity, we raise the root of
     * unity in result to all the powers that are co-prime to the
     * cyclotomic order. In power-of-two cyclotomics, this will be the
     * set of all odd powers, but here we use a more general routine
     * to support arbitrary cyclotomics.
     *
     */
  
      IntType mu = modulo.ComputeMu();
      IntType x(1);
      x.ModMulEq(result, modulo, mu);
      IntType minRU(x);
      IntType curPowIdx(1);
      std::vector<IntType> coprimes = GetTotientList<IntType>(m);
      for (usint i = 0; i < coprimes.size(); i++) {
          auto nextPowIdx = coprimes[i];
          IntType diffPow(nextPowIdx - curPowIdx);
          for (IntType j(0); j < diffPow; j += IntType(1)) {
              x.ModMulEq(result, modulo, mu);
          }
          if (x < minRU && x != IntType(1)) {
              minRU = x;
          }
          curPowIdx = nextPowIdx;
      }
      return minRU;
  }
  
  template <typename IntType>
  std::vector<IntType> RootsOfUnity(usint m, const std::vector<IntType> moduli) {
      std::vector<IntType> rootsOfUnity(moduli.size());
      for (usint i = 0; i < moduli.size(); i++) {
          rootsOfUnity[i] = RootOfUnity(m, moduli[i]);
      }
      return rootsOfUnity;
  }
  
  template <typename IntType>
  IntType GreatestCommonDivisor(const IntType& a, const IntType& b) {
      OPENFHE_DEBUG_FLAG(false);
      IntType m_a, m_b, m_t;
      m_a = a;
      m_b = b;
      OPENFHE_DEBUG("GCD a " << a << " b " << b);
      while (m_b != IntType(0)) {
          m_t = m_b;
          OPENFHE_DEBUG("GCD m_a.Mod(b) " << m_a << "( " << m_b << ")");
          m_b = m_a % m_b;
  
          m_a = m_t;
          OPENFHE_DEBUG("GCD m_a " << m_b << " m_b " << m_b);
      }
      OPENFHE_DEBUG("GCD ret " << m_a);
      return m_a;
  }
  
  /*
   The Miller-Rabin Primality Test
   Input: p the number to be tested for primality.
   Output: true if p is prime,
   false if p is not prime
   */
  template <typename IntType>
  bool MillerRabinPrimalityTest(const IntType& p, const usint niter) {
      OPENFHE_DEBUG_FLAG(false);
      if (p < IntType(2) || ((p != IntType(2)) && (p.Mod(2) == IntType(0))))
          return false;
      if (p == IntType(2) || p == IntType(3) || p == IntType(5))
          return true;
  
      IntType d = p - IntType(1);
      usint s   = 0;
      OPENFHE_DEBUG("start while d " << d);
      while (d.Mod(2) == IntType(0)) {
          d.DividedByEq(2);
          s++;
      }
      OPENFHE_DEBUG("end while s " << s);
      bool composite = true;
      for (usint i = 0; i < niter; i++) {
          OPENFHE_DEBUG(".1");
          IntType a = RNG(p - IntType(3)).ModAdd(2, p);
          OPENFHE_DEBUG(".2");
          composite = (WitnessFunction(a, d, s, p));
          if (composite)
              break;
      }
      OPENFHE_DEBUG("done composite " << composite);
      return (!composite);
  }
  
  /*
   The Pollard Rho factorization of a number n.
   Input: n the number to be factorized.
   Output: a factor of n.
   */
  template <typename IntType>
  const IntType PollardRhoFactorization(const IntType& n) {
      OPENFHE_DEBUG_FLAG(false);
      IntType divisor(1);
  
      IntType c(RNG(n));
      IntType x(RNG(n));
      IntType xx(x);
  
      // check divisibility by 2
      if (n.Mod(2) == IntType(0))
          return IntType(2);
  
      // Precompute the Barrett mu parameter
      IntType mu = n.ComputeMu();
  
      do {
          x  = x.ModMul(x, n, mu).ModAdd(c, n, mu);
          xx = xx.ModMul(xx, n, mu).ModAdd(c, n, mu);
          xx = xx.ModMul(xx, n, mu).ModAdd(c, n, mu);
  
          divisor = GreatestCommonDivisor((x > xx) ? x - xx : xx - x, n);
  
          OPENFHE_DEBUG("PRF divisor " << divisor.ToString());
      } while (divisor == IntType(1));
  
      return divisor;
  }
  
  /*
   Recursively factorizes and find the distinct primefactors of a number
   Input: n is the number to be prime factorized,
   primeFactors is a set of prime factors of n.
   */
  template <typename IntType>
  void PrimeFactorize(IntType n, std::set<IntType>& primeFactors) {
      OPENFHE_DEBUG_FLAG(false);
      OPENFHE_DEBUG("PrimeFactorize " << n);
  
      // primeFactors.clear();
      OPENFHE_DEBUG("In PrimeFactorize n " << n);
      OPENFHE_DEBUG("set size " << primeFactors.size());
  
      if (n == IntType(0) || n == IntType(1))
          return;
      OPENFHE_DEBUG("calling MillerRabinPrimalityTest(" << n << ")");
      if (MillerRabinPrimalityTest(n)) {
          OPENFHE_DEBUG("Miller true");
          primeFactors.insert(n);
          return;
      }
  
      OPENFHE_DEBUG("calling PrFact n " << n);
      IntType divisor(PollardRhoFactorization(n));
  
      OPENFHE_DEBUG("calling PF " << divisor);
      PrimeFactorize(divisor, primeFactors);
  
      OPENFHE_DEBUG("calling div " << divisor);
      // IntType reducedN = n.DividedBy(divisor);
      n /= divisor;
  
      OPENFHE_DEBUG("calling PF reduced n " << n);
      PrimeFactorize(n, primeFactors);
  }
  
  // issue-881: make sure we don't overflow an IntType
  template <typename IntType>
  IntType FirstPrime(uint64_t nBits, uint64_t m) {
  #if (NATIVEINT == 64 && !defined(HAVE_INT128)) || (NATIVEINT == 128)
      if ((typeid(IntType) == typeid(NativeInteger)) && (nBits >= MAX_MODULUS_SIZE)) {
          OPENFHE_THROW(math_error, "Requested modulus size " + std::to_string(nBits + 1) +
                                         " exceeds maximum allowed size " + std::to_string(MAX_MODULUS_SIZE));
      }
  #endif
      try {
          OPENFHE_DEBUG_FLAG(false);
          IntType r = IntType(2).ModExp(nBits, m);
          OPENFHE_DEBUG("r " << r);
          IntType qNew = (IntType(1) << nBits);
          double d1    = qNew.ConvertToDouble();
          double d2    = std::pow(2, static_cast<double>(nBits));
          if ((qNew == IntType(0)) || (std::llround(std::log2(d1)) < std::llround(std::log2(d2)))) {
              // too big for IntType
              OPENFHE_THROW(math_error,
                             "FirstPrime parameters are too large for this integer "
                             "implementation");
          }
          IntType qNew2 = (r > IntType(0)) ? (qNew + (IntType(m) - r) + IntType(1)) : (qNew + IntType(1));  // if r == 0
          if (qNew2 < qNew) {
              OPENFHE_THROW(math_error, "FirstPrime parameters overflow this integer implementation");
          }
          qNew = qNew2;
  
          while (!MillerRabinPrimalityTest(qNew)) {
              qNew2 = qNew + IntType(m);
              if (qNew2 < qNew) {
                  OPENFHE_THROW(math_error, "FirstPrime overflow growing candidate");
              }
              qNew = qNew2;
          }
  
          return qNew;
      }
      catch (...) {
          OPENFHE_THROW(math_error, "FirstPrime math exception");
      }
      //  try {
      //    OPENFHE_DEBUG_FLAG(false);
      //    IntType r = IntType(2).ModExp(nBits, m);
      //    OPENFHE_DEBUG("r "<<r);
      //    IntType qNew = (IntType(1) << nBits);
      //    double d1 = qNew.ConvertToDouble();
      //    double d2 = std::pow(2, double(nBits));
      //    if ((qNew == IntType(0)) || (std::llround(std::log2(d1)) <
      // std::llround(std::log2(d2)))) {
      //      // too big for IntType
      //      OPENFHE_THROW(math_error, "FirstPrime parameters are
      // too large for this integer implementation");
      //    }
      //    IntType mint = IntType(m);
      //    qNew -= mint + r - IntType(1);
      //
      //    while (!MillerRabinPrimalityTest(qNew)) {
      //      qNew -= mint;
      //    }
      //    return qNew;
      //  } catch (...) {
      //    OPENFHE_THROW(math_error, "FirstPrime math exception");
      //  }
  }
  
  template <typename IntType>
  IntType NextPrime(const IntType& q, uint64_t m) {
      IntType M(m);
      IntType qOld(q), qNew(q);
  
      do {
          qNew += M;
          if (qNew < qOld)
              OPENFHE_THROW(math_error, "NextPrime overflow growing candidate");
      } while (!MillerRabinPrimalityTest(qNew));
  
      return qNew;
  
      return qNew;
  }
  
  template <typename IntType>
  IntType PreviousPrime(const IntType& q, uint64_t m) {
      IntType M(m);
      IntType qNew = q - M;
  
      while (!MillerRabinPrimalityTest(qNew)) {
          qNew -= M;
      }
  
      return qNew;
  }
  
  template <typename IntType>
  IntType NextPowerOfTwo(const IntType& n) {
      usint result = ceil(log2(n));
      return result;
  }
  
  /*Naive Loop to find coprimes to n*/
  template <typename IntType>
  std::vector<IntType> GetTotientList(const IntType& n) {
      std::vector<IntType> result;
      IntType one(1);
      for (IntType i = IntType(1); i < n; i = i + IntType(1)) {
          if (GreatestCommonDivisor(i, n) == one)
              result.push_back(i);
      }
  
      return result;
  }
  
  /* Calculate the remainder from polynomial division */
  template <typename IntVector>
  IntVector PolyMod(const IntVector& dividend, const IntVector& divisor, const typename IntVector::Integer& modulus) {
      usint divisorLength  = divisor.GetLength();
      usint dividendLength = dividend.GetLength();
  
      IntVector result(divisorLength - 1, modulus);
      usint runs = dividendLength - divisorLength + 1;  // no. of iterations
  
      // Precompute the Barrett mu parameter
      auto mu = modulus.ComputeMu();
  
      auto mat = [](const typename IntVector::Integer& x, const typename IntVector::Integer& y,
                    const typename IntVector::Integer& z, const typename IntVector::Integer& mod,
                    const typename IntVector::Integer& muBarrett) {
          typename IntVector::Integer result(z.ModSub(x * y, mod, muBarrett));
          return result;
      };
  
      IntVector runningDividend(dividend);
  
      usint divisorPtr;
      for (usint i = 0; i < runs; i++) {
          typename IntVector::Integer divConst(runningDividend.at(dividendLength - 1));  // get the highest degree coeff
          divisorPtr = divisorLength - 1;
          for (usint j = 0; j < dividendLength - i - 1; j++) {
              if (divisorPtr > j) {
                  runningDividend.at(dividendLength - 1 - j) = mat(
                      divisor.at(divisorPtr - 1 - j), divConst, runningDividend.at(dividendLength - 2 - j), modulus, mu);
              }
              else {
                  runningDividend.at(dividendLength - 1 - j) = runningDividend.at(dividendLength - 2 - j);
              }
          }
      }
  
      for (usint i = 0, j = runs; i < divisorLength - 1; i++, j++) {
          result.at(i) = runningDividend.at(j);
      }
  
      return result;
  }
  
  template <typename IntVector>
  IntVector PolynomialMultiplication(const IntVector& a, const IntVector& b) {
      usint degreeA = a.GetLength() - 1;
      usint degreeB = b.GetLength() - 1;
  
      usint degreeResultant = degreeA + degreeB;
  
      const auto& modulus = a.GetModulus();
  
      IntVector result(degreeResultant + 1, modulus);
  
      for (usint i = 0; i < a.GetLength(); i++) {
          for (usint j = 0; j < b.GetLength(); j++) {
              const auto& valResult = result.at(i + j);
              const auto& valMult   = a.at(i) * b.at(j);
              result.at(i + j)      = (valMult + valResult).Mod(modulus);
          }
      }
  
      return result;
  }
  
  template <typename IntVector>
  IntVector GetCyclotomicPolynomial(usint m, const typename IntVector::Integer& modulus) {
      auto intCP = GetCyclotomicPolynomialRecursive(m);
      IntVector result(intCP.size(), modulus);
      for (usint i = 0; i < intCP.size(); i++) {
          auto val = intCP.at(i);
          if (intCP.at(i) > -1) {
              result.at(i) = typename IntVector::Integer(val);
          }
          else {
              val *= -1;
              result.at(i) = modulus - typename IntVector::Integer(val);
          }
      }
  
      return result;
  }
  
  template <typename IntVector>
  typename IntVector::Integer SyntheticRemainder(const IntVector& dividend, const typename IntVector::Integer& a,
                                                 const typename IntVector::Integer& modulus) {
      auto val = dividend.at(dividend.GetLength() - 1);
  
      // Precompute the Barrett mu parameter
      auto mu = modulus.ComputeMu();
  
      for (int i = dividend.GetLength() - 2; i > -1; i--) {
          val = dividend.at(i) + a * val;
          val = val.Mod(modulus, mu);
      }
  
      return val;
  }
  
  template <typename IntVector>
  IntVector SyntheticPolyRemainder(const IntVector& dividend, const IntVector& aList,
                                   const typename IntVector::Integer& modulus) {
      IntVector result(aList.GetLength(), modulus);
      for (usint i = 0; i < aList.GetLength(); i++) {
          result.at(i) = SyntheticRemainder(dividend, aList.at(i), modulus);
      }
  
      return result;
  }
  
  template <typename IntVector>
  IntVector PolynomialPower(const IntVector& input, usint power) {
      usint finalDegree = (input.GetLength() - 1) * power;
      IntVector finalPoly(finalDegree + 1, input.GetModulus());
      finalPoly.at(0) = input.at(0);
      for (usint i = 1; i < input.GetLength(); i++) {
          finalPoly.at(i * power) = input.at(i);
      }
      return finalPoly;
  }
  
  template <typename IntVector>
  IntVector SyntheticPolynomialDivision(const IntVector& dividend, const typename IntVector::Integer& a,
                                        const typename IntVector::Integer& modulus) {
      usint n = dividend.GetLength() - 1;
      IntVector result(n, modulus);
  
      auto mu = modulus.ComputeMu();
  
      result.at(n - 1) = dividend.at(n);
      auto val(dividend.at(n));
      for (int i = n - 1; i > 0; i--) {
          val              = val * a + dividend.at(i);
          val              = val.Mod(modulus, mu);
          result.at(i - 1) = val;
      }
  
      return result;
  }
  
  }  // namespace lbcrypto
  