//? cpmalgorithms.h
//? 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.


#ifndef CPM_ALGORITHMS_H
#define CPM_ALGORITHMS_H

#include <cpmx.h>
#include <cpmvo.h>
#include <cpminterval.h>
#include <cpmf.h>
#include <cpmrvector.h>
#include <cpmrecordhandler.h>
#include <concepts>
#include <functional>

/*
   Description: This namespace will be populated by functions and classes which
    express algorithms of general interest, which don't fit well into
    one of the classes already introduced.
*/
namespace CpmAlgorithms{

   using namespace CpmStd;
   using namespace CpmRoot;

   using CpmArrays::V;
   using CpmArrays::VV;
   using CpmArrays::Vo;
   using CpmArrays::X2;
   using CpmArrays::X3;
   using CpmArrays::R_Vector;
   using CpmFunctions::FncObj;
   using CpmFunctions::F;
#ifdef CPM_Fn
   using CpmFunctions::F2;
   using CpmFunctions::F3;
#endif
   using CpmGeo::Iv;
   using CpmLinAlg::Z2;
   using CpmRootX::RecordHandler;

//////////// functions gcd, lcm, lesDiv, prmFac, isPrm /////////////////
// Division properties of natural numbers
   N gcd(N const& a, N const& b, bool verify=false);
      //! greatest common divisor
      // No assumption concerning the arguments. If a or b is 0,
      // we return 0, which never occurs for meaningful arguments.
      // If verify==true, it is tested whether the result actually
      // divides a and b; if not an error is issued.
   N lcm(N const& a, N const& b);
      //! least common multiple
      // Implementation based on a*b=gcd(a,b)*lcm(a,b) and
      // ise careful in avoiding temporaries that are larger
      // than the result. If any of the arguments is 0, we return 0.
   N lesDiv(N const& n);
      //: least divisor
      // returns the least (true, i.e. !=1) divisor of n.
   V<N> prmFac(N const& n, bool verify=false);
      //: prime factors
      // Returns the prime factors of n as an array in
      // increasing order. If verify==true, it is tested whether the
      // components of the result actually multiply to n;
      // if not an error is issued.
   bool isPrm(N const& n);
      //: is prime
      // returns true iff n is prime. Although the algoritm is trivial
      // it takes less than 10 s to write a list of primes from 2
      // to 99991 (the range covered in Abramowitz, Stegun pp. 870-873)
      // to file.

/////////// functions ZfromBits, bitsFromZ, next_1 ///////////////////////

Z ZfromBits(const V<B>& vb);
   // for n:=vb.dim() we return
   // sum(i=1,...n) vb[i]*2^(n-i)
   // i.e. the positive integer the binary expansion is given by
   // vb where vb[1] is the most significant bit
   // and vb[n] the least significant one

V<B> bitsFromZ(Z n, Z p);
   // we assume and test 0<=n<=(2^p)-1 and return vb such that
   // ZfromBits(vb)==n; 2002-10-02: tested for (n,p)=(673,10),(972,10),
   // (0,12),(4095,12) always OK

bool next_1(std::vector<B>& v);
   //: next
   // Notice that the trailing '_1' says that the argument
   // gets changed by a function call.
   // After the call, v is the lexicographic successor of
   // the input v (order is false < true, and as successor
   // only lists are considered which have the same number
   // of components as v. If in v already all components
   // have the value true, there is no successsor and
   // the return value is false (i.e. no success in finding
   // a successor). Otherwise, the return value is true.
   // If v.size()==0 nothing is done on v and 'false' is returned.
   // I use std::vector here since I need it in a context were
   // high efficiency is desirable. Since the template parameter
   // <bool> is a very light one, the reference counting advantage
   // of V does not exist here.

V<L> digFromR(R const& x);
   //: (list of) digits from R
   // The digital seperation symbol (here the period) treated as a special
   // digit.

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

class LatticeEnumerator{ // enumerating the points in cubic lattices
   // consider the lattice {1,2,...b[1]}x{1,2,...b[2]}x...x{1,2,...b[n]}.
   // the class provides a way to enumerate the points and to
   // construct a (i[1],...i[n]) from the enumerating integer.

// indpendent data
   V<Z> b;
// dependent data
   Z n; // n=b.dim()
   V<Z> c;
      // c.dim()=n
      // example n=4
      // c[1]=b[2]*b[3]*b[4]
      // c[2]=     b[3]*b[4]
      // c[3]=          b[4]
      // c[4]= 1

   Z nLattice; // c[1]*b[1]=b[1]*b[2]*...*b[n]
   void update();
      // initializes dependent data
   typedef LatticeEnumerator Type;
public:
   CPM_IO
   LatticeEnumerator(const V<Z>& bounds):b(bounds){update();}
      // obvious
   LatticeEnumerator(){}
   Z operator()(const V<Z>& x)const;
      // We assume (and check) x.dim()==n
      // 1<=x[i]<=b[i] for all 1<=i<=n.
      // for given bounds, there are b[1]*...*b[n]=nLattice
      // points x which satisfy these conditions. The function
      // provides a suitable enumeration of these points by
      // assigning the enumeration number to x as a function value.
   V<Z> operator()(Z k)const;
      // inverse mapping of the previous one
};

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

class LatticeBox{
 // putting points into cells defined by a cartesian product of intervals

