/*******************************************
 * Mathematik II - Praktikum 02	- Grp MAMF *
 * FH Aachen - Fachbereich 05 - M.Cuvelier *
 * Thema: POLYNOM - NULLSTELLEN	           *
 * Ersterstellungsdatum: 28.04.2003        *
 * Letzte Ueberarbeitung: 29.05.2003       *
 *******************************************/

#include <stdio.h>
#include <math.h>
#include <conio.h>						// wird fuer getch() benoetigt
#include <stdlib.h>						// wird fuer malloc benoetigt


void Horner_W2(int n,double c[], double x_quer, double *pn, double *pn_strich);

double Horner_Abdiv(int n,double c[], double x_quer, double *d);

int Newton_Poly_Nst(int n,double c[], double x_quer,double eps_x, double eps_f,int ITMAX,double *x_nst);


int main (void)
{
	char auswahl;						// Variable fuer die Menustruktur
	int n, I;							// (Lauf-)Variable
	double *c,x_quer;					// Uebergabevariablen

	do
	{
		system("cls");
		printf("\n\n POLYNOM-NULLSTELLEN - Praktikum Aufgabe 2 - M.Cuvelier - Gruppe MAMF");
		printf("\n ====================================================================");
		printf("\n\n <1> --> Berechnung eines Polynoms und seine Ableitung");
		printf("\n <2> --> Berechnung des Deflationspolynoms");
		printf("\n <3> --> Berechnung einer Polynom-Nullstelle mit dem Newtonverfahren");
		printf("\n <4> --> Berechnung und Ausdruck saemtlicher reeler Nullstellen");
		printf("\n <5> --> Berechnung und Ausdruck der Nullstellen von p7*(x)");
		printf("\n <6> --> Berechnung der Nullstellen fuer ein beliebiges Polynom");
		printf("\n <9> --> Beenden");
		printf("\n\n Bitte treffen Sie eine Auswahl: ");

		fflush(stdout);
		auswahl= getch(); 
		printf("%c\n", auswahl);

		switch (auswahl)
		{
			case '1':							// Polynomberechnung und seine Ableitung	
			{
				double pn,pn_strich;

				system("cls");
				printf("\n\n <1> BERECHNUNG EINES POLYNOMS UND SEINE ABLEITUNG");
				printf("\n -------------------------------------------------");
				printf("\n\n Geben Sie bitte Anzahl der Polynomkoeffizienten n ein: ");
				scanf("%i", &n);
				rewind(stdin);
				printf("\n");

				c=(double*)malloc(sizeof (double)*(n+1));	// Dynam. Speicherallokation

				for (I=n;I>=0;I--)
				{
					printf("\n Geben Sie nun den %i.Polynomkoeffizienten ein: ", n-I+1);
					scanf("%lf", &c[I]);
					rewind(stdin);
				}
				printf("\n\n Geben Sie nun die Auswertestelle x_quer ein: ");
				scanf("%lf", &x_quer);
				rewind(stdin);

				system("cls");
				printf("\n\n <1> ANZEIGEN DER WERTE und BERECHNUNG");
				printf("\n -------------------------------------");
				printf("\n\n Von Ihnen wurden folgende Werte eingegeben:\n");
				for(I=n; I>=0; I--) printf("\n     c%i = %lf", n-I+1, c[I]);
				printf("\n x_quer = %lf\n", x_quer);

				Horner_W2(n,c,x_quer,&pn,&pn_strich);		// Unterprg-Aufruf

				printf("\n\n Berechnet wurde folgendes Polynom: pn  = %lf", pn);
				printf("\n und seine Ableitung:               pn' = %lf\n\n ", pn_strich);

				system("pause");
				break;
			}


			case '2':										// Deflationspolynomberechnung
			{
				double *d,pn;								// d  -> Vektor - Koeffizienten des Deflationspolynoms
															// pn -> Y-Wert für die Stelle x_quer 
				system("cls");
				printf("\n\n <2> BERECHNUNG DES DEFLATIONSPOLYNOMS");
				printf("\n -------------------------------------");
				printf("\n\n Geben Sie nun bitte den Grad n des Polynoms ein: ");
				scanf("%i",&n);
				rewind(stdin);
		
				c=(double*)malloc(sizeof (double)*(n+1));	// Speicherplatzallokation - Polynomkoeffizienten
				d=(double*)malloc(sizeof (double)*(n+1));	// Speicherplatzallokation Deflationspolynomkoeff.
			
				for (I=n;I>=0;I--)							// Rueckwaertsrekursion
				{
					printf("\n Geben Sie nun den %i. Polynomkoeffizienten ein: ",n-I+1);
					scanf("%lf",&c[I]);
					rewind(stdin);
				}
			
				printf("\n Bitte geben Sie die Nullstelle x_quer ein: ");
				scanf("%lf",&x_quer);
				rewind(stdin);
		
				pn=Horner_Abdiv(n,c,x_quer,d);				// Unterprg-Aufruf

				printf("\n\n Y-Wert fuer die Nullstelle x_quer (pn):     pn = %lf",pn);
				printf("\n Das dazugehoerige Deflationspolynom lautet: pn-1(x)= ");
			
				for (I=n;I>0;I--)
				{
					if(I==n);
					else if(I>0) printf(" + ");
					printf("%0.2lfx^%i",d[I],I-1);
				}
				printf("\n\n ");
				system("pause");
				break;
			}


			case '3':										// Polynom-Nullstellenbestimmung mit Newtonverfahren
			{
				double eps_x, eps_f, pn, x_nst;
				int ITMAX;									// Maximale Anzahl der Iterationen

				system("cls");
				printf("\n\n <3> BERECHNUNG DER POLYNOM-NULLSTELLEN MIT DEM NEWTONVERFAHREN");
				printf("\n --------------------------------------------------------------");
				printf("\n\n Geben Sie nun bitte den Grad n des Polynoms ein: ");
				scanf("%i",&n);
				rewind(stdin);
			
				c=(double*)malloc(sizeof (double)*(n+1));
			
				for (I=n;I>=0;I--)							// Rueckwaertsrekursion
				{
					printf("\n Geben Sie nun den %i. Polynomkoeffizienten ein: ",n-I+1);
					scanf("%lf",&c[I]);
					rewind(stdin);
				}
				printf("\n Geben Sie nun den Startwert x_quer ein: ");
				scanf("%lf",&x_quer);
				rewind(stdin);

				printf("\n Bitte geben Sie epsylon_x ein: ");
				scanf("%lf",&eps_x);
				rewind(stdin);

				printf("\n Bitte geben Sie epsylon_f ein: ");
				scanf("%lf",&eps_f);
				rewind(stdin);

				printf("\n Bitte geben Sie die Anzahl der Iterationen ein: ");
				scanf("%d",&ITMAX);
				rewind(stdin);

				pn=Newton_Poly_Nst(n, c, x_quer, eps_x, eps_f, ITMAX, &x_nst);

				if (pn==1) printf("\n Die Nullstelle konnte mit ITMAX Iterationen nicht berechnet werden.");
				else printf("\n\n Berechnete Nullstelle: %1.20lf\n", x_nst);
				printf("\n\n ");
				system("pause");
				break;
			}


			case '4':									// Berechnung und Ausdruck saemtlicher reeler Nullstellen
			{
				int ITMAX=50,N=7;						// ITMAX - Anzahl der Iterationen,  N -> Polynomgrad
				double eps_x,eps_f,x_nst;				// eps_x u. eps_f -> Fehler   x_nst -> Nullstelle
				double c[8];
				double d[8];
				FILE * Pointer1;						// 1.Dateipointer fuer Ausdruck in Datei

				eps_x=eps_f=0.00000000000001;			// Epsylon x und f auf 10^-14 setzen
				c[7]=1;		
				c[6]=-28;				
				c[5]=322;
				c[4]=-1960;
				c[3]=6769;
				c[2]=-13132;
				c[1]=13068;
				c[0]=-5040;

				Pointer1 = fopen("Nullstellen1.sav","w");	// Definition der Datei
				
				system("cls");
				printf("\n\n <4> BERECHNUNG UND AUSDRUCK SAEMTLICHER REELLER NULLSTELLEN");
				printf("\n -----------------------------------------------------------");
				printf("\n\n Folgendes Polynom ist gegeben:\n p7(x) = x^7 - 28x^6 + 322x^5 - 1960x^4");
				printf("  + 6769x^3 - 13132x^2 + 13068x - 5040\n\n\n");

				Newton_Poly_Nst(N, c, 0, eps_x, eps_f, ITMAX, &x_nst);	// Berechnung der ersten Nullstelle

				printf(" Berechnet wurden folgende Nullstellen:\n\n");
				printf(" Nullstelle Original Polynom  X0:  %#.14f\n",x_nst);
				fprintf(Pointer1, "Mathematik II Praktikum Gruppe MAMF - Marcel Cuvelier\n\n\n");
				fprintf(Pointer1, "Berechnung saemtlicher reeller Nullstellen:\n\n");
				fprintf(Pointer1, "Original Polynom  X0 = %#.14f\n",x_nst);

				for(I=1;I<N;I++)										// Weitere Nullstellen berechnen
				{
					Horner_Abdiv(N, c, x_nst, d);						// Deflationspolynom berechnen

					Newton_Poly_Nst(N, d, x_nst, eps_x, eps_f, ITMAX, &x_nst);	// Deflationspolynomnullst.
			
					printf(" Nullstelle Deflationspolynom X%i:  %#.14f\n",I,x_nst);			// Ausgabe der Nullstellen
					fprintf(Pointer1, "Deflationspolynom X%i = %#.14f\n",I,x_nst);			// Ausgabe der NSt in die Datei
					
					Newton_Poly_Nst(N, c, x_nst, eps_x, eps_f, ITMAX, &x_nst);	
					
					printf(" Nullstelle Original Polynom  X%i:  %#.14f\n",I,x_nst);			// Ausgabe der Orig. Polyn. Nullstellen
					fprintf(Pointer1, "Original Polynom  X%i = %#.14f\n",I,x_nst);			// Ausgabe der Orig. Polyn. NSt in die Datei

				}
				printf("\n ");
				fprintf(Pointer1, "\nDatei 'Nullstellen1.sav'");
				fclose(Pointer1);										// Datei wird geschlossen
				printf("\n Die berechneten Werte wurden in die Datei 'Nullstellen1.sav' gespeichert.\n ");
				system("pause");
				break;
			}


			case '5':										// Berechnung reelle Nullstellen von p7*(x)
			{
				int ITMAX=50,N=7;							// ITMAX -> Anzahl der Iterationen
				double eps_x,eps_f,x_nst;					// Epsylon x und f auf 10^-14 setzen
				double c[8];
				double d[8];
				FILE * Pointer2;							// 2.Dateipointer fuer Ausdruck in Datei

				eps_x=eps_f=0.00000000000001;	
				c[7]=1;
				c[6]=-28;
				c[5]=322;						
				c[4]=-1960;
				c[3]=6769;
				c[2]=-13133;
				c[1]=13068;
				c[0]=-5040;

				Pointer2 = fopen("Nullstellen2.sav","w");		// Definition der 2.Datei

				system("cls");
				printf("\n\n <5> BERECHNUNG UND AUSDRUCK DER NULLSTELLEN VON P7*(X)");
				printf("\n ------------------------------------------------------");	

				Newton_Poly_Nst(N, c, 0., eps_x, eps_f, ITMAX, &x_nst);

				printf("\n\n Nullstelle Original Polynom  X0: %#.14f\n",x_nst);
				fprintf(Pointer2, "Mathematik II Praktikum Gruppe MAMF - Marcel Cuvelier\n\n\n");
				fprintf(Pointer2, "Berechnung der Nullstellen von P7*(X):\n\n");
				fprintf(Pointer2, "Nullstelle Original Polynom  X0 = %#.14f\n",x_nst);

				for(I=1;I<5;I++)	
				{

					Horner_Abdiv(N, c, x_nst, d);

					Newton_Poly_Nst(N, d, x_nst, eps_x, eps_f, ITMAX, &x_nst);
			
					printf(" Nullstelle Deflationspolynom X%i: %#.14lf\n",I,x_nst);
					fprintf(Pointer2, "Nullstelle Deflationspolynom X%i = %#.14f\n",I,x_nst);

					Newton_Poly_Nst(N, c, x_nst, eps_x, eps_f, ITMAX, &x_nst);	
					
					printf(" Nullstelle Original Polynom  X%i: %#.14f\n",I,x_nst);				// Ausgabe der Orig. Polyn. Nullstellen
					fprintf(Pointer2, "Nullstelle Original Polynom  X%i = %#.14f\n",I,x_nst);	// Ausgabe der Orig. Polyn. NSt in die Datei
				} 
				printf("\n\n ");
				fprintf(Pointer2, "\nDatei 'Nullstellen2.sav'");
				fclose(Pointer2);																// Datei wird geschlossen
				printf("\n Die berechneten Werte wurden in die Datei 'Nullstellen2.sav' gespeichert.\n ");
				system("pause");
				break;
			}


			case '6':								// Berechnung der Nullstellen fuer ein beliebiges Polynom
			{
				double eps_x, eps_f, x_nst,*d;
				int ITMAX;							//Anzahl der Iterationen

				system("cls");
				printf("\n\n <6> BERECHNUNG DER NULLSTELLEN FUER EIN BELIEBIGES POLYNOM");
				printf("\n ----------------------------------------------------------");
			
				printf("\n\n Geben Sie bitte den Grad n des Polynoms ein: ");
				scanf("%i",&n);
				rewind(stdin);
			
				c=(double*)malloc(sizeof (double)*(n+1));
				d=(double*)malloc(sizeof (double)*(n+1));
			
				for (I=n;I>=0;I--)					// Rückwärtsrekursion
				{
					printf("\n Geben Sie bitte den %i. Polynomkoeffizienten ein: ",n-I+1);
					scanf("%lf",&c[I]);
					rewind(stdin);
				}
				
				printf("\n Geben Sie nun den Startwert x_quer ein: ");
				scanf("%lf",&x_quer);
				rewind(stdin);

				printf("\n Geben Sie nun epsylon_x ein: ");
				scanf("%lf",&eps_x);
				rewind(stdin);

				printf("\n Geben Sie nun epsylon_f ein: ");
				scanf("%lf",&eps_f);
				rewind(stdin);

				printf("\n Geben Sie die Anzahl der Iterationen ein: ");
				scanf("%d",&ITMAX);

				Newton_Poly_Nst(n, c, x_quer, eps_x, eps_f, ITMAX, &x_nst);

				system("cls");
				printf("\n\n <6> ANZEIGEN DER BERECHNETEN WERTE");
				printf("\n ----------------------------------");
				printf("\n\n Nullstelle X0: %lf\n",x_nst);

				for(I=1;I<n;I++)
				{
					Horner_Abdiv(n, c, x_nst, d);

					Newton_Poly_Nst(n, d, x_nst, eps_x, eps_f, ITMAX, &x_nst);
			
					printf("\n Nullstelle X%i: %lf\n",I,x_nst);
				}
				printf("\n\n ");
				system("pause");
				break;
			}


			case '9':							// BEENDEN
			{
				system("cls");
				printf("\n\n\n\n\n >>> Das Programm wird beendet\n >>> ");
				return 0;
			}
		}
	}while(auswahl!='9');
return 0;
}


