1/* { dg-do run } */ 2/* { dg-options "-O2 -lm" } */ 3/* { dg-options "-O2 -msse2 -mfpmath=sse" { target { { i?86-*-* x86_64-*-* } && ia32 } } } */ 4/* { dg-require-effective-target sse2_runtime { target { { i?86-*-* x86_64-*-* } && ia32 } } } */ 5 6extern double fabs (double); 7extern void abort (void); 8 9const int MAX_ITERATIONS = 50; 10const double SMALL_ENOUGH = 1.0e-10; 11const double RELERROR = 1.0e-12; 12 13typedef struct p 14{ 15 int ord; 16 double coef[7]; 17} 18polynomial; 19 20static double 21polyeval (double x, int n, double *Coeffs) 22{ 23 register int i; 24 double val; 25 26 val = Coeffs[n]; 27 for (i = n - 1; i >= 0; i--) 28 val = val * x + Coeffs[i]; 29 30 return (val); 31} 32 33static int 34regula_falsa (int order, double *coef, double a, double b, double *val) 35{ 36 int its; 37 double fa, fb, x, fx, lfx; 38 39 fa = polyeval (a, order, coef); 40 fb = polyeval (b, order, coef); 41 42 if (fa * fb > 0.0) 43 return 0; 44 45 if (fabs (fa) < SMALL_ENOUGH) 46 { 47 *val = a; 48 return 1; 49 } 50 51 if (fabs (fb) < SMALL_ENOUGH) 52 { 53 *val = b; 54 return 1; 55 } 56 57 lfx = fa; 58 59 for (its = 0; its < MAX_ITERATIONS; its++) 60 { 61 x = (fb * a - fa * b) / (fb - fa); 62 fx = polyeval (x, order, coef); 63 if (fabs (x) > RELERROR) 64 { 65 if (fabs (fx / x) < RELERROR) 66 { 67 *val = x; 68 return 1; 69 } 70 } 71 else 72 { 73 if (fabs (fx) < RELERROR) 74 { 75 *val = x; 76 return 1; 77 } 78 } 79 80 if (fa < 0) 81 { 82 if (fx < 0) 83 { 84 a = x; 85 fa = fx; 86 if ((lfx * fx) > 0) 87 fb /= 2; 88 } 89 else 90 { 91 b = x; 92 fb = fx; 93 if ((lfx * fx) > 0) 94 fa /= 2; 95 } 96 } 97 else 98 { 99 if (fx < 0) 100 { 101 b = x; 102 fb = fx; 103 if ((lfx * fx) > 0) 104 fa /= 2; 105 } 106 else 107 { 108 a = x; 109 fa = fx; 110 if ((lfx * fx) > 0) 111 fb /= 2; 112 } 113 } 114 115 if (fabs (b - a) < RELERROR) 116 { 117 *val = x; 118 return 1; 119 } 120 121 lfx = fx; 122 } 123 124 return 0; 125} 126 127static int 128numchanges (int np, polynomial * sseq, double a) 129{ 130 int changes; 131 double f, lf; 132 polynomial *s; 133 changes = 0; 134 135 lf = polyeval (a, sseq[0].ord, sseq[0].coef); 136 137 for (s = sseq + 1; s <= sseq + np; s++) 138 { 139 f = polyeval (a, s->ord, s->coef); 140 if (lf == 0.0 || lf * f < 0) 141 changes++; 142 143 lf = f; 144 } 145 146 return changes; 147} 148 149int 150sbisect (int np, polynomial * sseq, double min_value, double max_value, 151 int atmin, int atmax, double *roots) 152{ 153 double mid; 154 int n1, n2, its, atmid; 155 156 if ((atmin - atmax) == 1) 157 { 158 if (regula_falsa (sseq->ord, sseq->coef, min_value, max_value, roots)) 159 return 1; 160 else 161 { 162 for (its = 0; its < MAX_ITERATIONS; its++) 163 { 164 mid = (min_value + max_value) / 2; 165 atmid = numchanges (np, sseq, mid); 166 if ((atmid < atmax) || (atmid > atmin)) 167 return 0; 168 169 if (fabs (mid) > RELERROR) 170 { 171 if (fabs ((max_value - min_value) / mid) < RELERROR) 172 { 173 roots[0] = mid; 174 return 1; 175 } 176 } 177 else 178 { 179 if (fabs (max_value - min_value) < RELERROR) 180 { 181 roots[0] = mid; 182 return 1; 183 } 184 } 185 186 if ((atmin - atmid) == 0) 187 min_value = mid; 188 else 189 max_value = mid; 190 } 191 192 roots[0] = mid; 193 return 1; 194 } 195 } 196 197 for (its = 0; its < MAX_ITERATIONS; its++) 198 { 199 mid = (min_value + max_value) / 2; 200 atmid = numchanges (np, sseq, mid); 201 if ((atmid < atmax) || (atmid > atmin)) 202 return 0; 203 204 if (fabs (mid) > RELERROR) 205 { 206 if (fabs ((max_value - min_value) / mid) < RELERROR) 207 { 208 roots[0] = mid; 209 return 1; 210 } 211 } 212 else 213 { 214 if (fabs (max_value - min_value) < RELERROR) 215 { 216 roots[0] = mid; 217 return 1; 218 } 219 } 220 221 n1 = atmin - atmid; 222 n2 = atmid - atmax; 223 224 if ((n1 != 0) && (n2 != 0)) 225 { 226 n1 = sbisect (np, sseq, min_value, mid, atmin, atmid, roots); 227 n2 = sbisect (np, sseq, mid, max_value, atmid, atmax, &roots[n1]); 228 229 return (n1 + n2); 230 } 231 232 if (n1 == 0) 233 min_value = mid; 234 else 235 max_value = mid; 236 } 237 238 roots[0] = mid; 239 return 1; 240} 241 242int 243main () 244{ 245 polynomial sseq[7] = { 246 {6, {0.15735259075109281, -5.1185263411378736, 1.8516070705868664, 247 7.348009172322695, -2.2152395279161343, -2.7543325329350692, 1.0}}, 248 {5, {-0.8530877235229789, 0.61720235686228875, 3.6740045861613475, 249 -1.4768263519440896, -2.2952771107792245, 1.0}}, 250 {4, {0.13072124257049417, 2.2220687798791126, -1.6299431586726509, 251 -1.6718404582408546, 1.0}}, 252 {3, {0.86776597575462633, -2.1051099695282511, -0.49008580100694688, 253 1.0}}, 254 {2, {-11.117984175064155, 10.89886635045883, 1.0}}, 255 {1, {0.94453099602191237, -1.0}}, 256 {0, {-0.068471716890574186}} 257 }; 258 259 double roots[7]; 260 int nroots; 261 262 nroots = sbisect (6, sseq, 0.0, 10000000.0, 5, 1, roots); 263 if (nroots != 4) 264 abort (); 265 266 return 0; 267} 268