multishift conjugate gradient added and a strong test: take a diagonal

but non-identity matrix
l1 0  0  0 ....
0  l2 0  0 ....
0  0  l3 0 ...
.  .   .
.  .   .
.  .   .

And apply the multishift CG to it. Sum the poles and residues.
Insist that this be the same as the exactly taken square root
where l1,l2,l3 >= 0.

lib/algorithms/approx/MultiShiftFunction.cc

@@ -0,0 +1,29 @@
+#include <Grid.h>
+
+namespace Grid {
+double MultiShiftFunction::approx(double x)
+{
+  double a = norm;
+  for(int n=0;n<poles.size();n++){
+    a = a + residues[n]/(x+poles[n]);
+  }
+  return a;
+}
+void MultiShiftFunction::gnuplot(std::ostream &out)
+{
+  out<<"f(x) = "<<norm<<"";
+  for(int n=0;n<poles.size();n++){
+    out<<"+("<<residues[n]<<"/(x+"<<poles[n]<<"))";
+  }
+  out<<";"<<std::endl;
+}
+void MultiShiftFunction::csv(std::ostream &out)
+{
+  for (double x=lo; x<hi; x*=1.05) {
+    double f = approx(x);
+    double r = sqrt(x);
+    out<< x<<","<<r<<","<<f<<","<<r-f<<std::endl;
+  }
+  return;
+}
+}

lib/algorithms/approx/MultiShiftFunction.h

@@ -0,0 +1,28 @@
+#ifndef MULTI_SHIFT_FUNCTION
+#define MULTI_SHIFT_FUNCTION
+namespace Grid {
+class MultiShiftFunction {
+public:
+  int order;
+  std::vector<RealD> poles;
+  std::vector<RealD> residues;
+  std::vector<RealD> tolerances;
+  RealD norm;
+  RealD lo,hi;
+  MultiShiftFunction(int n,RealD _lo,RealD _hi): poles(n), residues(n), lo(_lo), hi(_hi) {;};
+  RealD approx(RealD x);
+  void csv(std::ostream &out);
+  void gnuplot(std::ostream &out);
+  MultiShiftFunction(AlgRemez & remez,double tol,bool inverse) :
+      order(remez.getDegree()),
+      tolerances(remez.getDegree(),tol),
+      poles(remez.getDegree()),
+      residues(remez.getDegree())
+  {
+    remez.getBounds(lo,hi);
+    if ( inverse ) remez.getIPFE (&residues[0],&poles[0],&norm);
+    else remez.getPFE (&residues[0],&poles[0],&norm);
+  }
+};
+}
+#endif

lib/algorithms/approx/Remez.h

@@ -125,8 +125,17 @@ class AlgRemez
   // Destructor
   virtual ~AlgRemez();
 
+  int getDegree(void){ 
+    assert(n==d);
+    return n;
+  }
   // Reset the bounds of the approximation
   void setBounds(double lower, double upper);
+  // Reset the bounds of the approximation
+  void getBounds(double &lower, double &upper) { 
+    lower=(double)apstrt;
+    upper=(double)apend;
+  }
 
   // Generate the rational approximation x^(pnum/pden)
   double generateApprox(int num_degree, int den_degree,