//? cpmalgorithms.cpp
//? C+- by Ulrich Mutze. Status of work 2022-05-06.
//? Copyright (c) 2022 Ulrich Mutze
//? contact: see contact-info at www.ulrichmutze.de
//?
//? This program is free software: you can redistribute it and/or
//? modify it under the terms of the GNU General Public License as
//? published by the Free Software Foundation, either version 3 of
//? the License, or (at your option) any later version.
//?
//? This program is distributed in the hope that it will be useful,
//? but WITHOUT ANY WARRANTY; without even the implied warranty of
//? MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//? GNU General Public License <http://www.gnu.org/licenses/> for
//? more details.


// notice that some algorithms are templates and are defined not here
// but in cpmalgorithms.h !!!!!!!!!!!!!!!!!!!!

#include <cpmva.h>
#include <cpmalgorithms.h>
#include <cpmmacros.h>

   // the names from CpmRootX will be used via using declarations in
   // function blocks.

using namespace std;

using CpmRoot::R;
using CpmRoot::Z;
using CpmRoot::N;
using CpmRoot::L;
using CpmRoot::Word;
using CpmRoot::B;
using CpmArrays::V;
using CpmArrays::VV;
using CpmArrays::Vo;
using CpmArrays::Va;
using CpmArrays::X2;
using CpmArrays::X3;
using CpmFunctions::F;
#ifdef CPM_Fn
using CpmFunctions::F1;
using CpmFunctions::F2;
using CpmFunctions::F3;
using CpmFunctions::F4;
#endif
using CpmFunctions::FncObj;
using CpmGeo::Iv;
using CpmAlgorithms::LatticeBox;
using CpmLinAlg::Z2;
   // these using declarations are active till the end of this file, also
   // within the anonymus namespace to be defined later.
   // The directive #include <cpmalgorithms.h> brings also some using
   // declarations in action. These are only valid within the namespace
   // CpmAlgorithms and  n o t  in the anonymus namespace mentioned
   // above.
   // They are also not valid as return values of CpmAlgorithms member
   // functions (only in argument lists and definition blocks)


N CpmAlgorithms::lesDiv(N const& n)
{
   if (n<=2) return n;
   N res=2;
   while (n%res) res++;
   return res;
}

V<N> CpmAlgorithms::prmFac(N const& n, bool verify)
{
   V<N> res;
   if (n<=2) res=V<N>(1,n);
   else{
      N na=n;
      N f=CpmAlgorithms::lesDiv(na);
      res<<f;
      while(f<na){
         na/=f;
         f=CpmAlgorithms::lesDiv(na);
         res<<f;
      }
   }
   if (verify){
      N nt=1;
      for (Z i=res.b();i<=res.e();++i) nt*=res[i];
      if (nt!=n) cpmerror("CpmAlgorithms::prmFac("&cpm(n)&") failed");
   }
   return res;
}

N CpmAlgorithms::gcd(N const& n1, N const& n2, bool verify)
{
   N a,b;
   if (n1>=n2){ a=n1; b=n2; }
   else{ a=n2; b=n1;}
   if (b==0) return 0;
   N r=a%b;
   while(r){
      a=b;
      b=r;
      r=a%b;
   }
   if (verify){
      if (n1%b || n2%b){
         cpmerror("CpmAlgorithms::gcd("&cpm(n1)&","&cpm(n2)&") failed");
      }
   }
   return b;
}

N CpmAlgorithms::lcm(N const& a, N const& b)
{
   N g=CpmAlgorithms::gcd(a,b);
   return g==0 ? 0 : ((a/g)*(b/g))*g;
   // avoids forming temporaries that are larger than the result
}

bool CpmAlgorithms::isPrm(N const& n)
{
   if (n<2) return false;
   return CpmAlgorithms::lesDiv(n)==n;
}

namespace{
   Z pow2(Z n)
   {
      if (n<0) cpmerror("pow2(Z): argument <0");
      Z res=1;
      for (Z i=1;i<=n;i++) res*=2;
      return res;
   }
}

Z CpmAlgorithms::ZfromBits(V<B> const& vb)
   // res=sum i=1,...vb.dim vb[i]*2^(i-1)
{
   Z mL=3;
   Word loc("CpmAlgorithms::ZfromBits(V<bool>)");
   cpmmessage(mL,loc&" started");
   Z n=vb.dim();
   // we need 2^(n-1) as the first sepatator
   cpmassert(n>0,loc);
   Z d=pow2(n-1);
   Z res=0;
   for (Z i=1;i<=n;i++){
      if (vb[i]==true) res+=d;
      d/=2;
   }
   cpmmessage(mL,loc&" done");
   return res;
}

