source: trunk/source/global/HEPNumerics/src/G4AnalyticalPolSolver.cc @ 1202

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: G4AnalyticalPolSolver.cc,v 1.7 2007/11/13 17:35:06 gcosmo Exp $
// GEANT4 tag $Name: geant4-09-02-ref-02 $
//

#include  "globals.hh"
#include  <complex>

#include "G4AnalyticalPolSolver.hh"

//////////////////////////////////////////////////////////////////////////////

G4AnalyticalPolSolver::G4AnalyticalPolSolver() {;}

//////////////////////////////////////////////////////////////////////////////

G4AnalyticalPolSolver::~G4AnalyticalPolSolver() {;}

//////////////////////////////////////////////////////////////////////////////
//
// Array r[3][5]  p[5]
// Roots of poly p[0] x^2 + p[1] x+p[2]=0
//
// x = r[1][k] + i r[2][k];  k = 1, 2

G4int G4AnalyticalPolSolver::QuadRoots( G4double p[5], G4double r[3][5] )
{
  G4double b, c, d2, d;

  b  = -p[1]/p[0]/2.;
  c  =  p[2]/p[0];
  d2 =  b*b - c;
  
  if( d2 >= 0. )
  {
    d       = std::sqrt(d2);
    r[1][1] = b - d;   
    r[1][2] = b + d;   
    r[2][1] = 0.; 
    r[2][2] = 0.;
  }
  else
  {
    d       = std::sqrt(-d2);
    r[2][1] =  d; 
    r[2][2] = -d;
    r[1][1] =  b; 
    r[1][2] =  b;
  }

  return 2;
}

//////////////////////////////////////////////////////////////////////////////
//
// Array r[3][5]  p[5]
// Roots of poly p[0] x^3 + p[1] x^2...+p[3]=0
// x=r[1][k] + i r[2][k]  k=1,...,3
// Assumes 0<arctan(x)<pi/2 for x>0

G4int G4AnalyticalPolSolver::CubicRoots( G4double p[5], G4double r[3][5] )
{
  G4double x,t,b,c,d;
  G4int k;

  if( p[0] != 1. )
  {
    for(k = 1; k < 4; k++ ) { p[k] = p[k]/p[0]; }
    p[0] = 1.;
  }
  x = p[1]/3.0; 
  t = x*p[1];
  b = 0.5*( x*( t/1.5 - p[2] ) + p[3] ); 
  t = ( t - p[2] )/3.0;
  c = t*t*t; 
  d = b*b - c;

  if( d >= 0. )
  {
    d = std::pow( (std::sqrt(d) + std::fabs(b) ), 1.0/3.0 );
    
    if( d != 0. )
    {
       if( b > 0. ) { b = -d; }
       else         { b =  d; }
       c =  t/b;
    }
    d       =  std::sqrt(0.75)*(b - c); 
    r[2][2] =  d; 
    b       =  b + c;
    c       = -0.5*b-x;
    r[1][2] =  c;

    if( ( b > 0. &&  x <= 0. ) || ( b < 0. && x > 0. ) )
    {
       r[1][1] =  c; 
       r[2][1] = -d; 
       r[1][3] =  b - x;
       r[2][3] =  0;
    }
    else
    {
       r[1][1] =  b - x; 
       r[2][1] =  0.; 
       r[1][3] =  c;
       r[2][3] = -d;
    }
  }              // end of 2 equal or complex roots 
  else
  {
    if( b == 0. ) { d =  std::atan(1.0)/1.5; }
    else          { d =  std::atan( std::sqrt(-d)/std::fabs(b) )/3.0; }

    if( b < 0. )  { b =  std::sqrt(t)*2.0; }
    else          { b = -2.0*std::sqrt(t); }

    c =  std::cos(d)*b; 
    t = -std::sqrt(0.75)*std::sin(d)*b - 0.5*c;
    d = -t - c - x; 
    c =  c - x; 
    t =  t - x;

    if( std::fabs(c) > std::fabs(t) ) { r[1][3] = c; }
    else
    {
      r[1][3] = t; 
      t       = c;
    }
    if( std::fabs(d) > std::fabs(t) ) { r[1][2] = d; }
    else
    {
      r[1][2] = t; 
      t       = d;
    }
    r[1][1] = t;

    for(k = 1; k < 4; k++ ) { r[2][k] = 0.; }
  }
  return 0;
}

