// ------------------------------------------------------------------
// qpatchmath.cpp
//
// This file contains routines for math with patches including:
// converting parametric to real coordinates, computing directional
// derivatives, computing normals, and intersecting patches with
// lines.
// ------------------------------------------------------------------
// Copyright (c) 1999 by Cornell University.  All rights reserved.
// 
// See the accompanying file 'Copyright' for authorship information,
// the terms of the license governing this software, and disclaimers
// concerning this software.
// ------------------------------------------------------------------
// This file is part of the QMG software.  
// Version 2.0 of QMG, release date September 3, 1999
// Version 2.0.1 -- rewritten to estimate generalized eigenvector
//  sensitivity, November 10, 1999.  Steve Vavasis
// ------------------------------------------------------------------
#ifdef _MSC_VER
#if _MSC_VER == 1200
#pragma warning (disable: 4786) 
#pragma warning (disable: 4788)
#endif
#endif





#include "qpatchmath.h"
#include "qmatvec.h"



#ifdef DEBUGGING
namespace QMG {
  namespace Local_to_PatchMath_CPP {
    using namespace QMG;
    ostream* debug_str;
    bool dump_everything = false;
  }
}
#endif


// Types used throughout this file.


// Types passed to Lapack routines

#ifdef LAPACK_DOUBLECOMPLEX_T
typedef LAPACK_DOUBLECOMPLEX_T Lapack_doublecomplex;
#else
typedef double Lapack_doublecomplex;
#endif

#ifdef LAPACK_INT_T
typedef LAPACK_INT_T Lapack_int;
#else
typedef long int Lapack_int;
#endif

typedef std::complex<double> Complex;

// LAPACK routine for complex standard nonsymmetric eigenvalues.


extern "C"
int zgeev_(char *jobvl, 
           char *jobvr, 
           Lapack_int *n,
           Lapack_doublecomplex *a, 
           Lapack_int *lda, 
           Lapack_doublecomplex *w, 
           Lapack_doublecomplex *vl,
           Lapack_int *ldvl, 
           Lapack_doublecomplex *vr, 
           Lapack_int *ldvr, 
           Lapack_doublecomplex* work,
           Lapack_int *lwork, 
           double *rwork, 
           Lapack_int *info);

// LAPACK routine for qr factorization

extern "C"
int zgeqrf_(Lapack_int *m, 
            Lapack_int *n, 
            Lapack_doublecomplex *a,
            Lapack_int *lda, 
            Lapack_doublecomplex *tau, 
            Lapack_doublecomplex *work, 
            Lapack_int *lwork,
            Lapack_int *info);

// LAPACK routine for applying Housholder reflections to a matrix.

extern "C"
int zunmqr_(char *side, 
            char *trans, 
            Lapack_int *m, 
            Lapack_int *n,
            Lapack_int *k, 
            Lapack_doublecomplex *a, 
            Lapack_int *lda, 
            Lapack_doublecomplex *tau,
            Lapack_doublecomplex *c, 
            Lapack_int *ldc, 
            Lapack_doublecomplex *work, 
            Lapack_int *lwork,
            Lapack_int *info);

// LAPACK routine for multiplying out a product of Householder reflections.

extern "C"
int zungqr_(Lapack_int *m, 
            Lapack_int *n, 
            Lapack_int *k,
            Lapack_doublecomplex *a, 
            Lapack_int *lda, 
            Lapack_doublecomplex *tau, 
            Lapack_doublecomplex *work, 
            Lapack_int *lwork, 
            Lapack_int *info);

// Lapack routine for pencil Hessenberg reduction.

extern "C" 
int zgghrd_(char *compq, 
            char *compz, 
            Lapack_int *n, 
            Lapack_int *ilo, 
            Lapack_int *ihi, 
            Lapack_doublecomplex *a, 
            Lapack_int *lda, 
            Lapack_doublecomplex *b,
            Lapack_int *ldb, 
            Lapack_doublecomplex *q, 
            Lapack_int *ldq, 
            Lapack_doublecomplex *z,
            Lapack_int *ldz, 
            Lapack_int *info);

// Lapack routine for QZ factorization.


extern "C" 
int zhgeqz_(char *job, 
            char *compq, 
            char *compz, 
            Lapack_int *n,
            Lapack_int *ilo, 
            Lapack_int *ihi, 
            Lapack_doublecomplex *a, 
            Lapack_int *lda, 
            Lapack_doublecomplex *b, 
            Lapack_int *ldb, 
            Lapack_doublecomplex *alpha, 
            Lapack_doublecomplex *beta, 
            Lapack_doublecomplex *q, 
            Lapack_int *ldq, 
            Lapack_doublecomplex *z, 
            Lapack_int *ldz,
            Lapack_doublecomplex *work, 
            Lapack_int *lwork, 
            double *rwork, 
            Lapack_int *info);


// Helper class that makes private routines public for this file only.

class QMG::PatchMath::Bezier_Helper {
public:
  static Real control_point_coord(const PatchMath& pm, int cpnum, int d) {
    return pm.control_point_coord_(cpnum,d);
  }
};

// Class Workspace is a workspace for math operations.
// Used so that we don't have to keep calling 'new' to
// get workspace for Lapack.
// Used as follows: Set a marker.  Then allocate vectors
// and matrices.  When the marker is destroyed, the space
// is deallocated.


class QMG::PatchMath::Workspace::RealVector {
private:
  Workspace& wkspa_;
  int base_;
  int sz_;
  // no copying, no assignment.
  RealVector(const RealVector& o) : wkspa_(o.wkspa_)  { }
  void operator=(const RealVector&) { }
public:
  RealVector(Workspace& wkspa, int sz) : wkspa_(wkspa),
    base_(wkspa.rwkspa_.size()),
    sz_(sz) {
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
    wkspa.rwkspa_.resize(base_ + sz_);
  }
  ~RealVector() { }
  Real& operator[](int i) {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.rwkspa_[base_ + i];
  }
  double* make_fortran_real() {
    return &(wkspa_.rwkspa_[base_]);
  }
};



// Matrices stored by column for fortran compatibility.

class QMG::PatchMath::Workspace::ComplexMatrix {
private:
  Workspace& wkspa_;
  int base_;
  int nr_;
  int nc_;
  // no copying, no assignment
  ComplexMatrix(const ComplexMatrix& o) : wkspa_(o.wkspa_) { }
  void operator=(const ComplexMatrix&) { }
public:
  ComplexMatrix(Workspace& wkspa, int nr, int nc) : wkspa_(wkspa),
    base_(wkspa.cwkspa_.size()),
    nr_(nr),
    nc_(nc) {
    wkspa.cwkspa_.resize(base_ + nr * nc);
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~ComplexMatrix() { }
  Complex& operator()(int i, int j) {
#ifdef RANGECHECK
    if (i < 0 || i >= nr_ || j < 0 || j >= nc_)
      throw_error("Range err");
#endif
    return wkspa_.cwkspa_[base_ + i + j * nr_];
  }
  Lapack_doublecomplex* make_fortran_complex() {
    return reinterpret_cast<Lapack_doublecomplex*>(
      static_cast<Complex*>(&(wkspa_.cwkspa_[base_])));
  }
  int numrow() const { return nr_; }
  int numcol() const { return nc_; }
};



class QMG::PatchMath::Workspace::ComplexVector {
private:
  Workspace& wkspa_;
  int base_;
  int sz_;
  // no copying, no assignment
  ComplexVector(const ComplexVector& o) : wkspa_(o.wkspa_) { }
  void operator=(const ComplexVector&) { }
public:
  ComplexVector(Workspace& wkspa, int sz) : wkspa_(wkspa),
    base_(wkspa.cwkspa_.size()),
    sz_(sz) {
    wkspa.cwkspa_.resize(base_ + sz_, Complex(0.0,0.0));
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~ComplexVector() { }
  Complex& operator[](int i) {
#ifdef RANGECHECK
    if (i < 0 || i >= sz_)
      throw_error("Range err");
#endif
    return wkspa_.cwkspa_[base_ + i];
  }
  Lapack_doublecomplex* make_fortran_complex() {
    return reinterpret_cast<Lapack_doublecomplex*>(
      static_cast<Complex*>(&(wkspa_.cwkspa_[base_])));
  }

};

class QMG::PatchMath::Workspace::RealMatrix {
private:
  Workspace& wkspa_;
  int base_;
  int nr_;
  int nc_;
  // no copying, no assignment
  RealMatrix(const RealMatrix& o) : wkspa_(o.wkspa_)  { }
  void operator=(const RealMatrix&) { }
public:
  RealMatrix(Workspace& wkspa, int nr, int nc) : wkspa_(wkspa),
    base_(wkspa.rwkspa_.size()),
    nr_(nr),
    nc_(nc) {
    wkspa.rwkspa_.resize(base_ + nr * nc);
#ifdef RANGECHECK
    if (!wkspa.marker_active_)
      throw_error("Attempt to allocate in inactive workspace");
#endif
  }
  ~RealMatrix() { }
  Real& operator()(int i, int j) {
#ifdef RANGECHECK
    if (i < 0 || i >= nr_ || j < 0 || j >= nc_)
      throw_error("Range err");
#endif
    return wkspa_.rwkspa_[base_ + i + j * nr_];
  }
  Real* make_fortran_real() {
    return &(wkspa_.rwkspa_[base_]);
  }
  int numrow() const { return nr_; }
  int numcol() const { return nc_; }
};


class QMG::PatchMath::Workspace::StartPosMarker {
private:
  Workspace& wkspa_;
  int rpos_;
  int cpos_;
  StartPosMarker(const StartPosMarker& o) : wkspa_(o.wkspa_) { }
  void operator=(const StartPosMarker&) { }
public:
  explicit StartPosMarker(Workspace& wkspa) : wkspa_(wkspa),
    rpos_(wkspa.rwkspa_.size()),
    cpos_(wkspa.cwkspa_.size()) {
    ++wkspa.marker_active_;
  }
  ~StartPosMarker() {
    wkspa_.rwkspa_.resize(rpos_);
    wkspa_.cwkspa_.resize(cpos_);
    --wkspa_.marker_active_;
  }
};



namespace {
  using namespace QMG;



  // Class for printing out a complex number.

  class Format_Complex;
  ostream& operator<<(ostream&, const Format_Complex&);

  class Format_Complex {
  private:
    Complex c;
  public:
    explicit Format_Complex(const Complex& c1) : c(c1) { }
    friend ostream& operator<<(ostream&, const Format_Complex&);
  };

  
  ostream& operator<<(ostream& os, const Format_Complex& fc) {
    if (fc.c.imag() >= 0)
      os << std::setprecision(17) << fc.c.real() << "+" << fc.c.imag() << "i";
    else
      os << std::setprecision(17) << fc.c.real() << "-" << -fc.c.imag() << "i";
    return os;
  }





  // ------------------------------------------------------------------
  // This routine carries out QR factorization on a matrix.
  // If sing_cutoff is nonzero, it also tries to find
  // dependent columns and quits as soon as a dependent column is found.
  // If it quits early then it returns the number of columns processed.


  void QR_Factor(PatchMath::Workspace::RealMatrix& A, //matrix to factor
    PatchMath::Workspace::RealMatrix& hhmat,
    PatchMath::Workspace::RealVector& betavec) {

    int m = A.numrow();
    int n = A.numcol();
    int nstep = (m < n)? m : n;
    for (int k = 0; k < nstep; ++k) {
      Real* hhmat_col = &hhmat(0,k);

      Real beta = 
        Matrix::compute_hh_transform(&hhmat_col[k], &A(k,k), 1, m - k);
      Real mu = (beta == 0.0)? 0 : (1.0 / fabs(beta * hhmat_col[k]));
      betavec[k] = beta;
      {
        for (int j = k; j < n; ++j) {
          Real ip1 = 0.0;
          {
            for (int i = k; i < m; ++i)
              ip1 += A(i,j) * hhmat_col[i];
          }
          ip1 *= beta;
          {
            for (int i = k; i < m; ++i)
              A(i,j) -= ip1 * hhmat_col[i];
          }
        }
      }
    }
  }


  void QR_Factor_with_col_pivoting(PatchMath::Workspace::RealMatrix& A, 
    PatchMath::Workspace::RealVector& select_col, 
    PatchMath::Workspace::RealMatrix& hhmat,
    PatchMath::Workspace::RealVector& betavec) {

    int m = A.numrow();
    int n = A.numcol();
    int nstep = (m < n)? m : n;

    for (int k = 0; k < nstep; ++k) {
      int select_col1;
      {
        Real mx1 = 0.0;
        for (int j = 0; j < n; ++j) {
          bool j_already_selected = false;
          {
            for (int i = 0; i < k; ++i) {
              if (static_cast<int>(select_col[i]) == j) {
                j_already_selected = true;
                break;
              }
            }
          }
          if (j_already_selected) continue;
          
          // Find the norm of column j residual entries.
          Real nr1 = 0.0;
          {
            for (int i = k; i < m; ++i) {
              nr1 += A(i,j) * A(i,j);
            }
          }
          if (nr1 >= mx1) {
            mx1 = nr1;
            select_col1 = j;
          }
        }
      }
      select_col[k] = static_cast<Real>(select_col1);
      if (k == m - 1)
        break;
      Real* hhmat_col = &hhmat(0,k);
      Real beta = 
        Matrix::compute_hh_transform(&hhmat_col[k], &A(k, select_col1), 1, m - k);
      betavec[k] = beta;
      {
        for (int j = 0; j < n; ++j) {
          Real ip1 = 0.0;
          {
            for (int i = k; i < m; ++i)
              ip1 += A(i,j) * hhmat_col[i];
          }
          ip1 *= beta;
          {
            for (int i = k; i < m; ++i)
              A(i,j) -= ip1 * hhmat_col[i];
          }
        }
      }
    }
  }

  // Multiply HH transforms together to get an explicit
  // orthogonal matrix Q.

  void formQ(PatchMath::Workspace::RealMatrix& hhmat,
    PatchMath::Workspace::RealVector& betavec,
    PatchMath::Workspace::RealMatrix& Q) {

    int m = hhmat.numrow();
    int n = hhmat.numcol();

    // Initialize Q to be the identity.

    {
      for (int i = 0; i < m; ++i)
        for (int j = 0; j < m; ++j)
          Q(i,j) = static_cast<Real>(static_cast<int>(i == j));
    }

    // Loop over HH transforms in reverse order 

    for (int k = n - 1; k >= 0; --k) {
      
      // Apply the kth HH transform to Q.
      // j loops over columns of Q.
      for (int j = k; j < m; ++j) {
        Real ip1 = 0.0;
        {
          for (int i = k; i < m; ++i)
            ip1 += hhmat(i,k) * Q(i,j);
        }
        ip1 *= betavec[k];
        {
          for (int i = k; i < m; ++i)
            Q(i,j) -= ip1 * hhmat(i,k);
        }
      }
    }
  }

        
  // This is a helper class for math on 
  // a one-variable Bezier.

  class Solve_OneVar_Bezier {
  public:
    virtual Real get_bez_coeff(int j) const = 0;
    Solve_OneVar_Bezier() { }
    virtual ~Solve_OneVar_Bezier() { }
    void solve(int degree,
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace, 
      double scfac, 
      double tol) const;
    Real eval(Real param,
      int degree) const;
    Real eval_deriv(Real param,
      int degree) const;
  };

  // Derived class for math on a one-variable Bezier that arises
  // from intersecting a line with a curve in 2D.

  class Solve_OneVar_Bezier1 : public Solve_OneVar_Bezier {
  private: 
    const Matrix& lineeq_;
    const Point& linerhs_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier1(const Matrix& lineeq,
      const Point& linerhs,
      const PatchMath& pm) : lineeq_(lineeq), linerhs_(linerhs), pm_(pm) { }
    Real get_bez_coeff(int cpnum) const {
      return lineeq_(0,0) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        lineeq_(0,1) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 1) -
        linerhs_[0];
    }
  };

  // Derived class for math on a one-variable Bezier that
  // arises from intersecting a curve with a plane in 3D.

