Context.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_CONTEXT_H
#define HELIB_CONTEXT_H
 
#include <helib/PAlgebra.h>
#include <helib/CModulus.h>
#include <helib/IndexSet.h>
#include <helib/recryption.h>
#include <helib/primeChain.h>
#include <helib/powerful.h>
#include <helib/apiAttributes.h>
#include <helib/range.h>
 
#include <NTL/Lazy.h>
 
namespace helib {
 
inline double lweEstimateSecurity(int n, double log2AlphaInv, int hwt)
{
  if (hwt < 0 || (hwt > 0 && hwt < 120)) {
    throw LogicError("Cannot estimate security for keys of weight<120");
  }
 
  // clang-format off
  constexpr double hwgts[] =
      {120, 150, 180, 210, 240, 270, 300, 330, 360, 390, 420, 450};
  constexpr double slopes[] =
      {2.4, 2.67, 2.83, 3.0, 3.1, 3.3, 3.3, 3.35, 3.4, 3.45, 3.5, 3.55};
  constexpr double cnstrms[] =
      {19, 13, 10, 6, 3, 1, -3, -4, -5, -7, -10, -12};
  // clang-format on
 
  constexpr size_t numWghts = sizeof(hwgts) / sizeof(hwgts[0]);
 
  const size_t idx = (hwt - 120) / 30; // index into the array above
  double slope = 0, consterm = 0;
  if (hwt == 0) { // dense keys
    slope = 3.8;
    consterm = -20;
  } else if (idx < numWghts - 1) { 
    // estimate prms on a line from prms[i] to prms[i+1]
    // how far into this interval
    double a = double(hwt - hwgts[idx]) / (hwgts[idx + 1] - hwgts[idx]);
    slope = slopes[idx] + a * (slopes[idx + 1] - slopes[idx]);
    consterm = cnstrms[idx] + a * (cnstrms[idx + 1] - cnstrms[idx]);
  } else { 
    // Use the params corresponding to largest weight (450 above)
    slope = slopes[numWghts - 1];
    consterm = cnstrms[numWghts - 1];
  }
 
  double x = n / log2AlphaInv;
  double ret = slope * x + consterm;
 
  return ret < 0.0 ? 0.0 : ret; // If ret is negative then return 0.0
}
 
long FindM(long k,
           long nBits,
           long c,
           long p,
           long d,
           long s,
           long chosen_m,
           bool verbose = false);
 
class EncryptedArray;
struct PolyModRing;
class Context
{
  std::vector<Cmodulus> moduli; // Cmodulus objects for the different primes
  // This is private since the implementation assumes that the list of
  // primes only grows and no prime is ever modified or removed.
 
public:
  // Context is meant for convenience, not encapsulation: Most data
  // members are public and can be initialized by the application program.
 
  PAlgebra zMStar;
 
  PAlgebraMod alMod;
 
  std::shared_ptr<const EncryptedArray> ea;
 
  std::shared_ptr<const PowerfulDCRT> pwfl_converter;
 
  std::shared_ptr<PolyModRing> slotRing;
 
  NTL::xdouble stdev;
 
  //======================= high probability bounds ================
  double scale; // default = 10
 
 
 
  //=======================================
 
 
  // NOTE: this is a bit heuristic: we assume that if we evaluate
  // f at a primitive root of unity, then we get something that well
  // approximates a normal random variable with the same variance,
  // which is equal to the sum of the variances of the individual
  // f_i's, which is (2*magBound)^2/12 = magBound^2/3.
  // We then multiply the sqrt of the variance by scale to get
  // the high probability bound.
 
  double noiseBoundForUniform(double magBound, long degBound) const
  {
    return scale * std::sqrt(double(degBound) / 3.0) * magBound;
  }
 
  NTL::xdouble noiseBoundForUniform(NTL::xdouble magBound, long degBound) const
  {
    return scale * std::sqrt(double(degBound) / 3.0) * magBound;
  }
 
  //=======================================
 
 
  // NOTE: for odd modulus, this means each f_i is uniformly distributed
  // over { -floor(modulus/2), ..., floor(modulus/2) }.
  // For even modulus, this means each f_i is uniformly distributed
  // over { modulus/2, ..., modulus/2 }, except that the two endpoints
  // (which represent the same residue class) occur with half the
  // probability of the others.
 
  // NOTE: this is a bit heuristic: we assume that if we evaluate
  // f at a primitive root of unity, then we get something that well
  // approximates a normal random variable with the same variance,
  // which is equal to the sum of the variances of the individual
  // f_i's, which is (modulus)^2/12 + 1/6 for even modulus,
  // and is at most (modulus^2)/12 for odd modulus.
  // We then multiply the sqrt of the variance by scale to get
  // the high probability bound.
 
