roots.c

Go to the documentation of this file.

/******************************************************************************/
#include "tpcclibConfig.h"
/******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/******************************************************************************/
#include "tpcextensions.h"
/******************************************************************************/
#include "tpclinopt.h"
/******************************************************************************/
 
/******************************************************************************/
int rootsCubic(
  const double a,
  const double b,
  const double c,
  double *x1,
  double *x2,
  double *x3
) {
  if(x1==NULL || x2==NULL || x3==NULL) return(0);
  *x1=*x2=*x3=nan("");
 
  if(a==0 && b==0 && c==0) return(0);
  if(!isfinite(a) || !isfinite(b) || !isfinite(c)) return(0);
 
  long double q=(a*a - 3*b);
  long double r=(2*a*a*a - 9*a*b + 27*c);
 
  long double Q=q/9;
  long double R=r/54;
 
  long double Q3=Q*Q*Q;
  long double R2=R*R;
 
  long double CR2=729*r*r;
  long double CQ3=2916*q*q*q;
 
  if(R==0 && Q==0) {
    *x1=-a/3; // *x2=-a/3; *x3=-a/3;
    return(1);
  } else if(CR2==CQ3) {
    long double sqrtQ=sqrt(Q);
    if(R>0) {
      *x1=-2*sqrtQ - a/3;
      *x2=sqrtQ - a/3; // *x3=sqrtQ - a/3;
    } else {
      *x1=-sqrtQ - a/3; // *x2=-sqrtQ - a/3; *x3=2*sqrtQ - a/3;
      *x2=2*sqrtQ - a/3;
    }
    return(2);
  } else if(R2<Q3) {
    double ratio=copysign(sqrt(R2/Q3), R);
    //if(!isfinite(ratio)) return(0);
    double theta=acos(ratio);
    double norm=-2*sqrt(Q);
    //if(!isfinite(norm)) return(0);
    *x1=norm*cos(theta/3) - a/3;
    *x2=norm*cos((theta + 2.0*M_PI)/3) - a/3;
    *x3=norm*cos((theta - 2.0*M_PI)/3) - a/3;
    if(*x1>*x2) {double s=*x1; *x1=*x2; *x2=s;}
    if(*x2>*x3) {
      double s=*x2; *x2=*x3; *x3=s;
      if(*x1>*x2) {double s=*x1; *x1=*x2; *x2=s;}
    }
    return(3);
  } else {
    long double A = -copysign(pow(fabsl(R) + sqrt(R2-Q3), 1.0/3.0), R);
    long double B=Q/A;
    *x1=A + B - a/3;
    return(1);
  }
}
/******************************************************************************/
 
/******************************************************************************/
int rootsQuadratic(
  const double a,
  const double b,
  const double c,
  double *x1,
  double *x2
) {
  if(x1==NULL || x2==NULL) return(0);
  *x1=*x2=nan("");
 
  if(!isfinite(a) || !isfinite(b) || !isfinite(c)) return(0);
 
  if(a==0) {
    if(b==0) return(0);
    else {*x1=-c/b; /* *x2=-c/b;*/ return(1);}
  }
  long double d=b*b - 4*a*c;
  if(d>0) {
    if(b==0) {
      long double r=fabs(sqrt(-c/a));
      *x1=-r; *x2=r;
    } else {
      long double t=-0.5*(b + copysign(sqrt(d), b));
      long double r1=t/a; 
      long double r2=c/t;
      if(r1<r2) {*x1=r1; *x2=r2;} else {*x1=r2; *x2=r1;}
    }
    return(2);
  } else if(d==0) {
    *x1=-0.5*b/a; // *x2=-0.5*b/a;
    return(1);
  } else {
    return(0);
  }
}
/******************************************************************************/
 
/******************************************************************************/