   V<Iv> viv;
   V<Z> vz;
   Z n;
   Z np;
   V<R> val;
   V< V<R> > centers;
   LatticeEnumerator le;
public:
   LatticeBox(const V<Iv>& vivIn, const V<Z>& vzIn);
   V< V<R> > points()const{ return centers;}
   V<R> values()const{ return val;}
   void set(const V<R>& pos, R valIn);
   Z dim()const{ return n;}
   Z nPoints()const{ return np;}
   VV<R> matrix()const;
      // values of the 1-2-plane in matrix-indexed arrangement
};

///////////// bracketing a root /////////////////////////////////////////

X3<R,R,Z> brack(const F<R,R>& f, R xMin, R xMax, R dxMin, R dxMax);
   // let f(x) be defined for xMin<=x<=xMax. We try to find x1 and x2
   // in [xMin,xMax] such that
   // (i) x1<x2 and sign(f(x1))!=sign(f(x2)) (then a root x0, f(x0)==0,
   //    is said to be bracketed, x1<=x0<=x2, see Press et al. p. 350)
   // and
   // (ii) x1 and x2-x1 as small as possible
   // The return values x1:=brack.c1(), x2:=brack.c2()  are guaranteed to
   // have property (i) if found>0 (found=brack.c3()) after call.
   // For found==0 the method found no indication for a sign change
   // of the function in the range [xMin,xMax].
   // For found==1, we got directly a root which is returned both as
   // x1 and x2.
   // For found=2, we found two bracketing points (normal behavior).
   //
   // dxMax is the step size of a first scan. The method ends with
   // found==0 if in this scan (testing a eqidistant chain of points)
   // no sign change and no  trend change (sign change of first
   // differences) was found. In case of detecting a trend change, we
   // performe a refined local scan of the interval which showed the
   // trend change. For scan length < dxMin no further refinement will
   // be done and a mere trend change will not cause further action.
   // Only sign changes will be considered.
   // Notice: We assume xMin<xMan, dxMin<dxMax, otherwise the behavior
   // is undefined.

X2<R,bool> root(const F<R,R>& f, R xMin, R xMax, R dxMin, R dxMax,
      bool silent=false);
   // Let the return value be written as (x,b).
   // For b=true, x is a root: f(x) close to 0.
   // For b=false, no root was found and x is the last value of x that
   // was treated. The arguments xMin,...dxMax are used as arguments in
   // an initial call to function brak. The final localization of the
   // root is done by bisection.

X2<R,bool> rootFast(const F<R,R>& f, R xMin, R xMax, R dx, R eps);
   // An experimental method to find a root of f in the interval
   // [xMin,xMax]. One starts investigating f at xMin and tries to
   // find the root closest to xMin (just as in ssqe). One does
   // successive parabolic interpolation. In order to avoid
   // instabilities we go not completely to the extrapolated position.
   // In the intended application to visualization of extended
   // bodies it did not result in a speed improvement. For details
   // see the code.

//////////// special solution of quadratic equation /////////////////////

R ssqe(R a, R b, R c, R xMin, R xMax);
// special solution of quadratic equation: returns the smallest
// solution x of the equation a*x*x+b*x+c=0 that satisfies xMin <= x <=
// xMax.
// If no such solution exists in this range, -1 is returned
// We assume xMin < xMax (otherwise behavior is not specified).
// This is a rather optimized algorithm which (for instance) helps to
// find effectively the point where a light ray hits a sufficiently
// regular body.

//////////////// solution of Kepler's equation (for 0<=eps<1) ////////////
R solKepEqu(R M, R e, R acc=1e-8);
//: Solution (of) Kepler's  equation.
// We return the solution E of the famous transcendental equation
//    E = M + eps * sin E
// In astronomical nomenclature M is the 'mean anamoly'
// E the 'excentric anomaly, and eps the numerical excentricity.

////////////////////  bisection /////////////////////////////////////////

   X2<R,R> bisection(const F<R,bool>& f, R xF, R xT,
      R xacc=1e-3, Z iterMax=40 );
      // Finding a boundary of a 1-dimensional domain characterized by an
      // 'indicator function'.
      // If f(xF)=false, f(xT)=true, xacc>0,  we return a pair (xf,xt)
      // such that
      // f(xf)= false and f(xt)=true where xt-xf = (xT-xF)/(2^n)
      // for some natural number n and |xt-xf|<xacc.
      // (Notice: if xF<xT then xt<xf; if xF>xT then xf>xt)
      // An unbiased guess for the boundary itself then is 0.5*(xf+xt).
      // The number of iterations is limited to iterMax, if the limit is
      // reached a warning is issued and the values reached so far are
      // returned (they, then, do not satisfy |xt-xf|<xacc).

   R rootByBisection(const F<R,R>& f, R x1, R x2, R xacc, Z iterMax=40 );
      // Finding the a root of f by bisection, assuming that
      // f(x1) and f(x2) differ in sign. The control of accuracy is just
      // as in the previous function. The return value is 0.5*(x1b+x2b);
      // where (x1b,x2b) resulted from the last bisection step.

//////////////////////  fixedPoint //////////////////////////////////////

