fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <complex.h>
  4. #include <tgmath.h>
  5.  
  6. /*
  7.  * Configuration
  8.  */
  9. #define DR double
  10. typedef complex DR D;
  11.  
  12. // Enable this to let you select which of the 3 complex cube roots to use
  13. // #define TEST_CUBEROOT
  14.  
  15. /*
  16.  * Output
  17.  */
  18. #define DUMP(var) dump(#var, var)
  19.  
  20. void dumpreal(char const *name, DR var)
  21. {
  22. 	printf("%s = %f\n", name, var);
  23. }
  24.  
  25. void dump(char const *name, D var)
  26. {
  27. 	if ( cabs(cimag(var)) < 0.00000001 )
  28. 		dumpreal(name, creal(var));
  29. 	else
  30. 		printf("%s = (%f,%f)\n", name, creal(var), cimag(var));
  31. }
  32.  
  33. /*
  34.  * Cube root calculation
  35.  */
  36. #ifdef TEST_CUBEROOT
  37. void ccbrt(D var, D roots[static 3])
  38. {
  39. 	DR mod = sqrt( creal(var) * creal(var) + cimag(var) * cimag(var) );
  40. 	DR exp = atan2(cimag(var), creal(var));
  41. 	DR pi = acos(-1.0);
  42.  
  43. 	//DUMP(var); DUMP(mod); DUMP(exp);
  44. 	//dump("pow result", pow(var, 1./3));
  45.  
  46. // atan2 computes in interval [-pi, pi]; convert to [0, 2pi]
  47. 	if ( exp < 0 )
  48. 		exp += 2*pi;
  49.  
  50. 	//DUMP(exp);
  51.  
  52. // Three cube roots
  53. 	for (int r = 0; r < 3; ++r)
  54. 	{
  55. 		DR cr_mod = cbrt(mod);
  56. 		DR cr_exp = (exp + r*2*pi) / 3;
  57.  
  58. 		roots[r] = cr_mod * (cos(cr_exp) + I*sin(cr_exp));
  59.  
  60. 		//DUMP(cr_mod);
  61. 		//DUMP(cr_exp);
  62. 		//DUMP(roots[r]);
  63. 	}
  64. 	//exit(0);
  65. }
  66. #endif
  67.  
  68. /*
  69.  * Quartic solver using Wikipedia algorithm
  70.  */
  71. int main()
  72. {
  73. // Inputs from http://m...content-available-to-author-only...e.com/questions/1796015/need-help-solving-x4-3x3-11x23x10-0/
  74. 	D a = 1, b = -3, c = -11, d = 3, e = 10;
  75.  
  76. // Solve
  77. 	D p = (8*a*c - 3*b*b) / (8*a*a);
  78. 	D q = (b*b*b - 4*a*b*c + 8*a*a*d) / (8*a*a*a);
  79.  
  80. 	D d0 = c*c - 3*b*d + 12*a*e;
  81. 	D d1 = 2*c*c*c - 9*b*c*d + 27*b*b*e + 27*a*d*d - 72*a*c*e;
  82.  
  83. 	D disc = d1*d1 - 4*d0*d0*d0;
  84. 	D discsq = sqrt(disc);		// can negate this
  85.  
  86. 	D Q3 = (d1 + discsq) / 2;
  87. 	D Q = pow(Q3, 1./3);
  88.  
  89. #ifdef TEST_CUBEROOT
  90. 	D Qx[3];
  91. 	ccbrt(Q3, Qx);
  92. 	Q = Qx[1];		// Choose 0,1,2 here to taste
  93. #endif
  94.  
  95. 	D S2 = -p*2/3 + (Q + d0/Q) / (3*a);
  96. 	D S = 0.5 * sqrt(S2);
  97.  
  98. 	D qdp = -4*S*S - 2*p + q/S;
  99. 	D qdpsq = sqrt(qdp);
  100.  
  101. 	D qdn = -4*S*S - 2*p - q/S;
  102. 	D qdnsq = sqrt(qdn);
  103.  
  104. 	D x1 = -b/(4*a) - S + qdpsq/2;
  105. 	D x2 = -b/(4*a) - S - qdpsq/2;
  106. 	D x3 = -b/(4*a) + S + qdnsq/2;
  107. 	D x4 = -b/(4*a) + S - qdnsq/2;
  108.  
  109. 	DUMP(p);
  110. 	DUMP(q);
  111. 	DUMP(d0);
  112. 	DUMP(d1);
  113. 	DUMP(disc);
  114. 	DUMP(discsq);
  115. 	DUMP(Q3);
  116. 	DUMP(Q);
  117. #ifdef TEST_CUBEROOT
  118. 	DUMP(Qx[0]);
  119. 	DUMP(Qx[1]);
  120. 	DUMP(Qx[2]);
  121. #endif
  122. 	DUMP(S2);
  123. 	DUMP(S);
  124. 	DUMP(qdp);
  125. 	DUMP(qdn);
  126. 	DUMP(x1);
  127. 	DUMP(x2);
  128. 	DUMP(x3);
  129. 	DUMP(x4);
  130. }
  131.  
Success #stdin #stdout 0s 2168KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
p = -14.375000
q = -16.875000
d0 = 268.000000
d1 = 7040.000000
disc = -27433728.000000
discsq = (0.000000,5237.721642)
Q3 = (3520.000000,2618.860821)
Q = (16.000000,3.464102)
S2 = 20.250000
S = 2.250000
qdp = 1.000000
qdn = 16.000000
x1 = -1.000000
x2 = -2.000000
x3 = 5.000000
x4 = 1.000000
https://ideone.com/6bigDf
language:
C99 (gcc 8.3)
created:
7 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!