roots.c File Reference

Finding the real roots of equations. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "tpcextensions.h"
#include "tpclinopt.h"

Go to the source code of this file.

Functions
int	rootsCubic (const double a, const double b, const double c, double *x1, double *x2, double *x3)
int	rootsQuadratic (const double a, const double b, const double c, double *x1, double *x2)

Detailed Description

Finding the real roots of equations.

Copyright

(c) Turku PET Centre

Author

Vesa Oikonen

Definition in file roots.c.

Function Documentation

rootsCubic()

int rootsCubic	(	const double	a,
		const double	b,
		const double	c,
		double *	x1,
		double *	x2,
		double *	x3 )

Find the real roots of cubic equation x^3 + a x^2 + b x + c = 0.

The number of real roots can be 1 or 3, but coincident roots are possible. Root values, excluding duplicates and triplicates, will be placed in pointers x1, x2, and x3 in ascending order, or set to NaN.

Based on GNU Scientific Library function gsl_poly_solve_cubic.

See also

rootsQuadratic, fitLine

Returns

Returns the number of real individual roots, and zero in case of an error.

Parameters

a	Input: Parameter a of the cubic equation.
b	Input: Parameter b of the cubic equation.
c	Input: Parameter c of the cubic equation.
x1	Output: Pointer for root 1 value. Will be set to NaN if there are no real roots.
x2	Output: Pointer for root 2 value. Will be set to NaN if there is only one individual real root.
x3	Output: Pointer for root 3 value. Will be set to NaN if there are only two individual real roots.

Definition at line 30 of file roots.c.

  {
  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);
  }
}

rootsQuadratic()

int rootsQuadratic	(	const double	a,
		const double	b,
		const double	c,
		double *	x1,
		double *	x2 )

Find the real roots of quadratic equation a x^2 + b x + c = 0.

The number of real roots can be 0, 1, or 2, but coincident roots are possible. Root values, excluding duplicate value, will be placed in pointers x1 and x2 in ascending order, or set to NaN.

Based on GNU Scientific Library function gsl_poly_solve_quadratic.

See also

rootsCubic, fitLine

Returns

Returns the number of real roots, excluding duplicate solutions.

Parameters

a	Input: Parameter a of the quadratic equation.
b	Input: Parameter b of the quadratic equation.
c	Input: Parameter c of the quadratic equation.
x1	Output: Pointer for root 1 value. Will be set to NaN if there are no real roots.
x2	Output: Pointer for root 2 value. Will be set to NaN if there is only one individual real root.

Definition at line 111 of file roots.c.

  {
  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);
  }
}