//////////////////////////////////////////////////////////////////////////////
//
// Array r[3][5]  p[5]
// Roots of poly p[0] x^4 + p[1] x^3...+p[4]=0
// x=r[1][k] + i r[2][k]  k=1,...,4

G4int G4AnalyticalPolSolver::BiquadRoots( G4double p[5], G4double r[3][5] )
{
  G4double a, b, c, d, e;
  G4int i, k, j, noRoots;

  if(p[0] != 1.0)
  {
    for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
    p[0] = 1.;
  }
  e = 0.25*p[1];
  b = 2*e;
  c = b*b;
  d = 0.75*c;
  b = p[3] + b*( c - p[2] );
  a = p[2] - d;
  c = p[4] + e*( e*a - p[3] );
  a = a - d;

  p[1] = 0.5*a;
  p[2] = (p[1]*p[1]-c)*0.25;
  p[3] = b*b/(-64.0);

  if( p[3] < 0. )
  {
    noRoots = CubicRoots(p,r);

    for( k = 1; k < 4; k++ )
    {
      if( r[2][k] == 0. && r[1][k] > 0 )
      {
        d = r[1][k]*4; 
        a = a + d;

        if     ( a >= 0. && b >= 0.) { p[1] =  std::sqrt(d); }
        else if( a <= 0. && b <= 0.) { p[1] =  std::sqrt(d); }
        else                         { p[1] = -std::sqrt(d); }

        b = 0.5*( a + b/p[1] );

        p[2]    = c/b; 
        noRoots = QuadRoots(p,r);

        for( i = 1; i < 3; i++ )
        {
          for( j = 1; j < 3; j++ ) { r[j][i+2] = r[j][i]; }
        }
        p[1]    = -p[1]; 
        p[2]    =  b; 
        noRoots = QuadRoots(p,r);

        for( i = 1; i < 5; i++ ) { r[1][i] = r[1][i] - e; }

        return 4;
      }
    }
  }
  if( p[2] < 0. )
  {
    b    = std::sqrt(c); 
    d    = b + b - a;
    p[1] = 0.; 
    
    if( d > 0. ) { p[1] = std::sqrt(d); }
  }
  else
  {
    if( p[1] > 0.) { b =  std::sqrt(p[2])*2.0 + p[1]; }
    else           { b = -std::sqrt(p[2])*2.0 + p[1]; }

    if( b != 0.) { p[1] = 0; }
    else
    {
      for(k = 1; k < 5; k++ )
      {
          r[1][k] = -e;
          r[2][k] =  0;
      }
      return 0;
    }
  }

  p[2]    = c/b; 
  noRoots = QuadRoots(p,r);

  for( k = 1; k < 3; k++ )
  {
    for( j = 1; j < 3; j++ ) { r[j][k+2] = r[j][k]; }
  }
  p[1]    = -p[1]; 
  p[2]    =  b; 
  noRoots = QuadRoots(p,r);

  for( k = 1; k < 5; k++ ) { r[1][k] = r[1][k] - e; }

  return 4;
}

//////////////////////////////////////////////////////////////////////////////