V<B> CpmAlgorithms::bitsFromZ(Z n, Z p)
   // we assume and test 0<=n<=2^(p-1) and return vb such that
   // ZfromBits(vb)==n
{
   Z mL=3;
   Word loc("CpmAlgorithms::bitsFromZ(Z n, Z p)");
   cpmmessage(mL,loc&" started");

   cpmassert(n>=0,loc);
   cpmassert(p>0,loc);
   cpmassert(p<32,loc);

   Z d=pow2(p);
      // e.g. p=8, d=256
   Z limit=d-1; // 255
   cpmassert(n<=limit,loc);
   V<B> res(p);
   Z separator=d/2; // 128
   d=separator; // 128
   for (Z i=1;i<=p;i++){
      d/=2; // 64
      if (n<separator){ //x<=127
         separator-=d;
      }
      else{
         separator+=d;
         res[i]=B(true);
      }
   }
   cpmmessage(mL,loc&" done");
   return res;
}

bool CpmAlgorithms::next_1(std::vector<B>& v)
// should be very efficient
{
   using std::vector;
   size_t n=v.size();
   if (n==0) return false;
   size_t i=n-1;
   while (v[i]==true){
      --i;
      if (i>n || i<0 ) return false; // then all entries were true
      // notice that size_t may not have negative values
   }
   // now v[i]==false
   v[i]=true; // so a change was made and return value
      // is rightfully 'true'
   for (size_t j=i+1; j<n; ++j) v[j]=false;
      // componets right to the point of action (if there are
      // such objects) get set to false.
   return true;
}

V<L> CpmAlgorithms::digFromR(R const& x)
{
   Word wx=Word::write(x);
   Z n= wx.dim();
   V<L> res(n);
   for (Z i=1;i<=n;++i) res[i]=wx[i];
   return res;
}

///////////////////// class LatticeEnumerator //////////////////////////

void CpmAlgorithms::LatticeEnumerator::update()
{
   Z mL=3;
   Word loc("CpmAlgorithms::LatticeEnumerator::update()");
   cpmmessage(mL,loc&" started");
   n=b.dim();
   cpmassert(n>0,loc);
   Z i,errorSum=0;
   for (i=1;i<=n;i++) errorSum+=(b[i]<1 ? 1 : 0);
   cpmassert(errorSum==0,loc);
   c=V<Z>(n,1);
   for (i=1;i<n;i++) c[n-i]=c[n-i+1]*b[n-i+1];
   nLattice=c[1]*b[1];
   cpmmessage(mL,loc&" done");
}

Z CpmAlgorithms::LatticeEnumerator::operator()(const V<Z>& x)const
{
   Z mL=3;
   Word loc("CpmAlgorithms::LatticeEnumerator::eval(V<Z>)");
   cpmmessage(mL,loc&" started");
   Z i;
   cpmassert(x.dim()==n,loc);
   Z errorSum=0;
   for (i=1;i<=n;i++) errorSum+=(x[i]<1 || x[i]>b[i] ? 1 : 0);
   cpmassert(errorSum==0,loc);
   Z res=1;
   for (i=1;i<=n;i++) res+=c[i]*(x[i]-1);
   cpmmessage(mL,loc&" done");
   return res;
}

V<Z> CpmAlgorithms::LatticeEnumerator::operator()(Z k)const
{
   Z mL=3;
   Word loc("CpmAlgorithms::LatticeEnumerator::eval(Z)");
   cpmmessage(mL,loc&" started");
   cpmassert(k>0 && k<=nLattice,loc);
   k--;
   V<Z> res(n);
   Z i=1;
   while (i<=n){
      Z ri=k/c[i];
      k-=ri*c[i];
      res[i]=ri+1;
      i++;
   }
   cpmmessage(mL,loc&" done");
   return res;
}

bool CpmAlgorithms::LatticeEnumerator::prnOn(ostream& str)const
{
   CpmRoot::writeTitle("CpmAlgorithms::LatticeEnumerator",str);
   V<Z> b_(b);
   cpmp(b_);
   return true;
}

bool CpmAlgorithms::LatticeEnumerator::scanFrom(istream& str)
{
   V<Z> b_;
   cpms(b_);
   b=b_;
   update();
   return true;
}


///////////////////// class LatticeBox /////////////////////////////////

