powerful.h

/* Copyright (C) 2012-2020 IBM Corp.
 * This program is Licensed under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance
 * with the License. You may obtain a copy of the License at
 *   http://www.apache.org/licenses/LICENSE-2.0
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License. See accompanying LICENSE file.
 */
#ifndef HELIB_POWERFUL_H
#define HELIB_POWERFUL_H
#include <helib/NumbTh.h>
#include <helib/ClonedPtr.h>
#include <helib/bluestein.h>
#include <helib/hypercube.h>
#include <helib/DoubleCRT.h>
#include <helib/Context.h>
 
namespace helib {
 
class PowerfulTranslationIndexes
{
public:
  long m;                // product of mi's
  long phim;             // phi(m) = product of phi(mi)'s
  NTL::Vec<long> mvec;   // mvec[i] = mi
  NTL::Vec<long> phivec; // phivec[i] = phi(mi)
  NTL::Vec<long> divvec; // divvec[i] = m/mi
  NTL::Vec<long> invvec; // invvec[i] = (m/mi)^{-1} mod mi
 
  CubeSignature longSig;  // hypercube of dimensions mi
  CubeSignature shortSig; // hypercube of dimensions phi(mi)
 
  NTL::Vec<long> polyToCubeMap; // index translation tables
  NTL::Vec<long> cubeToPolyMap;
  NTL::Vec<long> shortToLongMap;
 
  NTL::Vec<NTL::ZZX> cycVec; // cycvec[i] = Phi_mi(X)
  NTL::ZZX phimX;
 
  PowerfulTranslationIndexes(const NTL::Vec<long>& mv);
};
 
class PowerfulConversion
{
  const PowerfulTranslationIndexes* indexes;
  NTL::zz_pContext zzpContext;
  NTL::Vec<NTL::zz_pXModulus> cycVec_p; // cycvec[i] = Phi_mi(X)
  NTL::zz_pXModulus phimX_p;
 
public:
  PowerfulConversion() : indexes(nullptr) {}
 
  explicit PowerfulConversion(const PowerfulTranslationIndexes& ind) :
      indexes(nullptr)
  {
    initPConv(ind);
  }
 
  void initPConv(const PowerfulTranslationIndexes& ind)
  {
    if (indexes != nullptr)
      return; // cannot re-initialize a non-nullptr object
    indexes = &ind;
 
    cycVec_p.SetLength(ind.cycVec.length());
    zzpContext.save(); // store the current modulus
    for (long i = 0; i < ind.cycVec.length(); i++) {
      cycVec_p[i] =
          NTL::conv<NTL::zz_pX>(ind.cycVec[i]); // convert to zz_pXModulus
    }
    phimX_p = NTL::conv<NTL::zz_pX>(ind.phimX); // convert to zz_pXModulus
  }
 
  void restoreModulus() const { zzpContext.restore(); }
  const CubeSignature& getLongSig() const { return indexes->longSig; }
  const CubeSignature& getShortSig() const { return indexes->shortSig; }
 
  long powerfulToPoly(NTL::zz_pX& poly,
                      const HyperCube<NTL::zz_p>& powerful) const;
  long polyToPowerful(HyperCube<NTL::zz_p>& powerful,
                      const NTL::zz_pX& poly) const;
};
 
class PowerfulDCRT
{
  const Context& context; // points to the context for the DoubleCRT's
 
  PowerfulTranslationIndexes indexes; // modulus-independent tables
 
  // a vector of PowerfulConversion tables, one for each modulus
  NTL::Vec<PowerfulConversion> pConvVec;
 
  // product_bits[i] is the number of bits in the product of primes [0..i)
  NTL::Vec<long> product_bits;
 
  // number of excess bits needed to ensure correct conversion
  long to_pwfl_excess_bits;
  long to_poly_excess_bits;
 
  bool triv;
 
public:
  PowerfulDCRT(const Context& _context, const NTL::Vec<long>& mvec);
 
  const PowerfulTranslationIndexes& getIndexTranslation() const
  {
    return indexes;
  }
  const PowerfulConversion& getPConv(long i) const { return pConvVec.at(i); }
 
  // coefficients are reduced to the interval [-Q/2,Q/2], where
  // Q = product of primes in dcrt.getIndexSet();
  void dcrtToPowerful(NTL::Vec<NTL::ZZ>& powerful, const DoubleCRT& dcrt) const;
 
