inline void
MakeLegendre( Matrix<F>& A )
{
#ifndef RELEASE
    CallStackEntry entry("MakeLegendre");
#endif
    if( A.Height() != A.Width() )
        LogicError("Cannot make a non-square matrix Legendre");
    MakeZeros( A );

    const Int n = A.Width();
    for( Int j=0; j<n-1; ++j )
    {
        const F gamma = F(1) / Pow( F(2)*(j+1), F(2) );
        const F beta = F(1) / (2*Sqrt(F(1)-gamma));
        A.Set( j+1, j, beta );
        A.Set( j, j+1, beta );
    }
}

void
Ledger::rawTxInsert (uint256 const& key,
                     std::shared_ptr<Serializer const
                     > const& txn, std::shared_ptr<
                     Serializer const> const& metaData)
{
    assert (metaData);

    // low-level - just add to table
    Serializer s(txn->getDataLength () +
                 metaData->getDataLength () + 16);
    s.addVL (txn->peekData ());
    s.addVL (metaData->peekData ());
    auto item = std::make_shared<
                SHAMapItem const> (key, std::move(s));
    if (! txMap().addGiveItem
            (std::move(item), true, true))
        LogicError("duplicate_tx: " + to_string(key));
}

void CharPolyMod(ZZ_pX& g, const ZZ_pX& a, const ZZ_pX& ff)
{
   ZZ_pX f = ff;
   MakeMonic(f);
   long n = deg(f);

   if (n <= 0 || deg(a) >= n) 
      LogicError("CharPoly: bad args");

   if (IsZero(a)) {
      clear(g);
      SetCoeff(g, n);
      return;
   }

   if (n > 25) {
      ZZ_pX h;
      MinPolyMod(h, a, f);
      if (deg(h) == n) {
         g = h;
         return;
      }
   }

   if (ZZ_p::modulus() < n+1) {
      HessCharPoly(g, a, f);
      return;
   }

   vec_ZZ_p u(INIT_SIZE, n+1), v(INIT_SIZE, n+1);

   ZZ_pX h, h1;
   negate(h, a);
   long i;

   for (i = 0; i <= n; i++) {
      u[i] = i;
      add(h1, h, u[i]);
      resultant(v[i], f, h1);
   }

   interpolate(g, u, v);
}

inline LDLPivot
SelectFromPanel
( const Matrix<F>& A,
  const Matrix<F>& X,
  const Matrix<F>& Y, 
  LDLPivotType pivotType,
  Base<F> gamma )
{
    DEBUG_CSE
    LDLPivot pivot;
    switch( pivotType )
    {
    case BUNCH_KAUFMAN_A: 
    case BUNCH_KAUFMAN_C: pivot = PanelBunchKaufmanA( A, X, Y, gamma ); break;
    case BUNCH_KAUFMAN_D: pivot = PanelBunchKaufmanD( A, X, Y, gamma ); break;
    default: LogicError("This pivot type not yet supported");
    }
    return pivot;
}

// get path of current executable
/*static*/ wstring File::GetExecutablePath()
{
#ifdef WIN32
    wchar_t path[33000];
    if (GetModuleFileNameW(NULL, path, _countof(path)) == 0)
        LogicError("GetExecutablePath: GetModuleFileNameW() unexpectedly failed.");
    return path;
#else
    // from http://stackoverflow.com/questions/4025370/can-an-executable-discover-its-own-path-linux
    pid_t pid = getpid();
    char path[PATH_MAX + 1] = { 0 };
    sprintf(path, "/proc/%d/exe", pid);
    char dest[PATH_MAX + 1] = { 0 };
    if (readlink(path, dest, PATH_MAX) == -1)
        RuntimeError("GetExecutableDirectory: readlink() call failed.");
    else
        return msra::strfun::utf16(dest);
#endif
}

void EDF(vec_ZZ_pEX& factors, const ZZ_pEX& ff, const ZZ_pEX& bb,
         long d, long verbose)

{
   ZZ_pEX f = ff;
   ZZ_pEX b = bb;

   if (!IsOne(LeadCoeff(f)))
      LogicError("EDF: bad args");

   long n = deg(f);
   long r = n/d;

   if (r == 0) {
      factors.SetLength(0);
      return;
   }

   if (r == 1) {
      factors.SetLength(1);
      factors[0] = f;
      return;
   }

   if (d == 1) {
      RootEDF(factors, f, verbose);
      return;
   }

   
   double t;
   if (verbose) { 
      cerr << "computing EDF(" << d << "," << r << ")..."; 
      t = GetTime(); 
   }

   factors.SetLength(0);

   RecEDF(factors, f, b, d, verbose);

   if (verbose) cerr << (GetTime()-t) << "\n";
}