   template <class X>
   X fixedPoint(const F<X,X>& f, const X& x0,
      R eps=1e-3, Z iterMax=40 );
      // Finding a solution x of x=f(x). Start point for the iterative
      // solution is x0, and the condition for terminating the iteration
      // x_0, x_1, x_2, ... is (x_(n+1)).dis(x_n)<=eps.
      // We thus assume that class X defines a function X.dis(X)
      // such that x1.dis(x2) is small for x1 and x2 being close
      // together. The number of iterations is limited to iterMax, if
      // the limit is reached a warning is issued and the value reached
      // so far is returned. A simple and effective provision against
      // oscillatory behavior is implemented.

// notice that the MS compiler (linker?) expects all template definitions
// to be placed in header files (deviating from the practice suggested by
// Stroustrup)

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

   R fixedPoint(const F<R,R>& f, const R& x0,
      R eps, Z iterMax);
      // Specialization of the template version; needed since R can't
      // provide member functions (and the legal work-arround with a
      // non-member dis(X,X) does not work always with the MS
      // compiler.

//////////////////////////  Place ////////////////////////////////////////

   template <class X>
   class Place{//creating X-valued lists, satisfying specified conditions
      F<Z,X> create; // generator for X's
      F<X,bool> condition1; // condition which each x should satisfy
         // in order to be acceptable for a result list to be created
         // via opertor ()
      F< X2<X,X>, bool> condition2;
         // condition which each pair x[i], x[j]  should satisfy
         // in order to make x[i] and x[j] both be acceptable as
         // members of the result list to be created via opertor ()
      bool useCondition1;
         // If this is false, the loop for testing condition1
         // will be omitted (added 2009-02-28)
      bool useCondition2;
         // If this is false, the loop for testing condition2
         // will be omitted (added 2006-02-21)
     public:
     Word nameOf()const
     { return Word("Place< "&CpmRoot::Name<X>()(X())&" >"); }

     Place(const F<Z,X>& cr,
         const F<X,bool>& cond1,
         const F< X2<X,X>, bool>& cond2,
         bool useCon1=true, bool useCon2=true):
         create(cr),condition1(cond1),
         condition2(cond2),
         useCondition1(useCon1),
         useCondition2(useCon2){}
         // bringing the condition functions in place, that will be
         // utilized in a call to operator()

      Place(const F<Z,X>& cr,
         const F<X,bool>& cond1):
         create(cr),condition1(cond1),
         condition2(),
         useCondition2(true),
         useCondition2(false){}
         // Faster form of case in which condition2 is the 'always true
         // function'. Also more convenient since one has not to
         // give a definition of the 'always true function'.

     V<X> operator()(Z n, Z trialMax=1000 )const;
      // From the X-valued sequence create(1), create(2), create(3), ....
      // we select a finite sublist x[1],....x[n] such that
      // condition1(x[i])==true for all 1 <= i <= n
      // condition2(x[i],x[j])==true for all 1 <= i,j <= n, i!=j
      // If we would need higher sequence members than create(n+trialMax)
      // the function stops working and returns the array it has obtained
      // so far (this than has less than n components). The strategy is
      // simple: create a new X (via create()) test the conditions, and
      // take a new one if any of the conditions fails.
   };

template <class X>
V<X> CpmAlgorithms::Place<X>::operator()(Z n, Z trialMax )const
{
   Z mL=2;
   Word loc=nameOf()&"::operator()(Z,Z)";
   CPM_MA
   Z i=0,nValid=0;
   V<X> res(n);
   X xi;
   cpmassert(n>0,loc);
   R nInv=1./n;
   Z iLimit=n+trialMax;
   while (nValid<n){
LABEL:
      if (cpmverbose>mL) cpmprogress("Place-loop",cpmtod(R(nValid)*nInv));
      xi=create(++i);
      if (i>iLimit) break;
      if  (useCondition1){
         if (condition1(xi)==false) goto LABEL;
      }
      if (useCondition2){
         for (Z j=1;j<=nValid;j++){
            if (condition2(X2<X,X>(res[j],xi))==false) goto LABEL;
         }
      }
      res[++nValid]=xi;
   }
   if (nValid<n){
      cpmwarning(loc&": could create only "&
         cpm(nValid)&" items of "&cpm(n));
   }
   CPM_MZ
   return res.resize(nValid);
}

///////////////////// random selector /////////////////////////////////////

class RandomSelector{
   R_Vector w_; // List of weigths. negative values are treated as 0.
   R_Vector wp_; // negative components of w_ replaced by 0
   R s_; // sum of the components of wp_
   R_Vector wn_; // version of wp_ normalized such that the components add
      // up to 1.
public:
   RandomSelector(R_Vector const& w):w_(w),wp_(w_.cutToPosCopy()),
   s_{wp_.sum()},wn_{wp_*cpminv(s_)}{}

   Z operator()(void)const{
      R r=randomR();
      Z i=1;
again:
      if ( r< wn_[i]) return i;
      i++;
      goto again;
      return wn_.e();
   }
};

////////////////////  RandomIntegration /////////////////////////////////

