1 /**********
2 *
3 * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
4 * University of Manchester, UK. All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without modification,
7 * are permitted provided that the following conditions are met:
8 *
9 * . Redistributions of source code must retain the above copyright notice,
10 * this list of conditions and the following disclaimer.
11 *
12 * . Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 *
16 * . Neither the name of the University of Manchester nor the names of its
17 * contributors may be used to endorse or promote products derived from this
18 * software without specific prior written permission.
19 *
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 *
33 **********
34 *
35 * Program : TINA
36 * File : $Source: /home/tina/cvs/tina-libs/tina/math/mathNum_eqn.c,v $
37 * Date : $Date: 2003/09/23 11:32:20 $
38 * Version : $Revision: 1.4 $
39 * CVS Id : $Id: mathNum_eqn.c,v 1.4 2003/09/23 11:32:20 matts Exp $
40 *
41 * Author : Legacy TINA
42 *
43 * Notes : Quadratic & cubic roots
44 *
45 *********
46 */
47 /**
48 * @file
49 * @brief Routines for finding the quadratic and cubic roots of a function.
50 *
51 *
52 *
53 */
54
55 #include "mathNum_eqn.h"
56
57 #if HAVE_CONFIG_H
58 #include <config.h>
59 #endif
60
61 #include <math.h>
62 #include <tina/math/math_GenDef.h>
63
64 Bool quadratic_roots(double a, double b, double c, double *x1, double *x2)
65 {
66 double disc = b * b - 4.0 * a * c;
67
68 /** assumes non-degenerate quadratic **/
69
70 if (b < 0)
71 {
72 a *= -1.0;
73 b *= -1.0;
74 c *= -1.0;
75 }
76 if (disc >= 0)
77 { /**real roots **/
78 double q = -0.5 * (b + sqrt(disc));
79
80 *x1 = q / a;
81 *x2 = c / q;
82 return (true);
83 } else
84 {
85 *x1 = -b / (2.0 * a); /** real part **/
86 *x2 = sqrt(-disc) / (2.0 * a); /** imag part **/
87 return (false);
88 }
89 }
90
91 Bool cubic_roots(double a, double b, double c, double d, double *x1, double *x2, double *x3)
92 {
93 double disc, q, r;
94
95 /** assumes non-degenerate cubic **/
96 b /= a;
97 c /= a;
98 d /= a;
99
100 q = (b * b - 3.0 * c) / 9.0;
101 r = (2.0 * b * b * b - 9.0 * b * c + 27.0 * d) / 54.0;
102 disc = q * q * q - r * r;
103
104 if (disc >= 0)
105 {
106 double theta;
107 double sq = sqrt(q), b3 = b / 3.0, sq3;
108
109 sq3 = sq * sq * sq;
110 if (sq3 == 0.0)
111 theta = 0.0; /**arbitrary**/
112 else
113 theta = acos(r / sq3);
114
115 *x1 = -2.0 * sq * cos(theta / 3.0) - b3;
116 *x2 = -2.0 * sq * cos((theta + TWOPI) / 3.0) - b3;
117 *x3 = -2.0 * sq * cos((theta + 2.0 * TWOPI) / 3.0) - b3;
118 return (true);
119 } else
120 {
121 double s = pow(fabs(r) + sqrt(-disc), 1.0 / 3.0);
122 double b3 = b / 3.0;
123 double x;
124
125 if (r >= 0.0)
126 x = -(s + q / s);
127 else
128 x = (s + q / s);
129
130 x -= b3;
131
132 b += x;
133 c += x * b;
134
135 *x1 = x; /** real root **/
136 return (quadratic_roots(1.0, b, c, x2, x3)); /** complex roots **/
137 }
138 }
139
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.