LatticeBox::LatticeBox(const V<Iv>& vivIn, const V<Z>& vzIn):
viv(vivIn),vz(vzIn)
{
   Z mL=3;
   Word loc("LatticeBox(V<Iv>,V<Z>)");
   CPM_MA
   using CpmAlgorithms::LatticeEnumerator;
   n=viv.dim();
   cpmassert(n==vz.dim(),loc);
   Z i,j;
   np=1;
   R zero(0.);
   for (i=1;i<=n;i++) np*=vz[i];
   val=V<R>(np,zero);
   centers=V< V<R> >(np,V<R>(n,zero));
   le=LatticeEnumerator(vz);
   for (i=1;i<=np;i++){
      V<Z> lei=le(i);
      V<R> vali(n);
      for (j=1;j<=n;j++) vali[j]=viv[j].centers(vz[j])[lei[j]];
      centers[i]=vali;
   }
   CPM_MZ
}

void LatticeBox::set(const V<R>& pos, R valIn)
{
   V<Z> zdim(n);
   Z i;
   for (i=1;i<=n;i++){
      zdim[i]=viv[i].quant(pos[i],vz[i]);
   }
   Z k=le(zdim);
   val[k]+=valIn;
}

VV<R> LatticeBox::matrix()const
{
   Z mL=3;
   Word loc("LatticeBox::matrix()");
   CPM_MA
   cpmassert(n>=2,loc);
   Z mres=vz[2];
   Z nres=vz[1]; // sic

   VV<R> res(mres,nres);
   Z i,j,k;
   for (k=1;k<=np;k++){
      V<Z> vij=le(k);
      Z ix=vij[1];
      Z iy=vij[2];
      i=1+mres-iy;
      j=ix;
      res[i][j]=val[k];
   }
   CPM_MZ
   return res;
}

namespace{ // anonymous namespace for tools

   R xLimit;
      // a number for communication between class xf and function brack

   struct xf{ // for encapsulating some functionality
      R x1,x2,x3;
      R f1,f2,f3;
      xf(R x, R f):x1(x),x2(x),x3(x),f1(f),f2(f),f3(f){}
      void shift(R x, R f){x1=x2;f1=f2;x2=x3;f2=f3;x3=x,f3=f;}
      Z eval()const;
         // the return value encodes the meaning of the data of *this
         // with respect to the aims of brack()
   };

   inline bool signChange(R x, R y)
   {
      return  (x>=0 && y>=0)||(x<=0 && y<=0) ? false : true;
   }

   Z xf::eval()const
   {
      const R fac=2;
      using CpmRootX::inf;
      if (x3>xLimit) return 1;
      if (f3==0) return 2;
      if (signChange(f2,f3)) return 3;
      R d1=f2-f1;
      R d2=f3-f2;
      if (signChange(d1,d2)){
         R var=cpmabs(d1)+cpmabs(d2);
         R mini=inf(cpmabs(f1),cpmabs(f2),cpmabs(f3));
         if (mini>fac*var) return 0; else return 4;
      }
      return 0; // no special event, no action triggered
   }

} // anonymous namespace

X3<R,R,Z> CpmAlgorithms::brack(const F<R,R>& f, R xMin,
   R xMax, R dxMin, R dxMax)
{
   const Z mL=4;

   Z found;
   R x1=xMin;
   R length=xMax-xMin;
   R nh=cpmfloor(length/dxMax+0.5);
   R hStart=length/nh; // after doing some steps hStart one should
      // come very close to xMax
   R h=hStart;
   xLimit=xMax; // x1,h,xLimit are the iter ini data
   bool local=false;
iter:
   R f1=f(x1); // x1 needed
   if (cpmverbose>mL) cpmdata<<endl<<"f1="<<f1;
   xf state(x1,f1);
   Z res=0;
   while (res==0){
      R xn=state.x3+h; // h needed
      R fn=f(xn);
      state.shift(xn,fn);
      res=state.eval(); // xLimit needed (hidden)
   }
   // now res!=0
   if (res==1 ){
      if ((!local) || state.x3>=xMax){
         found=0;
         return X3<R,R,Z>(state.x2,state.x3,found);
      }
      else{
         x1=state.x3; // start where you ended, progress
         h=hStart;
         xLimit=xMax;
         local=false;
         goto iter;
      }
   }
   else if(res==2){
      found=1;
      return X3<R,R,Z>(state.x3,state.x3,found);

   }
   else if(res==3){
      found=2;
      return X3<R,R,Z>(state.x2,state.x3,found);
   }
   else{ // then res==4, we found a change in trend
      // which can be an indication that a local scan may find a
      // root between state.x1 and state.x3
      if (local){
         cpmmessage(mL,
            "CpmAlgorithms::brack(): trend change and no root");
         x1=state.x3; // start where you ended, progress
         h=hStart;
         xLimit=xMax;
         local=false;
      }
      else{
         cpmmessage(mL,
            "CpmAlgorithms::brack(): refinement on trend change");
         x1=state.x1; // start again at the same point, no progress
         xLimit=state.x3;
         h=dxMin;  // have to take a closer look to a limited range at
            // higher resolution
         local=true;
      }
       // now we have new values for the iter ini data
      goto iter;
   }
}