//----------------------------------- Unterprogramme --------------------------------------

void Horner_W2(int n,double c[], double x_quer, double *pn, double *pn_strich)
{
	double *C;									// c[]		--->	Polynomkoeffizienten
	int I;										// x_quer	--->	Auswertestelle
												// pn(x) = c0 + c1x + c2x^2 + ... + cnx^n,    ck,xEps

	C=(double*)malloc(sizeof (double)*(n+1)*2);	// Dynam. Speicherallok. des doppelten Vektors C

	C[n]=c[n];									// Setzen des neu angelegten Vektors (1.Haelfte - Polynom)
	C[(n*2)+1]=c[n];							// Setzen des neu angelegten Vektors (2.Haelfte - Abl.d.Polnoms)

	for (I=n-1;I>=0;I--)						// Polynomberechnung durch Rueckwaertsrekursion u. Horner Schema
	{								
		C[I]=C[I+1]*x_quer+c[I];
	}

	for(I=n*2;I>(n-1);I--)						// 1. Ableitung des Polynoms
	{
		C[I]=C[I+1]*x_quer+C[I-n-1];
	}

	*pn=C[0];									// Polynom in Rueckgabeparameter schreiben
	*pn_strich=C[n+2];							// Ableitung des Polynoms in Rueckgabeparameter schreiben
}	

