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
6 extern double fabs (double);
7 extern void abort (void);
8
9 const int MAX_ITERATIONS = 50;
10 const double SMALL_ENOUGH = 1.0e-10;
11 const double RELERROR = 1.0e-12;
12
13 typedef struct p
14 {
15 int ord;
16 double coef[7];
17 }
18 polynomial;
19
20 static double
polyeval(double x,int n,double * Coeffs)21 polyeval (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
33 static int
regula_falsa(int order,double * coef,double a,double b,double * val)34 regula_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
127 static int
numchanges(int np,polynomial * sseq,double a)128 numchanges (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
149 int
sbisect(int np,polynomial * sseq,double min_value,double max_value,int atmin,int atmax,double * roots)150 sbisect (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
242 int
main()243 main ()
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