PAlgebra.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_PALGEBRA_H
#define HELIB_PALGEBRA_H
#include <exception>
#include <utility>
#include <vector>
#include <complex>
 
#include <helib/NumbTh.h>
#include <helib/zzX.h>
#include <helib/hypercube.h>
#include <helib/PGFFT.h>
#include <helib/ClonedPtr.h>
#include <helib/apiAttributes.h>
 
namespace helib {
 
struct half_FFT
{
  PGFFT fft;
  std::vector<std::complex<double>> pow;
 
  half_FFT(long m);
};
 
struct quarter_FFT
{
  PGFFT fft;
  std::vector<std::complex<double>> pow1, pow2;
 
  quarter_FFT(long m);
};
 
class PAlgebra
{
  long m; // the integer m defines (Z/mZ)^*, Phi_m(X), etc.
  long p; // the prime base of the plaintext space
 
  long phiM;      // phi(m)
  long ordP;      // the order of p in (Z/mZ)^*
  long nfactors;  // number of distinct prime factors of m
  long radm;      // rad(m) = prod of distinct primes dividing m
  double normBnd; // max-norm-on-pwfl-basis <= normBnd * max-norm-canon-embed
  double polyNormBnd; // max-norm-on-poly-basis <= polyNormBnd *
                      // max-norm-canon-embed
 
  long pow2; // if m = 2^k, then pow2 == k; otherwise, pow2 == 0
 
  std::vector<long> gens; // Our generators for (Z/mZ)^* (other than p)
 
  //  native[i] is true iff gens[i] has the same order in the quotient
  //  group as its order in Zm*.
  NTL::Vec<bool> native;
 
  // frob_perturb[i] = j if gens[i] raised to its order equals p^j,
  // otherwise -1
  NTL::Vec<long> frob_perturb;
 
  CubeSignature cube; // the hypercube structure of Zm* /(p)
 
  NTL::ZZX PhimX; // Holds the integer polynomial Phi_m(X)
 
  double cM; // the "ring constant" c_m for Z[X]/Phi_m(X)
  // NOTE: cM is related to the ratio between the l_infinity norm of
  // a "random" ring element in different bases. For example, think of
  // choosing the power-basis coefficients of x uniformly at random in
  // [+-a/2] (for some parameter a), then the powerful basis norm of x
  // should be bounded whp by cM*a.
  //
  // More precisely, for an element x whose coefficients are chosen
  // uniformly in [+-a/2] (in either the powerful or the power basis)
  // we have a high-probability bound |x|_canonical < A*a for some
  // A = O(sqrt(phi(m)). Also for "random enough" x we have some bound
  //       |x|_powerful < |x|_canonical * B
  // where we "hope" that B = O(1/sqrt(phi(m)). The cM value is
  // supposed to be cM=A*B.
  //
  // The value cM is only used for bootstrapping, see more comments
  // for the method RecryptData::setAE in recryption.cpp. Also see
  // Appendix A of https://ia.cr/2014/873 (updated version from 2019)
 
  std::vector<long> T;    // The representatives for the quotient group Zm* /(p)
  std::vector<long> Tidx; // i=Tidx[t] is the index i s.t. T[i]=t.
                          // Tidx[t]==-1 if t notin T
 
  std::vector<long> zmsIdx; // if t is the i'th element in Zm* then zmsIdx[t]=i
                            // zmsIdx[t]==-1 if t notin Zm*
 
  std::vector<long> zmsRep; // inverse of zmsIdx
 
  std::shared_ptr<PGFFT> fftInfo; // info for computing m-point complex FFT's
                                  // shard_ptr allows delayed initialization
                                  // and lightweight copying
 
  std::shared_ptr<half_FFT> half_fftInfo;
  // an optimization for FFT's with even m
 
  std::shared_ptr<quarter_FFT> quarter_fftInfo;
  // an optimization for FFT's with m = 0 (mod 4)
 
public:
  PAlgebra& operator=(const PAlgebra&) = delete;
 
  PAlgebra(long mm,
           long pp = 2,
           const std::vector<long>& _gens = std::vector<long>(),
           const std::vector<long>& _ords = std::vector<long>()); // constructor
 
  bool operator==(const PAlgebra& other) const;
  bool operator!=(const PAlgebra& other) const { return !(*this == other); }
  // comparison
 
  /* I/O methods */
 
  void printout(std::ostream& out = std::cout) const;
  void printAll(std::ostream& out = std::cout) const; // print even more
 
