#include <stdio.h> 
#include <math.h>

int N; double A[100],X[100],Y[100];

double a,b,c,d,p,q,err,th,zx,zy,

       D,DD,BA3,BA4,CA,DA,AA,B,ApB,AsB,

       X1,X2,X3,X4,Y1,Y2,Y3,Y4;

double aa,bb,cc,dd,ee,r,pp,qq,rr,TX,TY;

double fzx,fzy;char stop;int flag;



FILE *fpA; FILE *fpB;

/***********************************************************
 
     　このセクションは　２次方程式の解を求めます。

            　　　  a*x*x + b*x + c = 0

       まず判別式　D=(b/a/2)*(b/a/2)-(c/a) を計算して、

　　　D の正負の値から、実数根か虚数根かを 判定します。

************************************************************/

int EQ2ji(void) { 

a=A[2];b=A[1];c=A[0];


D=(b/a/2)*(b/a/2)-(c/a);

printf("\n\n D = (b/a/2)*(b/a/2)-(c/a)=%f \n\n",D);

if (D<0.0) goto KYOSU;

X1=-b/a/2 + sqrt(D);Y1=0;
X2=-b/a/2 - sqrt(D);Y2=0;

goto NEXT;

KYOSU:

X1=-b/a/2 ; Y1=   sqrt(-D);
X2=-b/a/2 ; Y2= - sqrt(-D);


NEXT:

X[1]=X1;X[2]=X2;
Y[1]=Y1;Y[2]=Y2; return 0; }




/*******************************************************
 
     このセクションは　３次方程式の解を求めます。

            a*x*x*x + b*x*x + c*x + d  = 0

********************************************************/
/*                                                     */
/*    This section computes the roots of               */
/*                                                     */
/*              ax*x*x+bx*x+cx+d = 0                   */
/*                                                     */
/*******************************************************/


int shift(void) { X1=X1-BA3;X2=X2-BA3;X3=X3-BA3; return 0; }

double root3( double x ) { err=1;if ( x < 0 ) { return  - exp(log(-x)/3); }
                                 if ( x > 0 ) { return    exp(log(x)/3 ); }

                                 return err=-1; }

int case1(void) { printf("\n\n ** Case 1**１つ実根で残る２つは虚数根 ** \n\n");

                  X1=root3(2*q);Y1=0;
 
                  X2=-X1/2;  Y2=-X2*root3(3);

                  X3=X2;     Y3=-Y2; shift();return 0;}

int case2(void) { 

                  X1=0;Y1=0; 

      if (p>=0) { printf("\n\n ** Case 2-1 **３つとも実根** \n\n");

                  X2=sqrt(3*p); Y2=0;

                  X3=-X2;Y3=0; shift();return 0;}

                  printf("\n\n ** Case 2-2 **実根１つと虚数２根** \n\n");

                  X2=0;Y2=sqrt(-3*p);

                  X3=0;Y3=-Y2; shift();return 0;}

int case3(void) { printf("\n\n ** Case 3**３つの重根で実数根**\n\n");

                  X1=0;Y1=0;X2=0;Y2=0;X3=0;Y3=0;  shift();return 0;}

int case4(void) { printf("\n\n ** Case 4** ３つとも実数で２つが重根 ** \n\n");

                  X1=2*q/p; Y1=0;

                  X2=-q/p;Y2=0;

                  X3=X2;Y3=0; shift();return 0;}



int case5(void) { printf("\n\n ** Case 5** 実根１つ、虚根２つ** \n\n");

                  AA=root3(q+sqrt(DD));B=root3(q-sqrt(DD)); 

                  X1=AA+B;Y1=0;

                  X2=-X1/2;X3=X2;

                  Y2=sqrt(3)*(AA-B)/2;

                  Y3=-Y2; 

                  shift();return 0;}



int case6(void) { double pp=p,qq=1;

                  printf("\n\n ** Case 6** ３つとも実根 ** \n\n");

                  if ( q > 0 )   th=atan(sqrt(-DD)/q);
                  if ( q < 0 ) { th=atan(sqrt(-DD)/(-q));qq=-1; }

                  if (p<0) pp=-p;

                  ApB=2*sqrt(pp)*cos(th/3)*qq;AsB=2*sqrt(pp)*sin(th/3)*qq;

                  X1=ApB;
                  X2=-X1/2-sqrt(3)*AsB/2;
                  X3=-X1/2+sqrt(3)*AsB/2;

                  Y1=0;Y2=0;Y3=0;shift();return 0;}

