1 /* rquad.f -- translated by f2c (version 19980913).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
8 /* Table of constant values */
9
10 static integer c__2 = 2;
11
12 /* $Procedure RQUAD ( Roots of a quadratic equation ) */
rquad_(doublereal * a,doublereal * b,doublereal * c__,doublereal * root1,doublereal * root2)13 /* Subroutine */ int rquad_(doublereal *a, doublereal *b, doublereal *c__,
14 doublereal *root1, doublereal *root2)
15 {
16 /* System generated locals */
17 doublereal d__1, d__2;
18
19 /* Builtin functions */
20 double sqrt(doublereal);
21
22 /* Local variables */
23 doublereal scale;
24 extern /* Subroutine */ int chkin_(char *, ftnlen), moved_(doublereal *,
25 integer *, doublereal *);
26 doublereal discrm;
27 logical zeroed;
28 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
29 ftnlen), setmsg_(char *, ftnlen);
30 extern logical return_(void);
31 doublereal con, lin, sqr;
32
33 /* $ Abstract */
34
35 /* Find the roots of a quadratic equation. */
36
37 /* $ Disclaimer */
38
39 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
40 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
41 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
42 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
43 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
44 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
45 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
46 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
47 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
48 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
49
50 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
51 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
52 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
53 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
54 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
55 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
56
57 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
58 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
59 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
60 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
61
62 /* $ Required_Reading */
63
64 /* None. */
65
66 /* $ Keywords */
67
68 /* MATH */
69 /* POLYNOMIAL */
70 /* ROOT */
71
72 /* $ Declarations */
73 /* $ Brief_I/O */
74
75 /* Variable I/O Description */
76 /* -------- --- -------------------------------------------------- */
77
78 /* A I Coefficient of quadratic term. */
79 /* B I Coefficient of linear term. */
80 /* C I Constant. */
81 /* ROOT1 O Root built from positive discriminant term. */
82 /* ROOT2 O Root built from negative discriminant term. */
83
84 /* $ Detailed_Input */
85
86 /* A, */
87 /* B, */
88 /* C are the coefficients of a quadratic polynomial */
89
90 /* 2 */
91 /* Ax + Bx + C. */
92
93 /* $ Detailed_Output */
94
95 /* ROOT1, */
96 /* ROOT2 are the roots of the equation, */
97
98 /* 2 */
99 /* Ax + Bx + C = 0. */
100
101
102 /* ROOT1 and ROOT2 are both arrays of length 2. The */
103 /* first element of each array is the real part of a */
104 /* root; the second element contains the complex part */
105 /* of the same root. */
106
107 /* When A is non-zero, ROOT1 represents the root */
108
109 /* _____________ */
110 /* / 2 */
111 /* - B + \/ B - 4AC */
112 /* --------------------------- */
113 /* 2A */
114
115
116 /* and ROOT2 represents the root */
117
118 /* _____________ */
119 /* / 2 */
120 /* - B - \/ B - 4AC */
121 /* --------------------------- . */
122 /* 2A */
123
124
125 /* When A is zero and B is non-zero, ROOT1 and ROOT2 */
126 /* both represent the root */
127
128 /* - C / B. */
129
130 /* $ Parameters */
131
132 /* None. */
133
134 /* $ Exceptions */
135
136 /* 1) If the input coefficients A and B are both zero, the error */
137 /* SPICE(DEGENERATECASE) is signalled. The output arguments */
138 /* are not modified. */
139
140 /* $ Files */
141
142 /* None. */
143
144 /* $ Particulars */
145
146 /* None. */
147
148 /* $ Examples */
149
150 /* 1) Humor us and suppose we want to compute the "golden ratio." */
151
152 /* The quantity r is defined by the equation */
153
154 /* 1/r = r/(1-r), */
155
156 /* which is equivalent to */
157
158 /* 2 */
159 /* r + r - 1 = 0. */
160
161 /* The following code frament does the job. */
162
163
164 /* C */
165 /* C Compute "golden ratio." The root we want, */
166 /* C */
167 /* C ___ */
168 /* C / */
169 /* C -1 + \/ 5 */
170 /* C -----------, */
171 /* C 2 */
172 /* C */
173 /* C */
174 /* C is contained in ROOT1. */
175 /* C */
176
177 /* CALL RQUAD ( 1.D0, 1.D0, -1.D0, ROOT1, ROOT2 ) */
178
179 /* PRINT *, 'The "golden ratio" is ', ROOT1(1) */
180
181
182 /* 2) The equation, */
183
184 /* 2 */
185 /* x + 1 = 0 */
186
187 /* can be solved by the code fragment */
188
189
190 /* C */
191 /* C Let's do one with imaginary roots just for fun. */
192 /* C */
193
194 /* CALL RQUAD ( 1.D0, 0.D0, 1.D0, ROOT1, ROOT2 ) */
195
196 /* PRINT *, 'ROOT1 is ', ROOT1 */
197 /* PRINT *, 'ROOT2 is ', ROOT2 */
198
199 /* The printed results will be something like: */
200
201
202 /* ROOT1 is 0.000000000000000 1.000000000000000 */
203 /* ROOT2 is 0.000000000000000 -1.000000000000000 */
204
205 /* $ Restrictions */
206
207 /* No checks for overflow of the roots are performed. */
208
209 /* $ Literature_References */
210
211 /* None. */
212
213 /* $ Author_and_Institution */
214
215 /* N.J. Bachman (JPL) */
216
217 /* $ Version */
218
219 /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */
220
221 /* Comment section for permuted index source lines was added */
222 /* following the header. */
223
224 /* - SPICELIB Version 1.0.0, 10-JUL-1990 (NJB) */
225
226 /* -& */
227 /* $ Index_Entries */
228
229 /* roots of a quadratic equation */
230
231 /* -& */
232
233 /* SPICELIB functions */
234
235
236 /* Local variables */
237
238
239 /* Standard SPICE error handling. */
240
241 if (return_()) {
242 return 0;
243 } else {
244 chkin_("RQUAD", (ftnlen)5);
245 }
246
247 /* The degree of the equation is zero unless at least one of the */
248 /* second or first degree coefficients is non-zero. */
249
250 if (*a == 0. && *b == 0.) {
251 setmsg_("Both 1st and 2nd degree coefficients are zero.", (ftnlen)46);
252 sigerr_("SPICE(DEGENERATECASE)", (ftnlen)21);
253 chkout_("RQUAD", (ftnlen)5);
254 return 0;
255 }
256
257 /* If we can scale the coefficients without zeroing any of them out, */
258 /* we will do so, to help prevent overflow. */
259
260 /* Computing MAX */
261 d__1 = abs(*a), d__2 = abs(*b), d__1 = max(d__1,d__2), d__2 = abs(*c__);
262 scale = max(d__1,d__2);
263 zeroed = *a != 0. && *a / scale == 0. || *b != 0. && *b / scale == 0. || *
264 c__ != 0. && *c__ / scale == 0.;
265 if (! zeroed) {
266 sqr = *a / scale;
267 lin = *b / scale;
268 con = *c__ / scale;
269 } else {
270 sqr = *a;
271 lin = *b;
272 con = *c__;
273 }
274
275 /* If the second-degree coefficient is non-zero, we have a bona fide */
276 /* quadratic equation, as opposed to a linear equation. */
277
278 if (sqr != 0.) {
279
280 /* Compute the discriminant. */
281
282 /* Computing 2nd power */
283 d__1 = lin;
284 discrm = d__1 * d__1 - sqr * 4. * con;
285
286 /* A non-negative discriminant indicates that the roots are */
287 /* real. */
288
289 if (discrm >= 0.) {
290
291 /* The imaginary parts of both roots are zero. */
292
293 root1[1] = 0.;
294 root2[1] = 0.;
295
296 /* We can take advantage of the fact that CON/SQR is the */
297 /* product of the roots to improve the accuracy of the root */
298 /* having the smaller magnitude. We compute the larger root */
299 /* first and then divide CON/SQR by it to obtain the smaller */
300 /* root. */
301
302 if (lin < 0.) {
303
304 /* ROOT1 will contain the root of larger magnitude. */
305
306 root1[0] = (-lin + sqrt(discrm)) / (sqr * 2.);
307 root2[0] = con / sqr / root1[0];
308 } else if (lin > 0.) {
309
310 /* ROOT2 will contain the root of larger magnitude. */
311
312 root2[0] = (-lin - sqrt(discrm)) / (sqr * 2.);
313 root1[0] = con / sqr / root2[0];
314 } else {
315
316 /* The roots have the same magnitude. */
317
318 root1[0] = sqrt(discrm) / (sqr * 2.);
319 root2[0] = -root1[0];
320 }
321
322 /* The only other possibility is that the roots are complex. */
323
324 } else {
325
326 /* The roots are complex conjugates, so they have equal */
327 /* magnitudes. */
328
329 root1[0] = -lin / (sqr * 2.);
330 root1[1] = sqrt(-discrm) / (sqr * 2.);
331 root2[0] = root1[0];
332 root2[1] = -root1[1];
333 }
334
335 /* If the second-degree coefficient is zero, we actually have a */
336 /* linear equation. */
337
338 } else if (lin != 0.) {
339 root1[0] = -con / lin;
340 root1[1] = 0.;
341
342 /* We set the second root equal to the first, rather than */
343 /* leaving it undefined. */
344
345 moved_(root1, &c__2, root2);
346 }
347 chkout_("RQUAD", (ftnlen)5);
348 return 0;
349 } /* rquad_ */
350
351