G4int G4AnalyticalPolSolver::QuarticRoots( G4double p[5], G4double r[3][5])
{
  G4double a0, a1, a2, a3, y1;
  G4double R2, D2, E2, D, E, R = 0.;
  G4double a, b, c, d, ds;

  G4double reRoot[4];
  G4int k, noRoots, noReRoots = 0;
  
  for( k = 0; k < 4; k++ ) { reRoot[k] = DBL_MAX; }

  if( p[0] != 1.0 )
  {
    for( k = 1; k < 5; k++) { p[k] = p[k]/p[0]; }
    p[0] = 1.;
  }
  a3 = p[1];
  a2 = p[2];
  a1 = p[3];
  a0 = p[4];

  // resolvent cubic equation cofs:

  p[1] = -a2;
  p[2] = a1*a3 - 4*a0;
  p[3] = 4*a2*a0 - a1*a1 - a3*a3*a0;

  noRoots = CubicRoots(p,r);

  for( k = 1; k < 4; k++ )
  {
    if( r[2][k] == 0. ) // find a real root
    {
      noReRoots++;
      reRoot[k] = r[1][k];
    }
    else reRoot[k] = DBL_MAX; // kInfinity;
  }
  y1 = DBL_MAX; // kInfinity;  
  for( k = 1; k < 4; k++ )
  {
    if ( reRoot[k] < y1 ) { y1 = reRoot[k]; }
  }

  R2 = 0.25*a3*a3 - a2 + y1;
  b  = 0.25*(4*a3*a2 - 8*a1 - a3*a3*a3);
  c  = 0.75*a3*a3 - 2*a2;
  a  = c - R2;
  d  = 4*y1*y1 - 16*a0;

  if( R2 > 0.)
  {
    R = std::sqrt(R2);
    D2 = a + b/R;
    E2 = a - b/R;

    if( D2 >= 0. )
    {
      D       = std::sqrt(D2);
      r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
      r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
      r[2][1] = 0.;
      r[2][2] = 0.;
    }
    else
    {
      D       = std::sqrt(-D2);
      r[1][1] = -0.25*a3 + 0.5*R;
      r[1][2] = -0.25*a3 + 0.5*R;
      r[2][1] =  0.5*D;
      r[2][2] = -0.5*D;
    }
    if( E2 >= 0. )
    {
      E       = std::sqrt(E2);
      r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
      r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
      r[2][3] = 0.;
      r[2][4] = 0.;
    }
    else
    {
      E       = std::sqrt(-E2);
      r[1][3] = -0.25*a3 - 0.5*R;
      r[1][4] = -0.25*a3 - 0.5*R;
      r[2][3] =  0.5*E;
      r[2][4] = -0.5*E;
    }
  }
  else if( R2 < 0.)
  {
    R = std::sqrt(-R2);
    G4complex CD2(a,-b/R);
    G4complex CD = std::sqrt(CD2);

    r[1][1] = -0.25*a3 + 0.5*real(CD);
    r[1][2] = -0.25*a3 - 0.5*real(CD);
    r[2][1] =  0.5*R + 0.5*imag(CD);
    r[2][2] =  0.5*R - 0.5*imag(CD);
    G4complex CE2(a,b/R);
    G4complex CE = std::sqrt(CE2);

    r[1][3] = -0.25*a3 + 0.5*real(CE);
    r[1][4] = -0.25*a3 - 0.5*real(CE);
    r[2][3] =  -0.5*R + 0.5*imag(CE);
    r[2][4] =  -0.5*R - 0.5*imag(CE);
  }
  else // R2=0 case
  {
    if(d >= 0.)
    {
      D2 = c + std::sqrt(d);
      E2 = c - std::sqrt(d);

      if( D2 >= 0. )
      {
        D       = std::sqrt(D2);
        r[1][1] = -0.25*a3 + 0.5*R + 0.5*D;
        r[1][2] = -0.25*a3 + 0.5*R - 0.5*D;
        r[2][1] = 0.;
        r[2][2] = 0.;
      }
      else
      {
        D       = std::sqrt(-D2);
        r[1][1] = -0.25*a3 + 0.5*R;
        r[1][2] = -0.25*a3 + 0.5*R;
        r[2][1] =  0.5*D;
        r[2][2] = -0.5*D;
      }
      if( E2 >= 0. )
      {
        E       = std::sqrt(E2);
        r[1][3] = -0.25*a3 - 0.5*R + 0.5*E;
        r[1][4] = -0.25*a3 - 0.5*R - 0.5*E;
        r[2][3] = 0.;
        r[2][4] = 0.;
      }
      else
      {
        E       = std::sqrt(-E2);
        r[1][3] = -0.25*a3 - 0.5*R;
        r[1][4] = -0.25*a3 - 0.5*R;
        r[2][3] =  0.5*E;
        r[2][4] = -0.5*E;
      }
    }
    else
    {
      ds = std::sqrt(-d);
      G4complex CD2(c,ds);
      G4complex CD = std::sqrt(CD2);

      r[1][1] = -0.25*a3 + 0.5*real(CD);
      r[1][2] = -0.25*a3 - 0.5*real(CD);
      r[2][1] =  0.5*R + 0.5*imag(CD);
      r[2][2] =  0.5*R - 0.5*imag(CD);

      G4complex CE2(c,-ds);
      G4complex CE = std::sqrt(CE2);

      r[1][3] = -0.25*a3 + 0.5*real(CE);
      r[1][4] = -0.25*a3 - 0.5*real(CE);
      r[2][3] =  -0.5*R + 0.5*imag(CE);
      r[2][4] =  -0.5*R - 0.5*imag(CE);
    }  
  }
  return 4;
}

//
//
//////////////////////////////////////////////////////////////////////////////