/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  src/1DSLF/1DSLF.cc

Purpose:  defines functions for 1d spinless fermion classes.

***********************************************************/

#include "ABACUS.h"

using namespace std;

namespace ABACUS {

  //***************************************************************************************************

  // Function definitions:  class ODSLF_Base

  ODSLF_Base::ODSLF_Base () : Mdown(0), Nrap(Vect<int>()), Nraptot(0), Ix2_infty(Vect<DP>()), Ix2_max(Vect<int>()) {}

  ODSLF_Base::ODSLF_Base (const ODSLF_Base& RefBase)  // copy constructor
    : Mdown(RefBase.Mdown), Nrap(Vect<int>(RefBase.Nrap.size())), Nraptot(RefBase.Nraptot),
      Ix2_infty(Vect<DP>(RefBase.Nrap.size())), Ix2_max(Vect<int>(RefBase.Nrap.size())),
      id(RefBase.id)
  {
    for (int i = 0; i < Nrap.size(); ++i) {
      Nrap[i] = RefBase.Nrap[i];
      Ix2_infty[i] = RefBase.Ix2_infty[i];
      Ix2_max[i] = RefBase.Ix2_max[i];
    }
  }

  ODSLF_Base::ODSLF_Base (const Heis_Chain& RefChain, int M)
    : Mdown(M), Nrap(Vect<int>(RefChain.Nstrings)), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
      Ix2_max(Vect<int>(RefChain.Nstrings))
  {
    for (int i = 0; i < RefChain.Nstrings; ++i) Nrap[i] = 0;
    Nrap[0] = M;

    Nraptot = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];

    // The id of this is zero by definition
    id = 0LL;

    // Now compute the Ix2_infty numbers

