#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <dos.h>
#include <alloc.h>
#include <float.h>
#include <conio.h>
#define rnd ((rand()/32767.0)*2-1)
#define EPS 0.001
//int kk;
//typedef struct polinom STR;


typedef struct polinom{ float a[40];
			  int n[40];
			  int k;
		      }pln;

void inpol(pln *el,char[]),outpol(pln *el,char[]);
void copypol(pln *el,pln *el1),normpol(pln *el),delpol(pln *el,pln *el1);


typedef float * FU ;
void Mat2Pol(pln *el,FU a,int n);
void Proizvodnaia95(pln *el,pln *el1);
void PolumnN(pln *el,float n);
float Verdict_Newtona(pln *el,int n);
float point(pln *el,float num);
void PolSetMin(pln *el);
int System_Shturma(pln *el);
void SimplCreate(pln *el,float n);

int zero(float a, float b, float eps,pln *el, float *x, float *fx);
void Korni(pln *el);
void Master(pln *el);

void prrivet(FU,FU,int,int);
void pechat(FU, int, int, char *);
void handinit(FU, int);
void umn(FU, FU, FU, int, int, int);
void umnn(FU, float, int, int);
void sub(FU, FU, FU, int);
void inition(FU, int, int);
void prisv(FU, FU, int, int);
void InitTest_4x4(FU a,int m,int n);
void Init_Ed(FU a,int m,int n);
float Sled(FU a,int m,int n);
FU memory(int, int);
  pln *pol[2],*pold,*pst[2],*work[16];
  pln *polm,*ish,*ism;
  float hi,lo;
  int NK;

void main(void)
{
  int t,n,w,i,j,g,v1,v2,dv;char q; long *v;
  float prom1,prom2;
  FU Ann[10],Bnn[10],Enn[10],Pnn[10];

  polm=(pln *)farmalloc(sizeof(pln));
  ish=(pln *)farmalloc(sizeof(pln));
  ism=(pln *)farmalloc(sizeof(pln));
  pol[0]=(pln *)farmalloc(sizeof(pln));
  pol[1]=(pln *)farmalloc(sizeof(pln));
  pst[0]=(pln *)farmalloc(sizeof(pln));
  pst[1]=(pln *)farmalloc(sizeof(pln));
  pold=(pln *)farmalloc(sizeof(pln));
//  pold[1]=(pln *)farmalloc(sizeof(pln));
  for(i=0;i<16;i++)
  {
   work[i]=(pln *)farmalloc(sizeof(pln));
  }


do {
   clrscr();
   n=4;
	 Ann[9]=memory(n,n);
   InitTest_4x4(Ann[9],n,n);
	 pechat(Ann[9],n,n,"Исходная :");
   printf("\n Sled Ann: %f",Sled(Ann[9],n,n));
   printf("\n\n\n\n\n\n\n \t\t Press Any Key............\n");
   getch();

	 Pnn[0]=memory(1,n);
	 Enn[0]=memory(n,n);
	 Ann[0]=memory(n,n);
	 Bnn[0]=memory(n,n);

    Init_Ed(Enn[0],n,n);
    prisv(Ann[9],Ann[0],n,n);        //Ann[0]=Ann[9]


   for(i=0;i<n;i++)
  {
   if(i) {umn(Ann[9],Bnn[0],Ann[0],n,n,n);}
   pechat(Ann[0],n,n,"Ann0 :");getch();

   *(Pnn[0]+i)=Sled(Ann[0],n,n)/(i+1);
    printf("\t Pnn %d=%f \n",i,*(Pnn[0]+i));
    umnn(Enn[0],*(Pnn[0]+i),n,n);
    sub(Ann[0],Enn[0],Bnn[0],n);
    pechat(Bnn[0],n,n,"Bnn0 :");
    Init_Ed(Enn[0],n,n);
  }
 pechat(Pnn[0],n,1,"Pnn :");
 inpol(pol[0],"FUCK");
// Mat2Pol(pol[0],Pnn[0],n);
 normpol(pol[0]);
 clrscr();
 outpol(pol[0],"Hello");
 getch();
 copypol(pol[0],ism);
 hi=Verdict_Newtona(pol[0],1);
 getch();
 PolSetMin(pol[0]);
 lo=-Verdict_Newtona(pol[0],1);
 printf("\nA:%f B:%f",hi,lo);
 getch();
 NK=System_Shturma(ism);

 //Master(ism);


 Korni(ism);
 getch();
// Master(ism);

 free(Ann[9]);
 free(Pnn);
 free(Ann[0]);
 free(Bnn[0]);
 free(Enn[0]);

  farfree(polm);
  farfree(pol[0]);
  farfree(pol[1]);
  farfree(pst[0]);
  farfree(pst[1]);
  farfree(pold);
  farfree(ish);
  farfree(ish);
  for(i=0;i<16;i++)
  {
   farfree(work[i]);
  }




// for(i=3;i<n;i++)
// {
// }
  getch();
  q='n';
//   printf("\n Снова ? :Y/N  ");
//   scanf("%s",&q);
} while( q!='n');
}
/*--Умножение матриц--*/
void umn(FU a,FU b,FU c,int n,int m,int q)
{
  int i,j,k;
  float s;

  for(i=0;i<n;i++)
    for(j=0;j<q;j++)
      {
	for(s=k=0;k<m;k++)
	  s+=*(a+i*m+k)**(b+k*q+j);
	*(c+i*q+j)=s;
      }
 }