  /* Access methods */
 
  long getM() const { return m; }
 
  long getP() const { return p; }
 
  long getPhiM() const { return phiM; }
 
  long getOrdP() const { return ordP; }
 
  long getNFactors() const { return nfactors; }
 
  long getRadM() const { return radm; }
 
  double getNormBnd() const { return normBnd; }
 
  double getPolyNormBnd() const { return polyNormBnd; }
 
  long getNSlots() const { return cube.getSize(); }
 
  long getPow2() const { return pow2; }
 
  const NTL::ZZX& getPhimX() const { return PhimX; }
 
  double get_cM() const { return cM; }
 
  //  const std::vector<long> getMfactors() const { return mFactors; }
 
  long numOfGens() const { return gens.size(); }
 
  long ZmStarGen(long i) const { return (i < long(gens.size())) ? gens[i] : 0; }
 
  // VJS: I'm moving away from all of this unsigned stuff...
  // Also, note that j really may be negative
  // NOTE: i == -1 means Frobenius
  long genToPow(long i, long j) const;
 
  // p to the power j mod m
  long frobeniusPow(long j) const;
 
  long OrderOf(long i) const { return cube.getDim(i); }
 
  long ProdOrdsFrom(long i) const { return cube.getProd(i); }
 
  bool SameOrd(long i) const { return native[i]; }
 
  // FrobPerturb[i] = j if gens[i] raised to its order equals p^j,
  // where j in [0..ordP), otherwise -1
  long FrobPerturb(long i) const { return frob_perturb[i]; }
 
 
  long ith_rep(long i) const { return (i < getNSlots()) ? T[i] : 0; }
 
  long indexOfRep(long t) const { return (t > 0 && t < m) ? Tidx[t] : -1; }
 
  bool isRep(long t) const { return (t > 0 && t < m && Tidx[t] > -1); }
 
  long indexInZmstar(long t) const { return (t > 0 && t < m) ? zmsIdx[t] : -1; }
 
  long indexInZmstar_unchecked(long t) const { return zmsIdx[t]; }
 
  long repInZmstar_unchecked(long idx) const { return zmsRep[idx]; }
 
  bool inZmStar(long t) const { return (t > 0 && t < m && zmsIdx[t] > -1); }
 
  long exponentiate(const std::vector<long>& exps,
                    bool onlySameOrd = false) const;
 
  long coordinate(long i, long k) const { return cube.getCoord(k, i); }
 
  std::pair<long, long> breakIndexByDim(long idx, long dim) const
  {
    return cube.breakIndexByDim(idx, dim);
  }
 
  long assembleIndexByDim(std::pair<long, long> idx, long dim) const
  {
    return cube.assembleIndexByDim(idx, dim);
  }
 
  long addCoord(long i, long k, long offset) const
  {
    return cube.addCoord(k, i, offset);
  }
 
  /* Miscellaneous */
 
  bool nextExpVector(std::vector<long>& exps) const
  {
    return cube.incrementCoords(exps);
  }
 
  long fftSizeNeeded() const { return NTL::NextPowerOfTwo(getM()) + 1; }
  // TODO: should have a special case when m is power of two
 