X2<R,bool> CpmAlgorithms::root(const F<R,R>& f, R xMin,
                              R xMax, R dxMin, R dxMax, bool silent)
{
   Z mL=3;
   Word loc("CpmAlgorithms::root(...)");
   cpmmessage(mL,loc&" started");
   const Z iterMax=40;
   const R redFac=0.05;
   X3<R,R,Z> prr=brack(f,xMin,xMax,dxMin,dxMax);
   R x0=prr.c2();
   bool k=true;
   Z found=prr.c3();
   if (found==0){
      if (!silent) cpmwarning(loc&" could not bracket the root");
      k=false;
   }
   else if (found==1){
      ; // x0 and k are OK as already defined
   }
   else if (found==2){ // bisection
      R x1=prr.c1();
      R x2=prr.c2();
      x0=rootByBisection(f,x1,x2,dxMin*redFac,iterMax);
   }
   else{
      cpmerror("CpmAlgorithms::root() unvalid exit, found="
         &cpmwrite(found));
   }
   return X2<R,bool>(x0,k);
}

// Newton iteration with initialization by averaged iteration.
//  Best method.
//“Iteration Methods for Solving Kepler’s Equation”
//http://demonstrations.wolfram.com/IterationMethodsForSolvingKeplersEquation/
//Wolfram Demonstrations Project
//Published: June 24, 2015
R CpmAlgorithms::solKepEqu(R M, R eps, R acc)
{
   R xLast=M+1000_R, xNew=M;
  // cpmdata<<" data from solKepEqu(R M, R eps, R acc)"<<endl;
   Z co=0;
   Z coMax=2000;
   while ( cpmabs(xLast-xNew) > acc && co<coMax){
      xLast=xNew;
      if (co<40){
         R x1=M+eps*cpmsin(xNew);
         R x2=M+eps*cpmsin(x1);
         xNew=(x1+x2)*0.5;   // My standard provision against oscillations.
         // Works extremely well. Never seen this elsewhere.
        // cpmdata<<" co = "<<co<<endl;
        // cpmdata<<" xNew = "<<xNew<<endl;
      }
      else{
         xNew += (M - xNew + eps*cpmsin(xNew))/(1. - eps*cpmcos(xNew));
        // cpmdata<<" co = "<<co<<endl;
        // cpmdata<<" xNew = "<<xNew<<endl;
      }
      co++;
   }
   R err=xNew-M-eps*cpmsin(xNew);
   cpmassert(err<=2.*acc,"Solution of solKepEqu() not accurate");
  // cpmdata<<" data from solKepEqu(R M, R eps, R acc) done"<<endl;
   return xNew;
}

R CpmAlgorithms::ssqe(R a, R b, R c, R xMin, R xMax)
{

   R x1,x2,w1,w2,n,xmin,xmax;
   if (a==0){   // linear equation
      if (b==0) return -1; else x1=-c/b;
      if (x1>=xMin && x1<=xMax) return x1; else return -1;
   }
   else{   // quadratic equation

      w2=b*b-a*c*4;

      if (w2>0){   // two real solutions
         n=R(0.5)/a;
         w1=cpmsqrt(w2);
     //    x1=( b<0. ? (-b+w1)*n : (-b-w1)*n );   // this is the
         // numerically
            // uncritical root (never a difference of large quantities)
                // also never 0
         if (b<R{0.}){
            x1=(-b+w1)*n;
         }
         else { x1=(-b-w1)*n;}

         x2= c/(a*x1);   // second root from Vieta, c=a*x1*x2
         if (x1<x2) { xmin=x1; xmax=x2;} else { xmin=x2; xmax=x1;}
         if (xmin>=xMin && xmin<=xMax){
            return xmin;
         }
         else if (xmax>=xMin && xmax<=xMax){
            return xmax;
         }
         else{
            return -1;
         }
      }
      else if (w2<0){    // two complex solutions
         return -1;
      }
      else{    // w2==0, this gives the case of coinciding solutions
              // of a quadratic equation
         x1=-b/(a*2);    // notice a!=0
         if (x1>=xMin && x1<=xMax) return x1; else -1;
        }
   }

   cpmerror("psdqe(): never reach this point");

   return x1; // never reached, to avoid compiler warning
}

