Source code: com/port80/graph/dot/impl/CubicSolver.java
1 //
2 // Copyright(c) 2002, Chris Leung
3 //
4
5 package com.port80.graph.dot.impl;
6
7 /**
8 */
9 public class CubicSolver {
10
11 /** @return number of roots.
12 */
13 public static int solve3(double[] coeff, double[] roots) {
14
15 if (almost0(coeff[3])) return solve2(coeff, roots);
16
17 int nroot;
18 double a = coeff[3];
19 double b = coeff[2];
20 double c = coeff[1];
21 double d = coeff[0];
22 double b3a = b / (3 * a);
23 double ca = c / a;
24 double da = d / a;
25 double p = b3a * b3a;
26 double q = 2 * b3a * p - b3a * ca + da;
27 p = ca / 3 - p;
28 double disc = q * q + 4 * p * p * p;
29
30 if (disc < 0) {
31 double r = 0.5 * Math.sqrt( -disc + q * q);
32 double theta = Math.atan2(Math.sqrt( -disc), -q);
33 double temp = 2 * cubicRoot(r);
34 roots[0] = temp * Math.cos(theta / 3);
35 roots[1] = temp * Math.cos((theta + 2 * Math.PI) / 3);
36 roots[2] = temp * Math.cos((theta - 2 * Math.PI) / 3);
37 nroot = 3;
38 } else {
39 double alpha = 0.5 * (Math.sqrt(disc) - q);
40 double beta = -q - alpha;
41 roots[0] = cubicRoot(alpha) + cubicRoot(beta);
42 if (disc > 0) nroot = 1;
43 else {
44 roots[2] = roots[1] = -0.5 * roots[0];
45 nroot = 3;
46 }
47 }
48 for (int i = 0; i < nroot; i++) roots[i] -= b3a;
49 return nroot;
50 }
51
52 private static int solve2(double[] coeff, double[] roots) {
53
54 if (almost0(coeff[2])) return solve1(coeff, roots);
55
56 double a = coeff[2];
57 double b = coeff[1];
58 double c = coeff[0];
59 double b2a = b / (2 * a);
60 double disc = b2a * b2a - c / a;
61
62 if (disc < 0) return 0;
63 else if (disc == 0) {
64 roots[0] = -b2a;
65 return 1;
66 } else {
67 roots[0] = -b2a + Math.sqrt(disc);
68 roots[1] = -2 * b2a - roots[0];
69 return 2;
70 }
71 }
72
73 private static int solve1(double[] coeff, double[] roots) {
74
75 double a = coeff[1];
76 double b = coeff[0];
77 if (almost0(a)) {
78 if (almost0(b)) return 4;
79 else return 0;
80 }
81 roots[0] = -b / a;
82 return 1;
83 }
84
85 ////////////////////////////////////////////////////////////////////////
86
87 private static double cubicRoot(double x) {
88 if (x < 0) return -1*Math.pow( -x, 1.0 / 3.0);
89 else return Math.pow(x, 1.0 / 3.0);
90 }
91
92 private static boolean almost0(double x) {
93 return (x < 1e-7) && (x > -1e-7);
94 }
95
96 ////////////////////////////////////////////////////////////////////////
97
98 }