  class Solve_OneVar_Bezier2 : public Solve_OneVar_Bezier {
  private:
    const Point& plane_normal_;
    Real plane_rhs_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier2(const Point& plane_normal,
      Real plane_rhs,
      const PatchMath& pm) : plane_normal_(plane_normal), plane_rhs_(plane_rhs), 
      pm_(pm) { }
    Real get_bez_coeff(int cpnum) const {
      return plane_normal_[0] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        plane_normal_[1] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum,1) +
        plane_normal_[2] * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum,2) -
        plane_rhs_;
    }
  };

  // Derived class for evaluating one dimension of a Bezier curve.

  class Solve_OneVar_Bezier3 : public Solve_OneVar_Bezier {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };


  // Derived class for evaluating the derivative of a Bezier curve.


  class Solve_OneVar_Bezier4 : public Solve_OneVar_Bezier {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
  public:
    Solve_OneVar_Bezier4(const Point& direc, int di,
      const PatchMath& pm) : direc_(direc), di_(di), pm_(pm) 
    { }
    Real get_bez_coeff(int cpnum) const {
      Real t = 0.0;
      for (int j = 0; j < di_; ++j) {
        t += (PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum + 1, j) -
          PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, j)) *
          direc_[j];
      }
      return t;
    }
  };


  // Evaluate a one-variable Bezier curve using Horner algorithm.
  

  Real Solve_OneVar_Bezier::eval(Real param,
    int degree) const {
    Real u = 1.0 - param;
    Real pwrv = 1.0; 
    Real result = 0.0;
    int binom = 1;
    for (int j = 0; j < degree + 1; ++j) {
      result *= u;
      result += binom * pwrv * get_bez_coeff(j);      
      pwrv *= param;
      binom *= (degree - j);
      binom /= (j + 1);
    }
    return result;
  }
  
  // Evaluate derivative a one-variable Bezier curve using Horner algorithm.
  

  Real Solve_OneVar_Bezier::eval_deriv(Real param,
    int degree) const {
    Real u = 1.0 - param;
    Real pwrv = 1.0; 
    Real result = 0.0;
    int binom = 1;
    Real prevbez = get_bez_coeff(0);
    for (int j = 0; j < degree; ++j) {
      result *= u;
      Real nextbez = get_bez_coeff(j + 1);
      result += binom * pwrv * (nextbez - prevbez);
      prevbez = nextbez;
      pwrv *= param;
      binom *= (degree - j - 1);
      binom /= (j + 1);
    }
    return result * degree;
  }


  // solve: find zeros of a one-variable Bezier.
  // Transforms to an eigenvalue problem and invokes LAPACK.
  // Doesn't set output_stack real coord fields -- caller must do this.
 
  void Solve_OneVar_Bezier::solve(int degree,
    vector<PatchMath::IntersectRecord>& output_stack,
    PatchMath::Workspace& workspace, 
    double scfac, 
    double tol) const {


#ifdef DEBUGGING
    using QMG::Local_to_PatchMath_CPP::debug_str;
    using QMG::Local_to_PatchMath_CPP::dump_everything;

    if (dump_everything) {
      *debug_str << " In solve bezier degree = " << degree 
        << " scfac = " << scfac << " tol = " << tol << "\n bezcoeff = [";
      for (int j = 0; j < degree + 1; ++j) 
        *debug_str << std::setprecision(17) << get_bez_coeff(j) << " ";
      *debug_str << "]\n";
    }
#endif

    // If the bezier control points miss 0, then early return

    {
      bool pos_seen = false;
      bool neg_seen = false;
      for (int i = 0; i <= degree; ++i) {
        Real bz = get_bez_coeff(i);
        if (bz <= scfac * tol)
          neg_seen = true;
        if (bz >= -scfac * tol)
          pos_seen = true;
      }
      if (!pos_seen || !neg_seen) {
#ifdef DEBUGGING
        if (dump_everything)
          *debug_str << "origin missed; returning early";
#endif
        return;
      }
    }


    if (degree == 1) {

      // In the linear case just solve the one-variable equation.

      Real denom = get_bez_coeff(0) - get_bez_coeff(1);
      if (denom == 0) {
        return;
      }
      Real sol = get_bez_coeff(0) / denom;
      Real uncer = tol * scfac / fabs(denom);
      if (sol < -uncer || sol > 1.0 + uncer)
        return;
      PatchMath::IntersectRecord result;
      result.paramcoord[0] = sol;
      result.param_uncertainty = uncer;
      result.degen = fabs(sol) <= uncer || fabs(sol - 1.0) <= uncer;
      output_stack.push_back(result);
      return;
    }
    else {
      using std::abs;
      using std::setprecision;
      PatchMath::Workspace::StartPosMarker startposmarker_(workspace);
      
      int degree1 = degree + 1;

      // Multiply by binomial coefs.

      PatchMath::Workspace::RealVector bvec(workspace, degree1);
     
      Real bezmax = 0;
      {
        int binom = 1;
        for (int i = 0; i < degree + 1; ++i) {
          bvec[i] = get_bez_coeff(i) * binom;
          binom *= (degree - i);
          binom /= (i + 1);
          if (fabs(bvec[i]) > bezmax)
            bezmax = fabs(bvec[i]);
        }
      }

      if (bezmax == 0) return;
      
      // grab space in the workspace.
      
      Lapack_int lapack_cworkspace_size = degree * 8 + 
        ((degree < 64)? 64 : degree) * degree + 256;
      
      PatchMath::Workspace::ComplexVector powerv(workspace, degree1);
      PatchMath::Workspace::ComplexVector test_shiftpoly(workspace, degree1);
      PatchMath::Workspace::ComplexMatrix A(workspace, degree, degree);
      PatchMath::Workspace::ComplexMatrix A0(workspace, degree, degree);
      PatchMath::Workspace::ComplexVector test_eigs(workspace, degree);
      PatchMath::Workspace::ComplexVector best_eigs(workspace, degree);
      PatchMath::Workspace::RealVector 
        lapack_rworkspace(workspace, degree * 2);
      PatchMath::Workspace::ComplexVector lapack_cworkspace(workspace, 
        lapack_cworkspace_size);

      // Make a substitution [u;v] = [1, zs ; -conj(zs), 1] * [u'; v']
      // where a is a random complex number.  This is to make sure
      // the leading coef does not vanish.

      Complex best_zs;
      Complex best_conjzs;
      Real best_lead = -1.0;

      for (int trycount = 0; trycount < 5; ++trycount) {

        // Choose the shift.

        Complex zs = Complex(random_real() + 0.5, random_real() + 0.5);
        Complex conjzs = Complex(zs.real(), -zs.imag());

#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "trycount = " << trycount 
                     << "\n zs = " << Format_Complex(zs)
                     << "\n conjzs = " << Format_Complex(conjzs)
                     << '\n';
        }
#endif

        // To make the substitution, use a variant of Horner's rule
        //   powerv is a polynomial that holds ( -conj(zs)*u'+ v') ^k
        // for increasing powers of k.
        //   shifpoly is the result of the substitution.

        {
          for (int j = 0; j < degree1; ++j)
            powerv[j] = Complex(0.0, 0.0);
          powerv[0] = Complex(1.0, 0.0);
        }

        {
          for (int j = 0; j < degree1; ++j)
            test_shiftpoly[j] = Complex(0.0, 0.0);
          test_shiftpoly[0] = Complex(bvec[degree], 0.0);
        }
        // Now repeatedly multiply shiftpoly by u, and then add v^k, for
        // increasing values of k.

        {
          for (int k = 0; k < degree; ++k) {
            
            // Multiply shiftpoly by u, which is u' + zs*v'
            {
              for (int l = k + 1; l > 0; --l) {
                test_shiftpoly[l] = test_shiftpoly[l - 1] + 
                  test_shiftpoly[l] * zs;
              }
              test_shiftpoly[0] = test_shiftpoly[0] * zs;
            }
          

            // Multiply powerv by v, which is  -conj(zs)  * u' + v'

            {
              for (int l = k + 1; l > 0; --l) {
                powerv[l] -= powerv[l - 1] * conjzs;
              }
            }
            // Add bez[k+1] * powerv to shiftpoly.

            {
              for (int l = 0; l <= k + 1; ++l) {
                test_shiftpoly[l] += bvec[degree - k - 1] * powerv[l];
              }
            }
          }
        }

        // If the leading coef of shiftpoly is too small, must try again.

        Complex ldgcoef = test_shiftpoly[degree];
        Real test_lead = abs(ldgcoef);

#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "after loop\npowerv = [\n";
          for (int l = 0; l <= degree; ++l)
            *debug_str << " " << Format_Complex(powerv[l]) << '\n';
          *debug_str << "]\ntest_shiftpoly = [\n";
          for (int l2 = 0; l2 <= degree; ++l2)
            *debug_str << " " << Format_Complex(test_shiftpoly[l2]) << '\n';
          *debug_str << "]\n bezmax = " << bezmax 
                     << "\n test_lead = " << test_lead
                     << '\n';
        }
#endif


        // Set up the companion matrix for shiftpoly.
        {
          for (int i = 0; i < degree - 1; ++i) {
            for (int j = 0; j < degree; ++j) {
              A(i,j) = Complex(0.0,0.0);
              A0(i,j) = Complex(0.0,0.0);
            }
            A(i,i+1) = Complex(1.0,0.0);
            A0(i,i+1) = Complex(1.0,0.0);
          }
        }
        {
          for (int j = 0; j < degree; ++j) {
            A(degree - 1, j) = -test_shiftpoly[j] / ldgcoef;
            A0(degree - 1, j) = -test_shiftpoly[j] / ldgcoef;
          }
        }



#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "A = [";
          for (int i = 0; i < degree; ++i) {
            for (int j = 0; j < degree; ++j) {
              *debug_str << "  " << Format_Complex(A(i,j));
            }
            *debug_str << '\n';
          }
          *debug_str << "]\n";
        }
#endif


        // Call lapack.

        Lapack_int lapack_info = 100000;
        Lapack_int degree_l = degree;
        Lapack_int lda = degree;
        Lapack_int ldvl = degree;
        zgeev_("N", "N", &degree_l, A.make_fortran_complex(), &lda,
          test_eigs.make_fortran_complex(), 0, &degree_l, 0, &ldvl, 
          lapack_cworkspace.make_fortran_complex(), 
          &lapack_cworkspace_size, 
          lapack_rworkspace.make_fortran_real(), 
          &lapack_info);
        // done with loop if lapack successful.

        if (lapack_info != 0) {

          // If lapack failure, print some info.
          
#ifdef DEBUGGING
          *debug_str << "Lapack zgeev failure code = " << lapack_info << '\n';
          *debug_str << "initial A =\n";
          {
            for (int i = 0; i < degree; ++i) {
              for (int j = 0; j < degree; ++j) {
                *debug_str << Format_Complex(A0(i,j)) << '\n';
              }
            }
          }
          *debug_str << "final A =\n";
          {
            for (int i = 0; i < degree; ++i) {
              for (int j = 0; j < degree; ++j) {
                *debug_str << Format_Complex(A(i,j)) << '\n';
              }
            }
          }
#endif
          continue;
        }

        // Store the eigenvalues if this is the first try
        // or the best try.

        if (best_lead < 0 || test_lead > best_lead) {
          best_lead = test_lead;
          best_zs = zs;
          best_conjzs = conjzs;
          for (int i = 0; i < degree; ++i)
            best_eigs[i] = test_eigs[i];
          if (test_lead > bezmax / 10)
            break;
        }
      }
      if (best_lead < 0)
        throw_error("Repeated LAPACK zgeev failure");


      // Eigenvals successfully computed and stored in eig.  Now
      // output the roots in [0,1].


      for (int rootidx = 0; rootidx < degree; ++rootidx) {
        Complex sol0 = best_eigs[rootidx];
        // Apply i transformation to [sol0;1] to get u,v.
        Complex u = sol0 + best_zs;
        Complex v = -best_conjzs * sol0 + Complex(1.0,0.0);

        // If way outside the unit interval, skip it.

        if (abs(u) >= 3 * abs(u + v)) continue;

        Complex sol = u / (u + v);
        Complex sol2 = Complex(1.0,0.0) - sol;

#ifdef DEBUGGING
        if (dump_everything)
          *debug_str << " sol0 = " << Format_Complex(sol0)
                     << "\n u = " << Format_Complex(u) 
                     <<  "\n v = " << Format_Complex(v) 
                     << "\n sol = " << Format_Complex(sol) 
                     << "\n sol2 = " << Format_Complex(sol2) 
                     << '\n';
#endif



        // Compute the residual, derivative, and
        // absolute evaluation of the polynomial at
        // sol using Horners rule.  This is for
        // estimating the uncertainty in the ansewr.

        Complex pwrv(1.0,0.0);
        Complex resid(0.0,0.0);
        Complex deriv(0.0,0.0);
        {
          Real prevbez = get_bez_coeff(0);
          Real abssol2 = abs(sol2);
          int binom = 1;
          int binom2 = 1;
          for (int j = 0; j < degree + 1; ++j) {
            resid *= sol2;
            resid += pwrv * get_bez_coeff(j) * static_cast<Real>(binom);
            if (j < degree) {
              deriv *= sol2;
              Real nextbez = get_bez_coeff(j + 1);
              deriv += pwrv * (nextbez - prevbez) * static_cast<Real>(binom2);
              prevbez = nextbez;
              binom2 *= (degree - j - 1);
              binom2 /= (j + 1);
            }
            pwrv *= sol;
            binom *= (degree - j);
            binom /= (j + 1);
          }
          deriv *= static_cast<Real>(degree);
        }
        Real trueres_bound = abs(resid) + scfac * tol;
        Real absderiv = abs(deriv);
        Real param_uncertainty = (absderiv == 0.0)? 1.0 :
          trueres_bound / absderiv;

        // Check if sol is within param_uncertainty_width of the unit interval.

        Real solreal = sol.real();
        Real solimag = sol.imag();
        if (solreal < -param_uncertainty || solreal > 1 + param_uncertainty
          || fabs(solimag) > param_uncertainty) continue;
        PatchMath::IntersectRecord result;
        result.paramcoord[0] = solreal;
        result.degen = (fabs(solreal) <= param_uncertainty || 
          fabs(solreal - 1.0) <= param_uncertainty);
#ifdef DEBUGGING
        if (dump_everything) 
          *debug_str << "keeping, degen = " << result.degen << "\n";
#endif
        result.param_uncertainty = param_uncertainty;
        output_stack.push_back(result);
      }
    }
  }

  // A point in R^2, used for representing Bezier parameters. 
  
  struct R2Point {
    Real coord[2];
    R2Point() { }
    R2Point(Real x, Real y) {
      coord[0] = x;
      coord[1] = y;
    }
    Real l2norm() const {
      return sqrt(coord[0] * coord[0] + coord[1] * coord[1]);
    }
  };

  // A point in C^2, used for representing complex Bezier parameters.

  struct C2Point {
    Complex coord[2];
    C2Point() { }
    C2Point(Complex x, Complex y) {
      coord[0] = x;
      coord[1] = y;
    }
    Real l2norm() const {
      return sqrt(std::norm(coord[0]) + std::norm(coord[1]));
    }
  };


  inline void compute_real_givens(Real a, Real b, Real& c, Real& s) {
    if (a == 0.0 && b == 0.0) {
      c = 1.0;
      s = 0.0;
      return;
    }
    Real denom = (fabs(a)>fabs(b))?  
      (a * sqrt(1.0+b*b/(a*a))) : (b * sqrt(1.0 + a*a/(b*b)));
    c = a / denom;
    s = b / denom;
  }


  // Generalized eigenvalue solver.  Solve A*x-lambda*B*x=0. 
  // Use LAPACK QZ algorithm.  Input variables are A and B.
  // (They are overwritten).  Output variable V is the eigenvectors.  
  // The return value
  // is the worst condition number among the eigenvectors whose
  // eigenvalue is between the specified bounds.  It
  // is -1.0 if the QZ iteration didn't converge.

  Real generalized_eigensolver(PatchMath::Workspace::ComplexMatrix& A,
    PatchMath::Workspace::ComplexMatrix& B,
    PatchMath::Workspace::ComplexMatrix& V,
    PatchMath::Workspace::RealVector& eveccond,
    Real eigenval_lb,
    Real eigenval_ub,

    // Workspace:

    PatchMath::Workspace& wksp) {




#ifdef DEBUGGING
    
    using QMG::Local_to_PatchMath_CPP::debug_str;
    using QMG::Local_to_PatchMath_CPP::dump_everything;
#endif

    PatchMath::Workspace::StartPosMarker startmark_(wksp);


    const int blksize = 64;

    Lapack_int n = static_cast<Lapack_int>(B.numrow());
    Lapack_int lwork = static_cast<Lapack_int>(n * (n + blksize));

    PatchMath::Workspace::ComplexVector cwksp(wksp, lwork);
    PatchMath::Workspace::ComplexVector tau(wksp, n);
    PatchMath::Workspace::ComplexVector alpha(wksp, n);
    PatchMath::Workspace::ComplexVector beta(wksp, n);
    PatchMath::Workspace::RealVector rwork(wksp, n);


    Lapack_int info = 89304;

    
    // Factor B = QR.

    zgeqrf_(&n, &n, B.make_fortran_complex(), &n, tau.make_fortran_complex(),
      cwksp.make_fortran_complex(), &lwork, &info);
    if (info != 0)
      throw_error("zgeqrf failure");

    // Apply the HH transforms to A also.

    zunmqr_("L", "C", &n, &n, &n, B.make_fortran_complex(),
      &n, tau.make_fortran_complex(), A.make_fortran_complex(),
      &n, cwksp.make_fortran_complex(), &lwork, &info);

    if (info != 0)
      throw_error("zgunmqr failure");


    // Zero out B below the diagonal.

    {
      for (int j = 0; j < n; ++j)
        for (int i = j + 1; i < n; ++i)
          B(i,j) = Complex(0.0,0.0);
    }


    // Reduce to Hessenberg form.  Save Z.


    Lapack_int one = 1;


    PatchMath::Workspace::ComplexMatrix Z(wksp, n, n);
    zgghrd_("N", "I", &n, &one, &n, A.make_fortran_complex(),
      &n, B.make_fortran_complex(), &n, 0, &one, Z.make_fortran_complex(),
      &n, &info);

    if (info != 0)
      throw_error("zgghrd failure");

    // QZ iteration on A and B.


    zhgeqz_("S", "N", "V", &n, &one, &n, A.make_fortran_complex(),
      &n, B.make_fortran_complex(), &n, alpha.make_fortran_complex(),
      beta.make_fortran_complex(), 0, &one, Z.make_fortran_complex(),
      &n, cwksp.make_fortran_complex(), &lwork, rwork.make_fortran_real(), &info);

    if (info != 0)
      return -fabs(static_cast<Real>(info));

    // Solve for right eigenvectors.  Estimate condition numbers
    // of eigenvectors (nonstandard definition of eigenvector
    // condition!) for eigenvectors whose eigenvalue is between
    // lb and ub.

    Real worstcond = 0.0;
    Real delta = eigenval_ub - eigenval_lb;
    Real ub1 = eigenval_ub + delta / 10.0;
    Real lb1 = eigenval_lb - delta / 10.0;
    Real cbd = delta / 10.0;

    PatchMath::Workspace::ComplexMatrix T(wksp, n, n);
    PatchMath::Workspace::ComplexVector x(wksp, n - 1);
    PatchMath::Workspace::ComplexVector y(wksp, n - 1);
    PatchMath::Workspace::ComplexVector xi(wksp, n - 1);
    PatchMath::Workspace::ComplexVector z(wksp, n - 1);
    PatchMath::Workspace::RealMatrix savegiv(wksp, n, 4);
    PatchMath::Workspace::ComplexVector eigenvec(wksp, n);


    using std::conj;
    using std::polar;
    using std::arg;
    using std::abs;

    for (int eigidx = 0; eigidx < n; ++eigidx) {
      Complex beta1 = B(eigidx, eigidx);
      Complex alpha1 = A(eigidx, eigidx);
      // Set T = B(i,i)*A - A(i,i) * B.
      {
        for (int j = 0; j < n; ++j)
          for (int i = 0; i <= j; ++i)
            T(i,j) = beta1 * A(i,j) - alpha1 * B(i,j);
      }
      // Zero out diagonal entries below the pivot using Givens rotations.

      bool T_singular = false;

      {
        for (int p = eigidx; p < n - 1; ++p) {
          // Scale rows p and p+1 so that leading entries are real.
          
          for (int i1 = p; i1 < p + 2; ++i1) {      
            if (T(i1, p+1) != Complex(0.0,0.0)) {
              Complex sc = polar(1.0, -arg(T(i1,p+1)));
              for (int j = p+1; j < n; ++j)
                T(i1,j) *= sc;
            }
          }
          
          // Compute real Givens rotation.
          Real c,s;
          compute_real_givens(T(p,p+1).real(),T(p+1,p+1).real(), c, s);
          // Apply Givens rotations.
          {
            for (int j = p+1; j < n; ++j) {
              Complex tmp = T(p,j);
              T(p,j) = T(p,j) * c + T(p+1,j) * s;
              T(p+1,j) = T(p+1,j) * c - tmp * s;
            }
          }
          if (T(p,p+1) == Complex(0.0,0.0))
            T_singular = true;
        }
      }
      
      // Zero out diagonal entries above the pivot using Givens rotations.
      // Save the rotations.
      
      {
        for (int p = eigidx - 1; p >= 0; --p) {
          // Scale cols p and p+1 so that leading coeffs are real.
          
          for (int j1 = p; j1 < p + 2; ++j1) {
            savegiv(p,j1-p) = 0.0;
            if (T(p,j1) != Complex(0.0,0.0)) {
              Real a = -arg(T(p,j1));
              savegiv(p,j1-p) = a;
              Complex sc = polar(1.0, a);
              for (int i = 0; i <= p; ++i)
                T(i,j1) *= sc;
            }
          }

          // Compute real Givens rotation.

          Real c,s;
          compute_real_givens(T(p,p+1).real(), T(p,p).real(), c, s);
          savegiv(p,2) = c;
          savegiv(p,3) = s;
          // Apply Givens rotations.
          {
            for (int i = 0; i <= p; ++i) {
              Complex tmp = T(i,p+1);
              T(i,p+1) = T(i,p+1) * c + T(i,p) * s;
              T(i,p) = T(i,p) * c - tmp * s;
            }
          }
          if (T(p,p+1) == Complex(0.0,0.0))
            T_singular = true;
        }
      }

      // Apply the saved Givens rotations to [1;0;...0] to get the eigenvector.
      { 
        for (int i = 1; i < n; ++i)
          eigenvec[i] = Complex(0.0,0.0);
        eigenvec[0] = Complex(1.0,0.0);
      }


      {
        for (int p = 0; p < eigidx; ++p) {
          Real c = savegiv(p,2);
          Real s = savegiv(p,3);
          Complex tmp = eigenvec[p+1];
          eigenvec[p+1] = (eigenvec[p+1] * c - eigenvec[p] * s) * polar(1.0, savegiv(p,1));
          eigenvec[p] = (eigenvec[p] * c + tmp * s) * polar(1.0, savegiv(p,0));
        }
      }

      // Multiply Z times the eigenvector; store the answer as a column of V.

      {
        for (int i = 0; i < n; ++i) {
          Complex t(0.0, 0.0);
          for (int j = 0; j < n; ++j) {
            t += Z(i,j) * eigenvec[j];
          }
          V(i,eigidx) = t;
        }
      }



#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "eigenvec(:," << eigidx+1 << ")= [\n";
        for (int j = 0;  j < n; ++j) {
          *debug_str << Format_Complex(V(j,eigidx)) << '\n';
        }
        *debug_str << "];\n";
        if (eigidx == 0) {
          *debug_str << "reducedT = [\n ";
          for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
              if (j < i)
                *debug_str << "0 ";
              else
                *debug_str << Format_Complex(T(i,j)) << " ";
            }
            *debug_str << '\n';
          }
          *debug_str << "];\n";
        }
      }
      