X2<R,bool> CpmAlgorithms::rootFast(F<R,R> const& f, R xMin, R xMax,
                                  R dx, R eps)
{
   const R red=0.667;
   const R redComp=R(1)-red;
   const Z iterMax=20;

   const Z mL=3;

   R range=xMax-xMin;
   if (range<=0.) cpmerror("rootFast(...): xMax<=xMin");

   R h=CpmRootX::inf(dx,range);
   R x1=xMin;
   R f1=f(x1);
   R err=cpmabs(f1);
   if (err<=eps){
      if (cpmverbose>=mL){
         Word mes="rootFast(...): root found immediately, err=";
         mes&=cpmwrite(err);
         cpmmessage(mes);
      }
      return X2<R,bool>(x1,true);
   }
   Z count=0;

Iter:
   count++;
   R x2=x1+h;
   R x3=x2+h;
   R f2=f(x2);
   R f3=f(x3);
   R hInv=R(1)/h;

   R a=(f1-f2*2+f3)*hInv*hInv*0.5;
   R b=(f1*3-f2*4+f3)*hInv*(-0.5);
   R t=ssqe(a,b,f1,0.,xMax-x1);

   if (t==-1){ //
      cpmmessage(mL,"rootFast(...): root called");
      return root(f,xMin,xMax,dx*0.1,dx*10);
   }

   R xRoot=x1+t;
   R fRoot=f(xRoot);
   err=cpmabs(fRoot);
   if (err<=eps){
      if (cpmverbose>=mL)
         cpmmessage("rootFast(...): reg. return, err"&cpmwrite(err));
      return X2<R,bool>(xRoot,true);
   }
   if (count==iterMax){
      if (cpmverbose>=mL)
         cpmmessage("rootFast(...): iterMax exhausted, count="&
            cpmwrite(count));
      return X2<R,bool>(xRoot,false);
   }
   x1+=red*t;
   f1=f(x1);
   h=CpmRootX::inf<R>(h,redComp*t);
   goto Iter;
}

X2<R,R> CpmAlgorithms::bisection(const F<R,bool>& f, R tF, R tT,
   R tacc, Z iterMax)
{
   if (tacc<=0) cpmerror("CpmAlgorithms::bisection(...): tacc<=0");
   Z count=0;
   R tM;
   while( cpmabs(tT-tF)>tacc && count++<iterMax ){
      tM=(tF+tT)*0.5;
      if ( f(tM) ) tT=tM; else tF=tM;
   }
   if (count>=iterMax)
      cpmwarning("CpmAlgorithms::bisection() exhausted iterMax");
   return X2<R,R>(tF,tT);
}

R CpmAlgorithms::rootByBisection(const F<R,R>& f, R x1, R x2,
                                 R xacc, Z iterMax)
{
   Z mL=3;
   Word loc("CpmAlgorithms::rootByBisection(...)");
   CPM_MA
   R f1=f(x1);
   R f2=f(x2);
   R xM,fM;
   if (f1==0){
      cpmwarning(loc&" x1 found to be a root");
      CPM_MZ
      return x1;
   }
   if (f2==0){
      cpmwarning(loc&" x2 found to be a root");
      CPM_MZ
      return x2;
   }
   if (!signChange(f1,f2)){
      cpmerror(loc&": no sign change for initial values");
      return 0; // never get here
   }
   Z count=1;
   while( cpmabs(x1-x2)>xacc && count<iterMax ){
      xM=(x1+x2)*0.5;
      fM=f(xM);
      if (fM==0) return xM;
      if ( f1*fM<=0) { // notice that sign change is what we need,
            // so f2 may be replaced by fM
         f2=fM;
         x2=xM;
      }
      else{ // then f1 and fM have the same sign , and this differs from
         // that of f2, so f2 and fM are a team, f1 gets replaced
         f1=fM;
         x1=xM;
      }
      count++;
   }
   if (count>=iterMax){
      cpmwarning(loc&": exhausted iterMax="&cpm(iterMax));
      cpmdebug(x1);
      cpmdebug(x2);
   }
   CPM_MZ
   return (x1+x2)*0.5;
}


