1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
/*
* Command-line numerical polynomial equation solver using the GNU Scientific Library (GSL).
* GSL is available from: "http://www.gnu.org/software/gsl/".
*
* Copyright (C) 2008-2011 George Gesslein II
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
The chief copyright holder can be contacted at gesslein@mathomatic.org, or
George Gesslein II, P.O. Box 224, Lansing, NY 14882-0224 USA.
*/
/*
* To compile:
./compile.roots
* or
gcc -O3 -Wall -o roots roots.c -lgsl -lgslcblas
* This program nicely solves any degree polynomial
* with real coefficients by calling the
* GSL function gsl_poly_complex_solve().
* This file is also useful for testing
* and as an example of using the GSL from C.
* The results are double precision floating point values
* that are sometimes accurate to 14 digits.
* Complex numbers are output if successful.
* Here is an example of it solving a cubic polynomial:
$ roots 1 1 1 1
The 3 approximate floating point solutions of:
+x^3 +x^2 +x +1 = 0
are:
x = +0 +1*i
x = +0 -1*i
x = -1
$
Try this for a large amount of error:
roots -1 -1 1 1
Proof the GSL isn't perfect.
*/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <float.h>
#include <math.h>
#include <errno.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_errno.h>
#define EPSILON 0.00000000000005 /* a good value for doubles */
#define HELP 1 /* display helpful text */
void display_root(int i);
char *prog_name = "roots"; /* name of this program */
double *a, *z; /* input and output arrays, respectively */
int precision = DBL_DIG - 1; /* display precision, it is not useful to set this higher than 15 */
void
usage(void)
{
printf("\n%s version 1.0 - numerical polynomial equation solver\n", prog_name);
printf("Uses the GNU Scientific Library (GSL).\n");
printf("\nSolves polynomial = 0 when given all real coefficients of the polynomial.\n");
printf("Double precision floating point math is used, accurate to about 14 digits.\n");
printf("\nUsage: %s highest-power-coefficient ... constant-term\n", prog_name);
printf("\nThe coefficients must be decimal, floating point, real numbers.\n");
printf("For example, if 4 real numbers are given, there will be 3 complex number\n");
printf("results or \"roots\" that are all valid solutions to polynomial = 0.\n");
}
int
main(int argc, char **argv)
{
int i, highest_power;
char *cp, *arg;
gsl_poly_complex_workspace *w;
if (argc <= 1) {
fprintf(stderr, "%s: The polynomial coefficients must be specified on the command line.\n", prog_name);
usage();
exit(2);
}
highest_power = argc - 2;
a = calloc(highest_power + 1, sizeof(double)); /* allocate real double input array */
z = calloc(2 * highest_power, sizeof(double)); /* allocate complex double output array */
/* parse the command line into the coefficient array a[] */
for (i = 0; i < argc - 1; i++) {
arg = argv[argc-i-1];
errno = 0;
a[i] = strtod(arg, &cp);
if (arg == cp || *cp) {
fprintf(stderr, "%s: Argument \"%s\" is not a floating point number.\n", prog_name, arg);
usage();
exit(2);
}
if (errno) {
fprintf(stderr, "%s: Argument \"%s\" is out of range.\n", prog_name, arg);
exit(2);
}
}
#if HELP
/* nicely display the actual polynomial equation we are solving */
printf("The %d approximate floating point solutions of:\n", highest_power);
for (i = highest_power; i >= 0; i--) {
if (a[i]) {
if (i && a[i] == 1.0) {
printf("+x");
} else {
printf("%+.*g", precision, a[i]);
if (i) {
printf("*x");
}
}
if (i > 1) {
printf("^%d", i);
}
printf(" ");
}
}
printf("= 0\nare:\n\n");
#endif
/* solve the polynomial equation */
w = gsl_poly_complex_workspace_alloc(highest_power + 1);
if (gsl_poly_complex_solve(a, highest_power + 1, w, z) != GSL_SUCCESS) {
fprintf(stderr, "%s: GSL approximation failed.\n", prog_name);
exit(1);
}
gsl_poly_complex_workspace_free(w);
/* display all solutions */
for (i = 0; i < highest_power; i++) {
#ifdef EPSILON /* zero out relatively very small values (which are floating point error) */
if (fabs(z[2*i] * EPSILON) > fabs(z[2*i+1]))
z[2*i+1] = 0.0;
else if (fabs(z[2*i+1] * EPSILON) > fabs(z[2*i]))
z[2*i] = 0.0;
#endif
#if HELP
printf("x = ");
#endif
display_root(i);
printf("\n");
}
exit(0);
}
void
display_root(int i)
{
printf("%+.*g", precision, z[2*i]); /* output real part */
if (z[2*i+1])
printf(" %+.*g*i", precision, z[2*i+1]); /* output imaginary part */
}
|