   template <class X>
   class RandomIntegration{
      // Calculating expectation values based on a X-valued list.
      // A X-valued list is the algorithmic substitute for the notion of
      // a probability measure on X:
      // Let S be a subset of X and let x[i]:=create(i) be the sequence
      // of X's resulting from the creation mechanism of the present
      // class. Then
      //    card{i | 1<=i<=n, x[i] element of S}/n
      // is a guess for (measure of S)/(measure of X).
      // If create() generates a random sequence, this will become
      // independent of n for large n 'in the sense of probabilistic
      // convergence'.
      // All random aspects ly in the function create; the algorithms
      // implemented here are strictly deterministic (and repeatable).
         F<Z,X> create; // generator for X's
         F<X,bool> chi; // this is a indicator function of a subset in X
            // over which all integrations will take place effectively.
            // chi(x)==true means, that x belongs to the 'effective
            // domain of integration'.
         Z norMet_; // normalization method
            // if 1, the identity function on X gets the integral 1
            // if 0, the identity function on X gets the the integral
            // (measure of {x \in X | chi(x)==true})/(measure of X)

      public:
        RandomIntegration(const F<Z,X>& cr,
           const F<X,bool>& ch, Z norMet=1):
           create(cr),chi(ch),norMet_(norMet){}
            // bringing the functions in place, that will be utilized
            // in a call to operator()

        V<R> operator()( V<F<X,R> > const&, Z n, Z trialMax=1000 )const;
      // For an array of functions the expectation value is calculated.
      // Actually, each of the functions is set to zero outside the set
      // defined by function chi. It is assured that n points falling
      // into this subset are used for calculation. If this set is so
      // small that more than trialMax points from function create() had
      // to be tested, an error is issued. It is essential for efficiency
      // to let one call to this function evaluate the expectation values
      // of a full list of functions: doing these functions in succession
      // would ask for creating the random sequence many times.
   };

template <class X>
V<R> RandomIntegration<X>::operator()
( V<F<X,R> > const& f, Z n, Z trialMax)const
{
   Word loc("V<R> RandomIntegration<X>::operator()");
   Z d=f.dim();
   V<R> res(d,0.);
   Z i=1;
   Z nDone=0;
   Z j;
   while (nDone<n){
      if (i>trialMax)
         cpmerror("Too many trials in RandomIntegration<X>::operator()");
      X x=create(++i);
      bool isValid=chi(x);
      if (norMet_ || isValid){
         nDone++;
      }
      if (isValid){
         for (j=1;j<=d;j++){
            res[j]+=(f[j])(x);
         }
      }
   }
   cpmassert(nDone>0,loc);
   R nInv=1./nDone;
   for (j=1;j<=d;j++) res[j]*=nInv;
   return res;
}

///////////////////// class AveragingInterpolator //////////////////////

template <class X, class Y>
// we assume R X::distSquare(X) , Y*=R, Y(), Y+=Y
// General interpolation method for vector-valued functions on
// metric spaces. Seems to work well.
// Method was invented (by me) in 1997 as a way to define
// a constructor for Fr<> from discrete data.
class AveragingInterpolator: public FncObj<X,Y>{
// interpolation method for vector-valued functions on metric spaces
   // AveragingInterpolator is a non-mutable class
   const V<X> vx;
   const V<Y> vy;
   const Z ns;
      // if this is different from 0, the ns nearest
      // neighbors are used for forming a weighted mean.
      // This helps to keep the scheme local and to have
      // only to compute a linear combination of ns terms.
      // (normally not a concern). The price is an vector ordering
      // operation (to be done for each evaluation of f)
   static Y avrIpl(const X& x, const V<X>& vx, const V<Y>& vy,
      const Z& nn);
   // algorithmic content of the class as a static function.
   // This seems to be the most direct way of injecting complicated
   // algorithms into my F<X,Y> class
public:
   AveragingInterpolator(const V<X>& x, const V<Y>& y, Z n=0);
      // default ns==0 guarantees a continuous interpolation
      // function. Be cautious when using a different value
   Y operator()(const X& x)const{ return avrIpl(x,vx,vy,ns);}
};

template <class X, class Y>
Y AveragingInterpolator<X,Y>::avrIpl(const X& x,
   const V<X>& vx, const V<Y>& vy, const Z& nn)
{
   const R p=-1.75;
   // p==-2 and smaller (sic, not in absolute value) gives something
   // close to step function interpolation (flat terrasses connected
   // by steep ridges). Has to be negative in order to give near points
   // heavy weight.

   Z mL=4;
   Word loc("AveragingInterpolator<>::avrIpl(...)");
   CPM_MA
   Z i,n=vx.dim(),ny=vy.dim();
   cpmassert(n==ny,loc);
   if (n==0) return Y();

   R val,sum=0.;
   Y res=Y();
   if (nn==0 || nn>n){ // weighted mean over all points
      for (i=1;i<=n;i++){
         R r2=x.distSquare(vx[i]); // distSquare more efficient than
         // distance, could be important for large n
         if (r2==0)  return vy[i];
         val=cpmpow(r2,p); // val ~ 1/(r^2p) for application to points
            // in 3-space the power of r^2/r^2p should be <1 for
            // so that integral from 1 to infinity of r^2/r^2p is
            // finite
         sum+=val;
         res+=(vy[i]*val);
      }
   }
   else{ // if the vy[i] are 'heavy' objects for which
         // computing linear combinations is expensive,
         // this method which works with short linear combinations at
         // the price of having to order an vector might be better
      cpmassert(nn>0,loc);
      Vo<R> vr2(n);
      for (i=1;i<=n;i++) vr2[i]=x.distSquare(vx[i]);
      V<Z> per=vr2.permutationForIncreasingOrder();
         // finding the nearest neighbors
      for (i=1;i<=nn;i++){  // short sum over nearest neighbors
         Z j=per[i];
         R r2=vr2[j];
         if (r2==0) return vy[j];
         val=cpmpow(r2,p);
         sum+=val;
         res+=(vy[j]*val);
      }
   }
   cpmassert(sum>0.,loc);
   R sumInv=R(1.)/sum;
   res*=sumInv;
   CPM_MZ
   return res;
}

template <class X, class Y>
AveragingInterpolator<X,Y>::AveragingInterpolator
(const V<X>& x, const V<Y>& y, Z n):vx(x),vy(y),ns(n)
{
   Z mL=3;
   Word loc("AveragingInterpolator(...)");
   CPM_MA
   cpmassert(vx.dim()<=vy.dim(),loc);
   CPM_MZ
}

//////////////////////// Quantizer /////////////////////////////////////

