polynomialRootSolution.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           polynomialRootSolution.cpp
// author(s):      Jan Novak
// language:       C, C++
// license:        This program 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 of the License, or
//                 (at your option) any later version.
//
//                 This program is distributed in the hope that it will be useful,
//                 but WITHOUT ANY WARRANTY; without even the implied warranty of
//                 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//                 GNU Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cmath>

#include "macros.h"
#include "polynomialRootSolution.h"
#include "mathlib.h"
#include "elementaryFunctions.h"


namespace mumech {

//********************************************************************************************************
//                                            PUBLIC FUNCTIONS
//********************************************************************************************************


//********************************************************************************************************
// description: Function gives the analitical solution of cubic
// furmula of the real coefficients a, b, c
// x^3 + ax^2 + bx + c = 0
// last edit:   30. 11. 2008
//-------------------------------------------------------------------------------------------------------
// Solution adopted from Numerical Recepies book
//********************************************************************************************************

std::complex<double> * polynomialRootSolution::GetCubicPolyRoots( double a, double b, double c, std::complex<double> * roots )
{
  double Q = ( SQR( a ) - 3. * b )/9. , R = ( 2. * CUB( a ) - 9. * a * b + 27. * c )/54., Theta = 0.;

  if( SQR( R ) < CUB( Q ) ){ //all roots real
    Theta = acos( R/ sqrt( CUB( Q ) ) );
    roots[0].real(-2. * sqrt( Q ) * cos( Theta/3. ) - a/3.);
    roots[1].real(-2. * sqrt( Q ) * cos( ( Theta + 2 * PI )/3. ) - a/3.);
    roots[2].real(-2. * sqrt( Q ) * cos( ( Theta - 2 * PI )/3. ) - a/3.);
    roots[0].imag(0.);
    roots[1].imag(0.);
    roots[2].imag(0.);
  }
  else{
    //    double A = -pow( ( R + sqrt( SQR( R ) - CUB( Q ) ) ), 1./3. );
    double A = -sgn( R ) * pow( ( ABS( R ) + sqrt( SQR( R ) - CUB( Q ) ) ), 1./3. ), B = ( A == 0 )? 0 : Q/A;
    //the only real root
    roots[0].real(( A + B ) - a/3.); roots[0].imag(0.);
    //complex roots
    roots[1].real(-.5 * ( A + B ) - a/3.); roots[2].real(-.5 * ( A + B ) - a/3.);
    roots[1].imag(sqrt( 3.)/2. * ( A - B ));
    roots[2].imag(-sqrt( 3.)/2. * ( A - B ));
  }
  return( roots );
}//end of function: GetCubicFormSolu()

//********************************************************************************************************
// description: Function gives the analitical solution of cubic furmula of generally
//              complex coefficients a, b, c; i.e. solution to: x^3 + ax^2 + bx + c = 0
// last edit:   02. 08. 2010
//-------------------------------------------------------------------------------------------------------
// note :       Solution adopted from Numerical Recepies book
//********************************************************************************************************
std::complex<double> * polynomialRootSolution::GetCubicPolyRoots( std::complex<double> a, std::complex<double> b,
                                  std::complex<double> c, std::complex<double>* roots )
{
  std::complex<double> Q = ( SQR( a ) - 3. * b )/9. , R = ( 2. * CUB( a ) - 9. * a * b + 27. * c )/54.;

  //initial cleansing
  roots[0] = roots[1] = roots[2] = (0. , 0.);

  //all roots real
  if( ABS( imag( Q ) ) <= _MACHINE_EPS_ &&
      ABS( imag( R ) ) <= _MACHINE_EPS_ &&
      ( SQR( real( R ) ) < CUB( real( Q ) ) ) )
    {
      double Theta = acos( real( R )/ sqrt( CUB( real( Q ) ) ) ),
          aux_mlt1 = -2. * sqrt( real( Q ) );
      std::complex<double> aux_mlt2 = (real( a )/3. , 0.);

      roots[0] = aux_mlt1 * cos( Theta/3. ) - aux_mlt2;
      roots[1] = aux_mlt1 * cos( ( Theta + 2 * PI )/3. ) - aux_mlt2;
      roots[2] = aux_mlt1 * cos( ( Theta - 2 * PI )/3. ) - aux_mlt2;
    }
  else{//generally complex roots
    std::complex<double> aux_mlt1 = sqrt( SQR( R ) - CUB( Q ) );
    int signum = ( real( aux_mlt1 ) < 0. )? -1 : 1;
    std::complex<double> A = -1. * pow( ( R + (double)signum * aux_mlt1 ), 1./3. ),
      B = ( ABS( real( A ) ) <= _MACHINE_EPS_ &&
        ABS( imag( A ) ) <= _MACHINE_EPS_ ) ? (0. , 0.) : Q/A,
      AplusB = A + B, AminusB = A - B, a_div_3 = a/3.,
      auxConjPart = (0. , AminusB * _sgrt3div2_);

    //generally complex roots
    roots[0] = AplusB - a_div_3;
    roots[1] = roots[2] = -.5 * AplusB - a/3.;
    roots[1] += auxConjPart;
    roots[2] -= auxConjPart;
  }

  return roots;
}//end of function: GetCubicFormSolu()


void
polynomialRootSolution::GauLegF(double *x, double *w, const int nPoints, const double minLim, const double maxLim)
{
  /* compute x(i) and w(i)  i=1,n  Legendre ordinates and weights           */
  /* on interval  -1.0 to 1.0 (length is 2.0) and scaled to given interval  */
  /* nPoints = number of base points                                        */
  /* use ordinates and weights for Gauss Legendre integration               */

  if( x == NULL )  _errorr( "Field for point positions must be allocated" );
  if( w == NULL )  _errorr( "Field for weights must be allocated" );

  int i, j, m;
  double jj, nn;
  double p1, p2, p3, pp, xl, xm, z, z1;

  // zero vectors
  CleanVector(x,nPoints);
  CleanVector(w,nPoints);

  m = (nPoints+1) / 2;
  xm = 0.5 * (maxLim + minLim);
  xl = 0.5 * (maxLim - minLim);

  nn = (double)nPoints;
  for(i = 0; i < m; i++)
  {
    z = cos(PI * ((double)i + 0.75) / ((double)nPoints + 0.5));
    do{
      p1 = 1.0;
      p2 = 0.0;
      for(j = 0; j < nPoints; j++)
      {
        jj = double(j);
        p3 = p2;
        p2 = p1;
        p1 = ((2.*jj + 1.)*z*p2 - jj*p3) / (jj + 1.);
      }
      pp = nn * (z*p1 - p2) / (z*z - 1.);
      z1 = z;
      z  = z1 - p1 / pp;
    } while( fabs(z - z1) >= 1.e-12);

    x[i]           = xm - xl * z;
    x[nPoints-1-i] = xm + xl * z;
    w[i]           = 2. * xl / ((1. - z*z) * pp*pp);
    w[nPoints-1-i] = w[i];
  }
}

void
polynomialRootSolution::GauHerQ(double *x, double *w, const int nPoints)   //-+ infinity
{
  /* nPoints = number of base points                                        */

  if( x == NULL )  _errorr( "Field for point positions must be allocated" );
  if( w == NULL )  _errorr( "Field for weights must be allocated" );

  int i, j, m;
  double p1, p2, p3, pp, z, z1;
  double nn, jj, PIM4;

  PIM4 = pow(PI, -0.25);
  z = 0.;
  m = (nPoints+1) / 2;

  nn = (double)nPoints;
  for( i = 1; i <= m; i++ )
  {
    if( i==1) {
      z = sqrt(2.*nn + 1.) - 1.85575 * pow(2.*nn + 1., -0.16667);

    } else if( i==2 ) {
      z -= 1.14*pow(nn, 0.426) / z;

    } else if( i==3 ) {
      z = 1.86*z - 0.86 * x[1];

    } else if( i==4 ) {
      z = 1.91*z - 0.91 * x[2];

    } else {
      z = 2.0*z - x[i-2];
    }

    do
    {
      p1 = PIM4;
      p2 = 0.0;
      for( j = 1; j <= nPoints; j++ )
      {
        jj = double(j);
        p3 = p2;
        p2 = p1;
        p1 = z*sqrt(2./jj)*p2 - sqrt((jj-1.)/jj)*p3;
      }
      pp = sqrt(2.*nn) * p2;
      z1 = z;
      z  = z1 - p1/pp;
    } while (fabs(z-z1) >= 1.e-12);

    x[i-1]       = z;
    x[nPoints-i] = -z;
    w[i-1]       = 2.0 / (pp * pp);
    w[nPoints-i] = w[i-1];
  }
}

} // end of namespace mumech

/*end of file*/