fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define M_PI					3.141592653589793238
  5. #define MATH_2PI                (2*M_PI)
  6. #define sgn_float(X)            ((X>0)-(X<0))   // returns the sign of X (1 or -1); for floats, doubles
  7.  
  8. double roots_cubic_constrained(double pc[]){
  9.     double a,b,c;               // coefficients of input polynomial
  10.     double Q, R, theta;         // intermediate components of the numerical solution
  11.     double A, B;                // intermediate value needed for the second case R^2>=Q^3
  12.     double x;                   // solution
  13.  
  14.     // applies to cubic equation x3 + ax2 + bx + c = 0
  15.     if(pc[0]){                  // ensure a cubic, so that we don't divide by 0
  16.         // pc[0] = 1;
  17.         a = pc[1] / pc[0];     // divide all coeffs by pc[0] to fit the form (making the coeff of x^3 to be 1)
  18.         b = pc[2] / pc[0];
  19.         c = pc[3] / pc[0];
  20.  
  21.         // Q = (pc[1] * pc[1] - 3 * pc[2])/9;              // component for the solution
  22.         // R = (2*pc[1]^3 - 9*px[1]*pc[2] + 27*pc[3])/54;
  23.         Q = (a*a - 3*b)/9;              // (5.6.10)
  24.         R = (2*a*a*a - 9*a*b + 27*c)/54;  // (5.6.10)
  25.  
  26.         if((R*R) < (Q*Q*Q)){                    // if true, then there are 3 real roots
  27.             theta = acos(R/sqrt(Q*Q*Q));
  28.  
  29.             x = -2*sqrt(Q)*cos(theta/3) - a/3;      // (5.6.12)
  30.             // if(0.99*(-8) < x & 1.01*(-4) > x){      // TODO use the fmcf floating pt comparison functions
  31.             if(-8 < x && -4 > x){      //  use float comparison becasue exact comparison does not matter
  32.                 return x;
  33.             }
  34.  
  35.             x = -2*sqrt(Q)*cos((theta + MATH_2PI)/3) - a/3; // (5.6.12)
  36.             // if(0.99*(-8) < x & 1.01*(-4) > x){
  37.             if(-8 < x && -4 > x){      //  use float comparison becasue exact comparison does not matter      
  38.                 return x;
  39.             }
  40.  
  41.             x = -2*sqrt(Q)*cos((theta + MATH_2PI)/3) - a/3; // (5.6.12)
  42.             return x;                   // Paolo assures that there will be 1 root in the acceptable range
  43.  
  44.         }
  45.         else{                   // otherwise compute A (since Q and R are both real (since a,b,c are Real))
  46.             // A = -sgn(R)*[|R| + sqrt(R^2-Q^3)]^(1/3)
  47.             // A = fabs(R);                            // in order to compute sgn(R) via (R/A)
  48.             // A = -(R/A)*(A+sqrt(R^2-Q^3))^(1/3);     // (5.6.15)
  49.             A = (R>0)-(R<0);                        // sign of R (sgn(R))
  50.             A = sgn_float(R);                       // sign of R
  51.             A = -(A)*cbrt(fabs(R)+sqrt(R*R-Q*Q*Q)); // (5.6.15)
  52.             if(0 == A){                             // (5.6.16)
  53.                 B = 0;
  54.             }
  55.             else{
  56.                 B = Q/A;
  57.             }
  58.             x = (A+B)-(a/3);                        // (5.6.17) the single real root.
  59.             return x;
  60.         }
  61.     }
  62.     else {                      // use the quadratic equation
  63.         a = pc[1];
  64.         b = pc[2];
  65.         c = pc[3];
  66.  
  67. 		printf ("%f, %f, %f\n", a,b,c);
  68.  
  69.         R = b*b-4*a*c;          // Use R as the discriminant
  70.         if(R >= 0){             // Should always be this case - otherwise there are no real roots
  71.             Q = -(b + sgn_float(b)*sqrt(R))/2;    // (5.6.4) with 'Q' substituted for 'q'
  72.             printf("Q = %f\n", Q);
  73.             x = Q/a;            // (5.6.5)
  74. 			printf ("first %f\n", x);
  75.             if(-8 < x && -4 > x){      //  use float comparison becasue exact comparison does not matter      
  76.                 return x;
  77.             }
  78.             x = c/Q;            // (5.6.5)
  79. 			printf ("second %f\n", x);
  80.         return x;
  81.         }
  82.         // else{
  83.         //     print("ERRGG!!!")
  84.         // }
  85.     }
  86.  
  87. }
  88.  
  89.  
  90. int main(void) {
  91. 	// your code goes here
  92. 	double temp[5];
  93. 	double res;
  94.  
  95. 	// temp[0] = 100.123;
  96. 	// temp[1] = .1234124;
  97. 	// temp[2] = .000234567891234;
  98. 	// temp[3] = .000000031415926535;
  99. 	// temp[4] = 492314.1;
  100.  
  101. 	temp[0] = 1;
  102. 	temp[1] = -5;
  103. 	temp[2] = -21;
  104. 	temp[3] = 145;
  105.  
  106. 	// res = pfit(&temp[2]);
  107. 	// res = pfit((double *)&temp[2]);
  108. 	res = roots_cubic_constrained(temp);
  109. 	printf("ans: %6.6f", res);
  110. 	// printf("main \n%3.3lg  %3.3g %3.3g \n", temp[0], temp[1], temp[2]);
  111.  
  112. 	return 0;
  113. }
  114.  
  115.  
Success #stdin #stdout 0s 9424KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
ans: -5.000000
https://ideone.com/M7NHuJ
language:
C (gcc 8.3)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!