   // Putting real numbers into boxes defined by subintervals of an
   // interval.

   // The interval [xL,xU] is divided into n equally sized disjoint
   // boxes.
   // These are numbered from 1 to n in the order of increasing x.
   // If x belongs to one of the boxes, we return the number of it. If x
   // belongs to the interval [xL,xU] it always belongs to exactly one
   // box by the definition that an x falling exactly on 'an internal
   // separation line' between boxes, belongs to the box with the lower
   // index. For x<xL we return 0 and for x>xU we return n+1.
   // The degenerate case is xL=xU. Here we return 1 for x=xL ( and,
   // again, n+1 for x>xU). The previous explanation describes the first
   // version of Quantizer considered. Meanwhile a second constructor,
   // taking as an argument a const V<R>& v, provides more flexible
   // quantisation services. Consider the set X of R's provided by the
   // components of v. (this may be void since v.dim() is allowed
   // to be 0. Notice card(X)<=v.dim().
   // Set card(X)=:n+1. We order the elements of X and get n ordered
   // boxes. Quantizing x means to find the index of the box to which it
   // belongs. Again we return 0 for x below every box and n+1(=card(X))
   // above every box.
   // For card(X)=0, it is thus consistent to return always 0, indicating
   // that every x is below every box and above every box.

   // This functionality of quantizing is here put into a class
   // (instead of a function) since this allows to separate the pre-
   // computations (which get done after constructor call) and the
   // actual quantization (which get done after call of operator()).
   // This could increase speed in a situation in which many numbers
   // have to be quantized with respect to one system of values xL,xU,n.

class Quantizer: public F<R,Z>{ // Putting real numbers into boxes
      // this derivation defines Z operator()(const R&)const and makes
      // available all member functions of F<R,Z> and all functions that
      // take F<R,Z>& as arguments. Notice that the only data member is
      // an address of a memory area that stores an instance of a class
      // that is known to be derived from FncObj<R,Z>. It thus defines Z
      // operator()(const R&)const; the actual definition of this virtual
      // operator may utilize data of types that vary from one value of
      // the adress to the next. So, the class itself does not determine
      // the information layout of its internal state. Very advanced
      // construct, indeed. Notice that this is an immutable class that
      // satisfies the value interface.

   Quantizer(const F<R,Z>& f):F<R,Z>(f){}
      // private tool for defining public constructors
public:
   Quantizer();
      // trivial round-Off, needs no data!
   Quantizer(R xL, R xU, Z n);
      // setting the limits of the interval and the number of boxes.
      // If n is <1, the actual number of boxes will be 1.
      // xL, xU needs not to be ordered on input.
   Quantizer(Iv, Z n);

   Quantizer(const V<R>& x, bool ordered=false);
      // here the limits of the boxes are obtained from an ordered and
      // unified version (multiple appearance of a component in x is
      // changed into single appearance) of the array x. Here the boxes
      // need not be of equal size. If there are no boxes (x.dim()==0),
      // operator() will always return 0, which is by far the most
      // logical prescription (see above). If the second argument is
      // true, we assume that x is already ordered and unified.
   F<R,Z> toBase()const{ return *this;}
};

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

   // Putting integer numbers into boxes defined by subintervals of an
   // integer interval.

   // The integer interval [xL,xU] i.e. the set {xL,xL+1,...xU}
   // is divided into n approximately equally sized disjoint
   // subintervals (equally sized they can be only if xU-xL+1
   // is a multiple of n)
   // These are numbered from 1 to n in the order of increasing x.
   // If x belongs to one of the boxes, we return the number of it. If x
   // belongs to the interval [xL,xU] it always belongs to exactly one
   // box by the definition.
   // For x<xL we return 1 and for x>xU we return n.
   // The degenerate case is xL=xU. Here we return 1 always.
   // Here the class has its own data and a function among them.
   // This seems to be the most powerful type of implementation.
   // The main point is here that in defining the constructor
   // one may pre-process the constructor arguments before the
   // data member f has to be initialized. Initializing f is
   // simply by assignment. Notice that it would be contraproductive
   // to derive InQuant from FuObj<Z,Z> to get operator()(Z).
   // This would make the result an immutable class wheras the
   // present definition gives a class that implements the strict
   // value interface.

class IntQuant{ // Putting integer numbers into boxes