    (*this).Compute_Ix2_limits(RefChain);

  }

  ODSLF_Base::ODSLF_Base (const Heis_Chain& RefChain, const Vect<int>& Nrapidities)
    : Mdown(0), Nrap(Nrapidities), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
      Ix2_max(Vect<int>(RefChain.Nstrings)), id (0LL)
  {

    // Check consistency of Nrapidities vector with RefChain

    if (RefChain.Nstrings != Nrapidities.size())
      ABACUSerror("Incompatible Nrapidities vector used in ODSLF_Base constructor.");

    int Mcheck = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += RefChain.Str_L[i] * Nrap[i];
    Mdown = Mcheck;

    Nraptot = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];

    // Compute id
    id += Nrapidities[0];
    long long int factor = 100000LL;
    for (int i = 1; i < RefChain.Nstrings; ++i) {
      id += factor * Nrapidities[i];
      factor *= 100LL;
    }

    // Now compute the Ix2_infty numbers
    (*this).Compute_Ix2_limits(RefChain);

  }

  ODSLF_Base::ODSLF_Base (const Heis_Chain& RefChain, long long int id_ref)
    : Mdown(0), Nrap(Vect<int>(RefChain.Nstrings)), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
      Ix2_max(Vect<int>(RefChain.Nstrings)), id (id_ref)
  {
    // Build Nrapidities vector from id_ref

    long long int factor = pow_ulli (10LL, 2* RefChain.Nstrings + 1);
    long long int id_eff = id_ref;
    for (int i = 0; i < RefChain.Nstrings - 1; ++i) {
      Nrap[RefChain.Nstrings - 1 - i] = id_eff/factor;
      id_eff -= factor * Nrap[RefChain.Nstrings - 1 - i];
      factor /= 100LL;
    }
    Nrap[0] = id_eff;

    int Mcheck = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += RefChain.Str_L[i] * Nrap[i];
    Mdown = Mcheck;

    Nraptot = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];

    // Now compute the Ix2_infty numbers
    (*this).Compute_Ix2_limits(RefChain);
  }


  ODSLF_Base& ODSLF_Base::operator= (const ODSLF_Base& RefBase)
  {
    if (this != & RefBase) {
      Mdown = RefBase.Mdown;
      Nrap = RefBase.Nrap;
      Nraptot = RefBase.Nraptot;
      Ix2_infty = RefBase.Ix2_infty;
      Ix2_max = RefBase.Ix2_max;
      id = RefBase.id;
    }
    return(*this);
  }

  bool ODSLF_Base::operator== (const ODSLF_Base& RefBase)
  {
    bool answer = (Nrap == RefBase.Nrap);

    return (answer);
  }

  bool ODSLF_Base::operator!= (const ODSLF_Base& RefBase)
  {
    bool answer = (Nrap != RefBase.Nrap);

    return (answer);
  }

  void ODSLF_Base::Compute_Ix2_limits (const Heis_Chain& RefChain)
  {

    if ((RefChain.Delta > 0.0) && (RefChain.Delta < 1.0)) {

      // Compute the Ix2_infty numbers

      DP sum1 = 0.0;
      DP sum2 = 0.0;

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	sum1 = 0.0;

	for (int k = 0; k < RefChain.Nstrings; ++k) {

	  sum2 = 0.0;

	  sum2 += (RefChain.Str_L[j] == RefChain.Str_L[k]) ? 0.0 :
	    2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
			   - 0.5 * fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) * RefChain.anis));
	  sum2 += 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
				 - 0.5 * (RefChain.Str_L[j] + RefChain.Str_L[k]) * RefChain.anis));

	  for (int a = 1; a < ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]); ++a)
	    sum2 += 2.0 * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
					 - 0.5 * (fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) + 2.0*a) * RefChain.anis));

	  sum1 += (Nrap[k] - ((j == k) ? 1 : 0)) * sum2;
	}

	Ix2_infty[j] = (1.0/PI) * fabs(RefChain.Nsites * 2.0
				       * atan(tan(0.25 * PI * (1.0 + RefChain.par[j])
						  - 0.5 * RefChain.Str_L[j] * RefChain.anis)) - sum1);

      }  // The Ix2_infty are now set.

      // Now compute the Ix2_max limits

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	Ix2_max[j] = int(floor(Ix2_infty[j]));  // sets basic integer

	// Reject formally infinite rapidities (i.e. if Delta is root of unity)


	// For ODSLF:  parity of quantum numbers is complicated:
	// if (N-1 + M_j + (M + 1) n_j + (1 - nu_j) N/2) is even, q# are integers,
	// so Ix2_max must be even (conversely, if odd, then odd).
	if ((Ix2_max[j] + RefChain.Nsites - 1 + Nrap[j] + (Nraptot + 1) * RefChain.Str_L[j]
	     + (RefChain.par[j] == 1 ? 0 : RefChain.Nsites)) % 2) Ix2_max[j] -= 1;

	while (Ix2_max[j] > RefChain.Nsites) {
	  Ix2_max[j] -= 2;
	}
      }
    } // if XXZ gapless

    else if (RefChain.Delta == 1.0) {

      // Compute the Ix2_infty numbers

      int sum1 = 0;

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	sum1 = 0;

	for (int k = 0; k < RefChain.Nstrings; ++k) {

	  sum1 += Nrap[k] * (2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]) - ((j == k) ? 1 : 0));
	}

	Ix2_infty[j] = (RefChain.Nsites + 1.0 - sum1); // to get counting right...

      }  // The Ix2_infty are now set.

      // Now compute the Ix2_max limits

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	Ix2_max[j] = int(floor(Ix2_infty[j]));  // sets basic integer

	// Give the correct parity to Ix2_max

	// If Nrap is even, Ix2_max must be odd.  If odd, then even.

	//if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;   DEACTIVATED FOR ODSLF !

	// If Ix2_max equals Ix2_infty, we reduce it by 2:

	if (Ix2_max[j] == int(Ix2_infty[j])) Ix2_max[j] -= 2;

	while (Ix2_max[j] > RefChain.Nsites) {
	  Ix2_max[j] -= 2;
	}
      }

    } // if XXX AFM

    else if (RefChain.Delta > 1.0) {

      // Compute the Ix2_infty numbers

      int sum1 = 0;

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	sum1 = 0;

	for (int k = 0; k < RefChain.Nstrings; ++k) {

	  sum1 += Nrap[k] * (2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]) - ((j == k) ? 1 : 0));
	}

	Ix2_infty[j] = (RefChain.Nsites - 1 + 2 * RefChain.Str_L[j] - sum1);

      }  // The Ix2_infty are now set.

      // Now compute the Ix2_max limits

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	Ix2_max[j] = int(floor(Ix2_infty[j]));  // sets basic integer

	// Give the correct parity to Ix2_max

	// If Nrap is even, Ix2_max must be odd.  If odd, then even.

	//if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;     DEACTIVATED FOR ODSLF !

	// If Ix2_max equals Ix2_infty, we reduce it by 2:

	//if (Ix2_max[j] == Ix2_infty[j]) Ix2_max[j] -= 2;

	while (Ix2_max[j] > RefChain.Nsites) {
	  Ix2_max[j] -= 2;
	}

      }

    } // if XXZ_gpd

  }

  void ODSLF_Base::Scan_for_Possible_Types (Vect<long long int>& possible_type_id, int& nfound,
					    int base_level, Vect<int>& Nexcitations)
  {
    if (base_level == 0) {
      // We set all sectors 0, 1, 2, 3

      int Nexc_max_0 = (Nrap[0] + 1)/2;  // holes on R
      int Nexc_max_1 = Nrap[0]/2;
      int Nexc_max_2 = (Ix2_max[0] - Nrap[0] + 1)/2;
      int Nexc_max_3 = Nexc_max_2;

      // Count the number of excitations in higher levels:
      int Nexc_above = 0;
      for (int i = 4; i < Nexcitations.size(); ++i) Nexc_above += Nexcitations[i];

      for (int nphpairs = 0; nphpairs <= NEXC_MAX_HEIS/2 - Nexc_above; ++nphpairs)
	for (int nexc0 = 0; nexc0 <= ABACUS::min(Nexc_max_0, nphpairs); ++nexc0) {
	  int nexc1 = nphpairs - nexc0;
	  for (int nexc2 = 0; nexc2 <= ABACUS::min(Nexc_max_2, nphpairs); ++nexc2) {
	    int nexc3 = nphpairs - nexc2;

	    if (nexc1 >= 0 && nexc1 <= Nexc_max_1 && nexc3 >= 0 && nexc3 <= Nexc_max_3) {
	      // acceptable type found
	      Nexcitations[0] = nexc0; Nexcitations[1] = nexc1; Nexcitations[2] = nexc2; Nexcitations[3] = nexc3;
	      possible_type_id[nfound] = Ix2_Offsets_type_id (Nexcitations);
	      nfound++;
	    }
	  }
	}
    }

    else { // base_level > 0
      // We scan over the possible distributions of the excitations on the L and R sides:

      // Count the number of slots available on R:
      int Nexc_R_max = (Ix2_max[base_level] + 2)/2;
      int Nexc_L_max = (Ix2_max[base_level] + 1)/2;
      for (int Nexcleft = 0; Nexcleft <= Nrap[base_level]; ++Nexcleft)
	if (Nexcleft <= Nexc_L_max && Nrap[base_level] - Nexcleft <= Nexc_R_max) {
	  Nexcitations[2*base_level + 2] = Nrap[base_level] - Nexcleft;
	  Nexcitations[2*base_level + 3] = Nexcleft;
	  Scan_for_Possible_Types (possible_type_id, nfound, base_level - 1, Nexcitations);
	}
    }

  }

  Vect<long long int> ODSLF_Base::Possible_Types ()
  {
    // Given a base, this returns a vector of possible types

    Vect<long long int> possible_type_id (0LL, 1000);
    int nfound = 0;

    Vect<int> Nexcitations (0, 2*Nrap.size() + 2);

    Scan_for_Possible_Types (possible_type_id, nfound, Nrap.size() - 1, Nexcitations);

    // Copy results into new vector:
    Vect<long long int> possible_type_id_found (nfound);
    for (int i = 0; i < nfound; ++i)
      possible_type_id_found[i] = possible_type_id[i];

    return(possible_type_id_found);
  }


  //***************************************************************************************************

  // Function definitions:  class Ix2_Config

  ODSLF_Ix2_Config::ODSLF_Ix2_Config () {}

  ODSLF_Ix2_Config::ODSLF_Ix2_Config (const Heis_Chain& RefChain, int M)
    : Nstrings(1), Nrap(Vect<int>(M,1)), Nraptot(M), Ix2(new int*[1]) // single type of string here
  {
    Ix2[0] = new int[M];

    for (int j = 0; j < M; ++j) {
      Ix2[0][j] = -(M - 1) + 2*j;
      // Use simplification:  Nrap[0] = M = Nraptot, Str_L[0] = 1, par = 1, so
      if ((Ix2[0][j] + RefChain.Nsites) % 2) Ix2[0][j] -= 1;
    }
  }

  ODSLF_Ix2_Config::ODSLF_Ix2_Config (const Heis_Chain& RefChain, const ODSLF_Base& base)
    : Nstrings(RefChain.Nstrings), Nrap(base.Nrap), Nraptot(base.Nraptot), Ix2(new int*[RefChain.Nstrings])
  {
    Ix2[0] = new int[base.Nraptot];
    for (int i = 1; i < RefChain.Nstrings; ++i) Ix2[i] = Ix2[i-1] + base[i-1];

    // And now for the other string types...
    for (int i = 0; i < RefChain.Nstrings; ++i) {
      for (int j = 0; j < base[i]; ++j) Ix2[i][j] = -(base[i] - 1) + 2*j + 1 - (base[i] % 2);
    }

    // ODSLF:  put back correct parity of quantum nr:
    for (int i = 0; i < RefChain.Nstrings; ++i)
      for (int j = 0; j < base[i]; ++j)
	if ((Ix2[i][j] + RefChain.Nsites - 1 + base.Nrap[i] + (base.Nraptot + 1) * RefChain.Str_L[i] + (RefChain.par[i] == 1 ? 0 : RefChain.Nsites)) % 2)
	  Ix2[i][j] -= 1;
  }

  ODSLF_Ix2_Config& ODSLF_Ix2_Config::operator= (const ODSLF_Ix2_Config& RefConfig)
  {
    if (this != &RefConfig) {
      Nstrings = RefConfig.Nstrings;
      Nrap = RefConfig.Nrap;
      Nraptot = RefConfig.Nraptot;
      if (Ix2 != 0) {
	delete[] (Ix2[0]);
	delete[] (Ix2);
      }
      Ix2 = new int*[Nstrings];
      Ix2[0] = new int[Nraptot];
      for (int i = 1; i < Nstrings; ++i) Ix2[i] = Ix2[i-1] + Nrap[i-1];
      for (int i = 0; i < Nstrings; ++i)
	for (int j = 0; j < Nrap[i]; ++j) Ix2[i][j] = RefConfig.Ix2[i][j];

    }
    return(*this);
  }

  ODSLF_Ix2_Config::~ODSLF_Ix2_Config()
  {
    if (Ix2 != 0) {
      delete[] (Ix2[0]);
      delete[] Ix2;
    }
  }

  std::ostream& operator<< (std::ostream& s, const ODSLF_Ix2_Config& RefConfig)
  {
    s << &RefConfig << "\t" << RefConfig.Nstrings << "\t" << RefConfig.Nraptot << endl;
    for (int i = 0; i < RefConfig.Nstrings; ++i) s << RefConfig.Nrap[i] << "\t";
    s << endl;
    for (int i = 0; i < RefConfig.Nstrings; ++i) {
      s << "i " << i << "\t";
      for (int j = 0; j < RefConfig.Nrap[i]; ++j) s << RefConfig.Ix2[i][j] << "\t";
      s << endl;
    }
    return(s);
  }

  //***************************************************************************************************

  // Function definitions:  class Lambda

  ODSLF_Lambda::ODSLF_Lambda () {}

  ODSLF_Lambda::ODSLF_Lambda (const Heis_Chain& RefChain, int M)
    : Nstrings(1), Nrap(Vect<int>(M,1)), Nraptot(M), lambda(new DP*[1]) // single type of string here
  {
    lambda[0] = new DP[M];

    for (int j = 0; j < M; ++j) lambda[0][j] = 0.0;
  }

  ODSLF_Lambda::ODSLF_Lambda (const Heis_Chain& RefChain, const ODSLF_Base& base)
    : Nstrings(RefChain.Nstrings), Nrap(base.Nrap), Nraptot(base.Nraptot), lambda(new DP*[RefChain.Nstrings])
  {
    lambda[0] = new DP[base.Nraptot];
    for (int i = 1; i < RefChain.Nstrings; ++i) lambda[i] = lambda[i-1] + base[i-1];

    for (int i = 0; i < RefChain.Nstrings; ++i) {
      for (int j = 0; j < base[i]; ++j) lambda[i][j] = 0.0;
    }
  }

  ODSLF_Lambda& ODSLF_Lambda::operator= (const ODSLF_Lambda& RefLambda)
  {
    if (this != &RefLambda) {
      Nstrings = RefLambda.Nstrings;
      Nrap = RefLambda.Nrap;
      Nraptot = RefLambda.Nraptot;
      if (lambda != 0) {
	delete[] (lambda[0]);
	delete[] (lambda);
      }
      lambda = new DP*[Nstrings];
      lambda[0] = new DP[Nraptot];
      for (int i = 1; i < Nstrings; ++i) lambda[i] = lambda[i-1] + Nrap[i-1];
      for (int i = 0; i < Nstrings; ++i)
	for (int j = 0; j < Nrap[i]; ++j) lambda[i][j] = RefLambda.lambda[i][j];

    }
    return(*this);
  }

  ODSLF_Lambda::~ODSLF_Lambda()
  {
    if (lambda != 0) {
      delete[] (lambda[0]);
      delete[] lambda;
    }
  }

  //***************************************************************************************************

  //  Function definitions:  class ODSLF_Ix2_Offsets

  ODSLF_Ix2_Offsets::ODSLF_Ix2_Offsets () : base(), Tableau(Vect<Young_Tableau>()), type_id(0LL), id(0LL), maxid(0LL) {};

  ODSLF_Ix2_Offsets::ODSLF_Ix2_Offsets (const ODSLF_Ix2_Offsets& RefOffset) // copy constructor
    : base(RefOffset.base), Tableau(Vect<Young_Tableau> (2 * base.Nrap.size() + 2)), type_id(RefOffset.type_id), id(RefOffset.id), maxid(RefOffset.maxid)
  {
    for (int i = 0; i < 2 * base.Nrap.size() + 2; ++i) Tableau[i] = RefOffset.Tableau[i];
  }

  ODSLF_Ix2_Offsets::ODSLF_Ix2_Offsets (const ODSLF_Base& RefBase, long long int req_type_id)
  // sets all tableaux to empty ones, with nparticles(req_type_id) at each level
  {
    // Build nparticles vector from req_type_id

    Vect<int> nparticles(0, 2* RefBase.Nrap.size() + 2);
    long long int factor = pow_ulli (10LL, nparticles.size() - 1);
    long long int id_eff = req_type_id;
    for (int i = 0; i < nparticles.size(); ++i) {
      nparticles[nparticles.size() - 1 - i] = id_eff/factor;
      id_eff -= factor * nparticles[nparticles.size() - 1 - i];
      factor /= 10LL;
    }

    // Check if we've got the right vector...
    long long int idcheck = Ix2_Offsets_type_id (nparticles);
    if (idcheck != req_type_id) ABACUSerror("idcheck != req_type_id in Ix2_Offsets constructor.");

    (*this) = ODSLF_Ix2_Offsets(RefBase, nparticles);
  }

  ODSLF_Ix2_Offsets::ODSLF_Ix2_Offsets (const ODSLF_Base& RefBase, Vect<int> nparticles)
    // sets all tableaux to empty ones, with nparticles at each level
    : base(RefBase), Tableau(Vect<Young_Tableau> (2 * base.Nrap.size() + 2)),
      type_id(Ix2_Offsets_type_id (nparticles)), id(0LL), maxid(0LL)
  {
    // Checks on nparticles vector:

    if (nparticles.size() != 2 * base.Nrap.size() + 2)
      ABACUSerror("Wrong nparticles.size in Ix2_Offsets constructor.");

    if (nparticles[3] + nparticles[2] != nparticles[0] + nparticles[1]) {
      cout << nparticles[0] << "\t" << nparticles[1] << "\t" << nparticles[2] << "\t" << nparticles[3] << endl;
      ABACUSerror("Wrong Npar[0-3] in Ix2_Offsets constructor.");
    }

    for (int base_level = 1; base_level < base.Nrap.size(); ++ base_level)
      if (base.Nrap[base_level] != nparticles[2*base_level + 2] + nparticles[2*base_level + 3]) {
	cout << base_level << "\t" << base.Nrap[base_level] << "\t" << nparticles[2*base_level + 2]
	     << "\t" <<  nparticles[2*base_level + 3] << endl;
	ABACUSerror("Wrong Nrap[] in Ix2_Offsets constructor.");
      }

    // nparticles[0,1]:  number of holes on R and L side in GS interval
    if (nparticles[0] > (base.Nrap[0] + 1)/2)
      ABACUSerror("nparticles[0] too large in Ix2_Offsets constructor.");
    if (nparticles[1] > base.Nrap[0]/2)
      ABACUSerror("nparticles[1] too large in Ix2_Offsets constructor.");

    // nparticles[2,3]:  number of particles of type 0 on R and L side out of GS interval
    if (nparticles[2] > (base.Ix2_max[0] - base.Nrap[0] + 1)/2)
      ABACUSerror("nparticles[2] too large in Ix2_Offsets constructor.");
    if (nparticles[3] > (base.Ix2_max[0] - base.Nrap[0] + 1)/2)
      ABACUSerror("nparticles[3] too large in Ix2_Offsets constructor.");

    for (int base_level = 1; base_level < base.Nrap.size(); ++ base_level)
      if ((nparticles[2*base_level + 2] > 0 && nparticles[2*base_level + 2]
	   > (base.Ix2_max[base_level] + 2)/2) // ODSLF modif
	  || (nparticles[2*base_level + 3] > 0
	      && nparticles[2*base_level + 3] > (base.Ix2_max[base_level] + 1)/2)) // ODSLF modif
	{
	  cout << base_level << "\t" << base.Ix2_max[base_level] << "\t" << base.Ix2_infty[base_level]
	       << "\t" << nparticles[2*base_level + 2] << "\t"
	       << (base.Ix2_max[base_level] - ((base.Nrap[base_level] + 1) % 2) + 2)/2
	       << "\t" << nparticles[2*base_level + 3] << "\t"
	       << (base.Ix2_max[base_level] - (base.Nrap[base_level] % 2) - 1)/2
	       << "\t" << (nparticles[2*base_level + 2] > 0) << "\t"
	       << (nparticles[2*base_level + 2] > (base.Ix2_max[base_level] - ((base.Nrap[base_level] + 1) % 2) + 2)/2)
	       << "\t" << (nparticles[2*base_level + 3] > 0) << "\t"
	       << (nparticles[2*base_level + 3] > base.Ix2_max[base_level] + 1
		   - (base.Ix2_max[base_level] - ((base.Nrap[base_level] + 1) % 2) + 2)/2)
	       << endl;
	  ABACUSerror("nparticles too large in Ix2_Offsets constructor.");
	}

    // Check sum of rapidities

    // Holes in GS interval
    Tableau[0] = Young_Tableau(nparticles[0], (base.Nrap[0] + 1)/2 - nparticles[0]);
    Tableau[1] = Young_Tableau(nparticles[1], base.Nrap[0]/2 - nparticles[1], Tableau[0]);

    // Particles of type 0 out of GS interval
    Tableau[2] = Young_Tableau(nparticles[2], (base.Ix2_max[0] - base.Nrap[0] + 1)/2 - nparticles[2], Tableau[0]);
    Tableau[3] = Young_Tableau(nparticles[3], (base.Ix2_max[0] - base.Nrap[0] + 1)/2 - nparticles[3], Tableau[2]);

    // Tableaux of index i = 2,...:  data about string type i/2-1.
    for (int base_level = 1; base_level < base.Nrap.size(); ++base_level) {
      Tableau[2*base_level + 2] =
	Young_Tableau(nparticles[2*base_level + 2],
		      (base.Ix2_max[base_level] - ((base.Nrap[base_level] + 1) % 2))/2 + 1
		      - nparticles[2*base_level + 2], Tableau[2]);
      Tableau[2*base_level + 3] =
	Young_Tableau(nparticles[2*base_level + 3],
		      (base.Ix2_max[base_level] - (base.Nrap[base_level] % 2) - 1)/2 + 1
		      - nparticles[2*base_level + 3], Tableau[3]);
    }

    maxid = 1LL;
    for (int i = 0; i < nparticles.size(); ++i) {
      maxid *= Tableau[i].maxid + 1LL;
    }
    maxid -= 1LL;

  }

  ODSLF_Ix2_Offsets& ODSLF_Ix2_Offsets::operator= (const ODSLF_Ix2_Offsets& RefOffset)
  {
    if (this != &RefOffset) {
      base = RefOffset.base;
      Tableau = RefOffset.Tableau;
      type_id = RefOffset.type_id;
      id = RefOffset.id;
      maxid = RefOffset.maxid;
    }
    return(*this);
  }

  bool ODSLF_Ix2_Offsets::operator<= (const ODSLF_Ix2_Offsets& RefOffsets)
  {
    // Check whether all nonzero tableau row lengths in RefOffsets
    // are <= than those in *this

    bool answer = true;
    for (int level = 0; level < 4; ++level) { // check fundamental level only
      // First check whether all rows which exist in both tableaux satisfy rule:
      for (int tableau_level = 0; tableau_level
	     < ABACUS::min(Tableau[level].Nrows, RefOffsets.Tableau[level].Nrows); ++tableau_level)
	if (Tableau[level].Row_L[tableau_level] > RefOffsets.Tableau[level].Row_L[tableau_level])
	  answer = false;
      // Now check whether there exist extra rows violating rule:
      for (int tableau_level = ABACUS::min(Tableau[level].Nrows, RefOffsets.Tableau[level].Nrows);
	   tableau_level < Tableau[level].Nrows; ++tableau_level)
	if (Tableau[level].Row_L[tableau_level] > 0) answer = false;
    }

    return(answer);
  }

  bool ODSLF_Ix2_Offsets::operator>= (const ODSLF_Ix2_Offsets& RefOffsets)
  {
    // Check whether all nonzero tableau row lengths in RefOffsets
    // are >= than those in *this

    bool answer = true;
    for (int level = 0; level < 4; ++level) { // check fundamental level only
      // First check whether all rows which exist in both tableaux satisfy rule:
      for (int tableau_level = 0; tableau_level
	     < ABACUS::min(Tableau[level].Nrows, RefOffsets.Tableau[level].Nrows); ++tableau_level)
	if (Tableau[level].Row_L[tableau_level] < RefOffsets.Tableau[level].Row_L[tableau_level])
	  answer = false;
      // Now check whether there exist extra rows violating rule:
      for (int tableau_level = ABACUS::min(Tableau[level].Nrows, RefOffsets.Tableau[level].Nrows);
	   tableau_level < RefOffsets.Tableau[level].Nrows; ++tableau_level)
	if (RefOffsets.Tableau[level].Row_L[tableau_level] > 0) answer = false;
    }

    return(answer);
  }

  void ODSLF_Ix2_Offsets::Set_to_id (long long int idnr)
  {
    // The idnr of the Offset is given by
    // sub_id[0] + (total number of tableaux of type 0) * (sub_id[1] + (total number of tableaux of type 1) * (sub_id[2] + ...
    // + total number of tableaux of type (2*base.Nrap.size()) * sub_id[2*base.Nrap.size() + 1]

    if (idnr > maxid) {
      cout << idnr << "\t" << maxid << endl;
      ABACUSerror("idnr too large in offsets.Set_to_id.");
    }

    id = idnr;

    Vect<long long int> sub_id(0LL, 2*base.Nrap.size() + 2);

    long long int idnr_eff = idnr;
    long long int temp_prod = 1LL;

    Vect<long long int> result_choose(2*base.Nrap.size() + 2);

    for (int i = 0; i <= 2*base.Nrap.size(); ++i) {
      result_choose[i] = Tableau[i].maxid + 1LL;
      temp_prod *= result_choose[i];
    }

    for (int i = 2*base.Nrap.size() + 1; i > 0; --i) {
      sub_id[i] = idnr_eff/temp_prod;
      idnr_eff -= sub_id[i] * temp_prod;
      temp_prod /= result_choose[i-1];
    }
    sub_id[0] = idnr_eff; // what's left goes to the bottom...

    for (int i = 0; i <= 2*base.Nrap.size() + 1; ++i) {
      if ((Tableau[i].Nrows * Tableau[i].Ncols == 0) && (sub_id[i] != 0))
	ABACUSerror("index too large in offset.Set_to_id.");
      if (Tableau[i].id != sub_id[i]) Tableau[i].Set_to_id(sub_id[i]);
    }

    Compute_type_id ();

    return;
  }

  void ODSLF_Ix2_Offsets::Compute_type_id ()
  {
    type_id = 0LL;

    for (int i = 0; i < 2*base.Nrap.size() + 2; ++i) {
      Tableau[i].Compute_id();
      type_id += Tableau[i].Nrows * pow_ulli(10LL, i);
    }
  }

  void ODSLF_Ix2_Offsets::Compute_id ()
  {
    long long int prod_maxid = 1LL;

    id = 0LL;

    for (int i = 0; i < 2*base.Nrap.size() + 2; ++i) {
      Tableau[i].Compute_id();
      id += Tableau[i].id * prod_maxid;
      prod_maxid *= Tableau[i].maxid + 1LL;
    }

  }

  bool ODSLF_Ix2_Offsets::Add_Boxes_From_Lowest (int Nboxes, bool odd_sectors)
  {
    // Tries to add Nboxes to Young tableau in given sector parity.
    // If Nboxes is too large to be accommodated, `false' is returned.

    ODSLF_Ix2_Offsets offsets_init(*this);  // keep a copy in case it fails

    if (Nboxes < 0) ABACUSerror("Can't use Nboxes < 0 in Add_Boxes_From_Lowest");
    else if (Nboxes == 0) return(true);  // nothing to do

    int Nboxes_put = 0;
    int working_sector = odd_sectors;
    // We start from the bottom up.
    while (Nboxes_put < Nboxes && working_sector < 2*base.Nrap.size() + 2) {
      Nboxes_put += Tableau[working_sector].Add_Boxes_From_Lowest(Nboxes - Nboxes_put);
      working_sector += 2;
    }
    Compute_id();

    if (Nboxes_put < Nboxes) (*this) = offsets_init;  // reset !
    else if (Nboxes_put > Nboxes) {
      cout << Nboxes_put << "\t" << Nboxes << endl;
      ABACUSerror("Putting too many boxes in Ix2_Offsets::Add_Boxes_From_Lowest.");
    }
    return(Nboxes_put == Nboxes);
  }

  //***************************************************************************************************

  // Function definitions:  class ODSLF_Ix2_Offsets_List

  ODSLF_Ix2_Offsets_List::ODSLF_Ix2_Offsets_List() : ndef(0), Offsets(Vect<ODSLF_Ix2_Offsets>(NBASESMAX)) {};

  ODSLF_Ix2_Offsets& ODSLF_Ix2_Offsets_List::Return_Offsets (ODSLF_Base& RefBase, Vect<int> nparticles)
  {

    long long int Ix2_Offsets_type_id_ref = Ix2_Offsets_type_id (nparticles);
    int n = 0;

    while (n < ndef && !(Ix2_Offsets_type_id_ref == Offsets[n].type_id && RefBase == Offsets[n].base)) n++;

    if (n == ndef) {
      Offsets[n] = ODSLF_Ix2_Offsets (RefBase, nparticles);
      ndef++;
    }

    if (ndef >= NBASESMAX) ABACUSerror("Ix2_Offsets_List:  need bigger Offsets vector.");

    return(Offsets[n]);
  }

  ODSLF_Ix2_Offsets& ODSLF_Ix2_Offsets_List::Return_Offsets (ODSLF_Base& RefBase, long long int req_type_id)
  {

    int n = 0;

    while (n < ndef && !(req_type_id == Offsets[n].type_id && RefBase == Offsets[n].base)) n++;

    if (n == ndef) {
      Offsets[n] = ODSLF_Ix2_Offsets (RefBase, req_type_id);
      ndef++;
    }

    if (ndef >= NBASESMAX) ABACUSerror("Ix2_Offsets_List:  need bigger Offsets vector.");

    return(Offsets[n]);
  }



  //***************************************************************************************************

  // Function definitions:  class ODSLF_Bethe_State

  ODSLF_Bethe_State::ODSLF_Bethe_State ()
    : chain(Heis_Chain()), base(ODSLF_Base()), offsets(ODSLF_Ix2_Offsets()), Ix2(ODSLF_Ix2_Config(chain, 1)),
      lambda(ODSLF_Lambda(chain, 1)), BE(ODSLF_Lambda(chain, 1)), diffsq(0.0), conv(0), iter(0), iter_Newton(0),
      E(0.0), iK(0), K(0.0), lnnorm(-100.0), base_id(0LL), type_id(0LL), id(0LL), maxid(0LL), nparticles(0)
  {
  };

  ODSLF_Bethe_State::ODSLF_Bethe_State (const ODSLF_Bethe_State& RefState) // copy constructor
    : chain(RefState.chain), base(RefState.base), offsets(RefState.offsets),
      Ix2(ODSLF_Ix2_Config(RefState.chain, RefState.base)),
      lambda(ODSLF_Lambda(RefState.chain, RefState.base)), BE(ODSLF_Lambda(RefState.chain, RefState.base)),
      diffsq(RefState.diffsq), conv(RefState.conv), iter(RefState.iter), iter_Newton(RefState.iter_Newton),
      E(RefState.E), iK(RefState.iK), K(RefState.K), lnnorm(RefState.lnnorm),
      base_id(RefState.base_id), type_id(RefState.type_id), id(RefState.id), maxid(RefState.maxid),
      nparticles(RefState.nparticles)
  {
    // copy arrays into new ones
    for (int j = 0; j < RefState.chain.Nstrings; ++j) {
      for (int alpha = 0; alpha < RefState.base[j]; ++j) {
	Ix2[j][alpha] = RefState.Ix2[j][alpha];
	lambda[j][alpha] = RefState.lambda[j][alpha];
      }
    }
  }

  ODSLF_Bethe_State::ODSLF_Bethe_State (const ODSLF_Bethe_State& RefState, long long int type_id_ref)
    : chain(RefState.chain), base(RefState.base), offsets(ODSLF_Ix2_Offsets(RefState.base, type_id_ref)),
      Ix2(ODSLF_Ix2_Config(RefState.chain, RefState.base)), lambda(ODSLF_Lambda(RefState.chain, RefState.base)),
      BE(ODSLF_Lambda(RefState.chain, RefState.base)), diffsq(1.0), conv(0), iter(0), iter_Newton(0),
      E(0.0), iK(0), K(0.0), lnnorm(-100.0),
      base_id(RefState.base_id), type_id(0LL), id(0LL), maxid(offsets.maxid), nparticles(0)
  {
    (*this).Set_Ix2_Offsets (offsets);
  }

  ODSLF_Bethe_State::ODSLF_Bethe_State (const Heis_Chain& RefChain, int M)
    : chain(RefChain), base(RefChain, M), offsets(base, 0LL), Ix2(ODSLF_Ix2_Config(RefChain, M)),
      lambda(ODSLF_Lambda(RefChain, M)),
      BE(ODSLF_Lambda(RefChain, M)), diffsq(1.0), conv(0), iter(0), iter_Newton(0),
      E(0.0), iK(0), K(0.0), lnnorm(-100.0),
      base_id(0LL), type_id(0LL), id(0LL), maxid(offsets.maxid), nparticles(0)
  {
  }

  ODSLF_Bethe_State::ODSLF_Bethe_State (const Heis_Chain& RefChain, const ODSLF_Base& RefBase)
    : chain(RefChain), base(RefBase), offsets(RefBase, 0LL),
      Ix2(ODSLF_Ix2_Config(RefChain, RefBase)), lambda(ODSLF_Lambda(RefChain, RefBase)),
      BE(ODSLF_Lambda(RefChain, RefBase)), diffsq(1.0), conv(0), iter(0), iter_Newton(0),
      E(0.0), iK(0), K(0.0), lnnorm(-100.0),
      base_id(RefBase.id), type_id(0LL), id(0LL), maxid(offsets.maxid), nparticles(0)
  {
    // Check that the number of rapidities is consistent with Mdown
    int Mcheck = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += base[i] * RefChain.Str_L[i];
    if (RefBase.Mdown != Mcheck)
      ABACUSerror("Wrong M from Nrapidities input vector, in ODSLF_Bethe_State constructor.");

  }

  ODSLF_Bethe_State::ODSLF_Bethe_State (const Heis_Chain& RefChain,
					long long int base_id_ref, long long int type_id_ref)
    : chain(RefChain), base(ODSLF_Base(RefChain, base_id_ref)), offsets(base, type_id_ref),
      Ix2(ODSLF_Ix2_Config(RefChain, base)), lambda(ODSLF_Lambda(RefChain, base)),
      BE(ODSLF_Lambda(RefChain, base)), diffsq(1.0), conv(0), iter(0), iter_Newton(0),
      E(0.0), iK(0), K(0.0), lnnorm(-100.0),
      base_id(base_id_ref), type_id(type_id_ref), id(0LL), maxid(offsets.maxid), nparticles(0)
  {
  }

  void ODSLF_Bethe_State::Set_Ix2_Offsets (const ODSLF_Ix2_Offsets& RefOffset)  // sets the Ix2 to given offsets
  {
    if (base != RefOffset.base)
      ABACUSerror("Offset with incompatible base used in ODSLF_Bethe_State::Set_Ix2_Offsets.");

    // For each base_level, we make sure that the Ix2 are ordered:  Ix2[0][j] < Ix2[0][k] if j < k.

    // Level 2:  particles in R part outside GS interval
    for (int j = 0; j < RefOffset.Tableau[2].Nrows; ++j) {
      Ix2[0][base.Nrap[0] - 1 - j]
	= base.Nrap[0] - 1 + 2* RefOffset.Tableau[2].Nrows - 2*j + 2*RefOffset.Tableau[2].Row_L[j];
    }
    nparticles = RefOffset.Tableau[2].Nrows;

    // Level 0:  holes in R part of GS interval
    for (int j = 0; j < RefOffset.Tableau[0].Ncols; ++j) {
      Ix2[0][base.Nrap[0] - 1 - RefOffset.Tableau[2].Nrows - j] = (base.Nrap[0] + 1) % 2
	+ 2*RefOffset.Tableau[0].Ncols - 2 - 2*j + 2*RefOffset.Tableau[0].Col_L[j];
    }
    nparticles += RefOffset.Tableau[0].Ncols;

    // Level 1:  holes in L part of GS interval
    for (int j = 0; j < RefOffset.Tableau[1].Ncols; ++j) {
      Ix2[0][base.Nrap[0] - 1 - RefOffset.Tableau[0].Ncols - RefOffset.Tableau[2].Nrows - j]
	= -1 - (base.Nrap[0] % 2) - 2* j - 2*RefOffset.Tableau[1].Col_L[RefOffset.Tableau[1].Ncols - 1 - j];
    }
    nparticles += RefOffset.Tableau[1].Ncols;

    // Level 3:  particles in L part outside GS interval
    for (int j = 0; j < RefOffset.Tableau[3].Nrows; ++j) {
      Ix2[0][base.Nrap[0] - 1 - RefOffset.Tableau[0].Ncols - RefOffset.Tableau[1].Ncols - RefOffset.Tableau[2].Nrows - j]
	= -(base.Nrap[0] + 1 + 2*j + 2*RefOffset.Tableau[3].Row_L[RefOffset.Tableau[3].Nrows - 1 - j]);
    }
    nparticles += RefOffset.Tableau[3].Nrows;

    // Subsequent levels:  particles on R and L
    for (int base_level = 1; base_level < base.Nrap.size(); ++base_level) {

      for (int j = 0; j < RefOffset.Tableau[2*base_level + 2].Nrows; ++j)
	Ix2[base_level][base.Nrap[base_level] - 1 - j]
	  = ((base.Nrap[base_level] + 1) % 2) + 2*(RefOffset.Tableau[2*base_level + 2].Nrows - 1) -2*j
	  + 2*RefOffset.Tableau[2*base_level + 2].Row_L[j];

      nparticles += RefOffset.Tableau[2*base_level + 2].Nrows;

      for (int j = 0; j < RefOffset.Tableau[2*base_level + 3].Nrows; ++j)
	Ix2[base_level][base.Nrap[base_level] - 1 - RefOffset.Tableau[2*base_level + 2].Nrows - j]
	  = -1 - ((base.Nrap[base_level]) % 2) -2*j
	  - 2*RefOffset.Tableau[2*base_level + 3].Row_L[RefOffset.Tableau[2*base_level + 3].Nrows - 1 - j];

      nparticles += RefOffset.Tableau[2*base_level + 3].Nrows;
    }

    // ODSLF:  put back correct parity of quantum nr:
    for (int j = 0; j < base.Nrap.size(); ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha)
	if ((Ix2[j][alpha] + chain.Nsites - 1 + base.Nrap[j] + (base.Nraptot + 1) * chain.Str_L[j]
	     + (chain.par[j] == 1 ? 0 : chain.Nsites)) % 2)
	  Ix2[j][alpha] -= 1;

    type_id = RefOffset.type_id;
    id = RefOffset.id;
    maxid = RefOffset.maxid;

    return;
  }

  void ODSLF_Bethe_State::Set_to_id (long long int id_ref)
  {
    offsets.Set_to_id (id_ref);
    (*this).Set_Ix2_Offsets (offsets);
  }

  void ODSLF_Bethe_State::Set_to_id (long long int id_ref, ODSLF_Bethe_State& RefState)
  {
    offsets.Set_to_id (id_ref);
    (*this).Set_Ix2_Offsets (offsets);
  }

  bool ODSLF_Bethe_State::Check_Symmetry ()
  {
    // Checks whether the I's are symmetrically distributed.
    bool symmetric_state = true;
    int arg, test1, test3;

    for (int j = 0; j < chain.Nstrings; ++j) {
      test1 = 0;
      test3 = 0;
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	// since Ix2 = N is same as Ix2 = -N by periodicity, this is symmetric.
	arg = (Ix2[j][alpha] != chain.Nsites) ? Ix2[j][alpha] : 0;
	test1 += arg;
	test3 += arg * arg * arg;  // to make sure that all I's are symmetrical...
      }
      if (!(symmetric_state && (test1 == 0) && (test3 == 0))) symmetric_state = false;
    }

    return symmetric_state;
  }

  void ODSLF_Bethe_State::Compute_diffsq ()
  {
    DP maxterm = 0.0;
    int jmax, alphamax;

    diffsq = 0.0;
    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	diffsq += pow(BE[j][alpha], 2.0);
	if (pow(BE[j][alpha], 2.0)/chain.Nsites > maxterm) {
	  jmax = j;
	  alphamax = alpha;
	  maxterm = pow(BE[j][alpha], 2.0)/chain.Nsites;
	}
      }
    diffsq /= DP(chain.Nsites);
  }

  void ODSLF_Bethe_State::Iterate_BAE ()
  {
    for (int j = 0; j < chain.Nstrings; ++j) {
      for (int alpha = 0; alpha < base[j]; ++alpha)
	{
	  lambda[j][alpha] = Iterate_BAE (j, alpha);
	}
    }

    iter++;

    (*this).Compute_BE();
    (*this).Compute_diffsq();
  }

  void ODSLF_Bethe_State::BAE_smackdown (DP max_allowed)
  {
    // Re-solves for all rapidities lambda[j][alpha] such that BE[j][alpha]^2/N > max_allowed.
    // Assumes that BE[][] is up-to-date.

    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha)
	if (pow(BE[j][alpha], 2.0)/chain.Nsites > max_allowed) (*this).Solve_BAE (j, alpha, max_allowed, 100);
  }

  void ODSLF_Bethe_State::Solve_BAE_smackdown (DP max_allowed, int maxruns)
  {
    int runs_done = 0;

    (*this).Compute_BE();
    (*this).Compute_diffsq();

    while (diffsq > chain.prec && diffsq > max_allowed && runs_done < maxruns) {
      (*this).BAE_smackdown (max_allowed);
      (*this).Compute_BE();
      (*this).Compute_diffsq();
      runs_done++;
    }
  }

  void ODSLF_Bethe_State::Iterate_BAE_Newton ()
  {
    // does one step of a Newton method on the rapidities...
    // Assumes that BE[j][alpha] have been computed

    Vect_CX dlambda (0.0, base.Nraptot);  // contains delta lambda computed from Newton's method
    SQMat_CX Gaudin (0.0, base.Nraptot);
    Vect_INT indx (base.Nraptot);

    ODSLF_Lambda lambda_ref(chain, base);
    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];

    (*this).Build_Reduced_Gaudin_Matrix (Gaudin);

    int index = 0;
    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	dlambda[index] = - BE[j][alpha] * chain.Nsites;
	index++;
      }

    DP d;

    try {
      ludcmp_CX (Gaudin, indx, d);
      lubksb_CX (Gaudin, indx, dlambda);
    }

    catch (Divide_by_zero) {
      diffsq = log(-1.0);  // reset to nan to stop Newton method
      return;
    }

    // Regularize dlambda:  max step is +-1.0 to prevent rapidity flying off into outer space.
    for (int i = 0; i < base.Nraptot; ++i) if (fabs(real(dlambda[i])) > 1.0) dlambda[i] = 0.0;

    index = 0;
    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	lambda[j][alpha] = lambda_ref[j][alpha] + real(dlambda[index]);
	index++;
      }

    (*this).Compute_BE();
    (*this).Iterate_BAE();
    (*this).Compute_diffsq();

    // if we've converged, calculate the norm here, since the work has been done...
    if (diffsq < chain.prec) {
      lnnorm = 0.0;
      for (int j = 0; j < base.Nraptot; j++) lnnorm += log(abs(Gaudin[j][j]));
    }

    iter_Newton++;

    return;
  }

  void ODSLF_Bethe_State::Solve_BAE (int j, int alpha, DP req_prec, int itermax)
  {
    // Finds the root lambda[j][alpha] to precision req_prec

    DP prec_obtained = 1.0;
    int niter_here = 0;

    while (prec_obtained > req_prec && niter_here < itermax) {
      lambda[j][alpha] = (*this).Iterate_BAE (j, alpha);
      (*this).Compute_BE (j, alpha);
      prec_obtained = pow(BE[j][alpha], 2.0)/chain.Nsites;
      niter_here++;
    }
  }

  void ODSLF_Bethe_State::Find_Rapidities (bool reset_rapidities)
  {
    // This function finds the rapidities of the eigenstate

    lnnorm = -100.0;  // sentinel value, recalculated if Newton method used in the last step of iteration.

    ODSLF_Lambda lambda_ref(chain, base);

    diffsq = 1.0;

    if (reset_rapidities)
      (*this).Set_Free_lambdas();
    (*this).Compute_BE();

    iter = 0;
    iter_Newton = 0;

    // Start with conventional iterations

    DP iter_prec = 1.0e-2/DP(chain.Nsites);

    clock_t interp_start;
    clock_t interp_stop;
    clock_t Newton_start;
    clock_t Newton_stop;
    int iter_interp_start, iter_interp_stop, iter_Newton_start, iter_Newton_stop;
    DP diffsq_interp_start, diffsq_interp_stop, diffsq_Newton_start, diffsq_Newton_stop;

    while (diffsq > chain.prec && iter < 300 && iter_Newton < 20) {

      do {
	interp_start = clock();
	iter_interp_start = iter;
	diffsq_interp_start = diffsq;

	if (diffsq > iter_prec) (*this).Solve_BAE_interp (iter_prec, 10);

	interp_stop = clock();
	iter_interp_stop = iter;
	diffsq_interp_stop = diffsq;

	iter_prec = diffsq * 0.001;

      } while (diffsq > chain.prec && !is_nan(diffsq) && diffsq_interp_stop/diffsq_interp_start < 0.001);

      // If we haven't even reached iter_prec, try normal iterations without extrapolations
      if (diffsq > iter_prec) {
	(*this).Solve_BAE_smackdown (0.1 * diffsq, 1);
      }

      // Now try Newton

      Newton_start = clock();
      iter_Newton_start = iter_Newton;
      diffsq_Newton_start = diffsq;

      if (diffsq > iter_prec) (*this).Solve_BAE_Newton (chain.prec, 5);

      Newton_stop = clock();
      iter_Newton_stop = iter_Newton;
      diffsq_Newton_stop = diffsq;

      if (diffsq > iter_prec) (*this).Solve_BAE_smackdown (0.1 * diffsq, 1);

      if (diffsq > iter_prec) (*this).Solve_BAE_straight_iter (0.1 * diffsq, 50);

      // If we haven't converged here, retry with more precision in interp method...

      if (is_nan(diffsq)) {
	(*this).Set_Free_lambdas();  // good start if we've messed up with Newton
	diffsq = 1.0;
      }

      iter_prec *= 1.0e-4;
      iter_prec = ABACUS::max(iter_prec, chain.prec);

    }

    // Check convergence:
    conv = (diffsq < chain.prec && (*this).Check_Rapidities()) ? 1 : 0;

    return;
  }

  void ODSLF_Bethe_State::Solve_BAE_Newton (DP Newton_prec, int max_iter_Newton)
  {
    // This function attempts to get convergence diffsq <= Newton_prec in at most max_iter_Newton steps.

    // The results are accepted if diffsq has decreased, otherwise the lambda's are reset to original values, defined as...

    ODSLF_Lambda lambda_ref(chain, base);
    DP diffsq_ref = 1.0;

    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
    diffsq_ref = diffsq;

    // Now begin solving...

    int iter_done_here = 0;

    while ((diffsq > Newton_prec) && (diffsq < 10.0) && (!is_nan(diffsq)) && (iter_done_here < max_iter_Newton)) {
      (*this).Iterate_BAE_Newton();
      iter_done_here++;
      (*this).Iterate_BAE();
    }

    if ((diffsq > diffsq_ref) || (is_nan(diffsq))) {
      // This procedure has failed.  We reset everything to begin values.
      for (int j = 0; j < chain.Nstrings; ++j)
	for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
      diffsq = diffsq_ref;
    }

    return;
  }

  void ODSLF_Bethe_State::Solve_BAE_straight_iter (DP interp_prec, int max_iter_interp)
  {
    // This function attempts to get convergence diffsq <= interp_prec in at most max_iter_interp steps.

    // If this fails, the lambda's are reset to original values, defined as...

    ODSLF_Lambda lambda_ref(chain, base);
    DP diffsq_ref = 1.0;

    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
    diffsq_ref = diffsq;

    // Now begin solving...

    diffsq = 1.0;
    int iter_done_here = 0;

    while ((diffsq > interp_prec) && (max_iter_interp > iter_done_here)) {
      (*this).Iterate_BAE();
      iter_done_here++;
    }

    if ((diffsq > diffsq_ref)) {
      // This procedure has failed.  We reset everything to begin values.
      for (int j = 0; j < chain.Nstrings; ++j)
	for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
      diffsq = diffsq_ref;
    }
  }

  void ODSLF_Bethe_State::Solve_BAE_interp (DP interp_prec, int max_iter_interp)
  {
    // This function attempts to get convergence diffsq <= interp_prec in at most max_iter_interp steps.

    // If this fails, the lambda's are reset to original values, defined as...

    ODSLF_Lambda lambda_ref(chain, base);
    DP diffsq_ref = 1.0;

    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
    diffsq_ref = diffsq;

    // Now begin solving...

    diffsq = 1.0;
    int iter_done_here = 0;

    ODSLF_Lambda lambda1(chain, base);
    ODSLF_Lambda lambda2(chain, base);
    ODSLF_Lambda lambda3(chain, base);
    ODSLF_Lambda lambda4(chain, base);
    ODSLF_Lambda lambda5(chain, base);

    while ((diffsq > interp_prec) && (max_iter_interp > iter_done_here)) {

      (*this).Iterate_BAE();
      for (int j = 0; j < chain.Nstrings; ++j)
	for (int alpha = 0; alpha < base[j]; ++alpha) lambda1[j][alpha] = lambda[j][alpha];
      (*this).Iterate_BAE();
      for (int j = 0; j < chain.Nstrings; ++j)
	for (int alpha = 0; alpha < base[j]; ++alpha) lambda2[j][alpha] = lambda[j][alpha];
      (*this).Iterate_BAE();
      for (int j = 0; j < chain.Nstrings; ++j)
	for (int alpha = 0; alpha < base[j]; ++alpha) lambda3[j][alpha] = lambda[j][alpha];
      (*this).Iterate_BAE();
      for (int j = 0; j < chain.Nstrings; ++j)
	for (int alpha = 0; alpha < base[j]; ++alpha) lambda4[j][alpha] = lambda[j][alpha];

      iter_done_here += 4;

      // now extrapolate the result to infinite number of iterations...

      Vect_DP rap(0.0, 4);
      Vect_DP oneoverP(0.0, 4);
      DP deltalambda = 0.0;

      for (int i = 0; i < 4; ++i) oneoverP[i] = 1.0/(1.0 + i*i);

      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) {
	  rap[0] = lambda1[j][alpha];
	  rap[1] = lambda2[j][alpha];
	  rap[2] = lambda3[j][alpha];
	  rap[3] = lambda4[j][alpha];

	  polint (oneoverP, rap, 0.0, lambda[j][alpha], deltalambda);
	}

      (*this).Iterate_BAE();
    }

    if ((diffsq >= diffsq_ref)) {
      // This procedure has failed.  We reset everything to begin values.
      for (int j = 0; j < chain.Nstrings; ++j)
	for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
      diffsq = diffsq_ref;
    }

    return;
  }

  void ODSLF_Bethe_State::Compute_lnnorm ()
  {
    if (true || lnnorm == -100.0) {  // else norm already calculated by Newton method
      // Actually, compute anyway to increase accuracy !

      SQMat_CX Gaudin_Red(base.Nraptot);

      (*this).Build_Reduced_Gaudin_Matrix(Gaudin_Red);
      complex<DP> lnnorm_check = lndet_LU_CX_dstry(Gaudin_Red);
      lnnorm = real(lnnorm_check);
    }

    return;
  }

  void ODSLF_Bethe_State::Compute_Momentum ()
  {
    int sum_Ix2 = 0;
    DP sum_M = 0.0;

    for (int j = 0; j < chain.Nstrings; ++j) {
      sum_M += 0.5 * (1.0 + chain.par[j]) * base[j];
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	sum_Ix2 += Ix2[j][alpha];
      }
    }

    iK = (chain.Nsites/2) * int(sum_M + 0.1) - (sum_Ix2/2);  // + 0.1:  for safety...

    while (iK >= chain.Nsites) iK -= chain.Nsites;
    while (iK < 0) iK += chain.Nsites;

    K = PI * sum_M - PI * sum_Ix2/chain.Nsites;

    while (K >= 2.0*PI) K -= 2.0*PI;
    while (K < 0.0) K += 2.0*PI;

    return;
  }

  void ODSLF_Bethe_State::Compute_All (bool reset_rapidities)     // solves BAE, computes E, K and lnnorm
  {
    (*this).Find_Rapidities (reset_rapidities);
    if (conv == 1) {
      (*this).Compute_Energy ();
      (*this).Compute_Momentum ();
      (*this).Compute_lnnorm ();
    }
    return;
  }

  bool ODSLF_Bethe_State::Boost_Momentum (int iKboost)
  {
    if (iKboost == 0) return(true);  // done

    ODSLF_Ix2_Offsets offsets_here = offsets;

    bool success = false;

    if (iKboost < 0)
      success = offsets_here.Add_Boxes_From_Lowest(-iKboost, 0);  // add boxes in even sectors to decrease iK

    else if (iKboost > 0)
      success = offsets_here.Add_Boxes_From_Lowest(iKboost, 1);  // add boxes in odd sectors to increase iK

    if (success) (*this).Set_Ix2_Offsets(offsets_here);

    return(success);
  }


  std::ostream& operator<< (std::ostream& s, const ODSLF_Bethe_State& state)
  {
    // sends all the state data to output stream

    s << endl << "******** Chain with Nsites = " << state.chain.Nsites
      << " Mdown (nr fermions) = " << state.base.Mdown
      << ":  eigenstate with base_id " << state.base_id << ", type_id " << state.type_id
      << " id " << state.id << " maxid " << state.offsets.maxid << endl
      << "E = " << state.E << " K = " << state.K << " iK = " << state.iK
      << " lnnorm = " << state.lnnorm << endl
      << "conv = " << state.conv << " iter = " << state.iter
      << " iter_Newton = " << state.iter_Newton << "\tdiffsq " << state.diffsq << endl;

    for (int j = 0; j < state.chain.Nstrings; ++j) {
      if (state.base.Nrap[j] > 0) {
	s << "Type " << j << " Str_L = " << state.chain.Str_L[j]
	  << " par = " << state.chain.par[j] << " M_j = " << state.base.Nrap[j]
	  << " Ix2_infty = " << state.base.Ix2_infty[j] << " Ix2_max = " << state.base.Ix2_max[j] << endl;
	Vect_INT qnumbers(state.base.Nrap[j]);
	Vect_DP rapidities(state.base.Nrap[j]);
	for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) {
	  qnumbers[alpha] = state.Ix2[j][alpha];
	  rapidities[alpha] = state.lambda[j][alpha];
	}
	qnumbers.QuickSort();
	rapidities.QuickSort();

	s << "Ix2 quantum numbers: " << endl;
	for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) s << qnumbers[alpha] << " ";
	s << endl;
	s << "Rapidities: " << endl;
	for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) s << rapidities[alpha] << " ";
	s << endl;
      }
    }
    s << endl;

    return s;

  }

} // namespace ABACUS