Graphviz  2.41.20170921.2350
solvers.c
Go to the documentation of this file.
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 
33 int 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 
76 int 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 
100 int 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 }
#define M_PI
Definition: solvers.c:27
#define cbrt(x)
Definition: solvers.c:24
int solve2(double *coeff, double *roots)
Definition: solvers.c:76
int solve3(double *coeff, double *roots)
Definition: solvers.c:33
#define AEQ0(x)
Definition: solvers.c:31
#define alpha
Definition: shapes.c:3902
int solve1(double *coeff, double *roots)
Definition: solvers.c:100