R CpmAlgorithms::fixedPoint(const F<R,R>& f, R const& x0,
   R eps, Z iterMax)
{
   using CpmRoot::dis;
   R xOld=x0;
   R xNew=f(xOld);
   Z count=0;
   while(dis(xOld,xNew)>eps && count++ < iterMax ){
      xOld=xNew;
      R xAux=f(xOld);
      xNew=(f(xAux)+xAux)*0.5;
   }
   if (count>=iterMax)
      cpmwarning("CpmAlgorithms::fixedPoint() exhausted iterMax");
   return xNew;
}

//////////////////////// Quantizer /////////////////////////////////////
// auxiliaries:

namespace{

Z qant0(R const& x){ return cpmrnd(x);}

Z quantV(R const& x, const Vo<R>& v)
{
   return v.loc3(x).c1();
}

Z quantIv(R const& x, const Iv& iv, const Z& n)
{
   return iv.stair(x,n);
}

} // namespace

// constructors

CpmAlgorithms::Quantizer::Quantizer()
{
   *this=F<R,Z>(qant0);
}

CpmAlgorithms::Quantizer::Quantizer(R xL, R xU, Z n)
{
   Word loc("Quantizer::Quantizer(R xL, R xU, Z n)");
   cpmassert(n>0,loc);
#ifdef CPM_Fn
   *this=F2<R,Iv,Z,Z>(Iv(xL,xU),n)(quantIv);
#else
   *this=F<R,Z>(bind(quantIv,_1,Iv(xL,xU),n));
#endif
}

CpmAlgorithms::Quantizer::Quantizer(Iv iv, Z n)
{
   Word loc("Quantizer::Quantizer(Iv iv, Z n)");
   cpmassert(n>0,loc);
#ifdef CPM_Fn
   *this=F2<R,Iv,Z,Z>(iv,n)(quantIv);
#else
   *this=F<R,Z>(bind(quantIv,_1,iv,n));
#endif
      // this stores the data iv, n that characterize the
      // instant of Quantizer constructed by the given constructor.
      // We could define another constructor of the class which stores a
      // completely different data structure; for example some complex
      // numbers. This shows that a class needs not to have a single data
      // member'layout' if sufficiently flexible data structures are
      // available.
}

CpmAlgorithms::Quantizer::Quantizer(const V<R>& x, bool ordered)
{
   Z n=x.dim();
   if (n<=0){
      Z constVal=0;
      *this=F<R,Z>(constVal);
   }
   else{
      Vo<R> y(x);
      if (!ordered) y=y.strictlyIncreasingVersion();
#ifdef CPM_Fn
      *this=F1<R,Vo<R>,Z>(y)(quantV);
#else
      *this=F<R,Z>(bind(quantV,_1,y));
#endif
   }
}

////////////////////////// IntQuant ///////////////////////////////

namespace{

   Z intQuant(const Z& i, const Z& zL, const Z& nP, R const& spanInv)
   {
      R p=R(i-zL)*spanInv;
      bool restrict=true;
      return Iv(0,1).quant(p,nP,restrict);
   }
}

CpmAlgorithms::IntQuant::IntQuant(Z xL, Z xU, Z n)
{
   if (xL<=xU){
      zL=xL;
      zU=xU;
   }
   else{
      zL=xU;
      zU=xL;
   }
   nP=(n<1 ? 1 : n);
   spanInv=( zU==zL ? 0 : 1./(zU-zL) );
#ifdef CPM_Fn
   f=F3<Z,Z,Z,R,Z>(xL,nP,spanInv)(intQuant);
#else
   f=F<Z,Z>(bind(intQuant,_1,xL,nP,spanInv));
#endif
}

//////////////////////// Percentilizer //////////////////////////////////

namespace { // auxiliaries in anonymous namespace
   using CpmArrays::Vo;
   using CpmArrays::V;

class quantil: public FncObj<R,R>{
   Vo<R> b;
public:
explicit quantil(const V<R>& a):b(a){}
   R operator()(R const&)const;
};

R quantil::operator()(R const& p)const
// notice: function quantil will only be called for b in ascendend order
// quantil(0.5,b)=median(b);
// quantil(0.25,b)=Hinge_1(b);
// quantil(0,b)=minimum(b);
// quantil(0.75,b)=Hinge_2(b);
// quantil(1,b)=maximum(b);
// see Weisstein: CRC Concise Encyclopedia of Mathematics p. 844 for
// Hinge.
// we use an obvious method to define this in cases where n is not of the
// form
// 4*m+5
{
   R pp=p;
   if (pp<0) pp=0;
   if (pp>1) pp=1;
   Z n=b.dim();
   if (n==0) return 0;
   R ir=R(n-1)*pp+1;
   if (ir==1.) return b[1];
   if (ir==n) return b[n];
   Z i=cpmtoz(ir);
   R frac=ir-i; // notice 0<=frac<1
   return frac*b[i]+(R(1.)-frac)*b[i+1];
}

} // anonymus namespace for hiding this construct

