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