/***************************實驗題目:計算非線性方程f(x)=0的求根問題***********************


一:設計內容和期望
     本程序把函數f(x)在X(k) 完成二次泰勒展開,并得到一個新的迭代公式計算X(k+1),以此來計算非線性方程
     f(x)=0的根,從而改進牛頓法的計算公式,使計算精度從平方收斂達到3階收斂(估計,尚待證明)


二:遇到的問題
    1:對于任意f(x)的輸入;
	2:f(x)導數的實現;
	3:誤差精度的控制;


三:根據以上問題,提出下列解決方案
    1:采用數據結構的相關知識解決輸入問題,(由于受到算法限制先完成對f(x)為多項式函數的形式,其他
	    例如包含sin,cos,exp的形式有待解決);
    2:編寫通用多項式求導和求值問題,以便完成新的迭代公式;
	3:控制x(k+1)與x(k)的誤差(精度可手動輸入控制),保證迭代數列收斂,以便達到精度為3階收斂;
	

四:實驗效果預期
    1:可手動輸入多項式形式的非線性方程;
	2:可手動輸入精度,迭代初值;
	3:通過與牛頓法進行對比分析,從而從實驗中證明精度有所提高;
	  (例如同樣的多項式非線性方程,在同樣的精度下,新方法比牛頓法迭代次數少即可證明精度提高)*/






#include "stdio.h"
#include "iostream.h"
#include "stdlib.h"
#include "math.h"
#include "iomanip.h"
/**********************************
*  	常量定義
*********************************/


#define TRUE 1
#define FALSE 0
#define OK 1
#define ERROR 0


/**********************************
*  	類型定義
*********************************/

typedef int Status;  /* 函數返回類型   */

/*  鏈表結點的數據元素的數據類型定義 */

typedef struct{
  float coef;
  int expn;
} Term,LkLtElemtype;

/*  鏈表結點的數據類型定義 */

typedef struct LNode{
  LkLtElemtype data;         /* 說明結點的數據元素類型為上面定義的結構體類型 */
  struct LNode *next;
} LNode,*Link, *Position;


/* 鏈表的數據類型定義 */

typedef struct{
  Link head,tail;   /*  head指向鏈表的頭結點，tail指向最后一個結點 */
  int len;         /* 鏈表的長度*/
}LinkList;

/* 定義一個多項式類型為鏈表類型  */

typedef LinkList Polynomial;


/********************************** 
*   對鏈表的操作實現
*********************************/


/* 分配由*p指向的值為e的結點，并返回OK, 若分配失敗返回ERROR  */

Status MakeNode(Link *p,LkLtElemtype e){

  *p=(Link) malloc(sizeof(LNode));
  if(!(*p)) return ERROR;
  (*p)->data=e;
  (*p)->next=NULL;
  return OK;
}


/* 釋放*p所指向的結點  */

void FreeNode(Link *p){
  if(!(*p)){
    free(*p);
    *p=NULL;
  }
}


/*  構造一個空的線性鏈表,成功返回OK,失敗返回ERROR  */

Status InitList(LinkList * l){
  if(!l) return ERROR;
  l->head=(Link)malloc(sizeof(LNode));
  if(!(l->head)) return ERROR;
  l->head->next=NULL;
  l->tail=NULL;
  l->len=0;
  return OK;
}


/* 銷毀線性表鏈表*/

Status DestroyList(LinkList * l){
  Link q=NULL;
  Link p=NULL;
  if(!l) return ERROR;
  p=l->head->next;
  while(!p){
    q=p->next;
    free(p);
    p=q;
  }
  free(l->head);
  return OK;
}

/* 已知h指向線性鏈表的頭結點，將s所指向結點插入在第一個結點之前 */

Status InsFirst(Link head,Link s){
  if ((!head ) || (!s)) return ERROR;
  s->next=head->next;
  head->next=s;
  return OK;
}

/*  已知*p指向線性鏈表中的一個結點，用e更新*p所指結點中數據元素的值  */

Status SetCurElem(Link *p,LkLtElemtype e){
  if (!(*p)) return ERROR;
  (*p)->data=e;
  return OK;
}


/*  已知p指向線性鏈表中的一個結點，返回p所指結點中數據元素的值  */

LkLtElemtype GetCurElem(Link p){
  return p->data;
}


/* 返回線性鏈表中頭結點的位置 */

Position GetHead(LinkList l){
   return l.head;
}

/*************************************************************
*   對多項式鏈表所特有的操作實現(見教材P42-43)
*************************************************************/

/*  依a的指數值小于、等于、大于b的指數值，分別返回-1,0,1 */


int cmp(Term a,Term b){
  if(a.expn<b.expn) return -1;
  if(a.expn==b.expn) return 0;
  if(a.expn>b.expn) return 1;
  return OK;
}