   Z zL;
   Z zU;
   Z nP;
   R spanInv;
   F<Z,Z> f;

public:
   IntQuant(Z xL=1, Z xU=1, Z n=1);
      // setting the limits of the interval and the number of boxes.
      // If n is <1, the actual number of boxes will be 1.
      // xL, xU needs not to be ordered on input.
      // For xL==xU, operator()(Z) will always return 1
   Z operator()(Z i)const{ return f(i);}
   F<Z,Z> toF()const{ return f;}
};

//////////////////// Percentilizer /////////////////////////////////////
// Consider an array a[1],...., a[n] of real numbers. Here we consider
// this as a distribution of the a's over the real axis.
// An object Percentilizer(a) defines characteristic quantities
// of this distribution.
// We have the minimum, the median and the maximum as
// the most obvious landmarkes, and the quartiles (hinges) as similarly
// built constructs ( see Weisstein: CRC Concise Encyclopedia of
// Mathematics p. 844 for hinge).
// All these are special cases of percentiles, i.e. just the numbers
// which the 'Percentilizer' will return via operator(). Consider:
// V<R> a=...;
// Percentilizer per(a);
// R x=0.25;
// R y=per(x);
// y is a number such that 25% of all components of a have
// values below y. An important application of Percentilizer is to
// exclude from a set of numbers a 'few extreme values' that may be
// considered as irregular.

class Percentilizer: public F<R,R>{ // handling quartiles (hinges)

      // Notice that this is an immutable class that does not exactly
      // satisfies the value interface since it does not define a default
      // constructor. Such a constructor could not be of any value.

   Vo<R> arr;
      // ordered list of numbers that define a statistical distribution;
      // In contrast to class Quantizer the basic data elements are not
      // held by the base class. So it is readily available for
      // efficient implementation of further methods needing it.

 public:
   Percentilizer(const V<R>& );
      // setting the list of real values (no ordering implied) that
      // define the distribution, which in turns defines the quantiles

   F<R,R> toBase()const{ return *this;}
      // evident

   R mean(const Iv& iv)const;
      // iv is a subinterval of Iv(0,1) (or Iv(0,1) itself).
      // This selects the components of arr which lie
      // between (*this)(iv[1]) and (*this)(iv[2]). The mean value
      // of these selected components is returned.

   R modulation( const Iv& ivL, const Iv& ivU)const
      // a flexible descriptor of the variability (deviation from
      // constancy) of the values in arr. Allows to ignore
      // extreme (potentially irregular) values in arr.
      // Details clear from code:
      { return (mean(ivU)-mean(ivL))*cpminv(mean(ivU|ivL));}

   V<R> select(const Iv& iv)const;
      // returns the sub-list of arr over which function mean
      // performs the average
};

R gaussRandomValue(const R& mean, const R& sigma, Z j=0);
// 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.
// The last argument has the same meaning as in randomR(Z).

///////////////// some counter-like data types //////////////////////////

//////////////////////// PeriodicAlarm //////////////////////////////////

// Let per be an instance of class PeriodicAlarm. Then for some t of type
// R we may ask for the value of per(t) which is either true or false.
// This boolean return value does not only depend on t but also on the
// values of the t's for which per was called earlier. The situation for
// which PeriodicAlarm is designed and in which its functionality follows
// an easy to describe pattern is that calls to per(t) come in order of
// increasing t. This means that a code fragment
//
//    PeriodicAlarm per(...);
//    R t1=...;
//    bool b1=per(t1);
//       ...
//    R t2=...;
//    bool b2=per(t2);
//       ...
// should ensure that t2>=t1.

// Let the tasks be 'write the state variables of a dynamical system to
// disk' be scheduled for instants of time which form an equidistant
// lattice of width period (period>0) and let t be the time variable of
// the dynamical system that changes from one computational step to the
// next by varying amounts (adaptive time stepping) (normally small
// compared to the fixed step 'period' from task to task).
// Then per(t)==true says that we should call writeState().

class PeriodicAlarm{ // periodic counter
   bool active;
      // if false, the counter is deactivated so that calling operator()
      // gives always 'false' and no further action is taken
   R period;
      // distance between points of alarm
   R lastAlarm;
      // t walue for which (*this)(t) was found 'true' the last time
   Z counter;
      // counts the alarms
   R nextLatticePoint;
      // the time points for alarm form a equidistant lattice (chain)
      // of points on the t-axis. Thus the name
   R tol;
      // tolerance as fraction of period, so that not
      // the condition t>=nextLatticePoint triggers
      // operator()(t)==true but already
      // t>=nextLatticePoint-tol. We always assume period>=0.
      // In the frequently occurring case that t grows in steps of
      // period between calls of operator() this makes the trigger
      // condition independent of numerical chance events.
public:
   PeriodicAlarm(R period=0, R tIni=0, R phase=0.);
      // time for the first alarm is tIni+phase*period. Then the points
      // of alarm follow with a spacing given by period.
      // If period==0, the object is created in not active state.
   bool operator()(R t);
      // gives the alarm status for time t. Is 'true' for
      // (t+tol) >= nextLatticePoint
      // updates nextLatticePoint by adding period the smallest integer
      // multiple of period such that t<nextLatticePoint
   Z count()const{ return counter;}
      // returns the number of calls to operator() with a true result.
      // More precisely, a call can become counted with a value larger
      // than one if t-tLast is larger than period.
   R getLastAlarm()const{ return lastAlarm;}
   void setCounter(Z c){counter=c;}
   void switchOff(){ active=false;}
      // if this is the case, always false==(*this)(t)
   void switchOn(){ active=true;}
   R per()const{ return period;}
};

////////////////////////// SingleAlarm //////////////////////////////////

// alarm is issued only once

class SingleAlarm{ // 1-step counter
   R alarmPoint;
   bool done;


public:
   SingleAlarm(R tAlarm=0):alarmPoint(tAlarm),done(false){}
   bool operator()(R t);
      // gives the alarm status for time t. Is 'true' for t>=alarmPoint
      // sets done=true
};

////////////////////////// CyclicAlarm //////////////////////////////////

// version of a perodic alarm with a discrete periodic counter. Each
// call of the status (operator()) updates the counter

class CyclicAlarm{ // periodic counter

// notice that all attributes are class typed and thus initialized