  const PGFFT& getFFTInfo() const { return *fftInfo; }
  const half_FFT& getHalfFFTInfo() const { return *half_fftInfo; }
  const quarter_FFT& getQuarterFFTInfo() const { return *quarter_fftInfo; }
};
 
enum PA_tag
{
  PA_GF2_tag,
  PA_zz_p_tag,
  PA_cx_tag
};
 
class DummyBak
{
  // placeholder class used in GF2X impl
 
public:
  void save() {}
  void restore() const {}
};
 
class DummyContext
{
  // placeholder class used in GF2X impl
 
public:
  void save() {}
  void restore() const {}
  DummyContext() {}
  DummyContext(long) {}
};
 
class DummyModulus
{};
// placeholder class for CKKS
 
// some stuff to help with template code
template <typename R>
struct GenericModulus
{};
 
template <>
struct GenericModulus<NTL::zz_p>
{
  static void init(long p) { NTL::zz_p::init(p); }
};
 
template <>
struct GenericModulus<NTL::GF2>
{
  static void init(long p)
  {
    assertEq<InvalidArgument>(p, 2l, "Cannot init NTL::GF2 with p not 2");
  }
};
 
class PA_GF2
{
  // typedefs for algebraic structures built up from GF2
 
public:
  static const PA_tag tag = PA_GF2_tag;
  typedef NTL::GF2 R;
  typedef NTL::GF2X RX;
  typedef NTL::vec_GF2X vec_RX;
  typedef NTL::GF2XModulus RXModulus;
  typedef DummyBak RBak;
  typedef DummyContext RContext;
  typedef NTL::GF2E RE;
  typedef NTL::vec_GF2E vec_RE;
  typedef NTL::mat_GF2E mat_RE;
  typedef NTL::GF2EX REX;
  typedef NTL::GF2EBak REBak;
  typedef NTL::vec_GF2EX vec_REX;
  typedef NTL::GF2EContext REContext;
  typedef NTL::mat_GF2 mat_R;
  typedef NTL::vec_GF2 vec_R;
};
 
class PA_zz_p
{
  // typedefs for algebraic structures built up from zz_p
 
public:
  static const PA_tag tag = PA_zz_p_tag;
  typedef NTL::zz_p R;
  typedef NTL::zz_pX RX;
  typedef NTL::vec_zz_pX vec_RX;
  typedef NTL::zz_pXModulus RXModulus;
  typedef NTL::zz_pBak RBak;
  typedef NTL::zz_pContext RContext;
  typedef NTL::zz_pE RE;
  typedef NTL::vec_zz_pE vec_RE;
  typedef NTL::mat_zz_pE mat_RE;
  typedef NTL::zz_pEX REX;
  typedef NTL::zz_pEBak REBak;
  typedef NTL::vec_zz_pEX vec_REX;
  typedef NTL::zz_pEContext REContext;
  typedef NTL::mat_zz_p mat_R;
  typedef NTL::vec_zz_p vec_R;
};
 
class PA_cx
{
  // typedefs for algebraic structures built up from complex<double>
 
public:
  static const PA_tag tag = PA_cx_tag;
  typedef double R;
  typedef std::complex<double> RX;
  typedef NTL::Vec<RX> vec_RX;
  typedef DummyModulus RXModulus;
  typedef DummyBak RBak;
  typedef DummyContext RContext;
 
  // the other typedef's should not ever be used...they
  // are all defined as void, so that PA_INJECT still works
  typedef void RE;
  typedef void vec_RE;
  typedef void mat_RE;
  typedef void REX;
  typedef void REBak;
  typedef void vec_REX;
  typedef void REContext;
  typedef void mat_R;
  typedef void vec_R;
};
 
 
class PAlgebraModBase
{
 
public:
  virtual ~PAlgebraModBase() {}
  virtual PAlgebraModBase* clone() const = 0;
 
  virtual PA_tag getTag() const = 0;
 
  virtual const PAlgebra& getZMStar() const = 0;
 
  virtual const std::vector<NTL::ZZX>& getFactorsOverZZ() const = 0;
 
  virtual long getR() const = 0;
 
  virtual long getPPowR() const = 0;
 
  virtual void restoreContext() const = 0;
 
  virtual zzX getMask_zzX(long i, long j) const = 0;
};
 
#ifndef DOXYGEN_IGNORE
#define PA_INJECT(typ)                                                         \
  static const PA_tag tag = typ::tag;                                          \
  typedef typename typ::R R;                                                   \
  typedef typename typ::RX RX;                                                 \
  typedef typename typ::vec_RX vec_RX;                                         \
  typedef typename typ::RXModulus RXModulus;                                   \
  typedef typename typ::RBak RBak;                                             \
  typedef typename typ::RContext RContext;                                     \
  typedef typename typ::RE RE;                                                 \
  typedef typename typ::vec_RE vec_RE;                                         \
  typedef typename typ::mat_RE mat_RE;                                         \
  typedef typename typ::REX REX;                                               \
  typedef typename typ::REBak REBak;                                           \
  typedef typename typ::vec_REX vec_REX;                                       \
  typedef typename typ::REContext REContext;                                   \
  typedef typename typ::mat_R mat_R;                                           \
  typedef typename typ::vec_R vec_R;
 
#endif
 
template <typename type>
class PAlgebraModDerived;
// forward declaration
 
template <typename type>
class MappingData
{
 
public:
  PA_INJECT(type)
 
  friend class PAlgebraModDerived<type>;
 
private:
  RX G;      // the polynomial defining the field extension
  long degG; // the degree of the polynomial
 
  REContext contextForG;
 
