test.cpp

牛顿法解多项式的根 包含重根的判断和多项式的输入

#include <stdio.h>

#include <malloc.h>





double Polynomial_Root(double c[], int n, double a, double b, double eps);

 

int main()
{

    double eps = 0.00005;

    double *c;

    int n;

    double a, b;

    int i;

 

    while (scanf("%d", &n)!= EOF){

       c = (double *)calloc((n+1), sizeof(double));

       for (i=n; i>=0; i--) 

           scanf("%lf", &c[i]);

       scanf("%lf %lf", &a, &b);

       printf("%.4lf\n", Polynomial_Root(c, n, a, b, eps));

       free(c);
    }

 
    return 0;

}



/*
double ff(double c[], int n, double x)
{
       double sum = 0;
	   int i;
			
	   for ( i = n; i >= 0; i-- )
			sum = sum * x + c[i];

	   return sum;
}

double dff(double c[], int n, double x)
{
     double sum = 0;
	 int i;

     for ( i = n; i >= 0; i-- )
		 sum = sum * x + c[i] * i;
	     
	 return sum;
}
*/

#include <math.h>

void HornerMethod(int n, double a[], double x0, double *f, double *df)
{
	double y, z;
	int j;

	y = a[n];
	z = a[n];

    for(j = n-1; j >= 1; j--)
	{
		y = x0 * y + a[j];	
		z = x0 * z + y;
	}

	y = x0 *y + a[0];
	
	*f = y;
	*df = z;
}

/*5 1 -5 10 -10 5 -1
-100 20
二分法结果时是1.0001，牛顿法就是1.0000*/


int NewtonMethod(double c[], int n, double eps, double p0, double *p, double a, double b)
{
   double f, df;
   int i;
   int MAX_N = 1000;


   for(i = 0; i < MAX_N; i++)
   {
	   //f = ff(c, n, p0);
	   //df = dff(c, n, p0);
	   HornerMethod(n, c, p0, &f, &df);
	   
	   if((df==0)&&(f==0))
	   {
           if(p0 < a || p0 > b)
      		   return -1;

	       if(fabs(*p) <= (eps*1000))
		   {
			   *p = +0.0;
			   return 0;
		   }
			
		   *p = p0;
		   return 0;
	   }
	   else
		   if((df==0)&&(f!=0))
		      return -1;

       *p = p0 - (f) / (df);
	   
	   if(fabs(*p - p0) <= eps)
	   {
		   if(p0 < a || p0 > b)
      		   return -1;

	       if(fabs(*p) <= (eps*1000))
		   {
			   *p = +0.0;
			   return 0;
		   }
			
		   *p = p0;
		   return 0;
	   }

	   p0 = *p;
   }

   return -1;
}

double Polynomial_Root(double c[], int n, double a, double b, double eps)
{
     double p0, p, h;
	 int MAX_N = 20;
	 int i = 0, m, mark;
	 
     if(a > b)
	 {
	    p = a;
		a = b;
		b = p;
	 }

     m = 20;
	 h = (b - a) / (double)m;
	 for(i = 0; i <= m; i++)
	 {
		p0 = a + h * (double)i;
		mark = NewtonMethod(c, n, eps/1000, p0, &p, a, b);

		if(mark == 0)
		{
			  return p;
		}
		else 
			  continue;
	 }

	 printf("Newton's method failed!!\n");
	 return 0; 
}