   B dormant;
      // if this is true, the object does not react on calls
      // to operator(), especially the value of 'alarm'
      // remains unchanged

   B done;
      // if 'true' the one has alarm==false
      // and operator() always returns false

   Z period{0};
      // after period increments of cog we set alarm==true

   Z cog{0}; // 'Zahn','Zahnrad' in German
      // each call to operator() cyclically increases cog
      // cyclic structure given by period.
      // One should have available residue class Z/n at this place

   Z counter{0};
      // each setting alarm==true increases counter

   Z counterMax{0};
      // after counter has reached the value counterMax
      // we have done==true
      // since counter>=0 always, a negative value for
      // counterMax let it never be done (which is considered
      // normal behavior: 'work is never done').

   B alarm;
      // if this is true, alarm is set

public:
   typedef CyclicAlarm Type;
   virtual Type* clone()const{ return new CyclicAlarm(*this);}
   CyclicAlarm(Z per=1, Z cogIni=0, Z maxAlarms=-1);
   // per=0 gives return value false for all calls of
   // operator().
   // Effectively cogIni<period (since enforced).
   // Negative values of cogIni create an initial latency
   // of -cogIni calls to step or operator().
   virtual void step(void);
      // moves *this to the next state
   bool operator()(void)
   {
      if (period<1) return false;
      else{bool alNow=alarm; step(); return alNow;}
   }
      // most convenient usage of the class e.g. in the implementation
      // of function 'record'. Returns the status of alarm at present
      // state and moves *this to the next state
   Z count()const{ return counter;}
   void setAlarm(bool b){ alarm=B(b);}
   bool getAlarm()const{ return alarm;}
   void switchOff(){ dormant=true;}
   void switchOn(){ dormant=false;}
   bool ready()const{ return done;}
   void restart(){ done=B(false); cog=0, counter=0; alarm=B(true);}
   void restartWithNewPeriod(Z per){ period=per; done=false;
      cog=0, counter=0; alarm=true;}
   V<B> record(Z n);
   Z per()const{ return period;}
      // returns the sequence resulting from n-fold application of
      // operator(). Means for testing.
};

////////////////////// NonCyclicAlarm //////////////////////////////////

// cog-based state machine that gives alarm iff
// cog==points_[counter]. Each alarm increases not only
// cog but also counter.
//
// The behavior is thus just the following:
// V<Z> vz=...;
// NonCyclicAlarm nca(vz);
// bool b;
// b=nca();
//    ...
// b=nca();
// ...
// b=nca(); // n-th call of operator()() after constructor call
// results in b==true iff there if a k \in {1,2,..., vz.dim()} such that
// n==vz[k]

class NonCyclicAlarm{ // cog-based state machine

// notice that all attributes are class typed and thus initialized

   B done;
      // if 'true' the one has alarm==false
      // and operator() always returns false

   Z cog{0}; // 'Zahn','Zahnrad' in German
      // each call to operator() increases cog by 1

   Z counter{0};
      // each setting alarm==true increases counter

   V<Z> points_;
      // list of values such that cog==points_[counter]
      // is the point to trigger

   Z counterMax{0};
      // after counter has reached the value counterMax
      // we have done==true

   B alarm;
      // if this is true, alarm is set

public:
   typedef NonCyclicAlarm Type;
   virtual Type* clone()const{ return new NonCyclicAlarm(*this);}
   NonCyclicAlarm(const V<Z>& points=V<Z>(0));
      // condition for alarm==true is cog==points[counter]
      // alarm==true implies counter->counter+1
      // Der Bedeutung von 'points' entsprechend sollte dieses
      // eine strict monoton wachsende Liste sein. points_
      // wird aus points durch umordnen und Weglassen von
      // Doubletten erhalten. Komponenten von points_, die
      // kleiner als 1 sind, werden nicht eliminiert, führen
      // aber dazu daß immer alarm==false gilt.
   virtual void step(void);
      // moves *this to the next state: cog->cog+1
   bool operator()(void){ bool alNow=alarm; step(); return alNow;}
      // most convenient usage of the class.
      // Returns the status of alarm at present
      // state and moves *this to the next state
   Z count()const{ return counter;}
   bool getAlarm()const{ return alarm;}
   bool ready()const{ return done;}
   void restart(){ done=false; cog=0, counter=1; alarm=false;}
   Z last()const{ return points_.last();}
};

////////////////////////// BurstAlarm //////////////////////////////////

// Similar to Cyclic alarm but opertor() gives true for if the
// call number belongs to to a doubly periodic lattice
// 2004-05-11: Triggers in some cases also if per=0. This is not
// intended. Not completely understood!. Better don't use it!

class BurstAlarm: public CyclicAlarm{ // double periodic counter