  /* the remaining fields are visible only to PAlgebraModDerived */
 
  std::vector<RX> maps;
  std::vector<mat_R> matrix_maps;
  std::vector<REX> rmaps;
 
public:
  const RX& getG() const { return G; }
  long getDegG() const { return degG; }
  void restoreContextForG() const { contextForG.restore(); }
 
  // copy and assignment
};
 
template <typename T>
class TNode
{
public:
  std::shared_ptr<TNode<T>> left, right;
  T data;
 
  TNode(std::shared_ptr<TNode<T>> _left,
        std::shared_ptr<TNode<T>> _right,
        const T& _data) :
      left(_left), right(_right), data(_data)
  {}
};
 
template <typename T>
std::shared_ptr<TNode<T>> buildTNode(std::shared_ptr<TNode<T>> left,
                                     std::shared_ptr<TNode<T>> right,
                                     const T& data)
{
  return std::shared_ptr<TNode<T>>(new TNode<T>(left, right, data));
}
 
template <typename T>
std::shared_ptr<TNode<T>> nullTNode()
{
  return std::shared_ptr<TNode<T>>();
}
 
template <typename type>
class PAlgebraModDerived : public PAlgebraModBase
{
public:
  PA_INJECT(type)
 
private:
  const PAlgebra& zMStar;
  long r;
  long pPowR;
  RContext pPowRContext;
 
  RXModulus PhimXMod;
 
  vec_RX factors;
  std::vector<NTL::ZZX> factorsOverZZ;
  vec_RX crtCoeffs;
  std::vector<std::vector<RX>> maskTable;
  std::vector<RX> crtTable;
  std::shared_ptr<TNode<RX>> crtTree;
 
  void genMaskTable();
  void genCrtTable();
 
public:
  PAlgebraModDerived& operator=(const PAlgebraModDerived&) = delete;
 
  PAlgebraModDerived(const PAlgebra& zMStar, long r);
 
  PAlgebraModDerived(const PAlgebraModDerived& other) // copy constructor
      :
      zMStar(other.zMStar),
      r(other.r),
      pPowR(other.pPowR),
      pPowRContext(other.pPowRContext)
  {
    RBak bak;
    bak.save();
    restoreContext();
    PhimXMod = other.PhimXMod;
    factors = other.factors;
    maskTable = other.maskTable;
    crtTable = other.crtTable;
    crtTree = other.crtTree;
  }
 
  virtual PAlgebraModBase* clone() const override
  {
    return new PAlgebraModDerived(*this);
  }
 
  virtual PA_tag getTag() const override { return tag; }
 
  virtual const PAlgebra& getZMStar() const override { return zMStar; }
 
  virtual const std::vector<NTL::ZZX>& getFactorsOverZZ() const override
  {
    return factorsOverZZ;
  }
 
  virtual long getR() const override { return r; }
 
  virtual long getPPowR() const override { return pPowR; }
 
  virtual void restoreContext() const override { pPowRContext.restore(); }
 
  /* In all of the following functions, it is expected that the caller
     has already restored the relevant modulus (p^r), which
     can be done by invoking the method restoreContext()
   */
 
  const RXModulus& getPhimXMod() const { return PhimXMod; }
 
  const vec_RX& getFactors() const { return factors; }
 
  const vec_RX& getCrtCoeffs() const { return crtCoeffs; }
 
  // logically, but not really, const
  const std::vector<std::vector<RX>>& getMaskTable() const { return maskTable; }
 
  zzX getMask_zzX(long i, long j) const override
  {
    RBak bak;
    bak.save();
    restoreContext();
    return balanced_zzX(maskTable.at(i).at(j));
  }
 
 
  void CRT_decompose(std::vector<RX>& crt, const RX& H) const;
 
  void CRT_reconstruct(RX& H, std::vector<RX>& crt) const;
 
  void mapToSlots(MappingData<type>& mappingData, const RX& G) const;
 
  void embedInAllSlots(RX& H,
                       const RX& alpha,
                       const MappingData<type>& mappingData) const;
 
  void embedInSlots(RX& H,
                    const std::vector<RX>& alphas,
                    const MappingData<type>& mappingData) const;
 
  void decodePlaintext(std::vector<RX>& alphas,
                       const RX& ptxt,
                       const MappingData<type>& mappingData) const;
 