/* 若有序鏈表中存在與e滿足判定函數compare()取值為0的元素，則*q批示第一個值為e的結點的位置，并返回TRUE；
*  否則*q批示第一個與e滿足判定函數compare()取值>0的元素的前驅的位置，并返回FALSE
*/
Status LocateElem2(LinkList l,LkLtElemtype e,Position *q,int (*compare)(LkLtElemtype,LkLtElemtype)){
  Link p;
  Link pp;
  p=l.head;
  do{pp=p,p=p->next;}
   while(p&&(*compare)(p->data,e)<0);

  if  (!p||(*compare)(p->data,e)>0)
 {*q=pp;return FALSE;}
else  {*q=p;return TRUE;}
}

/* 輸入m項的系數和指數，建立表示一元多項式的按指數有序的鏈表 */

void CreatePolyn(Polynomial *p,int m)
{
  Term *t;
  int i;
  Link q;
  Link s;
  Link h;
  float f;
  InitList(p);
 /* printf("Ployn length:%d\n",(*p).len);*/
  h=GetHead(*p);
 /* printf("head:%o",h);*/
  t=(Term *)malloc(sizeof(Term));
  t->coef=0.0f;
  t->expn=-1;
  SetCurElem(&h,*t);
/* printf("head:expn:%d\n",(h->data).expn);*/
  printf("請輸入一元多項式f(x):(格式:系數 指數,順序可任意)\n");
  for(i=1;i<=m;++i)
  {
    //printf("please input the %d coef:",i);
    scanf("%f",&f);
    t->coef=f;
    //printf("please input the %d expn:",i);
    scanf("%d",&(t->expn));
    if(!LocateElem2(*p,*t,&q,(*cmp))){
      if(MakeNode(&s,*t)) InsFirst(q,s);
    }
  }
}



/* 顯示一個用鏈表表示的多項式 */


void PrintPolyn(Polynomial p){
  Link q;
  int i=1;
  float coef;
  int expn;
  q=p.head->next;
  while(q)
  {
      coef=(q->data).coef;
      expn=(q->data).expn;//cout<<coef<<"x^"<<expn<<endl;
      if(i%10==0) printf("\n");
      else
	  {
		  if(expn!=0)
		{
          if(coef>0)
		  {
            if(i==1)  printf("%3.2fx^%d",coef,expn);
            else      printf("+%3.2fx^%d",coef,expn);
          }
          else printf("%3.2fx^%d",coef,expn);
       }
       else
	   {
         if(coef>0)
		 {
            if(i==1)  printf("%3.2f",coef);
            else      printf("+%3.2f",coef);
          }
          else printf("%3.2f",coef);
		 }
	}
      
      q=q->next;
      i++;
  }
  printf("\n");
 
}


/*多項式求導函數*/

void qiudao(Polynomial *p,Polynomial pa,int m)
{
    Link q1;
	q1=pa.head->next;
	Term *t;
  //int i;
  Link q;
  Link s;
  Link h;
  
  InitList(p);
   h=GetHead(*p);
 /* printf("head:%o",h);*/
  t=(Term *)malloc(sizeof(Term));
  t->coef=0.0f;
  t->expn=-1;
  SetCurElem(&h,*t);
  while (q1)
  {
    t->coef=q1->data.coef*q1->data.expn;
    t->expn=q1->data.expn-1;
	if(t->coef!=0)
	{
       if(!LocateElem2(*p,*t,&q,(*cmp)))
	   {
         if(MakeNode(&s,*t)) InsFirst(q,s);
	   }
	}
	
	
	q1=q1->next;
  }
}

/*多項式求值*/

double qiuzhi(Polynomial p,double x)
{
	Link q;double sum=0.0;
    q=p.head->next;
	while(q)
	{
		sum=sum+q->data.coef*pow(x,q->data.expn);
		q=q->next;
	}
	return sum;
}

/*初值條件的判斷*/

int czpanduan(double *x,Polynomial pa,Polynomial pb,Polynomial pc)
{
	double f,f1,f2,sq;
        f=qiuzhi(pa,*x);f1=qiuzhi(pb,*x);f2=qiuzhi(pc,*x);
        sq=f1*f1-2*f2*f;	//cout<<"sq="<<sq<<"  "<<"f1="<<f1<<"  "<<"f2="<<f2<<endl;
	if(sq>=0&&f2!=0&&f1!=0) return OK;
	else return ERROR;
	
}

/*利用新的迭代公式計算收斂數列x(k+1)*/