#endif

      // Compute the condition number of the reduced matrix T to get a
      // condition estimation.  Use Hager's method.

      Real conde = 0.0;

      if (T_singular) {
        conde = BIG_REAL;
      }
      else {
        const int hager_trial = 2;
        const int max_hager_it = 6;
        
        for (int tr = 0; tr < hager_trial; ++tr) {
          // Select a random x with 1-norm equal to 1.
          {
            Real sm = 0.0;
            {
              for (int j = 0; j < n - 1; ++j) {
                x[j] = Complex(random_real() - 0.5, random_real() - 0.5);
                sm += abs(x[j]);
              }
            }
            {
              for (int j = 0; j < n - 1; ++j) {
                x[j] /= sm;
              }
            }
          }
          
          for (int it = 0;; ++it) {
            
            // Backsolve: y = T \ x.
            {
              for (int i = n - 2; i >= 0; --i) {
                Complex t = x[i];
                for (int j = i + 1; j < n - 1; ++j)
                  t -= T(i,j+1) * y[j];
                y[i] = t / T(i, i + 1);
              }
            }
            
            // xi = sgn(y)
            
            {
              for (int i = 0; i < n - 1; ++i)
                if (y[i]== Complex(0.0, 0.0))
                  xi[i] = Complex(1.0,0.0);
                else
                  xi[i] = polar(1.0, arg(y[i]));
            }
            // Forward sub: z = T' \ xi
            {
              for (int i = 0; i < n - 1; ++i) {
                Complex t = xi[i];
                for (int j = 0; j < i; ++j) 
                  t -= conj(T(j, i + 1)) * z[j];
                z[i] = t / conj(T(i, i + 1));
              }
            }
            Real normzinf = -1.0;
            Complex ipxz(0.0, 0.0);
            int maxp;
            {
              for (int i = 0; i < n - 1; ++i) {
                Real e1 = abs(z[i]);
                if (e1 > normzinf) {
                  normzinf = e1;
                  maxp = i;
                }
                ipxz += conj(z[i]) * x[i];
              }
            }
            if (abs(ipxz) >= normzinf || it == max_hager_it) {
              Real normy = 0.0;
              for (int i = 0; i < n - 1; ++i)
                normy += abs(y[i]);
#ifdef DEBUGGING
              if (dump_everything) {
                *debug_str << "normest = " << normy << '\n';
              }
#endif
              if (normy >= conde)
                conde = normy;
              break;
            }
            {
              for (int i = 0; i < n - 1; ++i)
                x[i] = Complex(0.0,0.0);
            }
            x[maxp] = Complex(1.0,0.0);
          }
        }
      }
        
      eveccond[eigidx] = conde;
      if (beta1 != Complex(0.0,0.0)) {
        Complex lambda = alpha1 / beta1;
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << " alpha1 = " << Format_Complex(alpha1) 
            << "\nbeta1 = " << Format_Complex(beta1)
            << "\nlambda = " << Format_Complex(lambda)
            << "\nconde = " << conde
            << "\nbounds = " << lb1 << "," << ub1 << "," << cbd <<'\n';
        }
#endif
        if (lambda.real() >= lb1 && lambda.real() <= ub1 && 
          fabs(lambda.imag()) <= cbd && conde > worstcond) {
          worstcond = conde;
        }
      }
    }
    return worstcond;
  }



  // Class for math with bezier triangles.

  class Solve_BezTri {
  public:
    virtual Real get_bez_coeff(int eqnum, int j) const = 0;
    Solve_BezTri() { }
    virtual ~Solve_BezTri() { }



    bool solve(int degree,
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace,
      double scfac, 
      double tol,
      bool newton_improve) const;
    Real eval(const R2Point& param,
      int eqnum,
      int degree) const;
    Real eval_direc_deriv(const R2Point& param,
      const R2Point& direc,
      int eqnum,
      int degree) const;
    Complex eval_complex(const C2Point& param,
      int eqnum,
      int degree) const;
    Complex eval_direc_deriv_complex(const C2Point& param,
      const C2Point& direc,
      int eqnum,
      int degree) const;
  };

  // Derived class for intersecting a line with a triangle.

  class Solve_BezTri1 : public Solve_BezTri {
  private: 
    const Matrix& lineeq_;
    const Point& linerhs_;
    const PatchMath& pm_;
  public:
    Solve_BezTri1(const Matrix& lineeq,
      const Point& linerhs,
      const PatchMath& pm) : lineeq_(lineeq), linerhs_(linerhs), pm_(pm) { }
    Real get_bez_coeff(int eqnum, int cpnum) const {
      return lineeq_(eqnum,0) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        lineeq_(eqnum,1) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 1) +
        lineeq_(eqnum,2) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 2) -
        linerhs_[eqnum];
    }
  };

  // Derived class for evaluating directional deriv.

  class Solve_BezTri2 : public Solve_BezTri {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
  public:
    Solve_BezTri2(const Point& direc, 
      int di,
      const PatchMath& pm) : direc_(direc), di_(di), pm_(pm) 
    { }
    Real get_bez_coeff(int eqnum, int cpnum) const;
  };
  
  

  // Derived class for evaluating one coordinate of a bezier triangle.

  class Solve_BezTri3 : public Solve_BezTri {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_BezTri3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int, int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };
  
  
  Real Solve_BezTri2::get_bez_coeff(int eqnum, int cpnum) const {
    Real t = 0.0;
    int offset = static_cast<int>(sqrt(2 * static_cast<double>(cpnum)));
    if (offset * (offset + 1) / 2 > cpnum)
      --offset;
    int bcp = cpnum + offset + 1;
    int ecp = (eqnum == 0)? bcp + 1 : cpnum;
    for (int j = 0; j < di_; ++j) {
      t += (PatchMath::Bezier_Helper::control_point_coord(pm_, ecp, j) -
        PatchMath::Bezier_Helper::control_point_coord(pm_, bcp, j)) *
        direc_[j];
    }
    return t;
  }
    
  // Compute real coord from param coord using Horner's algorithm.

  Real Solve_BezTri::eval(const R2Point& param,
    int eqnum,
    int degree) const {
    Real p0 = param.coord[0];
    Real p1 = param.coord[1];
    Real p2 = 1.0 - p0 - p1;
    int outer_binom = 1;
    Real pwrv = 1.0;
    Real result = 0.0;
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Real pwru = 1.0;
      Real row_result = 0.0;
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += inner_binom * 
          get_bez_coeff(eqnum, base_idx + col) * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * outer_binom;
      outer_binom *= row;
      outer_binom /= (degree - row + 1);
      pwrv *= p1;
    }
    return result;
  }

  // Compute complex coord from param coord using Horner's algorithm.

  Complex Solve_BezTri::eval_complex(const C2Point& param,
    int eqnum,
    int degree) const {
    Complex p0 = param.coord[0];
    Complex p1 = param.coord[1];
    Complex p2 = Complex(1.0,0.0) - p0 - p1;
    int outer_binom = 1;
    Complex pwrv(1.0,0.0);
    Complex result(0.0,0.0);
    for (int row = degree; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Complex pwru(1.0,0.0);
      Complex row_result(0.0,0.0);
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += Complex(static_cast<Real>(inner_binom), 0.0) * 
          Complex(get_bez_coeff(eqnum, base_idx + col),0.0) * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * Complex(static_cast<Real>(outer_binom), 0);
      outer_binom *= row;
      outer_binom /= (degree - row + 1);
      pwrv *= p1;
    }
    return result;
  }
    
  // Compute directional derivative using Horner's rule.
  
  Real Solve_BezTri::eval_direc_deriv(const R2Point& param,
    const R2Point& direc,
    int eqnum,
    int degree) const {   
    Real p0 = param.coord[0];
    Real p1 = param.coord[1];
    Real p2 = 1.0 - p0 - p1;
    Real d0 = direc.coord[0];
    Real d1 = direc.coord[1];
    Real d2 = -d0 - d1;    
    int outer_binom = 1;
    Real pwrv = 1.0;
    Real result = 0.0;
    for (int row = degree - 1; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Real pwru = 1.0;
      Real row_result = 0.0;
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += inner_binom * 
          (d1 * get_bez_coeff(eqnum, base_idx + col) +
          d2 * get_bez_coeff(eqnum, base_idx + col + row + 1) +
          d0 * get_bez_coeff(eqnum, base_idx + col + row + 2))
          * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * outer_binom;
      outer_binom *= row;
      outer_binom /= (degree - row);
      pwrv *= p1;
    }
    return degree * result;
  }

  // Compute directional derivative using Horner's rule.
  
  Complex Solve_BezTri::eval_direc_deriv_complex(const C2Point& param,
    const C2Point& direc,
    int eqnum,
    int degree) const {   
    Complex p0 = param.coord[0];
    Complex p1 = param.coord[1];
    Complex p2 = Complex(1.0,0.0) - p0 - p1;
    Complex d0 = direc.coord[0];
    Complex d1 = direc.coord[1];
    Complex d2 = -d0 - d1;    
    int outer_binom = 1;
    Complex pwrv(1.0,0.0);
    Complex result(0.0,0.0);
    for (int row = degree - 1; row >= 0; --row) {
      int base_idx = row * (row + 1) / 2;
      int inner_binom = 1;
      Complex pwru(1.0,0.0);
      Complex row_result(0.0,0.0);
      for (int col = 0; col <= row; ++col) {
        row_result *= p2;
        row_result += Complex(static_cast<Real>(inner_binom), 0.0) * 
          (d1 * get_bez_coeff(eqnum, base_idx + col) +
          d2 * get_bez_coeff(eqnum, base_idx + col + row + 1) +
          d0 * get_bez_coeff(eqnum, base_idx + col + row + 2))
          * pwru;
        pwru *= p0;
        inner_binom *= (row - col);
        inner_binom /= (col + 1);
      }
      result += row_result * pwrv * Complex(static_cast<Real>(outer_binom), 0.0);
      outer_binom *= row;
      outer_binom /= (degree - row);
      pwrv *= p1;
    }
    return result * static_cast<Real>(degree);
  }

  // Function to convert a monomial degree pair into a column number.

  inline int deg_pair_to_colnum(int u_deg, int v_deg, int total_deg) {
#ifdef RANGECHECK
    if (u_deg < 0 || v_deg < 0 || u_deg + v_deg > total_deg)
      throw_error("degrees out of range");
#endif
    return (total_deg - v_deg) * (total_deg - v_deg + 1) / 2 + u_deg;
  }

  // Inverse function -- convert a column number to a monomial degree pair.

  inline pair<int,int> colnum_to_deg_pair(int colnum, int total_deg) {
    int k = static_cast<int>(sqrt(2 * static_cast<double>(colnum)));
    if (k * (k + 1) / 2 > colnum)
      --k;
    return pair<int,int>(colnum -  k * (k + 1) / 2, total_deg - k);
  }

  
  
  Real linear_solve_C2(PatchMath::Workspace::ComplexMatrix& A,
    C2Point& rhs,
    C2Point& sol) {
    using std::abs;
    using std::norm;
    Complex det = A(0,0) * A(1,1) - A(1,0) * A(0,1);
    if (abs(det) == 0.0) {
      return BIG_REAL;
    }
    sol.coord[1] = (A(0,0) * rhs.coord[1] - A(1,0) * rhs.coord[0]) / det;
    if (abs(A(0,0)) >= abs(A(1,0))) {
      sol.coord[0] = (rhs.coord[0] - A(0,1) * sol.coord[1]) / A(0,0);
    }
    else {
      sol.coord[0] = (rhs.coord[1] - A(1,1) * sol.coord[1]) / A(1,0);
    }
    return sqrt(norm(A(0,0)) + norm(A(0,1)) + norm(A(1,0))
      + norm(A(1,1))) / abs(det);
  }


  // Solve a system of two polynomial equations
  // in Bezier triangle form.  Use Macaulay-type
  // result and random u-direction.  Reduce
  // to small system and invoke LAPACK for
  // generalized eigenvalues.
  // Returns true if there was a global degeneracy
  // during the computation.

  
  bool Solve_BezTri::solve(int degree,
    vector<PatchMath::IntersectRecord>& output_stack,
    PatchMath::Workspace& workspace, 
    double scfac, 
    double tol,
    bool newton_improve) const {




#ifdef DEBUGGING
    static int invocation_count = 0;
    ++invocation_count;
    using QMG::Local_to_PatchMath_CPP::debug_str;
    using QMG::Local_to_PatchMath_CPP::dump_everything;
    using std::setprecision;
    using std::setiosflags;
    using std::ios;


    // if (invocation_count > 962)
    //  dump_everything = true;
    // dump_everything = (invocation_count == 148);


    if (dump_everything) {
      *debug_str << " In solve beztri degree = " << degree 
        << " scfac = " << scfac << " tol = " << tol << "\n bez  = [\n";
      int ncp = (degree + 1) * (degree + 2) / 2; 
      for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < ncp; ++j) {
          *debug_str << std::setprecision(17) << get_bez_coeff(i,j) << " ";
        }
        *debug_str << '\n';
      }
      *debug_str << "];\n";
      *debug_str << "invocation_count = " << invocation_count << '\n';
    }
#endif

    bool is_degen = false;

    if (degree == 1) {
      Real x[3];
      Real rhs[2];
      Real matrixspace[4];
      Matrix m(2,2,matrixspace);
      for (int eqnum = 0; eqnum < 2; ++eqnum) {
        m(eqnum,0) = get_bez_coeff(eqnum,2) - get_bez_coeff(eqnum,1);
        m(eqnum,1) = get_bez_coeff(eqnum,0) - get_bez_coeff(eqnum,1);
        rhs[eqnum] = -get_bez_coeff(eqnum,1);
      }
      Real invnormest = m.linear_solve(rhs,x);
      if (invnormest >= BIG_REAL) {
        is_degen = true;
        return is_degen;
      }

      x[2] = 1.0 - x[0] - x[1];
      Real condest = scfac * invnormest;
      if (condest > sqrt(tol))
        is_degen = true;
      Real unc = condest * tol;
      if (x[0] < -unc || x[1] < -unc || x[2] < -unc)
        return is_degen;
      PatchMath::IntersectRecord result;
      result.paramcoord[0] = x[0];
      result.paramcoord[1] = x[1];
      result.param_uncertainty = unc;
      result.degen = fabs(x[0]) <= unc || fabs(x[1]) <= unc || fabs(x[2]) <= unc;
      output_stack.push_back(result);
      return is_degen;
    }
    else {    
      using std::abs;
      using std::norm;

      PatchMath::Workspace::StartPosMarker start_pos_marker_(workspace);
      
      // Check if (0,0) is in the bounding box; if not then early return.
      
      int num_control_pt = (degree + 1) * (degree + 2) / 2;

      {
        for (int eqnum = 0; eqnum < 2; ++eqnum) {
          bool pos_seen = false;
          bool neg_seen = false;
          for (int i = 0; i < num_control_pt; ++i) {
            Real bz = get_bez_coeff(eqnum, i);
            if (bz <= scfac * tol)
              neg_seen = true;
            if (bz >= -scfac * tol)
              pos_seen = true;
          }
          if (!neg_seen || !pos_seen) {
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << " early return; equation " << eqnum << " misses 0.\n";
            }
#endif
            return false;
          }
        }
      }
      
      int tmrow = degree * (degree + 1);
      int tmcol = degree * (2 * degree + 1);
      int evsize = degree * degree; // = tmcol-tmrow
      int monomcount = (2 * degree - 1) * degree;
      int num_row_monom = degree * (degree + 1) / 2;

      Real bezmax_poly[2];
      bezmax_poly[0] = 0.0;
      bezmax_poly[1] = 0.0;

      PatchMath::Workspace::RealMatrix coef_perturb(workspace, 2, num_control_pt);

      {
        for (int i = 0; i < 2; ++i) {
          for (int j = 0; j < num_control_pt; ++j) {
            coef_perturb(i,j) = 0.0;
            Real c1 = fabs(get_bez_coeff(i,j));
            if (c1 > bezmax_poly[i])
              bezmax_poly[i] = c1;
          }
        }
      }
      if (bezmax_poly[0] == 0 || bezmax_poly[1] == 0) {
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << " early return; one equation is all zeros\n";
        }
#endif
        is_degen = true;
        return is_degen;
      }

      // Find a shift so that none of the leading coefficients
      // vanish.  This shift is an orthogonal matrix applied to
      // u,v,w.

      // Make up a random 3-by-3 matrix, and QR factor it.

      PatchMath::Workspace::RealMatrix test_shiftQ(workspace, 3, 3);
      PatchMath::Workspace::RealMatrix best_shiftQ(workspace, 3, 3);
      PatchMath::Workspace::RealMatrix shift0(workspace, 3, 3);
      PatchMath::Workspace::RealMatrix sh_hhmat(workspace, 3, 3);
      PatchMath::Workspace::RealVector sh_betavec(workspace, 3);
      PatchMath::Workspace::RealVector pwrv(workspace, num_control_pt);
      PatchMath::Workspace::RealVector pwruv(workspace, num_control_pt);      
      PatchMath::Workspace::RealVector temp_poly(workspace, num_control_pt);    
      PatchMath::Workspace::RealMatrix test_shiftpoly(workspace, 2, num_control_pt);   
      PatchMath::Workspace::RealMatrix best_shiftpoly(workspace, 2, num_control_pt);
      PatchMath::Workspace::RealMatrix shiftpoly_thisrow(workspace, 2, num_control_pt);

      // Try four shifts; keep the one with the best leading coefficients.

      Real best_lead = -1.0;
      
      for (int shift_trycount = 0; shift_trycount < 4; ++shift_trycount) {

        for (int qevalid = 0;;++qevalid) {        
          // Make a matrix shift0 that is a perturbation of the identity.
          {
            for (int i = 0; i < 3; ++i) {
              for (int j = 0; j < 3; ++j) {
                shift0(i,j) = random_real();
              }
              shift0(i,i) += 1.0;
            } 
          }
          // Qr factor it
          QR_Factor(shift0, sh_hhmat, sh_betavec);
          formQ(sh_hhmat, sh_betavec, test_shiftQ);
          // Switch signs of rows so that diagonal entries are positive.
          {
            for (int i = 0; i < 3; ++i) {
              if (test_shiftQ(i,i) < 0) {             
                for (int j = 0; j < 3; ++j) {
                  test_shiftQ(i,j) *= -1.0;
                }
              }
            }
          }
          // Q must have the property that every entry of Q*[1;1;1] is
          // at least 0.1.   If this property
          // doesn't hold, try another time.
          bool Qe_valid = true;
          {
            for (int i = 0; i < 3; ++i) {
              Real t = 0.0;
              for (int j = 0; j < 3; ++j) {
                t += test_shiftQ(i,j);
              }
              if (t < 0.1)
                Qe_valid = false;
            }
          }
          if (Qe_valid) break;
        }


 #ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "test_shiftQ = [\n";
          for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
              *debug_str << " " << std::setprecision(17) << test_shiftQ(i,j);
            }
            *debug_str << '\n';
          }
          *debug_str << "];\n";
        }