  void ZZXtoPowerful(NTL::Vec<NTL::ZZ>& powerful, const NTL::ZZX& poly) const;
  void powerfulToZZX(NTL::ZZX& poly, const NTL::Vec<NTL::ZZ>& powerful) const;
};
 
/********************************************************************/
/****************    UNUSED CODE - COMMENTED OUT   ******************/
/********************************************************************/
 
#if 0
long computeProd(const Vec< Pair<long, long> >& vec);
 
inline long computePhi(const Pair<long, long>& x)
{
  long p = x.a;
  long e = x.b;
  return power_long(p, e - 1) * (p-1);
}
 
void computePhiVec(Vec<long>& phiVec,
                   const Vec< Pair<long, long> >& factors);
 
void computeDivVec(Vec<long>& divVec, long m,
                   const Vec<long>& powVec);
 
void computeInvVec(Vec<long>& invVec,
                   const Vec<long>& divVec, const Vec<long>& powVec);
 
 
void computeShortToLongMap(Vec<long>& shortToLongMap,
                           const CubeSignature& shortSig,
                           const CubeSignature& longSig);
 
void computeLongToShortMap(Vec<long>& longToShortMap,
                           long m,
                           const Vec<long>& shortToLongMap);
 
void recursiveReduce(const CubeSlice<zz_p>& s,
                     const Vec<zz_pXModulus>& cycVec,
                     long d,
                     zz_pX& tmp1,
                     zz_pX& tmp2);
 
void convertPolyToPowerful(HyperCube<zz_p>& cube,
                           HyperCube<zz_p>& tmpCube,
                           const zz_pX& poly,
                           const Vec<zz_pXModulus>& cycVec,
                           const Vec<long>& polyToCubeMap,
                           const Vec<long>& shortToLongMap);
 
void convertPowerfulToPoly(zz_pX& poly,
                           const HyperCube<zz_p>& cube,
                           long m,
                           const Vec<long>& shortToLongMap,
                           const Vec<long>& cubeToPolyMap,
                           const zz_pXModulus& phimX);
#endif
/********************************************************************/
#if 0
inline long computePow(const Pair<long, long>& x)
{
  long p = x.a;
  long e = x.b;
  return power_long(p, e);
}
 
void computePowVec(Vec<long>& powVec,
                   const Vec< Pair<long, long> >& factors);
 
void mapIndexToPowerful(Vec<long>& pow, long j, const Vec<long>& phiVec);
 
void computeCycVec(Vec<zz_pXModulus>& cycVec, const Vec<long>& powVec);
 
void computePowerToCubeMap(Vec<long>& polyToCubeMap,
                           Vec<long>& cubeToPolyMap,
                           long m,
                           const Vec<long>& powVec,
                           const Vec<long>& invVec,
                           const CubeSignature& longSig);
 
void computeMultiEvalPoints(Vec< Vec<zz_p> >& multiEvalPoints,
                            const zz_p& base,
                            long m,
                            const Vec<long>& powVec,
                            const Vec<long>& phiVec);
 
void computeLinearEvalPoints(Vec<zz_p>& linearEvalPoints,
                             const zz_p& base,
                             long m, long phim);
 
void computeCompressedIndex(Vec< Vec<long> >& compressedIndex,
                            const Vec<long>& powVec);
 
void computePowToCompressedIndexMap(Vec<long>& powToCompressedIndexMap,
                                    long m,
                                    const Vec<long>& powVec,
                                    const Vec< Vec<long> >& compressedIndex,
                                    const CubeSignature& shortSig);
 
void recursiveEval(const CubeSlice<zz_p>& s,
                   const Vec< Vec<zz_p> >& multiEvalPoints,
                   long d,
                   zz_pX& tmp1,
                   Vec<zz_p>& tmp2);
 
inline void eval(HyperCube<zz_p>& cube,
         const Vec< Vec<zz_p> >& multiEvalPoints)
{
   zz_pX tmp1;
   Vec<zz_p> tmp2;
 
   recursiveEval(CubeSlice<zz_p>(cube), multiEvalPoints, 0, tmp1, tmp2);
}
 
void mapPowerfulToPoly(ZZX& poly,
                       const Vec<long>& pow,
                       const Vec<long>& divVec,
                       long m,
                       const ZZX& phimX);
#endif
/********************************************************************/
#if 0
class FFTHelper {
 
private:
  long m;
  zz_p m_inv;
  zz_p root, iroot;
  zz_pXModulus phimx;
  Vec<bool> coprime;
  long phim;
 
  mutable zz_pX powers, ipowers;
  mutable Vec<mulmod_precon_t> powers_aux, ipowers_aux;
  mutable fftRep Rb, iRb;
  mutable fftrep_aux Rb_aux, iRb_aux;
  mutable fftRep Ra;
  mutable zz_pX tmp;
 
public:
  FFTHelper(long _m, zz_p x);
 
  void FFT(const zz_pX& f, Vec<zz_p>& v) const;
    // compute v = { f(x^i) }_{i in Z_m^*}
  void iFFT(zz_pX& f, const Vec<zz_p>& v, bool normalize = true) const;
    // computes inverse transform. If !normalize, result is scaled by m
 
  const zz_p& get_m_inv() const { return m_inv; }
  // useful for aggregate normalization
 
};
 
 
void computeMultiEvalPoints(Vec< CopiedPtr<FFTHelper> >& multiEvalPoints,
                            const zz_p& base,
                            long m,
                            const Vec<long>& powVec,
                            const Vec<long>& phiVec);
 
void recursiveEval(const CubeSlice<zz_p>& s,
                   const Vec< CopiedPtr<FFTHelper> >& multiEvalPoints,
                   long d,
                   zz_pX& tmp1,
                   Vec<zz_p>& tmp2);
 
inline void eval(HyperCube<zz_p>& cube,
          const Vec< CopiedPtr<FFTHelper> >& multiEvalPoints)
{
   zz_pX tmp1;
   Vec<zz_p> tmp2;
 
   recursiveEval(CubeSlice<zz_p>(cube), multiEvalPoints, 0, tmp1, tmp2);
}
 
void recursiveInterp(const CubeSlice<zz_p>& s,
                     const Vec< CopiedPtr<FFTHelper> >& multiEvalPoints,
                     long d,
                     zz_pX& tmp1,
                     Vec<zz_p>& tmp2);
 
void interp(HyperCube<zz_p>& cube,
        const Vec< CopiedPtr<FFTHelper> >& multiEvalPoints);
#endif
 
} // namespace helib
 
#endif // ifndef HELIB_POWERFUL_H