  CyclicAlarm burst;

public:
   typedef BurstAlarm Type;
   typedef CyclicAlarm Base;
   virtual CyclicAlarm* clone()const{ return new BurstAlarm(*this);}

   BurstAlarm(Z per=1, Z burstPer=0, Z burstLen=0);
      // We exemplify the meaning of the data by considering the
      // sequence created by function record (1=true, 0=false)
      // e.g. for
      //       per=10,  burstPer=2,    burstLen=3:
      // | 1 0 | 1 0 | 1 0 | 0 0 0 0 | 1 0 | 1 0 | 1 0 | 0 0 0 0
      // | 1 0 | 1 0 | 1 0 | 0 0 0 0 ...
      // one Burst is | 1 0 | 1 0 | 1 0 |. Period from burst beginning
      // to burst beginning is 10
   void step(void);
      // redefined
};

////////////////////////////// enum Switch ////////////////////////////

enum Switch{ //  to be used as an argument type in cases where the values
      // 'true' and 'false' of a boolean argument would not suggest the
      // meaning with sufficient clarity
   ON,
   OFF
};

//////////////////////////// getNerFrc_1(...) ////////////////////////////

Z getNerFrc_1(R& x, R xL);
   //: get nearest fraction
   // Notice that '_1' indicates that the first argument may be changed by
   // by a function call. Assumed is x<=xL and x!=0.
   // Consider  S(xL):={ y: there is n \in { 1,2,...} such that y*n=xL}
   // The function outputs the n and changes x into
   // x0 \in S(xL) such that |x0-x|=inf{ |s-x| : s \in S(xL)}.

//////////////////////////////// genSqrRoot(Z) ///////////////////////////

Z2 genSqrRoot(Z n);
//: generalized square root
// The square root in R solves the problem to  write a number x as a
// product of two factors which are equal.
// If we consider only natural numbers, a factorization with just
// this properties is possible only for the square numbers
// (1,4,9,16,...). It remains a meaningful problem to find two
// factors which 'differ as little as possible' and represent
// x as a product 'as accurate as possible'. It is this analogy
// to the square root to which the name refers, the generalization
// thus also comprises a specialization --- the specialization to
// integer or natural numbers.
// Let (z1,z2) be the return value of the function. Then z1*z2>=n, z2>=z1.
// Among the infinite number of solutions of these inequalities
// we choose one for which d1:=z1*z2-n and d2:=z2-z1 are both as small as
// they together can be. For instance, if n is a square number n=m*m, then
// z1=z2=m are oviously what we want since then d1=d2=0.
// For n=7 we have for z1=1, z2=7 the desirable property d1=0, but
// d2=7-1=6 is not small. Here we could consider the pairs  (2,4), (3,3).
// In the first case we have d1=8-7=1, d2=4-2=2, and
// in the second case d1=9-7=2, d2=3-3=0.
// d1+d2 is smaller in the second case and so (3,3) will be selected.
// See the implementation in cpmalgorithms.cpp for details.
// Defining this function was suggested by a graphical problem. Given
// n scalable figures; how to layout them on screen if only
// matrix subdivisions of the screen are under consideration.

Z appSqrRoot(Z n);
//: approximate square root
// simply round cpmsqr(R(n) to the nearest integer. n>=0 assumed

////////////////// operator norm of a linear mapping //////////////////////
/*
   X is assumed to define
      R X::abs()const
   and
      X operator*(R const&)const
   The intended application is that X is a linear space over R and f
   is a linear mapping X --> X. Then Norm<X>(f).run() yields the smalles N
   in R such that f(x).abs() <= N * x.abs() for all x in X.
*/

template<typename X>
class Norm{
   F<X,X> f_;
      // the function corresponding to f of the previous expanation
   R eps_; // accuracy measure which says when to finish the iteration
   Z c_; // 'complexity', mainly size of the starting element of X.
   Z j_; // controls the random selection.Clear from the first three
      // lines of code in run()

public:
   Norm(F<X,X> const& f, R eps=1.e-12, Z const& c=10, Z const& j=137):
      f_(f), eps_(eps), c_(c), j_(j){}
   R run()const{
      X x0;
      X x1=x0.test(c_);
      X x2=x1.ran(j_);
      X psiOld=x2;
      X chi=f_(psiOld);
      R dis1=1000_R;
      R dis2=dis1;
      R etaOld=-1_R;
      R etaNew=chi.abs();
      X psiNew=chi*inv(etaNew);
      while(cpmabs(dis1+dis2)>eps_){
         dis1=etaNew-etaOld;
         psiOld=psiNew;
         etaOld=etaNew;
         chi=f_(psiOld);
         etaNew=chi.abs();
         psiNew=chi*inv(etaNew);
         dis2=etaNew-etaOld;
         cout<<"cpmabs(dis1+dis2) = "<<cpmabs(dis1+dis2)<<endl;
      }
      return (etaOld+etaNew)*0.5_R;
   }
};


template<typename T>
concept disCpm =
requires (T x, T y)
{
   x.dis(y);
};

template<typename T>
concept nameOfCpm =
requires (CpmRoot::Word res, T x)
{
   res=x.nameOf();
};

} // namespace

#endif