  // NOTE: this is slightly more accurate that just calling
  // noiseBoundForUniform with magBound=modulus/2.
 
  double noiseBoundForMod(long modulus, long degBound) const
  {
    double var = fsquare(modulus) / 12.0;
    if (modulus % 2 == 0)
      var += 1.0 / 6.0;
 
    return scale * std::sqrt(degBound * var);
  }
 
  //=======================================
 
 
  // NOTE: if we evaluate f at a primitive root of unity,
  // then we get a normal random variable variance degBound * sigma^2.
  // We then multiply the sqrt of the variance by scale to get
  // the high probability bound.
 
  double noiseBoundForGaussian(double sigma, long degBound) const
  {
    return scale * std::sqrt(double(degBound)) * sigma;
  }
 
  //=======================================
 
 
  // NOTE: this is a bit heuristic: we assume that if we evaluate
  // f at a primitive root of unity, then we get something that
  // well approximates a normal random variable with the same variance,
  // which is equal to the sum of the individual variances,
  // which is degBound*prob.
  // We then multiply the sqrt of the variance by scale to get
  // the high probability bound.
 
  double noiseBoundForSmall(double prob, long degBound) const
  {
    return scale * std::sqrt(double(degBound)) * std::sqrt(prob);
  }
 
  //=======================================
 
 
  // NOTE: this is a bit heuristic: we assume that if we evaluate
  // f at a primitive root of unity, then we get something that
  // well approximates a normal random variable with the same variance,
  // which is hwt.
  // We then multiply the sqrt of the variance by scale to get
  // the high probability bound.
 
  // NOTE: degBound is not used here, but I include it
  // for consistency with the other noiseBound routines
 
  double noiseBoundForHWt(long hwt, UNUSED long degBound) const
  {
    return scale * std::sqrt(double(hwt));
  }
 
  //=======================================
 
 
 
 
  double stdDevForRecryption(long skHwt = 0) const
  {
    if (!skHwt)
      skHwt = rcData.skHwt;
    // the default reverts to rcData.skHwt, *not* rcData.defSkHwt
 
    long k = zMStar.getNFactors();
    // number of prime factors of m
 
    long m = zMStar.getM();
    long phim = zMStar.getPhiM();
 
    double mrat = double(phim) / double(m);
 
    return std::sqrt(mrat * double(skHwt) * double(1L << k) / 3.0) * 0.5;
  }
 
  double boundForRecryption(long skHwt = 0) const
  {
    double c_m = zMStar.get_cM();
    // multiply by this fudge factor
 
    return 0.5 + c_m * scale * stdDevForRecryption(skHwt);
  }
 
  IndexSet ctxtPrimes;
 
  IndexSet specialPrimes;
 
  IndexSet smallPrimes;
 
  ModuliSizes modSizes;
  void setModSizeTable() { modSizes.init(*this); }
 
  // Digits of ctxt/columns of key-switching matrix
  std::vector<IndexSet> digits;
 
  // includes both thin and thick
  ThinRecryptData rcData;
 
  /******************************************************************/
  // constructor
  Context(unsigned long m,
          unsigned long p,
          unsigned long r,
          const std::vector<long>& gens = std::vector<long>(),
          const std::vector<long>& ords = std::vector<long>());
 
  void makeBootstrappable(const NTL::Vec<long>& mvec,
                          long skWht = 0,
                          bool build_cache = false,
                          bool alsoThick = true)
  {
    rcData.init(*this, mvec, alsoThick, skWht, build_cache);
  }
 
  bool isBootstrappable() const { return rcData.alMod != nullptr; }
 
  IndexSet fullPrimes() const { return ctxtPrimes | specialPrimes; }
 
  IndexSet allPrimes() const
  {
    return smallPrimes | ctxtPrimes | specialPrimes;
  }
 
  // returns first nprimes ctxtPrimes
  IndexSet getCtxtPrimes(long nprimes) const
  {
    long first = ctxtPrimes.first();
    long last = std::min(ctxtPrimes.last(), first + nprimes - 1);
    return IndexSet(first, last);
  }
 
  // FIXME: replacement for bitsPerLevel...placeholder for now
  long BPL() const { return 30; }
 
  bool operator==(const Context& other) const;
  bool operator!=(const Context& other) const { return !(*this == other); }
 