double xindegongshi(double x,Polynomial pa,Polynomial pb,Polynomial pc)
{
     double f,f1,f2;
	 double t,t1;    //cout<<"x="<<x;
	f=qiuzhi(pa,x);
    //printf("f(%2f)=%f\n",x,f);
      f1=qiuzhi(pb,x);
    //printf("f`(%2f)=%f\n",x,f1);
     f2=qiuzhi(pc,x);
    //printf("f``(%2f)=%f\n",x,f2);
    t=x-(f1+sqrt(f1*f1-2*f2*f))/f2;
	t1=x-(f1-sqrt(f1*f1-2*f2*f))/f2;//新的迭代公式

	//cout<<"t="<<t;cout<<"t="<<t1;
    
	if(fabs(t-x)>fabs(t1-x)) return t1;
	else return t;//通過與上一次的值比較,選擇絕對值差小的作為下一次的迭代值,從而保證數列收斂
}




int xinfangfa(double x,double exp,Polynomial pa,Polynomial pb,Polynomial pc)
{ 
	
	double x1;int i=0;
	x1=xindegongshi(x,pa,pb,pc);
     if(!czpanduan(&x1,pa,pb,pc)) return ERROR;
	        cout<<"x["<<i++<<"]="<<setprecision(8)<<x<<endl;
			cout<<"x["<<i++<<"]="<<setprecision(8)<<x1<<endl;
     while (fabs(x1-x)>=exp)
	 {
		 x=x1;x1=xindegongshi(x,pa,pb,pc);
		 cout<<"x["<<i++<<"]="<<setprecision(8)<<x1<<endl;
	 }
	 return OK;
}

/*newton迭代公式*/

double gongshi(double x,Polynomial pa,Polynomial pb,Polynomial pc)
{
     double f,f1;
	 double t;   
	f=qiuzhi(pa,x);
    f1=qiuzhi(pb,x);
    t=x-f/f1;
    return t;
}


void newton(double x,double exp,Polynomial pa,Polynomial pb,Polynomial pc)
{ 
	
	double x1;int i=0;
	x1=gongshi(x,pa,pb,pc);
	        cout<<"x["<<i++<<"]="<<setprecision(8)<<x<<endl;
			cout<<"x["<<i++<<"]="<<setprecision(8)<<x1<<endl;
     while (fabs(x1-x)>=exp)
	 {
		 x=x1;x1=gongshi(x,pa,pb,pc);
		 cout<<"x["<<i++<<"]="<<setprecision(8)<<x1<<endl;
	 }
}





/*************************************************************
*               主函數
*************************************************************/

void main()
{

  Polynomial pa,pb,pc;/* 定義多項式鏈表pa,pa,pc分別保留原函數,一階導函數,二階導函數 */

  int ma; /*pa的項數*/

  double x;/*x作為初值*/
  
  double exp;/*定義精度*/
  
  /*輸入多項式函數項數*/

  printf("請輸入多項式函數項數:");
  scanf("%d",&ma);

  /*按指數有序構造多項式鏈表pa,pb,pc*/

  CreatePolyn(&pa,ma);qiudao(&pb,pa,ma);qiudao(&pc,pb,ma);

  /*顯示多項式pa*/

  printf("你輸入的多項式函數f(x)=");
  PrintPolyn(pa);
  
  /*顯示多項式一階和二階導函數*/

  printf("f`(x)=");
  PrintPolyn(pb);
  printf("f``(x)=");PrintPolyn(pc);cout<<endl;

  /*輸入迭代初值*/
 
  printf("現在計算f(x)=0的根,請輸入一個迭代初值:");
  cin>>x;

  /*判斷迭代初值是否滿足要求*/
  while(1)
  {
   if(czpanduan(&x,pa,pb,pc)) break;
   else 
   {
	  printf("初值不和要求,請重新輸入一個迭代初值:");
      cin>>x;
   }
  }
			  
  /*輸入精度*/

  printf("請輸入精度(例如1e-3,表示0.001)");
  cin>>exp;cout<<endl;

  /*利用新的迭代公式計算*/

  cout<<"利用新的迭代公式計算的收斂數列如下:"<<endl;
  while(1)
  {
	  if(xinfangfa(x,exp,pa,pb,pc)) break;
      printf("中間變量不和要求,請重新輸入一個迭代初值,重新進行計算:");
      cin>>x;
  }
  cout<<endl;
     
 /*利用牛頓法的迭代公式計算*/

  cout<<"利用牛頓法的迭代公式計算的收斂數列如下:"<<endl;
  newton(x,exp,pa,pb,pc);cout<<endl;

  /*計算完畢后銷毀多項式函數，釋放空間*/
  DestroyList(&pa);DestroyList(&pb);DestroyList(&pc);
  
  
  cout<<"通過數據請分析該算法的優缺點,以便加以改進,謝謝!!"<<endl;
}