/*--Умножение матрицы на число--*/
void umnn(FU a,float l,int m,int n)
{
  int i;

  for(i=0;i<m*n;i++)
    {
      *(a+i)*=l;
    }
}
/*--Вычитание матриц--*/
void sub(FU a,FU b,FU c,int n)
{
  int i,j;

  for(i=0;i<n;i++)
    for(j=0;j<n;j++)
      *(c+i*n+j)=*(a+i*n+j)-*(b+i*n+j);
}

void pechat(FU a,int m,int n,char *str)
{
  int i,j;

  puts(str);
  for(i=0;i<m;i++)
    {
      for(j=0;j<n;j++)
	printf("%6.3f ",*(a+i*n+j));
      printf("\n");
    }
}
/*--Выделение памяти--*/
FU memory(int m,int n)
{
  FU mas;

  mas=(FU)malloc(m*n*sizeof(float));
  return(mas);
}
/*--Инициализация матрицы--*/
void inition(FU a,int m,int n)
{
  int j;

  for(j=0;j<m*n;j++)
    *(a+j)=rnd;
}

void InitTest_4x4(FU a,int m,int n)
{
  int j,i;

      *(a)=2;  *(a+1)=-1;  *(a+2)=1;  *(a+3)=2;
    *(a+4)=0;   *(a+5)=1;  *(a+6)=1;  *(a+7)=0;
   *(a+8)=-1;   *(a+9)=1; *(a+10)=1; *(a+11)=1;
   *(a+12)=1;  *(a+13)=1; *(a+14)=1; *(a+15)=0;
}

float Sled(FU a,int m,int n)
{
  int i,j;float sl=0;

  for(i=0;i<m;i++)
    {
      for(j=0;j<n;j++)
      if(i==j) sl+=*(a+i*n+j);
    }
    return (sl);
}

void Init_Ed(FU a,int m,int n)
{
  int i,j;

  for(i=0;i<m;i++)
    {
      for(j=0;j<n;j++)
      if(i==j) *(a+i*n+j)=1;
      else *(a+i*n+j)=0;

    }
}






/*--Присвоение матрице B значения матрицы A--*/
void prisv(FU a,FU b,int m,int n)
{
  int i;

  for(i=0;i<m*n;i++)
    *(b+i)=*(a+i);
}
/*--Ручной ввод--*/
void handinit(FU a,int n)
{
  int x,y,i,j;
  float p;

  for(i=0,y=3,x=3,j=1;i<n*n;i++,j++)
    {
      x+=8;
      gotoxy(x,y); scanf("%f",&p);
      *(a+i)=p;
      if(j==n) { j=0; x=3; y+=2; }
    }
}
void prrivet(FU a,FU b,int n,int t)
{
  int i,j;
  for(i=0;i<t;i++)
    {
      for(j=0;j<t;j++)
	*(a+i*t+j)=*(b+i*n+j);
    }
}