#endif
        
        // Now apply the shift.  Use horner's rule 'symbolically';
        // Substitute Q(0,:)(u',v',w') for u, Q(1,:)(u',v',w') for v,
        // and Q(2,:)(u',v',w') for w.
        
        int outer_binom = 1;
        pwrv[0] = 1.0;
        {
          for (int i = 0; i < num_control_pt; ++i)
            for (int eqnum = 0; eqnum < 2; ++eqnum)
              test_shiftpoly(eqnum, i) = 0.0;
        }
        {
          for (int vdeg = 0; vdeg <= degree; ++vdeg) {
            int base_idx = (degree - vdeg) * (degree - vdeg + 1) / 2;
            int inner_binom = 1;
            {
              int vdeg_numcoef = (vdeg + 1) * (vdeg + 2) / 2;
              for (int i = 0; i < vdeg_numcoef; ++i)
                pwruv[i] = pwrv[i];
            }
            {
              for (int i = 0; i < num_control_pt; ++i)
                for (int eqnum = 0; eqnum < 2; ++eqnum)
                  shiftpoly_thisrow(eqnum, i) = 0.0;
            }
            for (int udeg = 0; udeg <= degree - vdeg; ++udeg) {
              int shiftpolydeg = udeg + vdeg;
              int wdeg = degree - shiftpolydeg;
              
              // temp_poly = shiftpoly_thisrow * w symbolically

              
              int shiftpoly_numcoef = 
                (shiftpolydeg + 1) * (shiftpolydeg + 2) / 2;
              for (int eqnum = 0; eqnum < 2; ++eqnum) {
                {
                  for (int j = 0; j < shiftpoly_numcoef; ++j) {
                    pair<int,int> rval = colnum_to_deg_pair(j, shiftpolydeg);
                    int ud1 = rval.first;
                    int vd1 = rval.second;
                    int wd1 = shiftpolydeg - ud1 - vd1;
                    Real t = 0.0;
                    if (ud1 > 0) {
                      int subp = deg_pair_to_colnum(ud1 - 1, vd1, shiftpolydeg - 1);
                      t += shiftpoly_thisrow(eqnum,subp) * test_shiftQ(2,0);
                    }
                    if (vd1 > 0) {
                      int subp = deg_pair_to_colnum(ud1, vd1 - 1, shiftpolydeg - 1);
                      t += shiftpoly_thisrow(eqnum,subp) * test_shiftQ(2,1);
                    }
                    if (wd1 > 0) {
                      int subp = deg_pair_to_colnum(ud1, vd1, shiftpolydeg - 1);
                      t += shiftpoly_thisrow(eqnum,subp) * test_shiftQ(2,2);
                    }
                    temp_poly[j] = t;
                  }
                }
                
                
                // Symbolic computation of 
                //  shiftpoly_thisrow = temp_poly + inner_binom * 
                //                get_bez_coef(eqnum, base_idx + udeg) * pwruv;
                
                Real g = get_bez_coeff(eqnum, base_idx + udeg) 
                  * inner_binom / bezmax_poly[eqnum];
                {
                  for (int j = 0; j < shiftpoly_numcoef; ++j) {
                    shiftpoly_thisrow(eqnum, j) = temp_poly[j] + g * pwruv[j];
                  }
                }
              }
              
              inner_binom *= (degree - vdeg - udeg);
              inner_binom /= (udeg + 1);

              // Symbolic computation of pwruv *= u;
              
              if (shiftpolydeg < degree) {
                int shiftpolydeg1 = shiftpolydeg + 1;
                int new_numcoef = (shiftpolydeg1 + 1) * (shiftpolydeg1 + 2) / 2;
                for (int j = 0; j < new_numcoef; ++j) {
                  pair<int,int> rval = colnum_to_deg_pair(j, shiftpolydeg1);
                  int ud1 = rval.first;
                  int vd1 = rval.second;
                  int wd1 = shiftpolydeg1 - ud1 - vd1;
                  Real t = 0.0;
                  if (ud1 > 0) {
                    int subp = deg_pair_to_colnum(ud1 - 1, vd1, shiftpolydeg);
                    t += pwruv[subp] * test_shiftQ(0,0);
                  }
                  if (vd1 > 0) {
                    int subp = deg_pair_to_colnum(ud1, vd1 - 1, shiftpolydeg);
                    t += pwruv[subp] * test_shiftQ(0,1);
                  }
                  if (wd1 > 0) {
                    int subp = deg_pair_to_colnum(ud1, vd1, shiftpolydeg);
                    t += pwruv[subp] * test_shiftQ(0,2);
                  }
                  temp_poly[j] = t;
                }
                // Transfer result back to pwruv.
                
                for (int j2 = 0; j2 < new_numcoef; ++j2) 
                  pwruv[j2] = temp_poly[j2];
              }
            }
            // Symbolic implementation of
            // shiftpoly += shiftpoly_thisrow * outer_binom

            
            {
              for (int j = 0; j < num_control_pt; ++j)
                for (int eqnum = 0; eqnum < 2; ++eqnum)
                  test_shiftpoly(eqnum,j) += shiftpoly_thisrow(eqnum,j) * outer_binom;
            }
          
            outer_binom *= (degree - vdeg);
            outer_binom /= (vdeg + 1);
            // Symbolic implementation of pwrv *= v;              
            
            if (vdeg < degree) {
              int vdeg1 = vdeg + 1;
              int new_numcoef = (vdeg1 + 1) * (vdeg1 + 2) / 2;
              for (int j = 0; j < new_numcoef; ++j) {
                pair<int,int> rval = colnum_to_deg_pair(j, vdeg1);
                int ud1 = rval.first;
                int vd1 = rval.second;
                int wd1 = vdeg1 - ud1 - vd1;
                Real t = 0.0;
                if (ud1 > 0) {
                  int subp = deg_pair_to_colnum(ud1 - 1, vd1, vdeg);
                  t += pwrv[subp] * test_shiftQ(1,0);
                }
                if (vd1 > 0) {
                  int subp = deg_pair_to_colnum(ud1, vd1 - 1, vdeg);
                  t += pwrv[subp] * test_shiftQ(1,1);
                }
                if (wd1 > 0) {
                  int subp = deg_pair_to_colnum(ud1, vd1, vdeg);
                  t += pwrv[subp] * test_shiftQ(1,2);
                }
                temp_poly[j] = t;
              }
              // Transfer result back to pwrv.
              
              for (int j2 = 0; j2 < new_numcoef; ++j2) 
                pwrv[j2] = temp_poly[j2];
            }
          }
        }

        
        // Now check if the leading coefs are suff large.
        
        bool size_ok = true;
        Real test_lead = BIG_REAL;
        for (int eqnum = 0; eqnum < 2; ++eqnum) {
          for (int p = 0; p < 3; ++p) {
            int ud = (p == 0)? degree : 0;
            int vd = (p == 1)? degree : 0;
            int co1 = deg_pair_to_colnum(ud, vd, degree);
            Real ldgco = fabs(test_shiftpoly(eqnum,co1));
            if (ldgco < test_lead) {
              test_lead = ldgco;
            }
          }
        }


#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << " test_shiftpoly = [\n";
          for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < num_control_pt; ++j)
              *debug_str << " " << std::setprecision(17)
                         << test_shiftpoly(i,j);
            *debug_str << '\n';
          }
          *debug_str << "]\ntest_lead = " << test_lead << "\n";
        }
#endif

        if (shift_trycount == 0 || test_lead > best_lead) {
#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << " keeping this one\n";
          }
#endif
          for (int eqnum2 = 0; eqnum2 < 2; ++eqnum2) {
            for (int j = 0; j < num_control_pt; ++j) {
              best_shiftpoly(eqnum2,j) = test_shiftpoly(eqnum2, j);
            }
          }
          for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
              best_shiftQ(i,j) = test_shiftQ(i,j);
            }
          }
        }
      }
      
      // Get here after shift is done.
      // Now we must factor the top part of the matrix.  If
      // factorization detects a singularity, add a perturbation
      // to the coefficients.
  
      // Compute the top part of the matrix,
      // which contains shifted copies of the bezier coefficients.

      // Note, in this code, the top part is transposed, so it's
      // actually the 'left part'.

      PatchMath::Workspace::RealVector row_u_deg(workspace,tmrow);
      PatchMath::Workspace::RealVector row_v_deg(workspace,tmrow);
      PatchMath::Workspace::RealVector rowtype(workspace,tmrow);
      PatchMath::Workspace::RealMatrix Tm(workspace, tmcol, tmrow);
      PatchMath::Workspace::RealMatrix Tmp(workspace, tmrow, tmcol);
      PatchMath::Workspace::RealMatrix hhmat_tmp(workspace, tmrow, tmrow);
      PatchMath::Workspace::RealVector betavec_tmp(workspace, tmrow);
      PatchMath::Workspace::RealVector select_col_tmp(workspace, tmrow);
      PatchMath::Workspace::RealVector fac(workspace, tmrow);
      PatchMath::Workspace::RealMatrix depmat(workspace, tmrow, tmrow);
      
      bool coefs_were_perturbed = false;
      Real expected_resd = scfac * tol;
      
      for (int top_part_trycount = 0;; ++top_part_trycount) {

        for (int row_monom_idx = 0; row_monom_idx < num_row_monom; ++row_monom_idx) {
          pair<int,int> row_idx = colnum_to_deg_pair(row_monom_idx, degree - 1);
          int u_row_deg = row_idx.first;
          int v_row_deg = row_idx.second;
          for (int eqnum = 0; eqnum < 2; ++eqnum) {
            int mtxrow = row_monom_idx * 2 + eqnum;
            row_u_deg[mtxrow] = static_cast<double>(u_row_deg);
            row_v_deg[mtxrow] = static_cast<double>(v_row_deg);
            rowtype[mtxrow] = static_cast<double>(eqnum);
            for (int j = 0; j < tmcol; ++j)
              Tm(j, mtxrow) = 0.0;
            for (int poly_monom_idx = 0; poly_monom_idx < num_control_pt; ++poly_monom_idx) {
              pair<int,int> poly_idx = colnum_to_deg_pair(poly_monom_idx, degree);
              int u_poly_deg = poly_idx.first;
              int v_poly_deg = poly_idx.second;
              int u_col_deg = u_poly_deg + u_row_deg;
              int v_col_deg = v_poly_deg + v_row_deg;
              int mtxcol = deg_pair_to_colnum(u_col_deg,v_col_deg, 2*degree - 1);
              Tm(mtxcol, mtxrow) = 
                best_shiftpoly(eqnum, poly_monom_idx) + coef_perturb(eqnum, poly_monom_idx);
            }
          }
        }
        
        if (top_part_trycount == 1) {
          break;
        }
        {
          for (int i = 0; i < tmcol; ++i) {
            for (int j = 0; j < tmrow; ++j) {
              Tmp(j,i) = Tm(i,j);
            }
          }
        }
        
        QR_Factor_with_col_pivoting(Tmp, select_col_tmp, hhmat_tmp, betavec_tmp);
        
        
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "after QR factorization with col pivoting, factors are\nTmp = [\n";
          {
            for (int i = 0; i < tmrow; ++i) {
              for (int j = 0; j < tmcol; ++j) {
                *debug_str << Tmp(i,j) << " ";
              }
              *debug_str << '\n';
            }
          }
          *debug_str << "];\nselect_col_tmp = [";
          {
            for (int i = 0; i < tmrow; ++i)
              *debug_str << static_cast<int>(select_col_tmp[i]) << " ";
          }
          *debug_str << "];\nhhmat_tmp = [";
          {
            for (int i = 0; i < tmrow; ++i) {
              for (int j = 0; j < i + 1; ++j) {
                *debug_str << hhmat_tmp(i,j) << " ";
              }
              for (int j2 = i + 1; j2 < tmrow; ++j2) {
                *debug_str << "0 ";
              }
              *debug_str << '\n';
            }
          }
          *debug_str << "];\nbetavec_tmp = [";
          {
            for (int i = 0; i < tmrow - 1; ++i)
              *debug_str << betavec_tmp[i] << " ";
          }
          *debug_str << "];\n";
        }
#endif
        
        int smallsvpos = tmrow - 1;
        Real bigsv = fabs(Tmp(0,static_cast<int>(select_col_tmp[0])));
        
        while (fabs(Tmp(smallsvpos,static_cast<int>(select_col_tmp[smallsvpos])))
          < bigsv * pow(MACHINE_EPS, 0.33333)) {
          --smallsvpos;
        }
        
        if (smallsvpos == tmrow - 1) {
          break;
        }
        
        
        // Make up approximate null vectors.  Initialize columns
        // of the identity matrix.
        
        {
          for (int j = smallsvpos + 1; j < tmrow; ++j) {
            for (int i = 0; i < tmrow; ++i) {
              depmat(i,j) = 0.0;
            }
            depmat(j,j) = 1.0;
          }
        }
        
        // Apply HH transforms.
        
        {
          for (int k = smallsvpos + 1; k < tmrow; ++k) {
            for (int j = tmrow - 2; j >= 0; --j) {
              Real ip = 0.0;
              {
                for (int i = j; i < tmrow; ++i)
                  ip += hhmat_tmp(i,j) * depmat(i,k);
              }
              ip *= betavec_tmp[j];
              {
                for (int i = j; i < tmrow; ++i)
                  depmat(i,k) -= ip * hhmat_tmp(i,j);
              }
            }
          }
        }
        
        // Compute the multipliers and update expected_resd.
        
        {
          for (int i = smallsvpos + 1; i < tmrow; ++i) {
            Real sv = fabs(Tmp(i,static_cast<int>(select_col_tmp[i])));
            if (sv < MACHINE_EPS * bigsv)
              sv = MACHINE_EPS * bigsv;
            Real scaledsv = sv / bigsv;
            fac[i] = bigsv * sqrt(MACHINE_EPS  / scaledsv) * (1.0 + random_real());
            Real er = scfac * sqrt(MACHINE_EPS * scaledsv);
            if (er > expected_resd)
              expected_resd = er;
          }
        }
        
        coefs_were_perturbed = true;
        
        {
          for (int i = 0; i < tmrow; ++i) {
            int eqnum = static_cast<int>(1 - rowtype[i]);
            int ud = static_cast<int>(row_u_deg[i]);
            int vd = static_cast<int>(row_v_deg[i]);
            int cpnum = deg_pair_to_colnum(ud, vd, degree);
            for (int j = smallsvpos + 1; j < tmrow; ++j) 
              coef_perturb(eqnum, cpnum) += depmat(i,j) * fac[j];
          }
        }
        
        
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "bigsv = " << bigsv << ";\ndepmat = [\n";
          {
            for (int i = 0; i < tmrow; ++i) {
              for (int j = smallsvpos + 1; j < tmrow; ++j) {
                *debug_str << depmat(i,j) << ' ';
              }
              *debug_str << '\n';
            }
          }
          *debug_str << "];\nsmallsv = [";
          {
            for (int j =  smallsvpos + 1; j < tmrow; ++j) {
              *debug_str << fabs(Tmp(j,static_cast<int>(select_col_tmp[j]))) << ' ';
            }
          }
          *debug_str << "];\nfac = [";
          {
            for (int j = smallsvpos + 1; j < tmrow; ++j) {
              *debug_str << fac[j] << " ";
            }
          }
          *debug_str << "];\ncoef_perturb = [\n";
          
          {
            for (int eqnum1 = 0; eqnum1 < 2; ++eqnum1) {
              for (int i = 0; i < num_control_pt; ++i) {
                *debug_str << coef_perturb(eqnum1,i) << ' ';
              }
              *debug_str << '\n';
            }
          }
          *debug_str << "];\n";
        }
#endif
      }
      
      
      PatchMath::Workspace::RealVector betavec(workspace, tmrow);
      PatchMath::Workspace::RealMatrix hhmat(workspace, tmcol, tmrow);
      
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "Before factoring, Tmi = [\n";
        for (int i = 0; i < tmcol; ++i) {
          for (int j = 0; j < tmrow; ++j) {
            *debug_str << Tm(i,j) << ' ';
          }
          *debug_str << '\n';
        }
        *debug_str << "];\n";
      }
#endif     
      
      QR_Factor(Tm, hhmat, betavec);
      
      
#ifdef DEBUGGING    
      if (dump_everything) {
        {
          *debug_str << "Tmi = [\n";
          for (int i = 0; i < tmcol; ++i) {
            for (int j = 0; j < tmrow; ++j) { 
              *debug_str << Tm(i,j) << " ";
            }
            *debug_str << '\n';
          }
        }
        *debug_str << "];\nhhmati = [\n";
        {
          for (int i = 0; i < tmcol; ++i) {
            int k = (i + 1 < tmrow)? (i + 1) : tmrow;
            for (int j = 0; j < k; ++j) {
              *debug_str << setprecision(17) << hhmat(i,j) << " ";
            }
            for (int j1 = k; j1 < tmrow; ++j1) {
              *debug_str << "0 ";
            }
            *debug_str << '\n';
          }
        }
        *debug_str << "];\n betavec = [\n";
        {
          for (int ii = 0; ii < tmrow; ++ii) {
            *debug_str << setprecision(17) << betavec[ii] << " ";
          }
        }
        *debug_str << "];\n";
      }
#endif
      
 
      // Figure out which rows from the bottom part should be dropped.
      // Make Tmsmall, which is like Tm except one lower degree monomials multiply the rows.
      // and not transposed.

      int tmsmall_row = degree * (degree - 1);
      int tmsmall_col = monomcount;

      PatchMath::Workspace::RealMatrix Tmsmall(workspace, tmsmall_row, tmsmall_col);
      {
        int num_row_monom = degree * (degree - 1) / 2;
        for (int row_monom_idx = 0; row_monom_idx < num_row_monom; ++row_monom_idx) {
          pair<int,int> row_idx = colnum_to_deg_pair(row_monom_idx, degree - 2);
          int u_row_deg = row_idx.first;
          int v_row_deg = row_idx.second;
          for (int eqnum = 0; eqnum < 2; ++eqnum) {
            int mtxrow = row_monom_idx * 2 + eqnum;
            for (int j = 0; j < tmsmall_col; ++j)
              Tmsmall(mtxrow,j) = 0.0;
            for (int poly_monom_idx = 0; poly_monom_idx < num_control_pt; ++poly_monom_idx) {
              pair<int,int> poly_idx = colnum_to_deg_pair(poly_monom_idx, degree);
              int u_poly_deg = poly_idx.first;
              int v_poly_deg = poly_idx.second;
              int u_col_deg = u_poly_deg + u_row_deg;
              int v_col_deg = v_poly_deg + v_row_deg;
              int mtxcol = deg_pair_to_colnum(u_col_deg,v_col_deg, 2*degree - 2);
              Tmsmall(mtxrow,mtxcol) = 
                best_shiftpoly(eqnum, poly_monom_idx) + coef_perturb(eqnum, poly_monom_idx);
            }
          }
        }
      }

#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "Tmsmall = [\n";
        for (int i = 0; i < tmsmall_row; ++i) {
          for (int j = 0; j < tmsmall_col; ++j)
            *debug_str << " " << std::setprecision(17) << Tmsmall(i,j);
          *debug_str << "\n";
        }
        *debug_str << "]\n";
      }
#endif
      
      //  save select_col.

      PatchMath::Workspace::RealMatrix hhmat_small(workspace, tmsmall_row, tmsmall_row);
      PatchMath::Workspace::RealVector betavec_small(workspace, tmsmall_row);
      PatchMath::Workspace::RealVector select_col(workspace, tmsmall_row);

      QR_Factor_with_col_pivoting(Tmsmall, select_col, hhmat_small, betavec_small);

#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "selectcol = [";
        for (int i = 0; i < tmsmall_row; ++i) 
          *debug_str << static_cast<int>(select_col[i]) << " ";
        *debug_str << "]\n";
      }