int sanji(void) { 

                  BA3=b/a/3;

                  CA=c/a;

                  DA=d/a;

                  p=BA3*BA3-CA/3;

                  q=-BA3*BA3*BA3+BA3*CA/2-DA/2;

                  DD=q*q-p*p*p;

                  printf("\n\n BA3=b/a/3= %f\n\n",BA3);
                  printf("\n\n CA =c/a  = %f\n\n",CA);
                  printf("\n\n DA =d/a  = %f\n\n",DA);


                  printf("\n\n p  = BA3*BA3-CA/3 = %f\n\n",p);
                  printf(    " q  =-BA3*BA3*BA3+BA3*CA/2-DA/2 = %f\n\n",q);
                  printf("     DD = %f\n",DD);
                  
      
if(p==0)  { if( q==0 ) { case3(); return 3;} case1();return 1 ;}
if(q==0)  { case2(); return 2; } 
if(DD==0) { case4(); return 4; }  
if(DD>0)  { case5(); return 5; }
if(DD<0)  { case6(); return 6; }  return -1; }       
          



int EQ3ji(void){  


a=A[3];b=A[2];c=A[1];d=A[0];
      
sanji();  

X[1]=X1;Y[1]=Y1;
X[2]=X2;Y[2]=Y2;
X[3]=X3;Y[3]=Y3;

return 0;}







/*******************************************************
 
     このセクションは　４次方程式の解を求めます。

         aa*x*x*x*x+bb*x*x*x+c*x*x+dd*x+ee = 0 

********************************************************/

/*******************************************************/
/*                                                     */
/*    This program computes the roots of               */
/*                                                     */
/*      aa*x*x*x*x+bb*x*x*x+c*x*x+dd*x+ee = 0          */
/*                                                     */
/*******************************************************/




double real( double AX, double AY, double BX, double BY ){ return AX*BX-AY*BY; }

double imag( double AX, double AY, double BX, double BY ){ return AX*BY+AY*BX; }


int shift4(void) { X1=X1-BA4;X2=X2-BA4;X3=X3-BA4;X4=X4-BA4; return 0; }


int ROOT2C( double x, double y) { double DDD=sqrt(x*x+y*y),zero=0.00000000001;

if(x>zero){DDD=sqrt(DDD+x);if(y> zero){TX=DDD/sqrt(2);TY=y/DDD/sqrt(2);goto ENDD;}
                           if(y<-zero){TX=DDD/sqrt(2);TY=y/DDD/sqrt(2);goto ENDD;}
                                                        TX=sqrt(x);TY=0;goto ENDD;}

if(x<-zero){DDD=sqrt(DDD-x);if(y> zero){TX=y/DDD/sqrt(2);TY=DDD/sqrt(2);goto ENDD;}
                            if(y<-zero){TX=y/DDD/sqrt(2);TY=DDD/sqrt(2);goto ENDD;}
                                                       TX=0;TY=sqrt(-x);goto ENDD;}         

if(y> zero) { TX=sqrt( y/2);TY=  TX;goto ENDD ;}
if(y<-zero) { TX=sqrt(-y/2);TY= -TX;goto ENDD ;} TX=0;TY=0;

ENDD:printf( "\n x=%f y=%f  root2( x+iy ) = %f  %f ",x,y,TX,TY );

fzx=TX*TX-TY*TY;fzy=2*TX*TY;

printf( "\n x=%f y=%f   checked\n",fzx,fzy);return;}




int checkf3( double AAA, double BBB){  double AA2,BB2,AA3,BB3;

AA2=real(AAA,BBB,AAA,BBB);BB2=imag(AAA,BBB,AAA,BBB);
AA3=real(AAA,BBB,AA2,BB2);BB3=imag(AAA,BBB,AA2,BB2);
fzx=a*AA3+b*AA2+c*AAA+d; fzy=a*BB3+b*BB2+c*BBB; 

printf("\n\n Checked AAA=%f  BBB=%f  fzx=%f fzy=%f\n",AAA,BBB,fzx,fzy); return; }



int checkf4( double AAA, double BBB){  double AA2,BB2,AA3,BB3,AA4,BB4;

AA2=real(AAA,BBB,AAA,BBB);BB2=imag(AAA,BBB,AAA,BBB);
AA3=real(AAA,BBB,AA2,BB2);BB3=imag(AAA,BBB,AA2,BB2);
AA4=real(AA2,BB2,AA2,BB2);BB4=imag(AA2,BB2,AA2,BB2);

fzx=aa*AA4+bb*AA3+cc*AA2+dd*AAA+ee; fzy=aa*BB4+bb*BB3+cc*BB2+dd*BBB; 
printf("\n\n Checked AAA=%f  BBB=%f  fzx=%f fzy=%f\n",AAA,BBB,fzx,fzy); return; }