void normpol(struct polinom *pol)
{int bufn,i,j,fn,k,kk,max,maxi,polk;
 float bufa,fa; polk=pol->k;
  if ((pol->k)<=1) return ;
  for(i=0;i<polk;i++)
  { for(j=i,max=pol->n[i],maxi=i;j<polk;j++)
     if (max<(pol->n[j])){maxi=j;max=pol->n[j];}
     bufn=pol->n[i];bufa=pol->a[i];
     pol->n[i]=pol->n[maxi];pol->a[i]=pol->a[maxi];
     pol->n[maxi]=bufn;pol->a[maxi]=bufa;
   }
   for(i=0;i<polk;i++)
    {while((pol->n[i])==(pol->n[i+1]))
	{pol->a[i]+=pol->a[i+1];
	 for(j=i+1;j<polk;j++) {pol->n[j]=pol->n[j+1];pol->a[j]=pol->a[j+1];}
	 polk--;pol->k=polk;
	 }
     }

   for(i=0;i<polk;i++)
    {if((fabs(pol->a[i]))<EPS)
	{polk--;pol->k=polk;
	 for(j=i;j<polk;j++)
	  {pol->n[j]=pol->n[j+1];
	   pol->a[j]=pol->a[j+1];
	  }
	 }
     }
}
void inpol(pln *el,char ch[30])
{ int i,fn;float fa;
  textbackground(9);el->k=0;
  cprintf("Input %s polinom until 20 elements,end [0] :\n\r",ch);
  for(i=0;i<20;i++)
  {scanf("%f",&fa);
   if (fa==0) break;
   scanf("%d",&fn);el->k+=1;;el->n[i]=fn;el->a[i]=fa;
  }
}
void Mat2Pol(pln *el,FU a,int n)
{ int i,fn;float fa;
  for(i=0;i<n;i++)
  {
  el->n[i]=i;el->a[i]=*(a+n-1-i);
  }
  el->n[n]=n;el->a[i]=1;
  el->k=n+1;
}

void PolumnN(pln *el,float n)
{ int i,j;
  j=el->k;
  for(i=0;i<j;i++)
  {
     el->a[i]*=n;
  }
}
void SimplCreate(pln *el,float n)
{ int i,j;
  el->k=2;
  el->a[0]=1.0;
  el->a[1]=-n;
}

void PolSetMin(pln *el)
{ int i,j;
  j=el->k;
  for(i=1;i<j;i++)
  {
    if(el->a[i]>0) el->a[i]*=-1.0;
  }
}


void Proizvodnaia95(pln *el,pln *el1)
{ int i,j,old_n,new_n;
  j=el->k;
  for(i=0;i<j;i++)
  {
  if(el->n[i]==0) {el1->k=el->k-1;break;}
 else {old_n=el->n[i];
      el1->n[i]=old_n-1;
      el1->a[i]=el->a[i]*old_n;}
  }
}


//void outpol(pln *el,char ch[30])
//{ int i,ii;
//  textbackground(12);cprintf("%s\n\r",ch);textbackground(1);
//  ii=el->k;
//  for(i=0;i<ii;i++)
//  {cprintf("%+.3f*x^(%d)",el->a[i],el->n[i]);} cprintf("\b\n\r");
//  textbackground(0);
//}

void outpol(pln *el,char ch[30])
{ int i,ii;
  printf("%s\n",ch);
  ii=el->k;
  for(i=0;i<ii;i++)
  {printf("%+.3f*x^(%d)",el->a[i],el->n[i]);} printf("\n");
  //textbackground(0);
}



void copypol(pln *el,pln *el1)
{int i,polk;polk=el->k;
 for(i=0;i<polk;i++)
  {el1->n[i]=el->n[i];el1->a[i]=el->a[i];}
   el1->k=el->k;
}

void delpol(pln *pol1,pln *pol2)
{
  int i,j,jj,fn;float fa;
  pold->k=0;
  for(i=0;(pol1->n[0])>=(pol2->n[0])&&(pol1->k)>0;i++)
  { pold->n[i]=(pol1->n[0])-(pol2->n[0]);
    pold->a[i]=(pol1->a[0])/(pol2->a[0]);
    pold->k+=1;
    for(j=(pol1->k),jj=0;jj<(pol2->k);j++,jj++)
     {fn=pold->n[i];fa=-1.0*(pold->a[i]);
      fn+=(pol2->n[jj]);
      pol1->n[j]=fn;
      fa*=(pol2->a[jj]);
      pol1->a[j]=fa;
      pol1->k+=1;
     }
    normpol(pol1);
    //outpol(pol1,"pol1");outpol(pold,"pold");getch();
  } copypol(pold,pol2);
}

