gcc/testsuite/gcc.dg/pr36584.c

   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
  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
  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
 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
 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
 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 }