int yoji(void) { double T1X,T1Y,T2X,T2Y,T3X,T3Y,Q2X,Q3X,Q2Y,Q3Y;

                        printf(" \n aa=%f",aa);printf(" \n bb=%f",bb);  
printf(" \n cc=%f",cc); printf(" \n dd=%f",dd);printf(" \n ee=%f",ee);  

BA4=bb/aa/4;

pp=3*BA4*BA4-cc/aa/2;printf("\n\n pp=%f\n",pp);
qq=-(BA4*BA4*BA4-BA4*cc/aa/4+dd/aa/8);printf(" qq=%f\n",qq);
rr=pp*pp/4+(3*BA4*BA4*BA4*BA4-BA4*BA4*cc/aa+BA4*dd/aa-ee/aa)/4;

printf(" rr=%f\n\n",rr);stop=getchar();stop=getchar();if(stop=='s')return;

a=1;b=-pp;c=rr;d=-qq*qq; sanji();  

printf("  \n   X1=%f  Y1=%f\n", X1,Y1); 
    printf("   X2=%f  Y2=%f\n", X2,Y2);
    printf("   X3=%f  Y3=%f\n", X3,Y3);printf("\n\n");

checkf3(X1,Y1);checkf3(X2,Y2);checkf3(X3,Y3);

printf("\n\n");stop=getchar(); stop=getchar();if( stop=='s' ) return;  

ROOT2C(X1,Y1); T1X=TX;T1Y=TY;  
ROOT2C(X2,Y2); T2X=TX;T2Y=TY;  
ROOT2C(X3,Y3); T3X=TX;T3Y=TY;

MORETX: Q2X=real(T1X,T1Y,T2X,T2Y);Q2Y=imag(T1X,T1Y,T2X,T2Y);
        Q3X=real(T3X,T3Y,Q2X,Q2Y);Q3Y=imag(T3X,T3Y,Q2X,Q2Y);

printf(" \n\nQ3X=%f  Q3Y=%f  qq=%f\n\n",Q3X,Q3Y,qq );

Q3X=Q3X*qq; if (Q3X<0) { T1X=-T1X;T1Y=-T1Y;goto MORETX; }


X1= T1X+T2X+T3X;Y1= T1Y+T2Y+T3Y;   X2= T1X-T2X-T3X;Y2= T1Y-T2Y-T3Y;
X3=-T1X+T2X-T3X;Y3=-T1Y+T2Y-T3Y;   X4=-T1X-T2X+T3X;Y4=-T1Y-T2Y+T3Y;

stop=getchar();stop=getchar(); if( stop=='s' ) return; shift4();return;}





int EQ4ji(void){  


aa=A[4];bb=A[3];cc=A[2];dd=A[1];ee=A[0];
      
yoji();  

X[1]=X1;Y[1]=Y1;
X[2]=X2;Y[2]=Y2;
X[3]=X3;Y[3]=Y3;
X[4]=X4;Y[4]=Y4;


checkf4(X1,Y1);checkf4(X2,Y2);checkf4(X3,Y3);checkf4(X4,Y4); return 0;}