#endif



      PatchMath::Workspace::RealVector use_monom(workspace, monomcount);      
      {
        for (int j = 0; j < monomcount; ++j)
          use_monom[j] = 1.0;
        for (int j2 = 0; j2 < tmsmall_row; ++j2)
          use_monom[static_cast<int>(select_col[j2])] = 0.0;
      }
      
      
      


      PatchMath::Workspace::RealVector expaArow(workspace, tmcol);
      PatchMath::Workspace::RealVector expaBrow(workspace, tmcol);
      PatchMath::Workspace::ComplexMatrix A(workspace, evsize, evsize);
      PatchMath::Workspace::ComplexMatrix Ao(workspace, evsize, evsize);
      PatchMath::Workspace::ComplexMatrix B(workspace, evsize, evsize);
      PatchMath::Workspace::ComplexMatrix Bo(workspace, evsize, evsize);
      PatchMath::Workspace::ComplexMatrix vr(workspace, evsize, evsize);
      PatchMath::Workspace::ComplexMatrix bestvr(workspace, evsize, evsize);
      PatchMath::Workspace::RealVector eveccond(workspace, evsize);
      PatchMath::Workspace::RealVector best_eveccond(workspace, evsize);
      
      const Real max_allowable_cond = 1.0e17;
      
      
      Real allowable_cond = 1.0e4;
      int lapack_failcount = 0;
      Real best_cond = -1.0;
      
      for (;;) {
        
        
        // Now compute the eigenval matrix.
        // Choose some random numbers.

        Real direc[3];
        Real angle = random_real() * 2 * PI;
        direc[0] = cos(angle);
        direc[1] = sin(angle);
        direc[2] = -cos(angle) - sin(angle); 
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "direc = [";
          for (int i = 0; i < 3; ++i) {
            *debug_str << setprecision(17) << direc[i] << " ";
          }
          *debug_str << "];\n";
        }
#endif
        // Compute upper and lower bounds on eigenvalues in range.

        Real eigenval_lb = BIG_REAL;
        Real eigenval_ub = -BIG_REAL;

        {
          Real extr[3];
          for (int k = 0; k < 3; ++k) {
            extr[0] = (k == 0)? 1.0 : 0.0;
            extr[1] = (k == 1)? 1.0 : 0.0;
            extr[2] = 1.0 - extr[0] - extr[1];
            Real numer = 0.0;
            Real denom = 0.0;
            {
              for (int i = 0; i < 3; ++i) {
                Real t = 0.0;
                for (int j = 0; j < 3; ++j) {
                  t += best_shiftQ(j,i) * extr[j];
                }
                numer += t * direc[i];
                denom += t;
              }
            }
            Real lambda = numer / denom;
            if (lambda > eigenval_ub)
              eigenval_ub = lambda;
            if (lambda < eigenval_lb)
              eigenval_lb = lambda;
          }
        }

        // Generate rows for the bottom part.
          
        {
          int rowcount = 0;
          for (int evrownum = 0; evrownum < monomcount; ++evrownum) {
            if (use_monom[evrownum] == 0.0) continue;

            for (int jj = 0; jj < tmcol; ++jj) {
              expaArow[jj] = 0.0;
              expaBrow[jj] = 0.0;
            }
            pair<int,int> rval = colnum_to_deg_pair(evrownum, 2 * degree - 2);
            int ud = rval.first;
            int vd = rval.second;
            int wd = 2 * degree - 2 - ud - vd;
            int colnum1 = deg_pair_to_colnum(ud + 1, vd, 2*degree - 1);
            expaArow[colnum1] = direc[0];
            expaBrow[colnum1] = 1.0;
            colnum1 = deg_pair_to_colnum(ud, vd + 1, 2*degree - 1);
            expaArow[colnum1] = direc[1];
            expaBrow[colnum1] = 1.0;
            colnum1 = deg_pair_to_colnum(ud,vd, 2*degree - 1);
            expaArow[colnum1] = direc[2];
            expaBrow[colnum1] = 1.0;
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << "expaArow(" << evrownum + 1 << ",:) = [";
              {
                for (int i = 0; i < tmcol; ++i) {
                  *debug_str << setprecision(17) << expaArow[i] << " ";
                }
              }
              *debug_str << "];\n expaBrow(" << evrownum + 1 << ",:) = [";
              {
                for (int i = 0; i < tmcol; ++i) {
                  *debug_str << setprecision(17) << expaBrow[i] << " ";
                }
              }
              *debug_str << "];\n";
            }
#endif          
            
            // Apply HH transforms 
            
            for (int j = 0; j < tmrow; ++j) {
              Real ip1 = 0.0;
              Real ip2 = 0.0;
              {
                for (int i = j; i < tmcol; ++i) {
                  ip1 += hhmat(i,j) * expaArow[i];
                  ip2 += hhmat(i,j) * expaBrow[i];
                }
              }
              ip1 *= betavec[j];
              ip2 *= betavec[j];
              {
                for (int i = j; i < tmcol; ++i) {
                  expaArow[i] -= ip1 * hhmat(i,j);
                  expaBrow[i] -= ip2 * hhmat(i,j);
                }
              }
            }
            
            for (int i = tmrow; i < tmcol; ++i) {
              A(rowcount, i - tmrow) = Complex(expaArow[i], 0.0);
              B(rowcount, i - tmrow) = Complex(expaBrow[i], 0.0);
            }
            ++rowcount;
          }
#ifdef DEBUGGING
          if (rowcount != evsize)
            throw_error("Counting problem");
#endif
        }

        {
          for (int i = 0; i < evsize; ++i)
            for (int j = 0; j < evsize; ++j) {
              Ao(i,j) = A(i,j);
              Bo(i,j) = B(i,j);
            }
        }

#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "];\nAi = [\n";
          {
            for (int i = 0; i < evsize; ++i) {
              for (int j = 0; j < evsize; ++j) {
                *debug_str << Format_Complex(A(i,j)) << " ";
              }
              *debug_str << '\n';
            }
          }
          *debug_str << "];\nBi = [\n";
          {
            for (int i = 0; i < evsize; ++i) {
              for (int j = 0; j < evsize; ++j) {
                *debug_str << Format_Complex(B(i,j)) << " ";
              }
              *debug_str << '\n';
            }
          }
          *debug_str << "];\n";
        }
#endif
        
        
        
        Real cond = generalized_eigensolver(A, B, vr, eveccond, eigenval_lb, 
          eigenval_ub, workspace);
        
        
        if (cond < 0) {
#ifdef DEBUGGING
          *debug_str << "Lapack failure in Solve_BezTri::solve\n";
          *debug_str << "A = ";
          {
            for (int i = 0; i < evsize; ++i) {
              for (int j = 0; j < evsize; ++j) {
                *debug_str << Format_Complex(Ao(i,j)) << " ";
              }
              *debug_str << '\n';
            } 
          }
          *debug_str << "B = ";
          {
            for (int i = 0; i < evsize; ++i) {
              for (int j = 0; j < evsize; ++j) {
                *debug_str << Format_Complex(Bo(i,j)) << " ";
              }
              *debug_str << '\n';
            }
          }  
          *debug_str << "Afinal = ";
          {
            for (int i = 0; i < evsize; ++i) {
              for (int j = 0; j < evsize; ++j) {
                *debug_str << Format_Complex(A(i,j)) << " ";
              }
              *debug_str << '\n';
            } 
          }
          *debug_str << "Bfinal = ";
          {
            for (int i = 0; i < evsize; ++i) {
              for (int j = 0; j < evsize; ++j) {
                *debug_str << Format_Complex(B(i,j)) << " ";
              }
              *debug_str << '\n';
            }
          }
#endif
          if (++lapack_failcount == 5) {
            throw_error("Repeated lapack failure");
          }
          continue;
        }
        
#ifdef DEBUGGING
        if (dump_everything)
          *debug_str << "cond = " << cond 
          << "\nbest_cond = " << best_cond 
          << "\nallowable_cond = " << allowable_cond << "\n";
#endif
        
        if (best_cond < 0 || cond < best_cond) {
          best_cond = cond; 
          {
            for (int i = 0; i < evsize; ++i) {
              for (int j = 0; j < evsize; ++j) {
                bestvr(i,j) = vr(i,j);
              }
            }
          }
          {
            for (int i = 0; i < evsize; ++i) {
              best_eveccond[i] = eveccond[i];
            }
          }
        }
        if (cond < allowable_cond || allowable_cond > max_allowable_cond)
          break;
        
        allowable_cond *= 10.0;
      }
    
        
      // Find the polynomial solutions from the eigenvectors.

      PatchMath::Workspace::ComplexVector transfevec(workspace, tmcol);
      PatchMath::Workspace::ComplexMatrix jac(workspace, 2, 2);
      
      {
        for (int vlrow = 0; vlrow < evsize; ++vlrow) {
#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << " vri = [\n";
            for (int i = 0; i < evsize; ++i) {
              *debug_str << Format_Complex(bestvr(i,vlrow)) << " ";
            }
            *debug_str << "\n];\n";
          }
#endif
          
          // Apply the HH transforms to get the transformed eigenvector.          
          
          {
            for (int j = 0; j < tmrow; ++j)
              transfevec[j] = Complex(0.0,0.0);
          }
          
          {
            for (int j = tmrow; j < tmcol; ++j) 
              transfevec[j] = bestvr(j - tmrow, vlrow);
          }
          {
            for (int j = tmrow - 1; j >= 0; --j) {
              Complex ip(0.0,0.0);
              {
                for (int i = j; i < tmcol; ++i) {
                  ip += Complex(hhmat(i,j),0.0) * transfevec[i];
                }
              }
              ip *= betavec[j];
              {
                for (int i = j; i < tmcol; ++i) {
                  transfevec[i] -= ip * hhmat(i,j);
                }
              }
            }
          }

          // Find whether u^d, v^d, or w^d is the largest coeff.
          // Then determine u,v,w from powers close to the largest power.

          Complex ud = transfevec[deg_pair_to_colnum(2 * degree - 1, 0, 2 * degree - 1)];
          Complex vd = transfevec[deg_pair_to_colnum(0, 2 * degree - 1, 2 * degree - 1)];
          Complex wd = transfevec[deg_pair_to_colnum(0, 0, 2 * degree - 1)];
          Complex u1, v1, w1;
          if (abs(ud) >= abs(vd) && abs(ud) >= abs(wd)) {
            u1 = ud;
            v1 = transfevec[deg_pair_to_colnum(2 * degree - 2, 1, 2 * degree - 1)];
            w1 = transfevec[deg_pair_to_colnum(2 * degree - 2, 0, 2 * degree - 1)];
          }
          else if (abs(vd) >= abs(wd)) {
            u1 = transfevec[deg_pair_to_colnum(1, 2 * degree - 2, 2 * degree - 1)];
            v1 = vd;
            w1 = transfevec[deg_pair_to_colnum(0, 2 * degree - 2, 2 * degree - 1)];
          }
          else {
            u1 = transfevec[deg_pair_to_colnum(1, 0, 2 * degree - 1)];
            v1 = transfevec[deg_pair_to_colnum(0, 1, 2 * degree - 1)];
            w1 = wd;
          }

          // Apply the inverse of shiftQ to (u1,v1,w1) to undo
          // undo effect of shift.

          Complex u2 = u1 * best_shiftQ(0,0) + v1 * best_shiftQ(0,1) + w1 * best_shiftQ(0,2);
          Complex v2 = u1 * best_shiftQ(1,0) + v1 * best_shiftQ(1,1) + w1 * best_shiftQ(1,2);
          Complex w2 = u1 * best_shiftQ(2,0) + v1 * best_shiftQ(2,1) + w1 * best_shiftQ(2,2);
          
          Complex denom = u2 + v2 + w2;
          if (abs(denom) == 0) {
#ifdef DEBUGGING
            if (dump_everything)
              *debug_str << " discarding, denom = 0\n";
#endif
            continue;
          }

          C2Point sol(u2 / denom, v2 / denom);

#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << "transfevec = [\n";
            for (int i = 0; i < tmcol; ++i)
              *debug_str << Format_Complex(transfevec[i]) << " ";
            *debug_str << "\n u1i = " 
                       << Format_Complex(u1) 
                       << ";\n v1i = " << Format_Complex(v1) 
                       << ";\n w1i = " << Format_Complex(w1)
                       << ";\n u2i = " << Format_Complex(u2) 
                       << ";\n v2i = " << Format_Complex(v2) 
                       << ";\nw2i = " << Format_Complex(w2)
                       << ";\n denomi = " << Format_Complex(denom) 
                       << ";\n sol = " << Format_Complex(sol.coord[0]) << " " << Format_Complex(sol.coord[1])
                       << ";\n";
          }
#endif
            

          if (sol.l2norm() > 2) {
#ifdef DEBUGGING
            if (dump_everything)
              *debug_str << "discarding sol; too large";
#endif
            continue;
          }

          C2Point fval, fvalcopy, direc, nstep, newsol, newfval;
          
          Real invnormest;
                   
          // Try to improve the solution with two steps of Newton's method
          // if newton_improve is true.
       
          
          for (int newtcount = 0; true; ++newtcount) {
            
            fval.coord[0] = eval_complex(sol, 0, degree);
            fval.coord[1] = eval_complex(sol, 1, degree);
            direc.coord[0] = 1.0;
            direc.coord[1] = 0.0;
            jac(0,0) = eval_direc_deriv_complex(sol, direc, 0, degree);
            jac(1,0) = eval_direc_deriv_complex(sol, direc, 1, degree);
            direc.coord[0] = 0.0;
            direc.coord[1] = 1.0;
            jac(0,1) = eval_direc_deriv_complex(sol, direc, 0, degree);
            jac(1,1) = eval_direc_deriv_complex(sol, direc, 1, degree);
            fvalcopy = fval;

            // solve 2-by-2 linear system
            
            invnormest = linear_solve_C2(jac, fval, nstep);

            
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << "soli = [" 
                         <<  Format_Complex(sol.coord[0]) << " ," 
                         << Format_Complex(sol.coord[1])
                         << "];\n fvali = [" << Format_Complex(fval.coord[0])
                         << "," << Format_Complex(fval.coord[1])
                         << "\n];\n jac1 = [" 
                         << Format_Complex(jac(0,0)) 
                         << "," << Format_Complex(jac(0,1)) 
                         << ";" << Format_Complex(jac(1,0) )
                         << " " << Format_Complex(jac(1,1));
              *debug_str << "];\n nstep = [" 
                         << Format_Complex(nstep.coord[0])
                         << "," << Format_Complex(nstep.coord[1])
                         << "];\n invnormest = " << invnormest;
              *debug_str << ";\n";
            }
            
#endif
            if (!newton_improve) break;
            if (newtcount == 2) break;
            if (invnormest >= BIG_REAL) break;

         
            newsol.coord[0] = sol.coord[0] - nstep.coord[0];
            newsol.coord[1] = sol.coord[1] - nstep.coord[1];
            newfval.coord[0] = eval_complex(newsol, 0, degree);
            newfval.coord[1] = eval_complex(newsol, 1, degree);
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << " newsoli = [" 
                         << Format_Complex(newsol.coord[0]) 
                         << ", " << Format_Complex(newsol.coord[1])
                         << "]; newfvali = [" 
                         << Format_Complex(newfval.coord[0]) 
                         << ", " << Format_Complex(newfval.coord[1]) 
                         << "];\n";
            }
#endif
            if (newfval.l2norm() < fval.l2norm()) {
#ifdef DEBUGGING
              if (dump_everything) {
                *debug_str << " changing\n";
              }
#endif
              sol = newsol;
            }
            else {
              break;
            }
          }
          
          Real fvalnorm = fval.l2norm();
          Real fvalcutoff = evsize * tol * scfac;
          Real param_uncertainty;
          
          if (invnormest >= BIG_REAL)
            param_uncertainty = 1.0;
          else 
            param_uncertainty = (fvalnorm + fvalcutoff) * invnormest;

#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << " fvalnorm = " << fvalnorm << " fvalcutoff = "
              << fvalcutoff << " param_uncertainty = "
              << param_uncertainty <<  "\n";
          }
#endif
          if (abs(sol.coord[0]) + abs(sol.coord[1]) > 3) {
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << " dropping root; norm too large";
            }
#endif
            continue;
          }
          if (coefs_were_perturbed && 
              fvalnorm > evsize * expected_resd * best_eveccond[vlrow]) {
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << "residual exceeds expected resd " 
                         << evsize << " * " << expected_resd << " * "
                         << best_eveccond[vlrow] << " = "
                         << evsize * expected_resd * best_eveccond[vlrow]
                         << " discarding root\n";
            }