void FastTraceVec(vec_ZZ_p& S, const ZZ_pX& f)
{
   long n = deg(f);

   if (n <= 0) 
      LogicError("FastTraceVec: bad args");

   if (n == 0) {
      S.SetLength(0);
      return;
   }

   if (n == 1) {
      S.SetLength(1);
      set(S[0]);
      return;
   }
   
   long i;
   ZZ_pX f1;

   f1.rep.SetLength(n-1);
   for (i = 0; i <= n-2; i++)
      f1.rep[i] = f.rep[n-i];
   f1.normalize();

   ZZ_pX f2;
   f2.rep.SetLength(n-1);
   for (i = 0; i <= n-2; i++)
      mul(f2.rep[i], f.rep[n-1-i], i+1);
   f2.normalize();

   ZZ_pX f3;
   InvTrunc(f3, f1, n-1);
   MulTrunc(f3, f3, f2, n-1);

   S.SetLength(n);

   S[0] = n;
   for (i = 1; i < n; i++)
      negate(S[i], coeff(f3, i-1));
}

void plain_mul_transpose_aux(mat_ZZ_p& X, const mat_ZZ_p& A, const mat_ZZ_p& B)  
{  
   long n = A.NumRows();  
   long l = A.NumCols();  
   long m = B.NumRows();  
  
   if (l != B.NumCols())  
      LogicError("matrix mul: dimension mismatch");  
  
   X.SetDims(n, m);  

   ZZ_pContext context;
   context.save();

   long sz = ZZ_p::ModulusSize();
   bool seq = (double(n)*double(l)*double(m)*double(sz)*double(sz) < PAR_THRESH);
  
   NTL_GEXEC_RANGE(seq, m, first, last)
   NTL_IMPORT(n)
   NTL_IMPORT(l)
   NTL_IMPORT(m)

   context.restore();

   long i, j, k;  
   ZZ acc, tmp;  

   for (j = first; j < last; j++) {
      const ZZ_p *B_col = B[j].elts();

      for (i = 0; i < n; i++) {
         clear(acc);
         for (k = 0; k < l; k++) {
            mul(tmp, rep(A[i][k]), rep(B_col[k]));
            add(acc, acc, tmp);
         }
         conv(X[i][j], acc);
      }
   }

   NTL_GEXEC_RANGE_END
}  