int main(void) { int i,j,k,n,s,up,dgt,M;char C[100],cc;double H;

fpA=fopen("A.txt","r");
fpB=fopen("B.txt","w");

printf("\n\n"); fprintf(fpB,"\n\n"); 

/**************************************************************************

        Input File (A.txt) の 最初の16行は コメントとして無視します。

***************************************************************************/

for (i=0;i<16;i++)  { fgets(C,99,fpA); printf("%s",C);fprintf(fpB,"%s",C);}

/**************************************************************************

       Input File (A.txt) の 17行目を 読み込み、N の値を決定します。

***************************************************************************/

fgets(C,99,fpA); 

i=-1;N=0;

NEXTA: i=i+1; printf("%c",C[i]);fprintf(fpB,"%c",C[i]); 

       if(C[i]=='=') goto NEXTB; goto NEXTA;

NEXTB:i=i+1;printf("%c",C[i]);fprintf(fpB,"%c",C[i]); dgt=10;

if (  C[i] == '1' ) dgt=1;if (  C[i] == '2' ) dgt=2;if (  C[i]== '3' ) dgt=3;
if (  C[i] == '4' ) dgt=4;if (  C[i] == '5' ) dgt=5;if (  C[i]== '6' ) dgt=6;
if (  C[i] == '7' ) dgt=7;if (  C[i] == '8' ) dgt=8;if (  C[i]== '9' ) dgt=9;

if(dgt>9){ if(N==0) goto NEXTB; if(N >0) goto NEXTD;}

if(N==0) { N=dgt;goto NEXTB;}

N=10*N+dgt; 

NEXTD:   /*** N の値は　99 以下としてここで決定します　***/


for (k=0;k<=N;k++) {    n = N - k ;

/**************************************************************************

       Input File (A.txt) の次の１行を読み込み、A[n] の値を決定します。

　　　　　　　　n の値は　n = N から　n = 0 までです。

***************************************************************************/

fgets(C,99,fpA);


i=-1;s=1;M=0;H=1;A[n]=0;

NEXTAA: i=i+1; printf("%c",C[i]);fprintf(fpB,"%c",C[i]); 

       if(C[i]=='=') goto NEXTBB; goto NEXTAA;

NEXTBB: i=i+1;printf("%c",C[i]);fprintf(fpB,"%c",C[i]); 

        if(C[i]=='.') { A[n]=M;  goto NEXTP; }

        if(C[i]=='-') { s=-1;  goto NEXTBB; }

        dgt=10;

if (  C[i] == '1' ) dgt=1;if (  C[i] == '2' ) dgt=2;if (  C[i]== '3' ) dgt=3;
if (  C[i] == '4' ) dgt=4;if (  C[i] == '5' ) dgt=5;if (  C[i]== '6' ) dgt=6;
if (  C[i] == '7' ) dgt=7;if (  C[i] == '8' ) dgt=8;if (  C[i]== '9' ) dgt=9;
if (  C[i] == '0' ) dgt=0;

if(dgt>9){ if(M==0) goto NEXTBB; if(M >0) { A[n]=M; goto NEXTDD;} }

M=M*10+dgt; goto NEXTBB;

NEXTP:i=i+1;printf("%c",C[i]);fprintf(fpB,"%c",C[i]); 

        dgt=10;

if (  C[i] == '1' ) dgt=1;if (  C[i] == '2' ) dgt=2;if (  C[i]== '3' ) dgt=3;
if (  C[i] == '4' ) dgt=4;if (  C[i] == '5' ) dgt=5;if (  C[i]== '6' ) dgt=6;
if (  C[i] == '7' ) dgt=7;if (  C[i] == '8' ) dgt=8;if (  C[i]== '9' ) dgt=9;
if (  C[i] == '0' ) dgt=0;

if(dgt>9)  goto NEXTDD; 

H=H/10; A[n]=A[n]+dgt*H; goto NEXTP;


NEXTDD: A[n]=s*A[n]; }

 fgets(C,99,fpA); printf("%s",C);fprintf(fpB,"%s",C);

printf("\n\n**********************************************\n");
printf("               %d 次方程式の解　　　　　　　　　　\n",N);
printf(    "**********************************************\n\n");


fprintf(fpB,"\n\n**********************************************\n");
fprintf(fpB,"               %d 次方程式の解　　　　　　　　　　\n",N);
fprintf(fpB,    "**********************************************\n\n");


for (i=0;i<=N;i++) { printf(     "A[%d]=%f  \n",N-i,A[N-i]) ;

                    fprintf(fpB, "A[%d]=%f  \n",N-i,A[N-i]) ;                                          }

if(N==2) EQ2ji( );
if(N==3) EQ3ji( );
if(N==4) EQ4ji( );

if(N>4) {  printf("\n\n**********************************************\n\n");
          fprintf(fpB,"\n\n**********************************************\n\n");


           printf(     " \n\n ５次方程式以上は解析的には求まりません。\n\n");
          fprintf(fpB, " \n\n ５次方程式以上は解析的には求まりません。\n\n");

          fprintf(fpB, " \n\n 現在工事中です。しばらくお待ちください。\n\n");
           printf(     " \n\n 現在工事中です。しばらくお待ちください。\n\n"); 

          goto FINAL;   }

printf("\n\n"); fprintf(fpB,"\n\n"); 

for (i=1;i<=N;i++) { printf(     "X[%d]=%f   Y[%d]=%f \n\n",i,X[i],i,Y[i]);

                    fprintf(fpB, "X[%d]=%f   Y[%d]=%f \n\n",i,X[i],i,Y[i]);}
                                                                         

FINAL:

 printf(    "\n**********************************************\n\n");
fprintf(fpB,"\n**********************************************\n\n");
printf("\n\n");fprintf(fpB,"\n\n");fclose(fpA);fclose(fpB);


cc=getchar( ); if(cc=='s') return 0;return 0; }