// constructors

CpmAlgorithms::Percentilizer::Percentilizer(const V<R>& a):arr(a)
{
   if (arr.dim()>0) arr.order_();
   p_=new quantil(arr); // no ordering done here !!!
   // p_ is a protected member base class P< FncObj<R,R> >
   // notice that the initialization of u_ by the default constructor is
   // just right.
   // This structure of a constructor seems to be the most flexible one.
}

V<R> CpmAlgorithms::Percentilizer::select(const Iv& iv)const
{
   Z mL=4;
   static Word loc("CpmAlgorithms::Percentilizer::select(Iv)");
   cpmmessage(mL,loc&" started");
   R iv1=iv[1];
   cpmassert(iv1>=0,loc);
   R iv2=iv[2];
   cpmassert(iv2<=1,loc);
   Z i,n=arr.dim();
   R x1=(*this)(iv1);
   R x2=(*this)(iv2);
   Iv ivOk(x1,x2);
   V<B> sel(n);
   for (i=1;i<=n;i++) sel[i]=ivOk.contains(arr[i]);
   cpmmessage(mL,loc&" done");
   return arr.select(sel);
}

R CpmAlgorithms::Percentilizer::mean(const Iv& iv)const
{
   Va<R> arrOk=select(iv);
   Z nOk=arrOk.dim();
   return nOk==0 ? R(0) : arrOk.sum()/nOk;
}

R CpmAlgorithms::gaussRandomValue(R const& mean, R const& sigma, Z j)
// Returns normally distributed random numbers with mean and sigma
// given by the respective arguments.
// Modified from gasdev(long*) of Press et al. by introducing another
// basic random generator and by not specializing the arguments.
{
    Z j2=2_Z*j;
    if (sigma==0_R) return mean;
    R fac,rsq,v1,v2,r1,r2;
    do{
         r1 = randomR(j2);
         r2 = j2==0_Z ? randomR(j2) : randomR(j2-1_Z);
         v1=2_R*r1-1_R;
         v2=2_R*r2-1_R;
         rsq=v1*v1+v2*v2;
         if (j2!=0_Z) j2+=2_Z;
     }
     while(rsq>=1_R || rsq==0_R);
        // these are the results we don't want
     fac=sigma*cpmsqrt(-2_R*cpmlog(rsq)/rsq);
     return mean+v1*fac;
}

///////////////////// class PeriodicAlarm //////////////////////////////

using CpmAlgorithms::PeriodicAlarm;

PeriodicAlarm::PeriodicAlarm(R per, R tIni, R phase):
active(true),
period(per),
lastAlarm(tIni),
counter(0)
{
   Z mL=3;
   static Word loc="PeriodicAlarm(R,R,R)";
   static R relTol=1e-8;
   CPM_MA
   if (period<=0) active=false;
   tol=relTol*period;
   nextLatticePoint=tIni+period*phase;
   CPM_MZ
}

bool PeriodicAlarm::operator()(R t)
{
   if (!active){ // period<=0 implies active==false
      return false;
   }
   R tDiff=t-nextLatticePoint+tol; // till 2016-09-11 tol was missing
      // here, and was placed as  'if (tDiff+tol>=0){...'
      // This allowed fr to start with -1 so that one could traverse a
      // alarm time without changing nextLatticePoint
      // (since in nextLatticePoint+=R(fr)*period; we had fr=0).
   if (tDiff>=0){
      Z fr=cpmtoz(tDiff/period); // starts with 0
      lastAlarm=t; // includes the case that t grows in steps larger
         // than period
      fr+=1;
         // fr should start with 1
      nextLatticePoint+=R(fr)*period;
     // nextLatticePoint+=period;
      counter+=fr;
     // counter++;
      return true;
   }
   else{
      return false;
   }
}

//////////////////////// class SingleAlarm /////////////////////////////

using CpmAlgorithms::SingleAlarm;

bool SingleAlarm::operator()(R t)
{
   if (done) return false;
   if (t>=alarmPoint){
      done=true;
      return true;
   }
   else{
      return false;
   }
}

//////////////////////// class CyclicAlarm /////////////////////////////

using CpmAlgorithms::CyclicAlarm;

