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