//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

double Horner_Abdiv(int n,double c[], double x_quer, double *d)
{
	double *C;		
	int I;		

	C=(double*)malloc(sizeof (double)*(n+1));	// Dynam. Speicherallokation

	C[n]=c[n];									// Setzen des neu angelegten Vektors 
	d[n]=c[n];

	for (I=n-1;I>=0;I--)						// Polynomberechnung durch Rueckwaertsrekursion u. Horner Schema
	{							
		C[I]=C[I+1]*x_quer+c[I];				// NUR fuer I>0!! Deflationspolynom besitzt selbe Nullstellen, 
		d[I]=C[I];								// wie das eigentliche Polynom ohne X0
	}							
	return C[0];								// Rueckgabe von C[0]
}

//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

int Newton_Poly_Nst(int n,double c[], double x_quer,double eps_x,double eps_f,int ITMAX,double *x_nst) 
{
	int I=0;
	double pn, pn_strich, x_hilf, x_hilf2;		// x_hilf  -> Hilfsvariable (vorherige Naeherung der Nullstelle)
												// x_hilf2 -> Hilfsvariable (neu berechnete Naeherung der Nullst.)
	x_hilf2=x_quer;								// Setzen des Startwertes durch x_quer

	do
	{	
		x_hilf=x_hilf2;							// vorherigen Wert auf neuen Wert setzen

		Horner_W2(n,c,x_hilf2,&pn,&pn_strich);	// Aufrufen des Horner-Schemas mit neuer Naeherung x_hilf2

		x_hilf2=x_hilf-(pn/pn_strich);			// Naeherung neu berechnen

	}while((I++<=ITMAX) && (fabs(x_hilf2-x_hilf)>eps_x*fabs(x_hilf2)) && (fabs(pn)>eps_f)); //Abbruchkriterium

	*x_nst=x_hilf2;								// Nullstelle bzw. Naeherung zurueckgeben

	if(I==ITMAX) return 1;						//Nullstelle konnte nicht berechnet werden
	else return 0;
}