#endif
            continue;
          }
            
            
          
          if (sol.coord[0].real() >= -param_uncertainty && 
            sol.coord[1].real() >= -param_uncertainty &&
            1.0 - sol.coord[0].real() - sol.coord[1].real() >= -param_uncertainty &&
            fabs(sol.coord[0].imag()) <= param_uncertainty &&
            fabs(sol.coord[1].imag()) <= param_uncertainty) {
            
            PatchMath::IntersectRecord result;
            result.paramcoord[0] = sol.coord[0].real();
            result.paramcoord[1] = sol.coord[1].real();
            result.degen = (fabs(sol.coord[0].real()) <= param_uncertainty || 
                fabs(sol.coord[1].real()) <= param_uncertainty ||
                fabs(1.0 - sol.coord[0].real() - sol.coord[1].real()) <= param_uncertainty);
            result.param_uncertainty = param_uncertainty;
#ifdef DEBUGGING
            if (dump_everything) {
              *debug_str << " passed second test, degen = " << result.degen << "\n";
            }
#endif
            output_stack.push_back(result);
          }
        }
      }
      return is_degen;
    }
  }

  // Helper class for math with bezier quad patches.

  class Solve_BezQuad {
  public:
    virtual Real get_bez_coeff(int eqnum, int j) const = 0;
    Solve_BezQuad() { }
    virtual ~Solve_BezQuad() { }
    bool solve(int degree1[2],
      int degree2[2],
      vector<PatchMath::IntersectRecord>& output_stack,
      PatchMath::Workspace& workspace, 
      double scfac, 
      double tol,
      bool newton_improve) const;
    Real eval(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Complex eval_complex(const C2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Real eval_partialu(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Complex eval_partialu_complex(const C2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Real eval_partialv(const R2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
    Complex eval_partialv_complex(const C2Point& param,
      int eqnum,
      int degree1,
      int degree2) const;
  };

  // Derived class for intersecting a patch with a line.
  
  class Solve_BezQuad1 : public Solve_BezQuad {
  private: 
    const Matrix& lineeq_;
    const Point& linerhs_;
    const PatchMath& pm_;
  public:
    Solve_BezQuad1(const Matrix& lineeq,
      const Point& linerhs,
      const PatchMath& pm) : lineeq_(lineeq), linerhs_(linerhs), pm_(pm) { }
    Real get_bez_coeff(int eqnum, int cpnum) const {
      return lineeq_(eqnum,0) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 0) +
        lineeq_(eqnum,1) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 1) +
        lineeq_(eqnum,2) * PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, 2) -
        linerhs_[eqnum];
    }
  };

  // derived class for directional derivative.

  class Solve_BezQuad2 : public Solve_BezQuad {
  private:
    const Point& direc_;
    int di_;
    const PatchMath& pm_;
    int degree1_;
    int degree2_;
  public:
    Solve_BezQuad2(const Point& direc, 
      int di,
      const PatchMath& pm,
      int degree1,
      int degree2) : direc_(direc), di_(di), pm_(pm), degree1_(degree1), degree2_(degree2)
    { }
    Real get_bez_coeff(int eqnum, int cpnum) const;
  };

  // derived class for evaluating one dimension.
    
  class Solve_BezQuad3 : public Solve_BezQuad {
  private:
    int select_dim_;
    const PatchMath& pm_;
  public:
    Solve_BezQuad3(int select_dim, 
      const PatchMath& pm) : select_dim_(select_dim), pm_(pm) 
    { }
    Real get_bez_coeff(int, int cpnum) const {
      return PatchMath::Bezier_Helper::control_point_coord(pm_, cpnum, select_dim_);
    }
  };

    
  Real Solve_BezQuad2::get_bez_coeff(int eqnum, int cpnum) const {
    Real t = 0.0;
    int ocp1, ocp2, fac;
    if (eqnum == 0) {
      int d2 = cpnum / degree1_;
      int d1 = cpnum % degree1_;
      ocp1 = d1 + 1 + d2 * (degree1_ + 1);
      ocp2 = d1 + d2 * (degree1_ + 1);
      fac = degree1_;
    }
    else {
      int d2 = cpnum / (degree1_ + 1);
      int d1 = cpnum % (degree1_ + 1);
      ocp1 = d1 + (d2 + 1) * (degree1_ + 1);
      ocp2 = d1 + d2 * (degree1_ + 1);
      fac = degree2_;
    }
    for (int j = 0; j < di_; ++j) {
        t += (PatchMath::Bezier_Helper::control_point_coord(pm_, ocp1, j) -
          PatchMath::Bezier_Helper::control_point_coord(pm_, ocp2, j)) * fac 
          * direc_[j];
    }
    return t;
  }

  // Horner algorithm for evaluating bezier quad (parametric
  // to real conversion)


  Real Solve_BezQuad::eval(const R2Point& param,
    int eqnum,
    int degree1,
    int degree2) const {

    Real u = param.coord[0];
    Real oppu = 1.0 - u;
    Real v = param.coord[1];
    Real oppv = 1.0 - v;

    int outer_binom = 1;
    Real outer_result = 0.0;
    Real pwrv = 1.0;
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Real pwru = 1.0;
      Real inner_result = 0.0;
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result *= oppu;
        inner_result += get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1) *
          inner_binom * pwru;
        pwru *= u;
        inner_binom *= (degree1 - d1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * outer_binom * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2);
      outer_binom /= (d2 + 1);
    }
    return outer_result;
  }

  
  // Horner algorithm for evaluating bezier quad (parametric
  // to real conversion)


  Complex Solve_BezQuad::eval_complex(const C2Point& param,
    int eqnum,
    int degree1,
    int degree2) const {

    Complex u = param.coord[0];
    Complex oppu = Complex(1.0,0.0) - u;
    Complex v = param.coord[1];
    Complex oppv = Complex(1.0,0.0) - v;

    int outer_binom = 1;
    Complex outer_result(0.0,0.0);
    Complex pwrv(1.0,0.0);
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Complex pwru(1.0,0.0);
      Complex inner_result(0.0,0.0);
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result *= oppu;
        inner_result += Complex(get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1), 0.0) *
          Complex(static_cast<Real>(inner_binom), 0.0) * pwru;
        pwru *= u;
        inner_binom *= (degree1 - d1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * Complex(static_cast<Real>(outer_binom), 0.0) * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2);
      outer_binom /= (d2 + 1);
    }
    return outer_result;
  }
  
  
  Real Solve_BezQuad::eval_partialu(const R2Point& param,
    int eqnum,
    int degree1,
    int degree2) const { 

    Real u = param.coord[0];
    Real oppu = 1.0 - u;
    Real v = param.coord[1];
    Real oppv = 1.0 - v;
    int outer_binom = 1;
    Real outer_result = 0.0;
    Real pwrv = 1.0;
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Real pwru = 1.0;
      Real inner_result = 0.0;
      Real prev_coef = get_bez_coeff(eqnum, d2 * (degree1 + 1));
      for (int d1 = 0; d1 < degree1; ++d1) {
        inner_result *= oppu;
        Real next_coef = get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1 + 1);
        inner_result += (next_coef - prev_coef) * inner_binom * pwru;
        prev_coef = next_coef;
        pwru *= u;
        inner_binom *= (degree1 - d1 - 1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * outer_binom * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2);
      outer_binom /= (d2 + 1);
    }
    return outer_result * degree1;
  }

  
  Complex Solve_BezQuad::eval_partialu_complex(const C2Point& param,
    int eqnum,
    int degree1,
    int degree2) const { 

    Complex u = param.coord[0];
    Complex oppu = Complex(1.0,0.0) - u;
    Complex v = param.coord[1];
    Complex oppv = Complex(1.0,0.0) - v;
    int outer_binom = 1;
    Complex outer_result(0.0,0.0);
    Complex pwrv(1.0,0.0);
    for (int d2 = 0; d2 <= degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Complex pwru(1.0,0.0);
      Complex inner_result(0.0,0.0);
      Real prev_coef = get_bez_coeff(eqnum, d2 * (degree1 + 1));
      for (int d1 = 0; d1 < degree1; ++d1) {
        inner_result *= oppu;
        Real next_coef = get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1 + 1);
        inner_result += Complex(next_coef - prev_coef, 0.0) * 
          Complex(static_cast<Real>(inner_binom), 0.0) * pwru;
        prev_coef = next_coef;
        pwru *= u;
        inner_binom *= (degree1 - d1 - 1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * 
        Complex(static_cast<Real>(outer_binom), 0.0) * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2);
      outer_binom /= (d2 + 1);
    }
    return outer_result * Complex(static_cast<Real>(degree1), 0.0);
  }


  Real Solve_BezQuad::eval_partialv(const R2Point& param,
    int eqnum,
    int degree1,
    int degree2) const { 

    Real u = param.coord[0];
    Real oppu = 1.0 - u;
    Real v = param.coord[1];
    Real oppv = 1.0 - v;
    int outer_binom = 1;
    Real outer_result = 0.0;
    Real pwrv = 1.0;
    for (int d2 = 0; d2 < degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Real pwru = 1.0;
      Real inner_result = 0.0;
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result *= oppu;
        Real coef =  get_bez_coeff(eqnum, (d2 + 1) * (degree1 + 1) + d1) -
          get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1);
        inner_result += coef * inner_binom * pwru;
        pwru *= u;
        inner_binom *= (degree1 - d1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * outer_binom * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2 - 1);
      outer_binom /= (d2 + 1);
    }
    return outer_result * degree2;
  }



  Complex Solve_BezQuad::eval_partialv_complex(const C2Point& param,
    int eqnum,
    int degree1,
    int degree2) const { 

    Complex u = param.coord[0];
    Complex oppu = Complex(1.0,0.0) - u;
    Complex v = param.coord[1];
    Complex oppv = Complex(1.0,0.0) - v;
    int outer_binom = 1;
    Complex outer_result(0.0,0.0);
    Complex pwrv(1.0,0.0);
    for (int d2 = 0; d2 < degree2; ++d2) {
      // Compute the inner coefficient.

      int inner_binom = 1;
      Complex pwru(1.0,0.0);
      Complex inner_result(0.0,0.0);
      for (int d1 = 0; d1 <= degree1; ++d1) {
        inner_result *= oppu;
        Real coef =  get_bez_coeff(eqnum, (d2 + 1) * (degree1 + 1) + d1) -
          get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1);
        inner_result += Complex(coef, 0.0) * 
          Complex(static_cast<Real>(inner_binom), 0.0) * pwru;
        pwru *= u;
        inner_binom *= (degree1 - d1);
        inner_binom /= (d1 + 1);
      }
      outer_result *= oppv;
      outer_result += inner_result * 
        Complex(static_cast<Real>(outer_binom), 0.0) * pwrv;
      pwrv *= v;
      outer_binom *= (degree2 - d2 - 1);
      outer_binom /= (d2 + 1);
    }
    return outer_result * Complex(static_cast<Real>(degree2), 0.0);
  }

  // Routine to solve a 2-var polynomial system in 
  // tensor product Bezier form.  Use Macaulay
  // resultant combined with u-resultant with randomized direction.
  // Return true if degenerate.
  
  bool Solve_BezQuad::solve(int udegree[2],
    int vdegree[2],
    vector<PatchMath::IntersectRecord>& output_stack,
    PatchMath::Workspace& workspace,
    double scfac, 
    double tol,
    bool newton_improve) const {
    
    PatchMath::Workspace::StartPosMarker startposmarker_(workspace);

    using std::abs;
    using std::norm;

#ifdef DEBUGGING
    static int invocation_count = 0;
    ++invocation_count;

    using QMG::Local_to_PatchMath_CPP::debug_str;
    using QMG::Local_to_PatchMath_CPP::dump_everything;
    using std::setprecision;
    using std::setiosflags;
    using std::ios;

    // *debug_str << "invoc count = " << invocation_count << '\n';

    // dump_everything = (invocation_count > 3558 && invocation_count < 3565);

    /*
    if (invocation_count == 587)
      dump_everything = true;
    if (invocation_count == 590)
      dump_everything = false;
*/


    if (dump_everything) {
      *debug_str << " In solve bez quad\n udegree = [" << udegree[0] << ',' << udegree[1]
        << "];\n vdegree = [" << vdegree[0] << ',' << vdegree[1] 
        << "];\n scfac = " << scfac << ";\n tol = " << tol << ";\n bez1 = [\n";
      int ncp1 = (udegree[0] + 1) * (vdegree[0] + 1); 
      {
        for (int j = 0; j < ncp1; ++j) 
          *debug_str << std::setprecision(17) << get_bez_coeff(0,j) << " ";
      }
      *debug_str << "];\n";
      int ncp2 = (udegree[1] + 1) * (vdegree[1] + 1); 
      *debug_str << "\n bez2 = [\n";
      {
        for (int j = 0; j < ncp2; ++j) 
          *debug_str << std::setprecision(17) << get_bez_coeff(1,j) << " ";
      }
      *debug_str << "];\n";    
      *debug_str << "invocation count = " << invocation_count << '\n';
    }
#endif

    // Check if (0,0) is in the bounding box; if not then early return.

    {
      for (int eqnum = 0; eqnum < 2; ++eqnum) {
        int ncp1 = (udegree[eqnum] + 1) * (vdegree[eqnum] + 1);
        bool pos_seen = false;
        bool neg_seen = false;
        for (int i = 0; i < ncp1; ++i) {
          Real bz = get_bez_coeff(eqnum, i);
          if (bz <= scfac * tol)
            neg_seen = true;
          if (bz >= -scfac * tol)
            pos_seen = true;
        }
        if (!neg_seen || !pos_seen) {
#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << " early return; equation " << eqnum << " misses 0.\n";
          }
#endif
          return false;
        }
      }
    }

   
    int tmrow = udegree[0] * vdegree[0] + udegree[1] * vdegree[1] + 
      udegree[0] + udegree[1] + vdegree[0] + vdegree[1];
    int umdegree = udegree[0] + udegree[1];
    int vmdegree = vdegree[0] + vdegree[1];
    int tmcol = (umdegree + 1) * (vmdegree + 1) - 1;
    int evsize = udegree[0] * vdegree[1] + udegree[1] * vdegree[0]; // = tmcol-tmrow
    int num_control_pt0 = (udegree[0] + 1) * (vdegree[0] + 1);
    int num_control_pt1 = (udegree[1] + 1) * (vdegree[1] + 1);
    int max_ncp = (num_control_pt1 > num_control_pt0)? num_control_pt1 : 
      num_control_pt0;
    int totdeg0 = udegree[0] + vdegree[0];
    int totdeg1 = udegree[1] + vdegree[1];
    int mxdegree = (totdeg1 > totdeg0)? totdeg1 : totdeg0;
    int monomcount = umdegree * vmdegree;


    PatchMath::Workspace::RealMatrix coef_perturb(workspace, 2, max_ncp);
    Real bezmax_poly[2];
    bezmax_poly[0] = 0.0;
    bezmax_poly[1] = 0.0;
    {
      for (int eqnum = 0; eqnum < 2; ++eqnum) {
        int ncp = (eqnum == 0)? num_control_pt0 : num_control_pt1;
        for (int i = 0; i < ncp; ++i) {
          coef_perturb(eqnum,i) = 0.0;
          Real c1 = fabs(get_bez_coeff(eqnum,i));
          if (c1 > bezmax_poly[eqnum])
            bezmax_poly[eqnum] = c1;
        }
      }
    }

    bool is_degen = false;
    if (bezmax_poly[0] == 0.0 || bezmax_poly[1] == 0.0) {
      is_degen = true;
      return is_degen;
    }

    // Compute two fractional linear transformations to ensure the
    // extreme coefficients are nonzero.

    PatchMath::Workspace::RealMatrix test_shiftpoly(workspace, 2, max_ncp);
    PatchMath::Workspace::RealMatrix best_shiftpoly(workspace, 2, max_ncp);
    PatchMath::Workspace::RealVector pwrv(workspace, max_ncp);
    PatchMath::Workspace::RealVector pwruv(workspace, max_ncp);
    PatchMath::Workspace::RealVector inner_result(workspace, max_ncp);
    PatchMath::Workspace::RealVector temp_poly(workspace, max_ncp);

    Real best_ushiftval;
    Real best_vshiftval;
    Real best_lead;

    for (int shift_trycount = 0; shift_trycount < 4; ++shift_trycount) {
      Real test_ushiftval = random_real() + 1.0;
      Real test_vshiftval = random_real() + 1.0;
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "shift trycount = " << shift_trycount 
                   << "\nshifts = [" << setprecision(17) << test_ushiftval
                   << ", " << setprecision(17) << test_vshiftval
                   << "];\n";
      }
#endif
      // Shift:  (u,oppu) = ushift * (u',oppu').
      //       (v,oppv) = vshift * (v',oppv')
      // Here ushift = [1,ushiftval;-ushiftval,1]
      // and similarly for vshift.
      // Use horner's rule symbolically to make the change of
      // variables.

      for (int eqnum = 0; eqnum < 2; ++eqnum) {
        int udegr = udegree[eqnum];
        int vdegr = vdegree[eqnum];
        int ncp = (udegr + 1) * (vdegr + 1);
        // shiftpoly = 0; pwrv = 1
        {
          for (int i = 0; i < ncp; ++i) {
            test_shiftpoly(eqnum,i) = 0.0;
            pwrv[i] = 0.0;
          }
        }
        pwrv[0] = 1.0;

        int outer_binom = 1;
        for (int vd = 0; vd <= vdegr; ++vd) {
          int inner_binom = 1;
          // pwruv = pwrv; inner_result = 0;
          {
            for (int i = 0; i < ncp; ++i) {
              pwruv[i] = pwrv[i];
              inner_result[i] = 0.0;
            }
          }
          for (int ud = 0; ud <= udegr; ++ud) {
            // temp_poly = inner_result * oppu;
            //           = inner_result - ushiftval * u' * inner_result
            {
              for (int udi = 0; udi <= udegr; ++udi)
                for (int vdi = 0; vdi <= vdegr; ++vdi) {
                  int cpn = udi + (vdi * (udegr + 1));
                  temp_poly[cpn] = inner_result[cpn] -
                    test_ushiftval *
                    ((udi > 0)? inner_result[cpn - 1] : 0.0);
                }
            }
            
            // inner_result = temp_poly + get_bez_coeff(eqnum, d2 * (degree1 + 1) + d1) *
            //  inner_binom * pwruv;

            {
              Real g = get_bez_coeff(eqnum, vd * (udegr + 1) + ud) 
                * inner_binom / bezmax_poly[eqnum];
              for (int udi = 0; udi <= udegr; ++udi) 
                for (int vdi = 0; vdi <= vdegr; ++vdi) {
                  int cpn = udi + (vdi * (udegr + 1));
                  inner_result[cpn] = temp_poly[cpn] + g * pwruv[cpn];
                }
            }

            // temp_poly = pwruv * u';
            {
              for (int udi = 0; udi <= udegr; ++udi)
                for (int vdi = 0; vdi <= vdegr; ++vdi) {
                  int cpn = udi + (vdi * (udegr + 1));
                  temp_poly[cpn] = ((udi > 0)? pwruv[cpn - 1] : 0.0) +
                    pwruv[cpn] * test_ushiftval;
                }
            }
            // pwruv = temp_poly;
            {
              for (int udi = 0; udi <= udegr; ++udi)
                for (int vdi = 0; vdi <= vdegr; ++vdi) {
                  int cpn = udi + (vdi * (udegr + 1));
                  pwruv[cpn] = temp_poly[cpn];
                }
            }         
            inner_binom *= (udegr - ud);
            inner_binom /= (ud + 1);
          }
          // temp_poly = shiftpoly * oppv;
          //           = shiftpoly - vshiftval * v' * shiftpoly
          {
            for (int udi = 0; udi <= udegr; ++udi)
              for (int vdi = 0; vdi <= vdegr; ++vdi) {
                int cpn = udi + (vdi * (udegr + 1));
                temp_poly[cpn] = test_shiftpoly(eqnum,cpn) -
                  test_vshiftval * 
                  ((vdi > 0)? test_shiftpoly(eqnum, cpn - udegr - 1) : 0.0);
              }
          }
          // shiftpoly = temp_poly + outer_binom * inner_result

          {
            for (int udi = 0; udi <= udegr; ++udi)
              for (int vdi = 0; vdi <= vdegr; ++vdi) {
                int cpn = udi + (vdi * (udegr + 1));
                test_shiftpoly(eqnum,cpn) = temp_poly[cpn] + 
                  outer_binom * inner_result[cpn];
              }
          }

          //temp_poly = pwrv * v;
          {
            for (int udi = 0; udi <= udegr; ++udi)
              for (int vdi = 0; vdi <= vdegr; ++vdi) {
                int cpn = udi + (vdi * (udegr + 1));
                temp_poly[cpn] = ((vdi > 0)? pwrv[cpn - udegr - 1] : 0.0) +
                test_vshiftval *  pwrv[cpn];
              }
          }
          //pwrv = temp_poly
          
          {
            for (int udi = 0; udi <= udegr; ++udi)
              for (int vdi = 0; vdi <= vdegr; ++vdi) {
                int cpn = udi + (vdi * (udegr + 1));
                pwrv[cpn] = temp_poly[cpn];
              }
          }
          outer_binom *= (vdegr - vd);
          outer_binom /= (vd + 1);
        }
      }
      // Now check if the shift is successful.


#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "shiftpoly0 = [\n";
        int ncp1 = (udegree[0] + 1) * (vdegree[0] + 1); 
        {
          for (int j = 0; j < ncp1; ++j) 
            *debug_str << std::setprecision(17) << test_shiftpoly(0,j) << " ";
        }
        *debug_str << "];\n";
        int ncp2 = (udegree[1] + 1) * (vdegree[1] + 1); 
        *debug_str << "shiftpoly1 = [\n";
        {
          for (int j = 0; j < ncp2; ++j) 
            *debug_str << std::setprecision(17) << test_shiftpoly(1,j) << " ";
        }
        *debug_str << "];\n";    
      }
#endif

      Real test_lead = BIG_REAL;
      
      for (int eqnum2 = 0; eqnum2 < 2; ++eqnum2) {
        int udegr = udegree[eqnum2];
        int vdegr = vdegree[eqnum2];
        for (int ud1 = 0; ud1 < udegr + 1; ud1 += udegr) {
          for (int vd1 = 0; vd1 < vdegr + 1; vd1 += vdegr) {
            Real c1 = fabs(test_shiftpoly(eqnum2, ud1 + vd1 * (udegr + 1)));
            if (c1 < test_lead) {
              test_lead = c1;
            }
          }
        }
      }
      
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << " testlead = " << test_lead << '\n';
      }
#endif
      
      if (shift_trycount == 0 || test_lead > best_lead) {

#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << " keeping this one\n";
        }
