/* { dg-do run } */ /* { dg-options "-O2 -lm" } */ /* { dg-options "-O2 -msse2 -mfpmath=sse" { target { { i?86-*-* x86_64-*-* } && ilp32 } } } */ #ifdef __i386__ #include "cpuid.h" #endif extern double fabs (double); extern void abort (void); const int MAX_ITERATIONS = 50; const double SMALL_ENOUGH = 1.0e-10; const double RELERROR = 1.0e-12; typedef struct p { int ord; double coef[7]; } polynomial; static double polyeval (double x, int n, double *Coeffs) { register int i; double val; val = Coeffs[n]; for (i = n - 1; i >= 0; i--) val = val * x + Coeffs[i]; return (val); } static int regula_falsa (int order, double *coef, double a, double b, double *val) { int its; double fa, fb, x, fx, lfx; fa = polyeval (a, order, coef); fb = polyeval (b, order, coef); if (fa * fb > 0.0) return 0; if (fabs (fa) < SMALL_ENOUGH) { *val = a; return 1; } if (fabs (fb) < SMALL_ENOUGH) { *val = b; return 1; } lfx = fa; for (its = 0; its < MAX_ITERATIONS; its++) { x = (fb * a - fa * b) / (fb - fa); fx = polyeval (x, order, coef); if (fabs (x) > RELERROR) { if (fabs (fx / x) < RELERROR) { *val = x; return 1; } } else { if (fabs (fx) < RELERROR) { *val = x; return 1; } } if (fa < 0) { if (fx < 0) { a = x; fa = fx; if ((lfx * fx) > 0) fb /= 2; } else { b = x; fb = fx; if ((lfx * fx) > 0) fa /= 2; } } else { if (fx < 0) { b = x; fb = fx; if ((lfx * fx) > 0) fa /= 2; } else { a = x; fa = fx; if ((lfx * fx) > 0) fb /= 2; } } if (fabs (b - a) < RELERROR) { *val = x; return 1; } lfx = fx; } return 0; } static int numchanges (int np, polynomial * sseq, double a) { int changes; double f, lf; polynomial *s; changes = 0; lf = polyeval (a, sseq[0].ord, sseq[0].coef); for (s = sseq + 1; s <= sseq + np; s++) { f = polyeval (a, s->ord, s->coef); if (lf == 0.0 || lf * f < 0) changes++; lf = f; } return changes; } int sbisect (int np, polynomial * sseq, double min_value, double max_value, int atmin, int atmax, double *roots) { double mid; int n1, n2, its, atmid; if ((atmin - atmax) == 1) { if (regula_falsa (sseq->ord, sseq->coef, min_value, max_value, roots)) return 1; else { for (its = 0; its < MAX_ITERATIONS; its++) { mid = (min_value + max_value) / 2; atmid = numchanges (np, sseq, mid); if ((atmid < atmax) || (atmid > atmin)) return 0; if (fabs (mid) > RELERROR) { if (fabs ((max_value - min_value) / mid) < RELERROR) { roots[0] = mid; return 1; } } else { if (fabs (max_value - min_value) < RELERROR) { roots[0] = mid; return 1; } } if ((atmin - atmid) == 0) min_value = mid; else max_value = mid; } roots[0] = mid; return 1; } } for (its = 0; its < MAX_ITERATIONS; its++) { mid = (min_value + max_value) / 2; atmid = numchanges (np, sseq, mid); if ((atmid < atmax) || (atmid > atmin)) return 0; if (fabs (mid) > RELERROR) { if (fabs ((max_value - min_value) / mid) < RELERROR) { roots[0] = mid; return 1; } } else { if (fabs (max_value - min_value) < RELERROR) { roots[0] = mid; return 1; } } n1 = atmin - atmid; n2 = atmid - atmax; if ((n1 != 0) && (n2 != 0)) { n1 = sbisect (np, sseq, min_value, mid, atmin, atmid, roots); n2 = sbisect (np, sseq, mid, max_value, atmid, atmax, &roots[n1]); return (n1 + n2); } if (n1 == 0) min_value = mid; else max_value = mid; } roots[0] = mid; return 1; } int main () { polynomial sseq[7] = { {6, {0.15735259075109281, -5.1185263411378736, 1.8516070705868664, 7.348009172322695, -2.2152395279161343, -2.7543325329350692, 1.0}}, {5, {-0.8530877235229789, 0.61720235686228875, 3.6740045861613475, -1.4768263519440896, -2.2952771107792245, 1.0}}, {4, {0.13072124257049417, 2.2220687798791126, -1.6299431586726509, -1.6718404582408546, 1.0}}, {3, {0.86776597575462633, -2.1051099695282511, -0.49008580100694688, 1.0}}, {2, {-11.117984175064155, 10.89886635045883, 1.0}}, {1, {0.94453099602191237, -1.0}}, {0, {-0.068471716890574186}} }; double roots[7]; int nroots; #ifdef __i386__ unsigned int eax, ebx, ecx, edx; if (!__get_cpuid (1, &eax, &ebx, &ecx, &edx)) return 0; if (!(edx & bit_SSE2)) return 0; #endif nroots = sbisect (6, sseq, 0.0, 10000000.0, 5, 1, roots); if (nroots != 4) abort (); return 0; }