Laguerre.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "Laguerre.hh"
#include "Base/Debug/Error.hh"
#include <Base/Common/BIASpragma.hh>

#include <limits>
using namespace std;
using namespace BIAS;

int LaguerreSolver::Solve(const std::vector< double >& coeffs,
           std::vector< std::complex<double> >& roots,
           const double eps)
{
  std::vector< std::complex<double> > complcoeffs;
  for(unsigned int i=0; i<coeffs.size(); i++)
  {
    //cout<<i<<"  "<<coeffs[i]<<endl;
    complcoeffs.push_back(complex<double>(coeffs[i],0));
  }
  return(Solve(complcoeffs,roots,eps));
}

int LaguerreSolver::Solve(const std::vector< std::complex<double> >& coeffs,
                          std::vector< std::complex<double> >& roots,
                          const double eps)
{
  
  int i,iter;
  complex<double> x,b,c;
  roots.clear();
  int m = coeffs.size()-1;
  vector< complex<double> > ad(m+1);
  for (int j=0;j<=m;j++)
  {
    ad[j]=coeffs[j];
  }
  for ( int j=m-1;j>=0;j--) 
  {
    x = 0;
    vector< complex<double> > ad_v(j+2);
    for ( int jj=0;jj<j+2;jj++)
    {
      ad_v[jj]=ad[jj];
    }
    iter = _Laguerre(ad_v,x);
    if(iter<0) return -1;
    if (abs(imag(x)) <= 2.0*eps*abs(real(x)))
    {
      x=complex<double>(real(x),0.0);
    }
    roots.push_back(x);
    b=ad[j+1];
    for (int jj=j;jj>=0;jj--) 
    {
      c=ad[jj];
      ad[jj]=b;
      b=x*b+c;
    }
  }
  if (_doPolish)
  {
    for(int j=0;j<m;j++)
    {
      iter = _Laguerre(coeffs,roots[j]);  //j is
      if(iter<0) return -1;
    }
  }
  for (int j=1;j<m;j++) 
  {
    x=roots[j];
    for(i=j-1;i>=0;i--) 
    {
      if(real(roots[i])<=real(x)) break;
      roots[i+1]=roots[i];
    }
    roots[i+1]=x;
  }
  return 0;
}

int LaguerreSolver::_Laguerre(const std::vector< std::complex<double> >& coeffs,
                              std::complex<double>& x)
{
  int iter = 0;
  const  int MR = 8;
  const  int MT = 10;
  const  int MAXIT = MT * MR;
  const double EPS = numeric_limits<double>::epsilon();
  static const double frac[MR+1] = {0.0,0.5,0.25,0.75,0.13,0.38,0.62,0.88,1.0};
  complex<double> dx,x1,b,d,f,g,h,sq,gp,gm,g2;
  unsigned int m=coeffs.size()-1;
  for (int it=1;it<=MAXIT;it++) 
  {
    iter=it;
    b=coeffs[m];
    double err=abs(b);
    d=f=0.0;
    double abx=abs(x);
    for (int j=m-1;j>=0;j--) 
    {
      f=x*f+d;
      d=x*d+b;
      b=x*b+coeffs[j];
      err=abs(b)+abx*err;
    }
    err *= EPS;
    if (abs(b) <= err) return iter;
    g=d/b;
    g2=g*g;
    h=g2-2.0*f/b;
    sq=sqrt(double(m-1)*(double(m)*h-g2));
    gp=g+sq;
    gm=g-sq;
    double abp=abs(gp);
    double abm=abs(gm);
    if(abp < abm)
    {
      gp=gm;
    }
    if(max(abp,abm)>0.0)
    {
      dx = ((double) m)/gp;
    }
    else
    {
      dx = std::polar(1+abx,((double)it)); 
    }   
    x1=x-dx;
    if(x==x1)
    { 
      return iter;
    }
    if((it%MT)!=0)
    { 
      x=x1;
    }
    else
    {
      x -= frac[it/MT]*dx;
    }
  }
  BIASWARN("too many iterations");
  return -1;
}