1 /* MIT License
2 *
3 * Copyright (c) 2016--2017 Felix Lenders
4 *
5 * Permission is hereby granted, free of charge, to any person obtaining a copy
6 * of this software and associated documentation files (the "Software"), to deal
7 * in the Software without restriction, including without limitation the rights
8 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9 * copies of the Software, and to permit persons to whom the Software is
10 * furnished to do so, subject to the following conditions:
11 *
12 * The above copyright notice and this permission notice shall be included in all
13 * copies or substantial portions of the Software.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21 * SOFTWARE.
22 *
23 */
24
25 #include "trlib.h"
26 #include "trlib_private.h"
27
28 #include "_c99compat.h"
29
trlib_quadratic_zero(trlib_flt_t c_abs,trlib_flt_t c_lin,trlib_flt_t tol,trlib_int_t verbose,trlib_int_t unicode,char * prefix,FILE * fout,trlib_flt_t * t1,trlib_flt_t * t2)30 trlib_int_t trlib_quadratic_zero(trlib_flt_t c_abs, trlib_flt_t c_lin, trlib_flt_t tol,
31 trlib_int_t verbose, trlib_int_t unicode, char *prefix, FILE *fout,
32 trlib_flt_t *t1, trlib_flt_t *t2) {
33 trlib_int_t n = 0; // number of roots
34 trlib_flt_t q = 0.0;
35 trlib_flt_t dq = 0.0;
36 trlib_flt_t lin_sq = c_lin*c_lin;
37 *t1 = 0.0; // first root
38 *t2 = 0.0; // second root
39
40 if (fabs(c_abs) > tol*lin_sq) {
41 // well-behaved non-degenerate quadratic
42 // compute discriminant
43 q = lin_sq - 4.0 * c_abs;
44 if ( fabs(q) <= (TRLIB_EPS*c_lin)*(TRLIB_EPS*c_lin) ) {
45 // two distinct zeros, but discrimant tiny --> numeric double zero
46 // initialize on same root obtained by standard formula with zero discrement, let newton refinement do the rest
47 n = 2;
48 *t1 = -.5*c_lin; *t2 = *t1;
49 }
50 else if ( q < 0.0 ) {
51 n = 2;
52 *t1 = 0.0; *t2 = 0.0;
53 return n;
54 }
55 else {
56 // discriminant large enough, two distinc zeros
57 n = 2;
58 // start with root according to plus sign to avoid cancellation
59 *t1 = -.5 * ( c_lin + copysign( sqrt(q), c_lin ) );
60 *t2 = c_abs/(*t1);
61 if (*t2 < *t1) { q = *t2; *t2 = *t1; *t1 = q; }
62 }
63 }
64 else {
65 n = 2;
66 if (c_lin < 0.0) { *t1 = 0.0; *t2 = - c_lin; }
67 else { *t1 = - c_lin; *t2 = 0.0; }
68 }
69
70 // newton correction
71 q = (*t1+c_lin)*(*t1)+c_abs; dq = 2.0*(*t1)+c_lin;
72 if (dq != 0.0) { *t1 = *t1 - q/dq; }
73 q = (*t2+c_lin)*(*t2)+c_abs; dq = 2.0*(*t2)+c_lin;
74 if (dq != 0.0) { *t2 = *t2 - q/dq; }
75 return n;
76 }
77