  void buildLinPolyCoeffs(std::vector<RX>& C,
                          const std::vector<RX>& L,
                          const MappingData<type>& mappingData) const;
private:
  /* internal functions, not for public consumption */
 
  static void SetModulus(long p)
  {
    RContext context(p);
    context.restore();
  }
 
  void mapToF1(RX& w, const RX& G) const { mapToFt(w, G, 1); }
 
  void mapToFt(RX& w, const RX& G, long t, const RX* rF1 = nullptr) const;
 
  void buildTree(std::shared_ptr<TNode<RX>>& res,
                 long offset,
                 long extent) const;
 
  void evalTree(RX& res,
                std::shared_ptr<TNode<RX>> tree,
                const std::vector<RX>& crt1,
                long offset,
                long extent) const;
};
 
template <>
class PAlgebraModDerived<PA_cx> : public PAlgebraModBase
{
  const PAlgebra& zMStar;
  long r; // counts bits of precision
 
public:
  PAlgebraModDerived(const PAlgebra& palg, long _r) : zMStar(palg), r(_r)
  {
    assertInRange<InvalidArgument>(r,
                                   1l,
                                   (long)NTL_SP_NBITS,
                                   "Invalid bit precision r");
  }
 
  PAlgebraModBase* clone() const override
  {
    return new PAlgebraModDerived(*this);
  }
  PA_tag getTag() const override { return PA_cx_tag; }
 
  const PAlgebra& getZMStar() const override { return zMStar; }
  long getR() const override { return r; }
  long getPPowR() const override { return 1L << r; }
  void restoreContext() const override {}
 
  // These function make no sense for PAlgebraModCx
  const std::vector<NTL::ZZX>& getFactorsOverZZ() const override
  {
    throw LogicError("PAlgebraModCx::getFactorsOverZZ undefined");
  }
  zzX getMask_zzX(UNUSED long i, UNUSED long j) const override
  {
    throw LogicError("PAlgebraModCx::getMask_zzX undefined");
  }
};
 
typedef PAlgebraModDerived<PA_cx> PAlgebraModCx;
 
PAlgebraModBase* buildPAlgebraMod(const PAlgebra& zMStar, long r);
 
// A simple wrapper for a pointer to an object of type PAlgebraModBase.
//
// Direct access to the virtual methods of PAlgebraModBase is provided,
// along with a "downcast" operator to get a reference to the object
// as a derived type, and == and != operators.
class PAlgebraMod
{
 
private:
  ClonedPtr<PAlgebraModBase> rep;
 
public:
  // copy constructor: default
  // assignment: deleted
  // destructor: default
  // NOTE: the use of ClonedPtr ensures that the default copy constructor
  // and destructor will work correctly.
 
  PAlgebraMod& operator=(const PAlgebraMod&) = delete;
 
  explicit PAlgebraMod(const PAlgebra& zMStar, long r) :
      rep(buildPAlgebraMod(zMStar, r))
  {}
  // constructor
 
  template <typename type>
  const PAlgebraModDerived<type>& getDerived(type) const
  {
    return dynamic_cast<const PAlgebraModDerived<type>&>(*rep);
  }
  const PAlgebraModCx& getCx() const
  {
    return dynamic_cast<const PAlgebraModCx&>(*rep);
  }
 
  bool operator==(const PAlgebraMod& other) const
  {
    return getZMStar() == other.getZMStar() && getR() == other.getR();
  }
  // comparison
 
  bool operator!=(const PAlgebraMod& other) const { return !(*this == other); }
  // comparison
 
  /* direct access to the PAlgebraModBase methods */
 
  PA_tag getTag() const { return rep->getTag(); }
  const PAlgebra& getZMStar() const { return rep->getZMStar(); }
  const std::vector<NTL::ZZX>& getFactorsOverZZ() const
  {
    return rep->getFactorsOverZZ();
  }
  long getR() const { return rep->getR(); }
  long getPPowR() const { return rep->getPPowR(); }
  void restoreContext() const { rep->restoreContext(); }
 
  zzX getMask_zzX(long i, long j) const { return rep->getMask_zzX(i, j); }
};
 
bool comparePAlgebra(const PAlgebra& palg,
                     unsigned long m,
                     unsigned long p,
                     unsigned long r,
                     const std::vector<long>& gens,
                     const std::vector<long>& ords);
 
// for internal consumption only
double calcPolyNormBnd(long m);
 
} // namespace helib
 
#endif // #ifndef HELIB_PALGEBRA_H