  long ithPrime(unsigned long i) const
  {
    return (i < moduli.size()) ? moduli[i].getQ() : 0;
  }
 
  const Cmodulus& ithModulus(unsigned long i) const { return moduli[i]; }
 
  long numPrimes() const { return moduli.size(); }
 
  bool isZeroDivisor(const NTL::ZZ& num) const
  {
    for (long i : range(moduli.size()))
      if (divide(num, moduli[i].getQ()))
        return true;
    return false;
  }
 
  bool inChain(long p) const
  {
    for (long i : range(moduli.size()))
      if (p == moduli[i].getQ())
        return true;
    return false;
  }
 
  void productOfPrimes(NTL::ZZ& p, const IndexSet& s) const;
  NTL::ZZ productOfPrimes(const IndexSet& s) const
  {
    NTL::ZZ p;
    productOfPrimes(p, s);
    return p;
  }
 
  // FIXME: run-time error when ithPrime(i) returns 0
  double logOfPrime(unsigned long i) const { return log(ithPrime(i)); }
 
  double logOfProduct(const IndexSet& s) const
  {
    if (s.last() >= numPrimes())
      throw RuntimeError("Context::logOfProduct: IndexSet has too many rows");
 
    double ans = 0.0;
    for (long i : s)
      ans += logOfPrime(i);
    return ans;
  }
 
  long bitSizeOfQ() const
  {
    IndexSet primes = ctxtPrimes | specialPrimes;
    return std::ceil(logOfProduct(primes) / log(2.0));
  }
 
  double securityLevel(int hwt = 0) const
  {
    IndexSet primes = ctxtPrimes | specialPrimes;
    if (primes.card() == 0) {
      throw LogicError(
          "Security level cannot be determined as modulus chain is empty.");
    }
 
    double s = to_double(stdev);
    if (zMStar.getPow2() == 0) { // not power of two
      s *= sqrt(zMStar.getM());
    }
    double log2AlphaInv = (logOfProduct(primes) - log(s)) / log(2.0);
    return lweEstimateSecurity(zMStar.getPhiM(), log2AlphaInv, hwt);
  }
 
  void printout(std::ostream& out = std::cout) const;
 
  void AddSmallPrime(long q);
  void AddCtxtPrime(long q);
  void AddSpecialPrime(long q);
 
 
  friend void writeContextBase(std::ostream& str, const Context& context);
 
  friend std::ostream& operator<<(std::ostream& str, const Context& context);
 
  friend void readContextBase(std::istream& str,
                              unsigned long& m,
                              unsigned long& p,
                              unsigned long& r,
                              std::vector<long>& gens,
                              std::vector<long>& ords);
 
  friend std::istream& operator>>(std::istream& str, Context& context);
 
  friend void writeContextBinary(std::ostream& str, const Context& context);
  friend void readContextBinary(std::istream& str, Context& context);
};
 
void writeContextBase(std::ostream& s, const Context& context);
void readContextBase(std::istream& s,
                     unsigned long& m,
                     unsigned long& p,
                     unsigned long& r,
                     std::vector<long>& gens,
                     std::vector<long>& ords);
std::unique_ptr<Context> buildContextFromAscii(std::istream& str);
 
void writeContextBaseBinary(std::ostream& str, const Context& context);
void writeContextBinary(std::ostream& str, const Context& context);
 
void readContextBaseBinary(std::istream& s,
                           unsigned long& m,
                           unsigned long& p,
                           unsigned long& r,
                           std::vector<long>& gens,
                           std::vector<long>& ords);
std::unique_ptr<Context> buildContextFromBinary(std::istream& str);
void readContextBinary(std::istream& str, Context& context);
 
// Build modulus chain with nBits worth of ctxt primes,
// using nDgts digits in key-switching.
// The willBeBootstrappable and skHwt parameters are needed to get around some
// some circularity when making the context boostrappable.
// If you later call context.makeBootstrappable with a given value
// of skHwt, you should first buildModChain with willBeBootstrappable
// set to true and the given value of skHwt.
// FIXME: We should really have a simpler way to do this.
// resolution ... FIXME
 
void buildModChain(Context& context,
                   long nBits,
                   long nDgts = 3,
                   bool willBeBootstrappable = false,
                   long skHwt = 0,
                   long resolution = 3,
                   long bitsInSpecialPrimes = 0);
// should be called if after you build the mod chain in some way
// *other* than calling buildModChain.
void endBuildModChain(Context& context);
 
// Should point to the "current" context
extern Context* activeContext;
 
} // namespace helib
 
#endif // ifndef HELIB_CONTEXT_H