1/* $Id$ $Revision$ */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14
15#include "config.h"
16#include <math.h>
17#include "solvers.h"
18
19#ifdef DMALLOC
20#include "dmalloc.h"
21#endif
22
23#ifndef HAVE_CBRT
24#define cbrt(x) ((x < 0) ? (-1*pow(-x, 1.0/3.0)) : pow (x, 1.0/3.0))
25#endif
26#ifndef M_PI
27#define M_PI 3.14159265358979323846
28#endif
29
30#define EPS 1E-7
31#define AEQ0(x) (((x) < EPS) && ((x) > -EPS))
32
33int solve3(double *coeff, double *roots)
34{
35 double a, b, c, d;
36 int rootn, i;
37 double p, q, disc, b_over_3a, c_over_a, d_over_a;
38 double r, theta, temp, alpha, beta;
39
40 a = coeff[3], b = coeff[2], c = coeff[1], d = coeff[0];
41 if (AEQ0(a))
42 return solve2(coeff, roots);
43 b_over_3a = b / (3 * a);
44 c_over_a = c / a;
45 d_over_a = d / a;
46
47 p = b_over_3a * b_over_3a;
48 q = 2 * b_over_3a * p - b_over_3a * c_over_a + d_over_a;
49 p = c_over_a / 3 - p;
50 disc = q * q + 4 * p * p * p;
51
52 if (disc < 0) {
53 r = .5 * sqrt(-disc + q * q);
54 theta = atan2(sqrt(-disc), -q);
55 temp = 2 * cbrt(r);
56 roots[0] = temp * cos(theta / 3);
57 roots[1] = temp * cos((theta + M_PI + M_PI) / 3);
58 roots[2] = temp * cos((theta - M_PI - M_PI) / 3);
59 rootn = 3;
60 } else {
61 alpha = .5 * (sqrt(disc) - q);
62 beta = -q - alpha;
63 roots[0] = cbrt(alpha) + cbrt(beta);
64 if (disc > 0)
65 rootn = 1;
66 else
67 roots[1] = roots[2] = -.5 * roots[0], rootn = 3;
68 }
69
70 for (i = 0; i < rootn; i++)
71 roots[i] -= b_over_3a;
72
73 return rootn;
74}
75
76int solve2(double *coeff, double *roots)
77{
78 double a, b, c;
79 double disc, b_over_2a, c_over_a;
80
81 a = coeff[2], b = coeff[1], c = coeff[0];
82 if (AEQ0(a))
83 return solve1(coeff, roots);
84 b_over_2a = b / (2 * a);
85 c_over_a = c / a;
86
87 disc = b_over_2a * b_over_2a - c_over_a;
88 if (disc < 0)
89 return 0;
90 else if (disc == 0) {
91 roots[0] = -b_over_2a;
92 return 1;
93 } else {
94 roots[0] = -b_over_2a + sqrt(disc);
95 roots[1] = -2 * b_over_2a - roots[0];
96 return 2;
97 }
98}
99
100int solve1(double *coeff, double *roots)
101{
102 double a, b;
103
104 a = coeff[1], b = coeff[0];
105 if (AEQ0(a)) {
106 if (AEQ0(b))
107 return 4;
108 else
109 return 0;
110 }
111 roots[0] = -b / a;
112 return 1;
113}
114