inline void
DistNodalMultiVec<F>::Pull
( const DistMap& inverseMap, const DistSymmInfo& info,
  const DistMultiVec<F>& X )
{
    DEBUG_ONLY(CallStackEntry cse("DistNodalMultiVec::Pull"))
    height_ = X.Height();
    width_ = X.Width();

    // Traverse our part of the elimination tree to see how many indices we need
    int numRecvInds=0;
    const int numLocal = info.localNodes.size();
    for( int s=0; s<numLocal; ++s )
        numRecvInds += info.localNodes[s].size;
    const int numDist = info.distNodes.size();
    for( int s=1; s<numDist; ++s )
        numRecvInds += info.distNodes[s].multiVecMeta.localSize;
    
    // Fill the set of indices that we need to map to the original ordering
    int off=0;
    std::vector<int> mappedInds( numRecvInds );
    for( int s=0; s<numLocal; ++s )
    {
        const SymmNodeInfo& nodeInfo = info.localNodes[s];
        for( int t=0; t<nodeInfo.size; ++t )
            mappedInds[off++] = nodeInfo.off+t;
    }
    for( int s=1; s<numDist; ++s )
    {
        const DistSymmNodeInfo& nodeInfo = info.distNodes[s];
        const Grid& grid = *nodeInfo.grid;
        const int gridSize = grid.Size();
        const int gridRank = grid.VCRank();
        const int alignment = 0;
        const int shift = Shift( gridRank, alignment, gridSize );
        for( int t=shift; t<nodeInfo.size; t+=gridSize )
            mappedInds[off++] = nodeInfo.off+t;
    }
    DEBUG_ONLY(
        if( off != numRecvInds )
            LogicError("mappedInds was filled incorrectly");
    )

Base<Field> LanczosDecomp
( const SparseMatrix<Field>& A,
        Matrix<Field>& V,
        Matrix<Base<Field>>& T,
        Matrix<Field>& v,
        Int basisSize )
{
    EL_DEBUG_CSE
    const Int n = A.Height();
    if( n != A.Width() )
        LogicError("A was not square");

    auto applyA =
      [&]( const Matrix<Field>& X, Matrix<Field>& Y )
      {
          Zeros( Y, n, X.Width() );
          Multiply( NORMAL, Field(1), A, X, Field(0), Y );
      };
    return LanczosDecomp( n, applyA, V, T, v, basisSize );
}

void BuildIrred(ZZ_pEX& f, long n)
{
   if (n <= 0)
      LogicError("BuildIrred: n must be positive");

   if (n == 1) {
      SetX(f);
      return;
   }

   ZZ_pEX g;

   do {
      random(g, n);
      SetCoeff(g, n);
   } while (!IterIrredTest(g));

   f = g;

}

XmlDomElement* FootprintPad::serializeToXmlDomElement() const throw (Exception)
{
    if (!checkAttributesValidity()) throw LogicError(__FILE__, __LINE__);

    QScopedPointer<XmlDomElement> root(new XmlDomElement("pad"));
    root->setAttribute("uuid", mUuid);
    root->setAttribute("type", typeToString(mType));
    root->setAttribute("x", mPosition.getX().toMmString());
    root->setAttribute("y", mPosition.getY().toMmString());
    root->setAttribute("rotation", mRotation);
    root->setAttribute("width", mWidth);
    root->setAttribute("height", mHeight);
    root->setAttribute("drill", mDrillDiameter);
    root->setAttribute("layer", mLayerId);
    foreach (const QString& locale, mNames.keys())
        root->appendTextChild("name", mNames.value(locale))->setAttribute("locale", locale);
    foreach (const QString& locale, mDescriptions.keys())
        root->appendTextChild("description", mDescriptions.value(locale))->setAttribute("locale", locale);
    return root.take();
}

void UndoStack::abortCommand() throw (Exception)
{
    Q_ASSERT(mCurrentIndex == mCommands.count());

    if (!mCommandActive)
        throw LogicError(__FILE__, __LINE__, QString(), tr("No command active!"));

    mCommands.last()->undo(); // throws an exception on error
    mCurrentIndex--;
    mCommandActive = false;
    delete mCommands.takeLast(); // delete and remove the aborted command from the stack

    // emit signals
    emit undoTextChanged(getUndoText());
    emit redoTextChanged(tr("Redo"));
    emit canUndoChanged(canUndo());
    emit canRedoChanged(false);
    emit cleanChanged(isClean());
    emit commandAborted(); // this is important!
}

Base<F> LanczosDecomp
( const DistSparseMatrix<F>& A,
        DistMultiVec<F>& V, 
        ElementalMatrix<Base<F>>& T,
        DistMultiVec<F>& v,
        Int basisSize )
{
    DEBUG_CSE
    const Int n = A.Height();
    if( n != A.Width() )
        LogicError("A was not square");

    auto applyA =
      [&]( const DistMultiVec<F>& X, DistMultiVec<F>& Y )
      {
          Zeros( Y, n, X.Width() );
          Multiply( NORMAL, F(1), A, X, F(0), Y );
      };
    return LanczosDecomp( n, applyA, V, T, v, basisSize );
}

void ComponentSignalInstance::init() throw (Exception)
{
    // create ERC messages
    mErcMsgUnconnectedRequiredSignal.reset(new ErcMsg(mCircuit.getProject(), *this,
        QString("%1/%2").arg(mComponentInstance.getUuid().toStr()).arg(mComponentSignal->getUuid().toStr()),
        "UnconnectedRequiredSignal", ErcMsg::ErcMsgType_t::CircuitError, QString()));
    mErcMsgForcedNetSignalNameConflict.reset(new ErcMsg(mCircuit.getProject(), *this,
        QString("%1/%2").arg(mComponentInstance.getUuid().toStr()).arg(mComponentSignal->getUuid().toStr()),
        "ForcedNetSignalNameConflict", ErcMsg::ErcMsgType_t::SchematicError, QString()));
    updateErcMessages();

    // register to component attributes changed
    connect(&mComponentInstance, &ComponentInstance::attributesChanged,
            this, &ComponentSignalInstance::updateErcMessages);

    // register to net signal name changed
    if (mNetSignal) connect(mNetSignal, &NetSignal::nameChanged, this, &ComponentSignalInstance::netSignalNameChanged);

    if (!checkAttributesValidity()) throw LogicError(__FILE__, __LINE__);
}

inline void
MakeGKS( Matrix<F>& A )
{
#ifndef RELEASE
    CallStackEntry entry("MakeGKS");
#endif
    const Int m = A.Height();
    const Int n = A.Width();
    if( m != n )
        LogicError("Cannot make a non-square matrix GKS");

    MakeZeros( A );
    for( Int j=0; j<n; ++j )
    {
        const F jDiag = F(1)/Sqrt(F(j));
        for( Int i=0; i<j; ++i )
            A.Set( i, j, -jDiag );
        A.Set( j, j, jDiag );
    }
}

inline 
Grid::Grid( mpi::Comm comm, int height, GridOrder order )
: haveViewers_(false), order_(order)
{
    DEBUG_ONLY(CallStackEntry cse("Grid::Grid"))

    // Extract our rank, the underlying group, and the number of processes
    mpi::Dup( comm, viewingComm_ );
    mpi::CommGroup( viewingComm_, viewingGroup_ );
    size_ = mpi::Size( viewingComm_ );

    // All processes own the grid, so we have to trivially split viewingGroup_
    owningGroup_ = viewingGroup_;

    height_ = height;
    if( height_ < 0 )
        LogicError("Process grid dimensions must be non-negative");

    SetUpGrid();
}

void InnerProduct(ZZ_p& x, const vec_ZZ_p& a, const vec_ZZ_p& b,
                  long offset)
{
    if (offset < 0) LogicError("InnerProduct: negative offset");
    if (NTL_OVERFLOW(offset, 1, 0))
        ResourceError("InnerProduct: offset too big");

    long n = min(a.length(), b.length()+offset);
    long i;
    NTL_ZZRegister(accum);
    NTL_ZZRegister(t);

    clear(accum);
    for (i = offset; i < n; i++) {
        mul(t, rep(a[i]), rep(b[i-offset]));
        add(accum, accum, t);
    }

    conv(x, accum);
}

void CompositeDataReader::StartEpoch(const EpochConfiguration& cfg)
{
    EpochConfiguration config = cfg;

    if (config.m_totalEpochSizeInSamples <= 0)
    {
        RuntimeError("Unsupported Epoch size '%d'.", (int)config.m_totalEpochSizeInSamples);
    }

    m_sequenceEnumerator->StartEpoch(config);

    // TODO: As the next step the packers should be moved into the network.
    switch (m_packingMode)
    {
    case PackingMode::sample:
        m_packer = std::make_shared<FramePacker>(
            m_provider,
            m_sequenceEnumerator,
            m_streams);
        break;
    case PackingMode::sequence:
        m_packer = std::make_shared<SequencePacker>(
            m_provider,
            m_sequenceEnumerator,
            m_streams);
        break;
    case PackingMode::truncated:
    {
        config.m_truncationSize = m_truncationLength;
        m_packer = std::make_shared<TruncatedBPTTPacker>(
            m_provider,
            m_sequenceEnumerator,
            m_streams);
        break;
    }
    default:
        LogicError("Unsupported type of packer '%d'.", (int)m_packingMode);
    }

    m_packer->StartEpoch(config);
}

void PowerCompose(ZZ_pEX& y, const ZZ_pEX& h, long q, const ZZ_pEXModulus& F)
{
   if (q < 0) LogicError("PowerCompose: bad args");

   ZZ_pEX z(INIT_SIZE, F.n);
   long sw;

   z = h;
   SetX(y);

   while (q) {
      sw = 0;

      if (q > 1) sw = 2;
      if (q & 1) {
         if (IsX(y))
            y = z;
         else
            sw = sw | 1;
      }

      switch (sw) {
      case 0:
         break;

      case 1:
         CompMod(y, y, z, F);
         break;

      case 2:
         CompMod(z, z, z, F);
         break;

      case 3:
         Comp2Mod(y, z, y, z, z, F);
         break;
      }

      q = q >> 1;
   }
}