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