#endif

        best_lead = test_lead;
        best_ushiftval = test_ushiftval;
        best_vshiftval = test_vshiftval;
        for (int eqnum3 = 0; eqnum3 < 2; ++eqnum3) {
          int ncp = (udegree[eqnum3] + 1) * (vdegree[eqnum3] + 1);
          for (int j = 0; j < ncp; ++j)
            best_shiftpoly(eqnum3,j) = test_shiftpoly(eqnum3,j);
        }
      }
    }



    // Compute the top part of the matrix.  Make sure it's full
    // rank.  If not, must perturb.


    PatchMath::Workspace::RealVector rowtype(workspace, tmrow);
    PatchMath::Workspace::RealVector row_u_deg(workspace, tmrow);
    PatchMath::Workspace::RealVector row_v_deg(workspace, tmrow);
    PatchMath::Workspace::RealMatrix Tm(workspace, tmcol, tmrow);    
    PatchMath::Workspace::RealMatrix Tmp(workspace, tmrow, tmcol);
    PatchMath::Workspace::RealMatrix hhmat_tmp(workspace, tmrow, tmrow);
    PatchMath::Workspace::RealVector betavec_tmp(workspace, tmrow);
    PatchMath::Workspace::RealVector select_col_tmp(workspace, tmrow);
    PatchMath::Workspace::RealVector fac(workspace, tmrow);
    PatchMath::Workspace::RealMatrix depmat(workspace, tmrow, tmrow);

    bool coefs_were_perturbed = false;
    Real expected_resd = scfac * tol;

    for (int top_part_trycount = 0;; ++top_part_trycount) {
    
      int rowind = 0;
      for (int eqnum = 0; eqnum < 2; ++eqnum) {
        for (int ud = 0; ud <= udegree[1-eqnum]; ++ud) {
          for (int vd = 0; vd <= vdegree[1-eqnum]; ++vd) {
            if (ud == udegree[1 - eqnum] && vd == vdegree[1 - eqnum]) 
              continue;
            rowtype[rowind] = eqnum;
            row_u_deg[rowind] = ud;
            row_v_deg[rowind] = vd;
            for (int j = 0; j < tmcol; ++j)
              Tm(j, rowind) = 0.0; 
            for (int ud1 = 0; ud1 <= udegree[eqnum]; ++ud1) {
              for (int vd1 = 0; vd1 <= vdegree[eqnum]; ++vd1) {
                int mtxcol = ud + ud1 + (vd + vd1) * (umdegree + 1);
                int cpnum = ud1 + vd1 * (udegree[eqnum] + 1);
                Real bz = best_shiftpoly(eqnum, cpnum);
                Real pt = coef_perturb(eqnum,cpnum);
                Tm(mtxcol, rowind) = bz + pt;
              }
            }
            ++rowind;
          }
        }
      }
      if (top_part_trycount == 1) {
        break;
      }
      {
        for (int i = 0; i < tmcol; ++i) {
          for (int j = 0; j < tmrow; ++j) {
            Tmp(j,i) = Tm(i,j);
          }
        }
      }

      QR_Factor_with_col_pivoting(Tmp, select_col_tmp, hhmat_tmp, betavec_tmp);


#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "after QR factorization with col pivoting, factors are\nTmp = [\n";
        {
          for (int i = 0; i < tmrow; ++i) {
            for (int j = 0; j < tmcol; ++j) {
              *debug_str << Tmp(i,j) << " ";
            }
            *debug_str << '\n';
          }
        }
        *debug_str << "];\nselect_col_tmp = [";
        {
          for (int i = 0; i < tmrow; ++i)
            *debug_str << static_cast<int>(select_col_tmp[i]) << " ";
        }
        *debug_str << "];\nhhmat_tmp = [";
        {
          for (int i = 0; i < tmrow; ++i) {
            for (int j = 0; j < i + 1; ++j) {
              *debug_str << hhmat_tmp(i,j) << " ";
            }
            for (int j2 = i + 1; j2 < tmrow; ++j2) {
              *debug_str << "0 ";
            }
            *debug_str << '\n';
          }
        }
        *debug_str << "];\nbetavec_tmp = [";
        {
          for (int i = 0; i < tmrow - 1; ++i)
            *debug_str << betavec_tmp[i] << " ";
        }
        *debug_str << "];\n";
      }
#endif



      int smallsvpos = tmrow - 1;
      Real bigsv = fabs(Tmp(0,static_cast<int>(select_col_tmp[0])));

      while (fabs(Tmp(smallsvpos,static_cast<int>(select_col_tmp[smallsvpos])))
        < bigsv * pow(MACHINE_EPS, 0.33333)) {
        --smallsvpos;
      }

      if (smallsvpos == tmrow - 1) {
        break;
      }


      // Make up approximate null vectors.  Initialize columns
      // of the identity matrix.

      {
        for (int j = smallsvpos + 1; j < tmrow; ++j) {
          for (int i = 0; i < tmrow; ++i) {
            depmat(i,j) = 0.0;
          }
          depmat(j,j) = 1.0;
        }
      }

      // Apply HH transforms.

      {
        for (int k = smallsvpos + 1; k < tmrow; ++k) {
          for (int j = tmrow - 2; j >= 0; --j) {
            Real ip = 0.0;
            {
              for (int i = j; i < tmrow; ++i)
                ip += hhmat_tmp(i,j) * depmat(i,k);
            }
            ip *= betavec_tmp[j];
            {
              for (int i = j; i < tmrow; ++i)
                depmat(i,k) -= ip * hhmat_tmp(i,j);
            }
          }
        }
      }

      // Compute the multipliers and update expected_resd.

      {
        for (int i = smallsvpos + 1; i < tmrow; ++i) {
          Real sv = fabs(Tmp(i,static_cast<int>(select_col_tmp[i])));
          if (sv < MACHINE_EPS * bigsv)
            sv = MACHINE_EPS * bigsv;
          Real scaledsv = sv / bigsv;
          fac[i] = bigsv * sqrt(MACHINE_EPS  / scaledsv) * (1.0 + random_real());
          Real er = scfac * sqrt(MACHINE_EPS * scaledsv);
          if (er > expected_resd)
            expected_resd = er;
        }
      }

      coefs_were_perturbed = true;

      {
        for (int i = 0; i < tmrow; ++i) {
          int eqnum = static_cast<int>(1 - rowtype[i]);
          int ud = static_cast<int>(row_u_deg[i]);
          int vd = static_cast<int>(row_v_deg[i]);
          int cpnum = ud + vd * (udegree[eqnum] + 1); 
          for (int j = smallsvpos + 1; j < tmrow; ++j) 
            coef_perturb(eqnum, cpnum) += depmat(i,j) * fac[j];
        }
      }


#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "bigsv = " << bigsv << ";\ndepmat = [\n";
        {
          for (int i = 0; i < tmrow; ++i) {
            for (int j = smallsvpos + 1; j < tmrow; ++j) {
              *debug_str << depmat(i,j) << ' ';
            }
            *debug_str << '\n';
          }
        }
        *debug_str << "];\nsmallsv = [";
        {
          for (int j =  smallsvpos + 1; j < tmrow; ++j) {
            *debug_str << fabs(Tmp(j,static_cast<int>(select_col_tmp[j]))) << ' ';
          }
        }
        *debug_str << "];\nfac = [";
        {
          for (int j = smallsvpos + 1; j < tmrow; ++j) {
            *debug_str << fac[j] << " ";
          }
        }
        *debug_str << "];\ncoef_perturb = [\n";
        
        {
          for (int i = 0; i < num_control_pt0; ++i)
            *debug_str << coef_perturb(0,i) << ' ';
          *debug_str << '\n';
        }
        {
          for (int i = 0; i < num_control_pt1; ++i)
            *debug_str << coef_perturb(1,i) << ' ';
          *debug_str << '\n';
        }
        *debug_str << "];\n";
      }
#endif
    }
    
    
    PatchMath::Workspace::RealVector betavec(workspace, tmrow);
    PatchMath::Workspace::RealMatrix hhmat(workspace, tmcol, tmrow);

#ifdef DEBUGGING
    if (dump_everything) {
      *debug_str << "Before factoring, Tmi = [\n";
      for (int i = 0; i < tmcol; ++i) {
        for (int j = 0; j < tmrow; ++j) {
          *debug_str << Tm(i,j) << ' ';
        }
        *debug_str << '\n';
      }
      *debug_str << "];\n";
    }
#endif

    


    QR_Factor(Tm, hhmat, betavec);
    

#ifdef DEBUGGING    
    if (dump_everything) {
      {
        *debug_str << "Tmi = [\n";
        for (int i = 0; i < tmcol; ++i) {
          for (int j = 0; j < tmrow; ++j) { 
            *debug_str << Tm(i,j) << " ";
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\nhhmati = [\n";
      {
        for (int i = 0; i < tmcol; ++i) {
          int k = (i + 1 < tmrow)? (i + 1) : tmrow;
          for (int j = 0; j < k; ++j) {
            *debug_str << setprecision(17) << hhmat(i,j) << " ";
          }
          for (int j1 = k; j1 < tmrow; ++j1) {
            *debug_str << "0 ";
          }
          *debug_str << '\n';
        }
      }
      *debug_str << "];\n betavec = [\n";
      {
        for (int ii = 0; ii < tmrow; ++ii) {
          *debug_str << setprecision(17) << betavec[ii] << " ";
        }
      }
      *debug_str << "];\n";
    }
#endif
     

    // Now select which rows for the bottom of the matrix.

    int tmsmall_row = udegree[0] * vdegree[0] + udegree[1] * vdegree[1];
    int tmsmall_col = monomcount;
    PatchMath::Workspace::RealMatrix Tmsmall(workspace, tmsmall_row, tmsmall_col);
    {
      int rowind = 0;
      for (int eqnum = 0; eqnum < 2; ++eqnum) {
        for (int ud = 0; ud < udegree[1-eqnum]; ++ud) {
          for (int vd = 0; vd < vdegree[1-eqnum]; ++vd) {
            for (int j = 0; j < tmsmall_col; ++j)
              Tmsmall(rowind, j) = 0.0; 
            for (int ud1 = 0; ud1 <= udegree[eqnum]; ++ud1) {
              for (int vd1 = 0; vd1 <= vdegree[eqnum]; ++vd1) {
                int mtxcol = ud + ud1 + (vd + vd1) * umdegree;
                int cpnum = ud1 + vd1 * (udegree[eqnum] + 1);
                Real bz = best_shiftpoly(eqnum, cpnum);
                Real pt = coef_perturb(eqnum,cpnum);
                Tmsmall(rowind, mtxcol) = bz + pt;
              }
            }
            ++rowind;
          }
        }
      }
      
#ifdef DEBUGGING
      if (rowind != tmsmall_row) {
        *debug_str << " rowind = " << rowind << " tmrow = " << tmrow << '\n';
        throw_error("Counting error B");
      }
#endif
    }

#ifdef DEBUGGING
    if (dump_everything) {
      *debug_str << "Tmsmall = [\n";
      for (int i = 0; i < tmsmall_row; ++i) {
        for (int j = 0; j < tmsmall_col; ++j)
          *debug_str << " " << std::setprecision(17) << Tmsmall(i,j);
        *debug_str << "\n";
      }
      *debug_str << "]\n";
    }
#endif


    PatchMath::Workspace::RealMatrix hhmat_small(workspace, tmsmall_row, tmsmall_row);
    PatchMath::Workspace::RealVector betavec_small(workspace, tmsmall_row);
    PatchMath::Workspace::RealVector select_col(workspace, tmsmall_row);
    QR_Factor_with_col_pivoting(Tmsmall, select_col, hhmat_small, betavec_small);

#ifdef DEBUGGING
    if (dump_everything) {
      *debug_str << "selectcol = [";
      for (int i = 0; i < tmsmall_row; ++i) 
        *debug_str << static_cast<int>(select_col[i]) << " ";
      *debug_str << "]\n";
    }
#endif

        
    PatchMath::Workspace::RealVector use_monom(workspace, monomcount);      
    {
      for (int j = 0; j < monomcount; ++j)
        use_monom[j] = 1.0;
      for (int j2 = 0; j2 < tmsmall_row; ++j2)
        use_monom[static_cast<int>(select_col[j2])] = 0.0;
    }

    // grab space in the workspace.


    PatchMath::Workspace::RealVector expaArow(workspace, tmcol);
    PatchMath::Workspace::RealVector expaBrow(workspace, tmcol);

    PatchMath::Workspace::ComplexMatrix A(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix B(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix Ao(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix Bo(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix vr(workspace, evsize, evsize);
    PatchMath::Workspace::ComplexMatrix bestvr(workspace, evsize, evsize);
    PatchMath::Workspace::RealVector eveccond(workspace, evsize);
    PatchMath::Workspace::RealVector best_eveccond(workspace, evsize);
    
    const Real max_allowable_cond = 1.0e17;


    Real allowable_cond = 1.0e4;
    int lapack_failcount = 0;
    Real best_cond = -1.0;

    for (;;) {

      
      Real direc[2];
      
      // Choose some random direction.
      Real angle = random_real() * 2 * PI;
      direc[0] = cos(angle);
      direc[1] = sin(angle);

#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "direc = [";
        for (int i = 0; i < 2; ++i) {
          *debug_str << setprecision(17) << direc[i] << " ";
        }
        *debug_str << "];\n";
      }
#endif
      

      // From the shift and direction, compute upper and lower
      // bounds on eigenvalues that land in the patch.

      Real umin1, umax1, vmin1, vmax1;
      if (direc[0] > 0) {
        umin1 = -best_ushiftval;
        umax1 = 1.0/best_ushiftval;
      }
      else {
        umin1 = 1.0/best_ushiftval;
        umax1 = -best_ushiftval;
      }
      if (direc[1] > 0) {
        vmin1 = -best_vshiftval;
        vmax1 = 1.0/best_vshiftval;
      }
      else {
        vmin1 = 1.0/best_vshiftval;
        vmax1 = -best_vshiftval;
      }

      Real eigenval_lb = direc[0] * umin1 + direc[1] * vmin1;
      Real eigenval_ub = direc[0] * umax1 + direc[1] * vmax1;

     
      // Compute candidate directions for the bottom half.
      Real ear1 = direc[0];
      Real ear2 = direc[1];
      Real ebr1 = 1.0;
      int evr = 0;
      for (int vd = 0; vd < vmdegree; ++vd) {
        for (int ud = 0; ud < umdegree; ++ud) {
          int exrownum = ud + vd * umdegree;
          if (!use_monom[exrownum]) continue;

          for (int jj = 0; jj < tmcol; ++jj) {
            expaArow[jj] = 0.0;
            expaBrow[jj] = 0.0;
          }
          expaArow[ud + 1 + vd * (umdegree + 1)] = ear1;
          expaArow[ud + (vd + 1) * (umdegree + 1)] = ear2;
          expaBrow[ud + vd * (umdegree + 1)] = ebr1;

          
#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str <<  " expaArow(" << exrownum << ",:) = [\n";
            {
              for (int i = 0; i < tmcol; ++i) {
                *debug_str << setprecision(17) << expaArow[i] << " ";
              }
            }
            *debug_str <<  "\n expaBrow(" << exrownum << ",:) = [\n";
            {
              for (int i = 0; i < tmcol; ++i) {
                *debug_str << setprecision(17) << expaBrow[i] << " ";
              }
            }
            *debug_str << "];\n";
          }
#endif
              
          
          // Apply HH transforms 
          
          for (int j = 0; j < tmrow; ++j) {
            Real ip1 = 0.0;
            Real ip2 = 0.0;
            {
              for (int i = j; i < tmcol; ++i) {
                ip1 += hhmat(i,j) * expaArow[i];
                ip2 += hhmat(i,j) * expaBrow[i];
              }
            }
            ip1 *= betavec[j];
            ip2 *= betavec[j];
            {
              for (int i = j; i < tmcol; ++i) {
                expaArow[i] -= ip1 * hhmat(i,j);
                expaBrow[i] -= ip2 * hhmat(i,j);
              }
            }
          }
          
          for (int i = tmrow; i < tmcol; ++i) {
            A(evr, i - tmrow) = Complex(expaArow[i], 0.0);
            B(evr, i - tmrow) = Complex(expaBrow[i], 0.0);
          }       
          ++evr;
        }
      }
      
#ifdef DEBUGGING
      if (evr != evsize) {
        *debug_str << " evr = " << evr << " evsize = " << evsize << '\n';
        throw_error("counting error C");
      }
      if (dump_everything) {
        *debug_str << "A = [\n";
        {
          for (int i = 0; i < evsize; ++i) {
            for (int j = 0; j < evsize; ++j) {
              *debug_str << Format_Complex(A(i,j)) << " ";
            }
            *debug_str << '\n';
          } 
        }
        *debug_str << "];\nB = [\n";
        {
          for (int i = 0; i < evsize; ++i) {
            for (int j = 0; j < evsize; ++j) {
              *debug_str << Format_Complex(B(i,j)) << " ";
            }
            *debug_str << '\n';
          }
        }          
        *debug_str << "];\n";
      }
#endif

      // Copy selected rows into Ao and Bo.
      {
        for (int i = 0; i < evsize; ++i) {
          for (int j = 0; j < evsize; ++j) {
            Ao(i,j) = A(i,j);
            Bo(i,j) = B(i,j);
          }
        }
      }


      Real cond = generalized_eigensolver(A, B, vr, eveccond, eigenval_lb, 
        eigenval_ub, workspace);
      
  
      if (cond < 0) {
#ifdef DEBUGGING
        *debug_str << "Lapack failure in Solve_BezQuad::solve\n";
        *debug_str << "A = ";
        {
          for (int i = 0; i < evsize; ++i) {
            for (int j = 0; j < evsize; ++j) {
              *debug_str << Format_Complex(Ao(i,j)) << " ";
            }
            *debug_str << '\n';
          } 
        }
        *debug_str << "B = ";
        {
          for (int i = 0; i < evsize; ++i) {
            for (int j = 0; j < evsize; ++j) {
              *debug_str << Format_Complex(Bo(i,j)) << " ";
            }
            *debug_str << '\n';
          }
        }  
        *debug_str << "Afinal = ";
        {
          for (int i = 0; i < evsize; ++i) {
            for (int j = 0; j < evsize; ++j) {
              *debug_str << Format_Complex(A(i,j)) << " ";
            }
            *debug_str << '\n';
          } 
        }
        *debug_str << "Bfinal = ";
        {
          for (int i = 0; i < evsize; ++i) {
            for (int j = 0; j < evsize; ++j) {
              *debug_str << Format_Complex(B(i,j)) << " ";
            }
            *debug_str << '\n';
          }
        }
#endif
        if (++lapack_failcount == 5) {
          throw_error("Repeated lapack failure");
        }
        continue;
      }

#ifdef DEBUGGING
      if (dump_everything)
        *debug_str << "cond = " << cond 
        << "\nbest_cond = " << best_cond 
        << "\nallowable_cond = " << allowable_cond << "\n";
#endif

      if (best_cond < 0 || cond < best_cond) {
        best_cond = cond; 
        {
          for (int i = 0; i < evsize; ++i) {
            for (int j = 0; j < evsize; ++j) {
              bestvr(i,j) = vr(i,j);
            }
          }
        }
        {
          for (int i = 0; i < evsize; ++i) {
            best_eveccond[i] = eveccond[i];
          }
        }
      }
      if (cond < allowable_cond || allowable_cond > max_allowable_cond)
        break;

      allowable_cond *= 10.0;
    }

      
    // Find the Bezier solutions.

    PatchMath::Workspace::ComplexVector transfevec(workspace, tmcol);
    PatchMath::Workspace::ComplexMatrix jac(workspace, 2, 2);
    
    
    for (int vlrow = 0; vlrow < evsize; ++vlrow) {
      
      
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << " - row " << vlrow << '\n';
        *debug_str << " bestvr row = [";
        for (int i = 0; i < evsize; ++i) {
          *debug_str << Format_Complex(bestvr(i, vlrow)) << " ";
        }
        *debug_str << "]\n";
      }
#endif
      
      // Apply the HH transforms to get the transformed eigenvector.
      
      {
        for (int j = 0; j < tmrow; ++j)
          transfevec[j] = Complex(0.0,0.0);
      }
      
      {
        for (int j = tmrow; j < tmcol; ++j) 
          transfevec[j] = bestvr(j - tmrow, vlrow);
      }
      {
        for (int j = tmrow - 1; j >= 0; --j) {
          Complex ip(0.0,0.0);
          {
            for (int i = j; i < tmcol; ++i) {
              ip += Complex(hhmat(i,j),0.0) * transfevec[i];
            }
          }
          ip *= Complex(betavec[j], 0.0);
          {
            for (int i = j; i < tmcol; ++i) {
              transfevec[i] -= ip * Complex(hhmat(i,j), 0.0);
            }
          }
        }
      }
      
      // Find the largest entry in the eigenvector.
      
      
      Complex u1, oppu1, v1, oppv1;
      {
        int maxud, maxvd;
        Real mx = 0.0;
        {
          for (int ud = 0; ud <= umdegree; ++ud) {
            for (int vd = 0; vd <= vmdegree; ++vd) {
              if (ud == umdegree && vd == vmdegree) continue;
              Complex en1 = transfevec[ud + vd * (umdegree + 1)];
              if (abs(en1) >= mx) {
                maxud = ud;
                maxvd = vd;
                mx = abs(en1);
              }
            }
          }
        }
        Complex c1 = transfevec[maxud + maxvd * (umdegree + 1)];
        if (maxud == 0) {
          oppu1 = c1;
          u1 = transfevec[maxud + 1 + maxvd * (umdegree + 1)];
        }
        else {
          u1 = c1;
          oppu1 = transfevec[maxud - 1 + maxvd * (umdegree + 1)];
        }
        if (maxvd == 0) {
          oppv1 = c1;
          v1 = transfevec[maxud + (maxvd + 1) * (umdegree + 1)];
        }
        else {
          v1 = c1;
          oppv1 = transfevec[maxud + (maxvd - 1) * (umdegree + 1)];
        }
      }
      
      // Apply the transformation to get u2,oppu2, v2, oppv2.
      
      
      Complex ushiftval_c(best_ushiftval, 0.0);
      Complex vshiftval_c(best_vshiftval, 0.0);
      Complex u2 = u1 + ushiftval_c * oppu1;
      Complex oppu2 = -ushiftval_c * u1 + oppu1;
      Complex v2 = v1 + vshiftval_c * oppv1;
      Complex oppv2 = -vshiftval_c * v1 + oppv1;
      
      Complex udenom = u2 + oppu2;
      Complex vdenom = v2 + oppv2;
      
      if (abs(udenom) == 0 || abs(vdenom) == 0) continue;
      
      C2Point sol(u2 / udenom, v2 / vdenom);
      if (sol.l2norm() >  pow(tol, -.6666)) continue;
      
      
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "transfevec = [";
        for (int i = 0; i < tmcol; ++i)
          *debug_str << Format_Complex(transfevec[i]) << " ";
        *debug_str << "]\n u1 = " << Format_Complex(u1) 
          << " u2 = " << Format_Complex(u2) 
          << " \n v1 = " << Format_Complex(v1)
          << " v2 = " << Format_Complex(v2)
          << " \n udenom = " << Format_Complex(udenom)
          << " vdenom = " << Format_Complex(vdenom) << '\n';
      }
#endif
      
      C2Point fval, newsol, newfval,  nstep;
      Real invnormest;
      
      // Try to improve the solution with two steps of Newton's method.
      
      for (int newtcount = 0;; ++newtcount) {  
        fval.coord[0] = eval_complex(sol, 0, udegree[0], vdegree[0]);
        fval.coord[1] = eval_complex(sol, 1, udegree[1], vdegree[1]);
        jac(0,0) = eval_partialu_complex(sol, 0, udegree[0], vdegree[0]);
        jac(1,0) = eval_partialu_complex(sol, 1, udegree[1], vdegree[1]);
        jac(0,1) = eval_partialv_complex(sol, 0, udegree[0], vdegree[0]);
        jac(1,1) = eval_partialv_complex(sol, 1, udegree[1], vdegree[1]);
        invnormest = linear_solve_C2(jac, fval, nstep);
        
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "sol = [" 
            <<    Format_Complex(sol.coord[0]) << " " << Format_Complex(sol.coord[1])
            << "]\n fval = [" << Format_Complex(fval.coord[0]) << " " 
            << Format_Complex(fval.coord[1])
            << "]\n jac = [" 
            << Format_Complex(jac(0,0)) << "," 
            << Format_Complex(jac(0,1)) << ";" 
            << Format_Complex(jac(1,0)) << "," 
            << Format_Complex(jac(1,1));
          *debug_str << "]\n nstep = [" 
            << Format_Complex(nstep.coord[0]) << " " 
            << Format_Complex(nstep.coord[1])
            << "]\n invnormest = " << invnormest;
          *debug_str << '\n';
        }
        
#endif
        if (!newton_improve) break;
        if (newtcount == 2) break;
        if (invnormest >= BIG_REAL) break;
        
        newsol.coord[0] = sol.coord[0] - nstep.coord[0];
        newsol.coord[1] = sol.coord[1] - nstep.coord[1];
        newfval.coord[0] = eval_complex(newsol, 0, udegree[0], vdegree[0]);
        newfval.coord[1] = eval_complex(newsol, 1, udegree[1], vdegree[1]);
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << " newsol = [" << Format_Complex(newsol.coord[0]) << 
            "," << Format_Complex(newsol.coord[1])
            << "]\nnewfval = [" << Format_Complex(newfval.coord[0]) << "," << 
            Format_Complex(newfval.coord[1]) << "]\n";
        }
#endif
        if (newfval.l2norm() < fval.l2norm()) {
#ifdef DEBUGGING
          if (dump_everything) {
            *debug_str << " changing\n";
          }
#endif
          sol = newsol;
        }
        else {
          break;
        }
      }
      
      Real fvalcutoff = evsize * tol * scfac;
      Real fval2 = fval.l2norm();
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << " cutoff = " << fvalcutoff << " l2norm = " << fval2
          << '\n';
      }
#endif
      Real fvalv = fval2 + fvalcutoff;
      Real param_uncertainty = (invnormest >= BIG_REAL)? 1.0 : (fvalv * invnormest);
      
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << " passed first test, param_uncertainty = "
          << param_uncertainty <<  "\n";
      }
#endif
      if (abs(sol.coord[0]) + abs(sol.coord[1]) > 3) {
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << " dropping root; norm too large\n";
        }
#endif
        continue;
      }
      
      if (coefs_were_perturbed &&
        fval2 > evsize * expected_resd * best_eveccond[vlrow]) {
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << "residual " << fval2 << " exceeds expected resd\n"
            << evsize << " * " << expected_resd << " * " << best_eveccond[vlrow]
            << " = " << evsize * expected_resd * best_eveccond[vlrow]
            << "\ndiscarding root\n";
        }
