Important alert: (current site time 4/18/2014 2:10:55 PM EDT)

VB icon

Solve cubic equations

Submitted on: 7/27/2000 8:43:30 PM
By: Bob Stout (republished under Open Content License)  
Level: Intermediate
User Rating: Unrated
Compatibility: C, C++ (general), Microsoft Visual C++, Borland C++, UNIX C++
Views: 4848
     CUBIC.C - Solve a cubic polynomial

INCLUDE files:

Can't Copy and Paste this?
Click here for a copy-and-paste friendly version of this code!
//INCLUDE files for :Solve cubic equations
/* +++Date last modified: 05-Jul-1997 */
** SNIPMATH.H - Header file for SNIPPETS math functions and macros
#ifndef SNIPMATH__H
#define SNIPMATH__H
#include <math.h>
#include "sniptype.h"
#include "round.h"
** Callable library functions begin here
voidSetBCDLen(int n); /* Bcdl.C */
longBCDtoLong(char *BCDNum); /* Bcdl.C */
voidLongtoBCD(long num, char BCDNum[]);/* Bcdl.C */
double bcd_to_double(void *buf, size_t len, /* Bcdd.C */
 int digits);
int double_to_bcd(double arg, char *buf, /* Bcdd.C */
 size_t length, size_t digits );
DWORDncomb1 (int n, int m);/* Combin.C*/
DWORDncomb2 (int n, int m);/* Combin.C*/
voidSolveCubic(double a, double b, double c, /* Cubic.C*/
 double d, int *solutions,
 double *x);
DWORDdbl2ulong(double t); /* Dbl2Long.C */
longdbl2long(double t);/* Dbl2Long.C */
double dround(double x); /* Dblround.C */
/* Use #defines for Permutations and Combinations -- Factoryl.C */
#define log10P(n,r) (log10factorial(n)-log10factorial((n)-(r)))
#define log10C(n,r) (log10P((n),(r))-log10factorial(r))
double log10factorial(double N); /* Factoryl.C */
double fibo(unsigned short term);/* Fibo.C */
double frandom(int n);/* Frand.C*/
double ipow(double x, int n);/* Ipow.C */
int ispow2(int x);/* Ispow2.C*/
longdouble ldfloor(long double a);/* Ldfloor.C */
int initlogscale(long dmax, long rmax);/* Logscale.C */
longlogscale(long d); /* Logscale.C */
floatMSBINToIEEE(float f); /* Msb2Ieee.C */
floatIEEEToMSBIN(float f); /* Msb2Ieee.C */
int perm_index (char pit[], int size);/* Perm_Idx.C */
int round_div(int n, int d); /* Rnd_Div.C */
longround_ldiv(long n, long d);/* Rnd_Div.C */
double rad2deg(double rad); /* Rad2Deg.C */
double deg2rad(double deg); /* Rad2Deg.C */
#include "pi.h"
#ifndef PHI
 #define PHI ((1.0+sqrt(5.0))/2.0) /* the golden number*/
 #define INV_PHI (1.0/PHI) /* the golden ratio */
** File: ISQRT.C
struct int_sqrt {
 unsigned sqrt,
void usqrt(unsigned long x, struct int_sqrt *q);
#endif /* SNIPMATH__H */
Can't Copy and Paste this?
Click here for a copy-and-paste friendly version of this code!
Terms of Agreement:   
By using this code, you agree to the following terms...   
  1. You may use this code in your own programs (and may compile it into a program and distribute it in compiled format for languages that allow it) freely and with no charge.
  2. You MAY NOT redistribute this code (for example to a web site) without written permission from the original author. Failure to do so is a violation of copyright laws.   
  3. You may link to this code from another website, but ONLY if it is not wrapped in a frame. 
  4. You will abide by any additional copyright restrictions which the author may have placed in the code or code's description.
// Name: Solve cubic equations
// Description:CUBIC.C - Solve a cubic polynomial
// By: Bob Stout (republished under Open Content License)
//This code is copyrighted and has// limited warranties.Please see details.//**************************************

/* +++Date last modified: 05-Jul-1997 */
** CUBIC.C - Solve a cubic polynomial
** public domain by Ross Cottrell
#include <math.h>
#include <stdlib.h>
#include "snipmath.h"
void SolveCubic(double a,
double b,
double c,
double d,
double *x)
 long doublea1 = b/a, a2 = c/a, a3 = d/a;
 long doubleQ = (a1*a1 - 3.0*a2)/9.0;
 long double R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0;
 doubleR2_Q3 = R*R - Q*Q*Q;
 if (R2_Q3 <= 0)
*solutions = 3;
theta = acos(R/sqrt(Q*Q*Q));
x[0] = -2.0*sqrt(Q)*cos(theta/3.0) - a1/3.0;
x[1] = -2.0*sqrt(Q)*cos((theta+2.0*PI)/3.0) - a1/3.0;
x[2] = -2.0*sqrt(Q)*cos((theta+4.0*PI)/3.0) - a1/3.0;
*solutions = 1;
x[0] = pow(sqrt(R2_Q3)+fabs(R), 1/3.0);
x[0] += Q/x[0];
x[0] *= (R < 0.0) ? 1 : -1;
x[0] -= a1/3.0;
#ifdef TEST
int main(void)
 double a1 = 1.0, b1 = -10.5, c1 = 32.0, d1 = -30.0;
 double a2 = 1.0, b2 = -4.5, c2 = 17.0, d2 = -30.0;
 double x[3];
 int solutions;
 SolveCubic(a1, b1, c1, d1, &solutions, x);
 /* should get 3 solutions: 2, 6 & 2.5*/
 SolveCubic(a2, b2, c2, d2, &solutions, x);
 /* should get 1 solution: 2.5*/
 return 0;
#endif /* TEST */

Other 155 submission(s) by this author


Report Bad Submission
Use this form to tell us if this entry should be deleted (i.e contains no code, is a virus, etc.).
This submission should be removed because:

Your Vote

What do you think of this code (in the Intermediate category)?
(The code with your highest vote will win this month's coding contest!)
Excellent  Good  Average  Below Average  Poor (See voting log ...)

Other User Comments

10/1/2000 6:02:55 AMHasan Alasadi

there is an error in file snipmath.h :
near the middle in the line :

longdouble ldfloor(long double a);/* Ldfloor.C */

(If this comment was disrespectful, please report it.)

5/1/2002 11:16:26 PMMitch Foral

I dont believe this solves for imaginary roots.....correct me if i'm wrong plz.
(If this comment was disrespectful, please report it.)


Add Your Feedback
Your feedback will be posted below and an email sent to the author. Please remember that the author was kind enough to share this with you, so any criticisms must be stated politely, or they will be deleted. (For feedback not related to this particular code, please click here instead.)

To post feedback, first please login.