Subversion Repositories openrisc

/* { dg-do run } */

/* { dg-options "-O2 -lm" } */

/* { dg-options "-O2 -msse2 -mfpmath=sse" { target { { i?86-*-* x86_64-*-* } && ilp32 } } } */

/* { dg-require-effective-target sse2 { target { { i?86-*-* x86_64-*-* } && ilp32 } } } */

/* { dg-require-effective-target sse2_runtime { target { { i?86-*-* x86_64-*-* } && ilp32 } } } */

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;

  nroots = sbisect (6, sseq, 0.0, 10000000.0, 5, 1, roots);

  if (nroots != 4)

    abort ();

  return 0;

}