CyclicAlarm::CyclicAlarm(Z per, Z cogIni, Z maxAlarms):
period(per),cog(cogIni),counterMax(maxAlarms)
{
   Z mL=3;
   static Word loc("CyclicAlarm(Z,Z,Z)");
   CPM_MA
   if (period<1){
      dormant=true; // alarm=false by initialization
      done=true;
   }
   else{
      while (cog>=period) cog-=period;
      alarm=true;
   }
   CPM_MZ
}

void CyclicAlarm::step(void)
{
   if (dormant){
      return;
      // don't change state
   }
   else if (done){
     alarm=false;
   }
   else{
     cog+=1;
     if (cog==period || cog==0){
         cog=0;
         alarm=true;
         counter+=1; // now this the number of alarms given
         if (counter==counterMax){
            done=true;
            alarm=false;
         }
      }
      else{
         alarm=false;
      }
   }
}

V<B> CyclicAlarm::record(Z n)
{
   V<B> res(n);
   for (Z i=1;i<=n;i++) res[i]=(*this)();
   return res;
}

////////////////////// class NonCyclicAlarm /////////////////////////////

using CpmAlgorithms::NonCyclicAlarm;

NonCyclicAlarm::NonCyclicAlarm(const V<Z>& points):
cog(1), counter(1), points_(points)
{
   Z mL=3;
   static Word loc("NonCyclicAlarm()");
   CPM_MA
   Vo<Z> p_=points_;
   points_=p_.strictlyIncreasingVersion();
   counterMax=points_.dim();

   if (counterMax<=0){
      alarm=false;
      done=true;
   }
   else{
      alarm=(cog==points_[counter]);
      if (alarm) counter+=1;
      done=(counter>counterMax);
   }
   CPM_MZ
}

void NonCyclicAlarm::step(void)
{
   if (done){
     alarm=false;
   }
   else{
     cog+=1;
     if (cog==points_[counter]){
         alarm=true;
         counter+=1; // now this the number of alarms given
         if (counter>counterMax){
            // never call points_[counterMax+1]
            done=true;
         }
      }
      else{
         alarm=false;
      }
   }
}

//////////////////////// class BurstAlarm /////////////////////////////

using CpmAlgorithms::BurstAlarm;

BurstAlarm::BurstAlarm(Z per, Z burPer, Z burLen):
Base(per,0,-1),burst(burPer,0,burLen)
{
   Z mL=3;
   static Word loc("BurstAlarm(Z,Z,Z)");
   CPM_MA
   CPM_MZ
}

void BurstAlarm::step(void)
{
   if (getAlarm() && burst.ready()) burst.restart();
   Base::step();
   burst.step();
   setAlarm(getAlarm()||burst.getAlarm());
}

////////////////////////////////////////////////////////////////////////
Z CpmAlgorithms::getNerFrc_1(R& x, R xL)
{
  Z mL=4;
  Word loc("CpmAlgorithms::getNerFrc_1(R&,R)");
  CPM_MA
  Z res=0;
  cpmassert(x<xL,loc);
  cpmassert(x!=0,loc);
  res=cpmrnd(xL/x);
  cpmassert(res!=0,loc);
  if (res!=0) x=xL/res;
  CPM_MZ
  return res;
}

namespace{
   Z penalty(Z n1, Z n2, Z n)
   {
      Z z1=n1, z2=n2;
      CpmRootX::order<Z>(z1,z2);
      Z p=z1*z2;
      Z a1=z2-z1;
      Z a2=p-n;
      if (a2<0) return 1+n*n;
      return a1+a2;
   }
}

Z2 CpmAlgorithms::genSqrRoot(Z n)
{
   Word loc("CpmAlgorithms::genSqrRoot(Z)");
   Z z1=cpmrnd(cpmsqrt(R(n)));
   Z z2=n/z1;
   Z n0=z1*z2;
   if (n0==n) return Z2(z1,z2);
   CpmRootX::order<Z>(z1,z2);
   if (n0>n){ // in principle OK, but try an improvement
      Z p1=penalty(z1,z2,n);
      Z p2=penalty(z1+1,z2-1,n);
      if (p2<p1){
         z1++; z2--;
      }
      return Z2(z1,z2);
   }
   while (n0<n){
      Z p1=penalty(z1+1,z2,n);
      Z p2=penalty(z1,z2+1,n);
      if (p1<=p2) z1++;
      else z2++;
      n0=z1*z2;
   }
   CpmRootX::order<Z>(z1,z2);
   cpmassert(z1*z2>=n,loc);
   return Z2(z1,z2);
}

Z CpmAlgorithms::appSqrRoot(Z n)
{
   Word loc("CpmAlgorithms::genSqrRoot(Z)");
   return cpmrnd(cpmsqrt(R(n)));
}