/* Конец */

float point(pln *el,float num)
{ float buf,pr;int i;
  for(i=0,buf=0;i<(el->k);i++)
 { if(num!=0) pr=pow(num,(el->n[i])*1.0);
   else { if(el->n[i]>0) pr=0.0;
	  else pr=1.0;
	}
 buf+=(el->a[i])*pr;
	}
 return buf;
}




float Verdict_Newtona(pln *el,int n)
{ int i,j,k,m;float x,koef,h1,h2;
 m=el->k;
// printf("\nM%d\n",m);
 k=el->n[0];
 copypol(el,work[0]);
 for(i=1;i<k;i++)
 {Proizvodnaia95(work[i-1],work[i]);
 }
 for(i=0;i<k;i++)
{ outpol(work[i],"..");
}

//Step 1 High terminator
j=work[k-1]->k;
//printf(".....K:%d",j);
if(j==2)
  { x=(-(work[k-1]->a[1]))/(work[k-1]->a[0]);
  //  printf(".....X:%f[%d]",x,j);
    koef=x;
   }
if(j==1)
  { x=work[k-1]->a[0];
    if(x>0) koef=0;
    else koef=-0;
   // printf("\n.....X:%f\n",koef);
   }
  for(i=k-1;i>=0;i--)
  {
    m1: h1=point(work[i],koef);
    printf("\n Pol:%d K:%d in point %f =%f",i,work[i]->k,koef,h1);
    if(n){    if(h1<0){ koef+=1.0;/*getch();*/ goto m1;
		       }
	 }
     else    if(h1<0){ koef-=1.0;/*getch();*/ goto m1;
		       }
  }
  //PolSetMin(el);

 //Verdict_Newtona(pln *el)
return ((float)koef);
}

int System_Shturma(pln *el)
{  int i,k,n,m,b1,b2;float temp[16],a,hh;
 copypol(el,work[0]);
 Proizvodnaia95(el,work[1]);
 //outpol(work[0],"..PR");
 k=el->n[0];
 printf("\nShturm\t K%d\n",k);
 for(i=2;i<k+1;i++)
 {
   copypol(work[i-1],pst[1]);
   copypol(work[i-2],pst[0]);
       delpol(pst[0],pst[1]);
     copypol(pst[0],work[i]);
     //outpol(work[i],"..");
     //  getch();
   //   copypol(pst[0],work[i]);
     PolumnN(work[i],-1.0);
 }
// outpol(work[2],"doirghn");
  for(i=0;i<k+1;i++)
{ outpol(work[i],"");
}
  hh=hi;
m2:;
  for(i=0;i<k+1;i++)
  {temp[i]=point(work[i],hh);
   printf("t:%f\n",temp[i]);
  }
     n=0;
      if(temp[0]>0) m=1;
      else m=0;
  for(i=1;i<k+1;i++)
  {  a=temp[i];
     if(a<0) {
	      if(m==1){ n++;m=0;
		      }
	     }
     else    {
	      if(m==0){ n++;m=1;
		      }
	     }
  }
 printf("Hi : %d",n);
 if(hh!=lo){b2=n;hh=lo;goto m2;}
 b1=n;
 b1-=b2;
 printf("Chislo korney:%d",b1);
 return (b1);
}

void Korni(pln *el)
{ float lok,hik; int i,j;
  int nu=0,kf;
  float a,b;
  float eps=0.0;
  float x,fx;
  int ir,kff=0;
      kf=0;
      copypol(el,ish);
      //printf("\nLO:%f ,HI:%f NK:%d",lo,hi,NK);
      //outpol(ish,""); getch();

      ir=zero(lo,hi,eps,el,&x,&fx);
      if(ir) printf("\n  %2d\t %3.1g\t %3.1g\t  - нет корня !!! ",
						     nu+1,lo,hi);
       else printf("\n %2d\t% 3.1f\t% 3.1f\t%2d\t% 8.7f\t% 8.3e ",
					     nu+1,lo,hi,kf,x,fx);
      kff+=kf;
      printf("\n Res:%f",point(el,x));

  copypol(ish,pst[0]);
  SimplCreate(pst[1],x);
  delpol(pst[0],pst[1]);
  printf("\nREAL");
  outpol(pst[1],"");
  getch();
  puts("Ostatok");
  getch();
  outpol(pst[0],"");
//  copypol(pst[1],ish);



}