#endif
        continue;
      }
      
      
      if (sol.coord[0].real() > -param_uncertainty && 
        sol.coord[1].real() > -param_uncertainty &&
        sol.coord[0].real() < 1.0 + param_uncertainty &&
        sol.coord[1].real() < 1.0 + param_uncertainty &&
        fabs(sol.coord[0].imag()) < param_uncertainty &&
        fabs(sol.coord[1].imag()) < param_uncertainty) {          
        PatchMath::IntersectRecord result;
        result.paramcoord[0] = sol.coord[0].real();
        result.paramcoord[1] = sol.coord[1].real();
        result.degen = (fabs(sol.coord[0].real()) <= param_uncertainty || 
          fabs(sol.coord[1].real()) <= param_uncertainty ||
          fabs(1.0 - sol.coord[0].real()) <= param_uncertainty ||
          fabs(1.0 - sol.coord[1].real()) <= param_uncertainty);
        result.param_uncertainty = param_uncertainty;
#ifdef DEBUGGING
        if (dump_everything) {
          *debug_str << " passed second test, degen = " << result.degen << "\n";
        }
#endif
        output_stack.push_back(result);
      }
    }
   
#ifdef DEBUGGING
    if (dump_everything)
      *debug_str << '\n';
#endif
    return is_degen;
  }
    
}


// Intersect a segment with a curve (2D) or patch (3D).

void QMG::PatchMath::intersect_seg(vector<IntersectRecord>& output_stack,
                                   Workspace& workspace,
                                   const Point& startpoint,
                                   const Point& endpoint,
                                   double scfac,
                                   double tol,
                                   bool newton_improve) const {
  
#ifdef DEBUGGING
  using Local_to_PatchMath_CPP::dump_everything;
  using Local_to_PatchMath_CPP::debug_str;
  if (dump_everything)
    *debug_str << " reached intersect_seg\n";
#endif

  if (gdim_ != di_ - 1)
    return;

  Point tangent = Point::subtract(endpoint,startpoint);
  Real lth = tangent.l2norm();
  if (lth == 0.0)
    throw_error("Zero-length segment in intersect_seg");


  // Convert the segment from implicit to parametric form.

  Real lineeqspace[MAXDIM * MAXDIM];
  Matrix lineeq(di_ - 1, di_, lineeqspace);

  if (di_ == 2) {
    lineeq(0,0) = tangent[1] / lth;
    lineeq(0,1) = -tangent[0] / lth;
  }
  else {
    Point hhtransform;
    Real beta = 
      Matrix::compute_hh_transform(&hhtransform[0], &tangent[0], 1, 3);
    for (int i = 0; i < 2; ++i) {
      Real mult = -beta * hhtransform[i + 1];
      for (int j = 0; j < 3; ++j) {
        lineeq(i,j) = mult * hhtransform[j];
      }
      lineeq(i, i + 1) += 1.0;
    }
  }

  Point linerhs;
  {
    for (int i = 0; i < di_ - 1; ++i) {
      Real t = 0.0;
      for (int j = 0; j < di_; ++j) {
        t += lineeq(i,j) * startpoint[j];
      }
      linerhs[i] = t;
    }
  }
  int init_stacksize = output_stack.size();
  intersect_line(output_stack, workspace, lineeq, linerhs, scfac, tol, newton_improve);
  int end_stacksize = output_stack.size();

  int output_pos = init_stacksize;
  int degree_fac;
  if (gdim_ == 1 || ptype_ == BEZIER_TRIANGLE) {
    degree_fac = (degree1_ + 1) * (degree1_ + 1) * 2;
  }
  else {
    degree_fac = (degree1_ + 1) * (degree2_ + 1) * 2;
  }



  Real ip0 = Point::inner_product(startpoint, tangent);
  Real ip1 = lth * lth;

  for (int input_pos = init_stacksize; input_pos < end_stacksize; ++input_pos) {
    Real ip2 = Point::inner_product(output_stack[input_pos].realcoord, tangent);
    Real t = (ip2 - ip0) / ip1;

    Real seg_unc = output_stack[input_pos].param_uncertainty * degree_fac *
      scfac / lth;
    
#ifdef DEBUGGING
    if (dump_everything) {
      *debug_str << "input_pos = " << input_pos << " output_pos = " << output_pos
                 << " ips = " << ip0 << " " << ip1 << " " << ip2
                 << "\n t = " << std::setprecision(17) << t <<  " segunc = " << seg_unc << '\n'
                 << "scfac = " << scfac << " lth = " << lth << " degree_fac ="
                 << degree_fac << '\n'
                 << "os uncertainty = " 
                 << output_stack[input_pos].param_uncertainty << '\n';
    }
#endif
    if (t > -seg_unc && t < 1.0 + seg_unc) {
      output_stack[output_pos] = output_stack[input_pos];
      output_stack[output_pos].dist = t;
      if (fabs(t) <= seg_unc || fabs(t - 1.0) <= seg_unc) {
        output_stack[output_pos].degen = true;
      }
#ifdef DEBUGGING
      if (dump_everything) {
        *debug_str << "keeping, degen = " << output_stack[output_pos].degen << '\n';
      }
#endif
      output_stack[output_pos].param_uncertainty_seg = seg_unc;
      ++output_pos;
    }
  }
  output_stack.resize(output_pos, IntersectRecord());
}


// Intersect a curve (in 2D) or patch (in 3D) with a line.


void QMG::PatchMath::intersect_line(vector<IntersectRecord>& output_stack,
                                    Workspace& workspace,
                                    const Matrix& lineeq,
                                    const Point& linerhs,
                                    double scfac,
                                    double tol,
                                    bool newton_improve) const {

  if (gdim_ != di_ - 1)
    return;

#ifdef DEBUGGING
  using Local_to_PatchMath_CPP::dump_everything;
  using Local_to_PatchMath_CPP::debug_str;
  if (dump_everything)
    *debug_str << " reached intersect_line\n";
#endif

  int initsize = output_stack.size();

  if (di_ == 2) {
    Solve_OneVar_Bezier1 helper(lineeq, linerhs, *this);
    helper.solve(degree1_, output_stack, workspace, scfac, tol);
  }
  else {
    if (ptype_ ==  BEZIER_TRIANGLE) {
      Solve_BezTri1 helper(lineeq, linerhs, *this);
      helper.solve(degree1_, output_stack, workspace, scfac, tol, newton_improve);
    }
    else {
      Solve_BezQuad1 helper(lineeq, linerhs, *this);
      int deg1[2];
      int deg2[2];
      deg1[0] = deg1[1] = degree1_;
      deg2[0] = deg2[1] = degree2_;
      helper.solve(deg1, deg2, output_stack, workspace, scfac, tol, newton_improve);
    }
  }
  for (int j = initsize; j < output_stack.size(); ++j) {
    output_stack[j].realcoord = 
      get_real_point(output_stack[j].paramcoord);
#ifdef DEBUGGING
    if (dump_everything)
      *debug_str << " real point = " << output_stack[j].realcoord << '\n';
#endif
    
  }
  return;
}


// Intersect a plane with a bezier curve.


void QMG::PatchMath::intersect_plane(vector<IntersectRecord>& output_stack,
                                     Workspace& workspace,
                                     const Point& plane_normal,
                                     Real rhs,
                                     double scfac,
                                     double tol,
                                     bool newton_improve) const {

#ifdef DEBUGGING
  if (di_ != 3 || gdim_ != 1)
    throw_error("Intersect_plane called on wrong dimensions");
#endif

  int initsize = output_stack.size();
  Solve_OneVar_Bezier2 helper(plane_normal, rhs, *this);
  helper.solve(degree1_, output_stack, workspace, scfac, tol);
  for (int j = initsize; j < output_stack.size(); ++j) {
    output_stack[j].realcoord = 
      get_real_point(output_stack[j].paramcoord);
  }
}

// Get extreme points of a Bezier patch or curve.
// This means extreme in a certain direc direc.
// Do this by finding roots of the directional derivative.


void QMG::PatchMath::get_extreme_points(vector<IntersectRecord>& output_stack,
                                        Workspace& workspace,
                                        const Point& direc,
                                        double scfac,
                                        double tol) const {
  int initsize = output_stack.size();
  if (gdim_ == 0) {
    IntersectRecord result;
    for (int j = 0; j < di_; ++j) 
      result.realcoord[j] = control_point_coord_(0,j);
    result.degen = false;
    result.param_uncertainty = tol;
    result.param_uncertainty_seg = 0.0;
    result.dist = 0.0;
    output_stack.push_back(result);
  }
  else {
    if (degree1_ == 1) return;
    if (gdim_ == 1) {
      Solve_OneVar_Bezier4 helper(direc, di_, *this);
      helper.solve(degree1_ - 1, output_stack, workspace, scfac, tol);
    }
    else { // get here if gdim == 2, di = 3.
      if (ptype_ == BEZIER_TRIANGLE) {
        Solve_BezTri2 helper(direc, di_, *this);
        helper.solve(degree1_ - 1, output_stack, workspace, scfac, tol, true);
      }
      else {
        if (degree2_ == 1) return;
        Solve_BezQuad2 helper(direc, di_, *this, degree1_, degree2_);
        int deg1[2];
        int deg2[2];
        deg1[0] = degree1_ - 1;
        deg2[0] = degree2_;
        deg1[1] = degree1_;
        deg2[1] = degree2_ - 1;
        helper.solve(deg1, deg2, output_stack, workspace, scfac, tol, true);
      }
    }
    for (int k = initsize; k < output_stack.size(); ++k) {
      output_stack[k].realcoord = 
        get_real_point(output_stack[k].paramcoord);
    }
  }
}
      

QMG::Point 
QMG::PatchMath::get_real_point(const Point& param_coord) const {
  Point real_coord;

  if (gdim_ == 0) {
    for (int j = 0; j < di_; ++j)
      real_coord[j] = control_point_coord_(0,j);  
  }
  else if (gdim_ == 1) {
    for (int k = 0; k < di_; ++k) {
      Solve_OneVar_Bezier3 e1(k, *this);
      real_coord[k] = e1.eval(param_coord[0], degree1_);
    }
  }
  else  {//  gdim == 2

    if (ptype_ == BEZIER_TRIANGLE) {
      
      // Use deCasteljau algorithm for a Bezier triangle

      for (int k = 0; k < di_; ++k) {
        Solve_BezTri3 e1(k, *this);
        real_coord[k] = e1.eval(R2Point(param_coord[0], param_coord[1]),
          0, degree1_);
      }
    }
    else {
      for (int k = 0; k < di_; ++k) {
        Solve_BezQuad3 e1(k, *this);
        real_coord[k] = e1.eval(R2Point(param_coord[0], param_coord[1]),
          0, degree1_, degree2_);
      }
    }
  }
  return real_coord;
}

QMG::Point 
QMG::PatchMath::direc_deriv(const Point& param_coord, 
                            const Point& direc) const {

  Point real_coord;

  if (gdim_ == 0) {
    throw_error("Directional derivative of 0D-patch not allowed.");
  }
  else if (gdim_ == 1) {
    for (int k = 0; k < di_; ++k) {
      Solve_OneVar_Bezier3 e1(k, *this);
      real_coord[k] = direc[0] * e1.eval_deriv(param_coord[0], 
        degree1_);
    }
  }
  else if (gdim_ == 2) {

    if (ptype_ == BEZIER_TRIANGLE) {
      
      for (int k = 0; k < di_; ++k) {
        Solve_BezTri3 e1(k, *this);
        real_coord[k] = e1.eval_direc_deriv(R2Point(param_coord[0], param_coord[1]),
          R2Point(direc[0], direc[1]), 0, degree1_);
      }
    }
    else {
      Real partial;
      for (int k = 0; k < di_; ++k) {
        real_coord[k] = 0.0;
        Solve_BezQuad3 e1(k, *this);
        if (direc[0] != 0) {
          partial = e1.eval_partialu(R2Point(param_coord[0], param_coord[1]),
            0, degree1_, degree2_);
          real_coord[k] += partial * direc[0];
        }
        if (direc[1] != 0) {
          partial = e1.eval_partialv(R2Point(param_coord[0], param_coord[1]),
            0, degree1_, degree2_);
          real_coord[k] += partial * direc[1];
        }
      }
    }
  }
  return real_coord;
}


QMG::Point 
QMG::PatchMath::normal(const Point& paramcoord) const {
  if (gdim_ != di_ - 1)
    throw_error("Normal called for non-full-dimensional patch");
  Point direc;
  Point returnval;
  if (gdim_ == 1) {
    direc[0] = 1;
    Point tan = direc_deriv(paramcoord, direc);
    returnval[0] = tan[1];
    returnval[1] = -tan[0];
  }
  else {
    direc[0] = 1;
    direc[1] = 0;
    Point tan1 = direc_deriv(paramcoord, direc);
    direc[0] = 0;
    direc[1] = 1;
    Point tan2 = direc_deriv(paramcoord, direc);
    returnval = Point::cross_product(tan1, tan2);
  }
  if (returnval.l2norm() == 0)
    throw_error("zero normal");
  returnval.normalize();
  if (!is_forward_) {
    for (int i = 0; i < di_; ++i) 
      returnval[i] = -returnval[i];
  }
  return returnval;
}

QMG::Point QMG::PatchMath::normal_at_barycenter() const {
  return normal(barycenter_param());
}

QMG::Point QMG::PatchMath::barycenter_param() const {
  Point returnval;
  if (gdim_ == 1)
    returnval[0] = 0.5;
  else if (gdim_ == 2) {
    if (ptype_ == BEZIER_TRIANGLE) {
      returnval[0] = 0.33333333333333333;
      returnval[1] = 0.33333333333333333;
    }
    else {
      returnval[0] = 0.5;
      returnval[1] = 0.5;
    }
  }
  return returnval;
}

 
int QMG::PatchMath::num_control_points() const { 
  if (gdim_ == 0) {
    return 1;
  }
  else if (gdim_ == 1) {
    return degree1_ + 1;
  }
  else if (ptype_ == BEZIER_TRIANGLE) {
    return (degree1_ + 1) * (degree1_ + 2) / 2;
  }
  else {
    return (degree1_ + 1) * (degree2_ + 1);
  }
}