void Master(pln *el)
{ float lok,hik; int i,j;
  int nu=0,kf;
  float a,b;
  float eps=0.0;
  float x,fx;
  int ir,kff=0;
      kf=0;
       copypol(el,ish);
       normpol(ish);
  do{
      printf("\nLO:%f ,HI:%f NK:%d",lo,hi,NK);
      outpol(ish,""); getch();
      ir=zero(lo,hi,eps,ish,&x,&fx);
      if(ir) printf("\n  %2d\t %3.1g\t %3.1g\t  - нет корня !!! ",
						     nu+1,lo,hi);
       else printf("\n %2d\t% 3.1f\t% 3.1f\t%2d\t% 8.7f\t% 8.3e ",
					     nu+1,lo,hi,kf,x,fx);
      kff+=kf;
      printf("\n Res:%f",point(ish,x));

  copypol(ish,pst[0]);
  SimplCreate(pst[1],x);
  delpol(pst[0],pst[1]);
  printf("\nREAL");
  outpol(pst[1],"");
  getch();
  puts("Ostatok");
  getch();
  outpol(pst[0],"");
  copypol(pst[1],ish);

 hi=Verdict_Newtona(pst[1],1);
 getch();
 PolSetMin(pst[1]);
 lo=-Verdict_Newtona(pst[1],1);
 printf("\nLO:%f HI:%f",lo,hi);
   NK--;
  }while(NK>0);
//  copypol(pst[1],ish);



}

int zero(float a, float b, float eps,pln *el, float *x, float *fx)
{
  float c,d,fa,fb,fc,fd,a0,a1,a2,ak,bk,ck,sa,d1,d2,buf;

// правильно задан интервал (a,b)?
  fa=point(el,a); fb=point(el,b);
  if(fa*fb>0) return 1;  // корня нет
  eps+=FLT_EPSILON; // защищаемся от eps=0
  c=FLT_MAX;
  do
    {
      if(c==FLT_MAX)  // запускаем метод хорд
       { c=a-fa*(b-a)/(fb-fa); fc=point(el,c); }
      if((c-a<eps)||(b-c<eps)||(fabs(fc)<1e-6)) { *x=c; *fx=fc; return 0; }
// квадратичная интерполяция
      a2=a*(fc-fb)+b*(fa-fc)+c*(fb-fa);
      ak=a*a; bk=b*b; ck=c*c;
      a1=ak*(fb-fc)+bk*(fc-fa)+ck*(fa-fb);
      a0=ak*(b*fc-c*fb)+bk*(c*fa-a*fc)+ck*(a*fb-b*fa);
      sa=sqrt(a1*a1-4*a2*a0);
      d1=-(a1+sa)/(2*a2); d2=-(a1-sa)/(2*a2);
      if((d1>a)&&(d1<b)) d=d1;    // ищем подходящий корень
       else if((d2>a)&&(d2<b)) d=d2;
	     else if(fa*fc<0) { a=c; fa=fc; }  // нет подходящего корня!
		   else { b=c; fb=fc; }
      fd=point(el,d);
      if((d-a<eps)||(b-d<eps)) { *x=d; *fx=fd; return 0; }
      if(c>d) { buf=c; c=d; d=buf; buf=fc; fc=fd; fd=buf; }
      if(fc*fd<0)  // разные знаки
	   { a=c; fa=fc; b=d; fb=fd; c=FLT_MAX; }
       else if(fa*fc<0) { b=d; fb=fd; }
	     else { a=c; fa=fc; c=d; fc=fd; }
    }
  while((b-a)>eps);
  if(fabs(fa)>fabs(fb)) { *x=b; *fx=fb; }
   else { *x=a; *fx=fa; }
  return 0;  // нашли корень!
}

