#define  LEGENDRE     0
#define  JACOBI       1
#define  CHEBYSHEV    2
#define  HERMITE      3 
#define  ALL          4


#define  COEFFICIENTS 0
#define  VALUE        1
#define  ZEROES       2
#define  SAME_DEGREE  3
#define  SAME_POL     4

#define  DEFAULT  0
#define  SHIFTED  1

#include "Polinom.h"
/**
* OrthogonalPolinom class encapsulates methods for computing the
* coefficients of Legendre's, Chebyshev's, Hermite's and Jacobi's
* orthogonal Polinoms on [-1,1] ( (-oo, +oo) for Hermite's) interval or
* on the shifted interval [0,1] using reccurence formulas. It also contains
* methods for calculating the value of the chosen Polinom at an isolated point in the given interval
* and methods for finding zeroes of the chosen Polinom ( using Bairstow method).
*/



class OrthogonalPolinom 
{
public:
 int PolinomType ;     /* type of the Polinom - Legendre's, Hermite's, Chebyshev's or the Jacobi's*/
 int degree ;     /* degree of the chosen Polinom*/
 int problemType ;     /* type of the problem to be solved */
 int interval ;     /* the interval that is chosen [-1,1] or the [0,1]*/

 double alpha;
 double beta;

Polinom *legendre; //[MAXORDER]
Polinom *hermite;
Polinom *chebyshev;
Polinom *jacobi;
Polinom tekuci; //za cuvanje nula i tekucih koeficijenata

/* for the computation of the value of the polinomial at the isolated point */
double isolatedPoint;
double PolinomValue;

OrthogonalPolinom(int polType,int deg,int interv,int prType,double isPoint=0.0,double al=0.0,double bt=0.0)
{
	
	PolinomType = polType;
	degree = deg;
	alpha = al;
	beta = bt;
	interval = interv;
	problemType = prType;
	isolatedPoint = isPoint;

        /* direct calculation of the results */

        doMethod();
}





void doMethod()
  {
     int i;
	
	if ((interval!=DEFAULT) && (interval != SHIFTED))
		throw " neispravan interval!";

	/* if the type of the Polinom is Hermite's the only appropriate interval is [-oo,+oo], not the shifted one*/
	if ((interval == SHIFTED) && (PolinomType == HERMITE))
		throw " ne postoji Hermite-ov polinom na [0,1] intervalu !!! ";

	/* degree of the Polinom has to be larger than 0 */
	if (degree < 0)
		throw " stepen manji od 0 ! ";


	/* the value of the Polinom can be calculated only if the point of calculation is in the selected interval*/
	if (PolinomType!= HERMITE) 
	{
		switch (interval) 
		{
			case DEFAULT:
				if ((isolatedPoint <-1.0) || (isolatedPoint>1.0))
					throw " Izolovana tacka je van intervala !";
				break;
			case SHIFTED:
				if ((isolatedPoint < 0.0) || (isolatedPoint>1.0))
					throw " Izolovana tacka je van intervala !";	
				break;
		}
	}


/** Legendre's Polinoms are orthogonal on the interval [-1,1] with weight function p(x)=1
* the reccurence formula for their calculation is:
* legendre[0]=1, legendre[1]=x,
* legendre[i] = ((2*i-1)*x*legendre[i-1]-(i-1)*legendre[i-2])/i
* Legendre's Polinoms on the shifted [0,1] interval are also orthogonal Polinoms with 
* the weigth function p(x)=1.The reccurence formula for their calculation is: 
* legendre[0]=1, legendre[1]=2*x-1
* legendre[i]=((2*i-1)*(2*x-1)*x*legendre[i-1]-(i-1)*legendre[i-2])/i. Since legendre[1]=2*x-1 on the shifted interval,
* reccurence formula for both intervals is:
* legendre[i] = ((2*i-1)*x*legendre[1]*legendre[i-1]-(i-1)*legendre[i-2])/i
* @returns an array of legendre Polinoms*/
	//delete[] legendre;
	legendre = new Polinom[degree+1];

	/*initialize legendre[0] and legendre[1] */
    legendre[0] = Polinom();
	legendre[0].coefs[0]=1;
        if (degree>0) 
		{
                legendre[1]= Polinom(1);
	
                if (interval== DEFAULT)
                        /* legendre[1]=x*/
                        legendre[1].coefs[1]=1;
                else {  /*legendre[1]= 2*x-1*/
                        legendre[1].coefs[0]=-1;
                        legendre[1].coefs[1]=2;
                      }
         }

	/*calculates the coefficients using the previously discussed reccurence formulas */
	for(i=2; i<= degree; i++)
	{
		legendre[i] = Polinom(i);
		legendre[i] = legendre[i] + (legendre[i-1]*(legendre[1]*(2*i-1)));
		legendre[i] = legendre[i] + (legendre[i-2] * (1-i));
		legendre[i] = legendre[i] * (1.0/i);
	}

/**
* Hermite's Polinoms are orthogonal on the interval (-oo,+oo) with the weight function
* p(x)=e^(-x^2). This method calculates the coefficients of the Hermite's Polinoms using the following reccurence formulas:
* hermite[0]=1, hermite[1]= 2*x
* hermite[i]= 2*x*hermite[i-1] - 2*(i-1)*hermite[i-2]= hermite[1]*hermite[i-1]-2*(i-1)*hermite[i-2]
* the method returns an array of Polinoms:
* hermite[i] represents Hermite's Polinom of the degree i
* @returns an array of hermite Polinoms
*/
 
	//delete[] hermite;
	hermite = new Polinom[degree+1]; 
	hermite[0] = Polinom();
	hermite[0].coefs[0]=1;

        if (degree > 0 )
		      {
                /* initialize hermite[1]=2*x */
                hermite[1] = Polinom(1); /*this constructor sets all coefficients of the Polinom to zero*/
                hermite[1].coefs[1]= 2;
               }
	/*calculating hermites Polinoms using the previously explained reccurence formulas*/
	for (i=2; i<= degree; i++) 
	    {
		hermite[i] = Polinom(i); /* creates a new Polinom of the degree i*/
		hermite[i] = hermite[i] + (hermite[i-1] * (hermite[1]));
		hermite[i] = hermite[i] + (hermite[i-2] * (2*(1-i)));
		}
/** Chebyshev's Polinoms are orthogonal Polinoms, the weigh function beeing: p(x)=(1-x^2)^(-1/2)
* This method calculates the coefficients of the  Chebyshev's Polinoms using the following reccurence formulas:
* on the interval [-1,1]:
* chebyshev[0]=1, chebyshev[1]=x, chebyshev[i]=2*x*chebyshev[i-1]-chebyshev[i-2]
* on the interval [0,1]:
*chebyshev[0]=1, chebyshev[1]=2*x-1, chebyshev[i]=2*(2*x-1)*chebyshev[i-1]-chebyshev[i-2];
* SInce chebyshev[1]=x in the first case, and chebyshev[1]=2*x-1 in the second case
*the reccurence formula for both intervals is:
* chebyshev[i]=2*chebyshev[1]*chebyshev[i-1]-chebyshev[i-2]
* the method returns an array of Polinoms:
* chebyshev[i] represents Chebyshev's Polinom of the degree i
* @returns an array of chebyshev Polinoms
*/
     //delete[] chebyshev;
	 chebyshev = new Polinom[degree+1];

	/* initialize chebyshev[0] */
    chebyshev[0] = Polinom();
	chebyshev[0].coefs[0]=1;

        if (degree > 0 ) 
		{

                /* initialize chebyshev[1] */
                chebyshev[1] = Polinom(1);
	
                /* depending on the type of the interval chebyshev[1]=x ([-1,1]) or chebyshev[1]=2*x-1 ([0,1]) */
                if(interval== DEFAULT)
                        chebyshev[1].coefs[1]=1;
                else 
				{
                        chebyshev[1].coefs[0]=-1;
                        chebyshev[1].coefs[1]=2;
		        }
        }
	/* compute the rest of the Polinoms using the previously explained reccurence formulas*/
	
	for(i=2; i<=degree; i++)
	{
		chebyshev[i] = Polinom(i); /* creates and initializes new array element */
		chebyshev[i] = chebyshev[i] + (chebyshev[i-1] * (chebyshev[1] * (2)));
		chebyshev[i] = chebyshev[i] + (chebyshev[i-2] * (-1.0));
	}
	
/** Jacobi Polinoms are orthogonal Polinoms on the [-1,1] interval 
* weight function beeing p(x)=(1-x)^alpha*(1+x)^beta 
* this method calculates the coefficients of the Jacobi's polynomilas using the following reccurence formulas:
* on the interval [-1,1]:
* jacobi[0]=1, jacobi[1]= (alpha-beta)/2 + ((alpha+beta +2)*x)/2;
* jacobi[i]=((a2i + a3i*x)*jacobi[i-1] - a4i*jacobi[i-2])/ a1i with
* a1i = 2*i*(i+alpha+beta)*(2*(i-1) + alpha + beta)
* a2i = (2*i - 1 + alpha + beta)*(alpha^2-beta^2)
* a3i = (2*i + alpha + beta)*(2*i + alpha + beta - 1)*(2*i + alpha + beta - 2)
* a4i = 2* ( i-1+alpha)*(i-1+beta)*(2*i + alpha+beta)
* for calculation of the Polinom's coefficients on the interval [0,1] use the following:
* jacobi[0]=1, jacobi[1]=-1*(beta+1) + (alpha+beta + 2)*x, the coefficients for the rest of the Polinoms can be
* calculated by substituting 2*x-1 instead of x.
* the method returns an array of Polinoms:
* jacobi[i] represents Jacobi's Polinom of the  ith degree
* @returns an array of jacobi Polinoms
*/

	//delete[] jacobi;
	jacobi = new Polinom[degree+1];/* allocates space in memory for an array of Polinoms*/
    
	//if ((alpha==0) && (beta==0))*jacobi = *legendre;
			//else
	//{	
	double a1, a2, a3, a4; /* coefficients needed for the calculation*/
	Polinom multiplyBy(1); /*multiplication factor x for [-1,1], 2*x-1 for the [0,1] interval*/


	/* initialization of the first two elements of the array jacobi[] */


	jacobi[0]= Polinom();
	jacobi[0].coefs[0]=1;

        if (degree > 0 )
		{
                jacobi[1] = Polinom(1);

                if (interval==DEFAULT) 
				{
                        multiplyBy.coefs[1] = 1; /* multiplyBy = x on the interval [-1,1] */
                        jacobi[1].coefs[0] = (alpha-beta)/2;
                        jacobi[1].coefs[1] = (alpha+beta+2)/2;
                }
                else {
                        multiplyBy.coefs[0] =- 1;
                        multiplyBy.coefs[1] = 2; /* multiplyBy = 2*x-1 on the interval [0,1] */
                        jacobi[1].coefs[0] = -1 * (1+beta);
                        jacobi[1].coefs[1] = alpha + beta +2;
                     }
         }

	/* computation of the rest of the Jacobi Polinoms from the 2nd degree to the selected degree*/
	

        for (i=2; i<=degree; i++) 
		{
 /* creates ith element of the array jacobi as a Polinom of the ith degree*/
                jacobi[i] = Polinom(i); 
	        
		/* compute the a1, a2 , a3, a4 coefficients of the reccurence formula */
		        a1 = 2*i*(i+alpha+beta)*(2*(i-1) + alpha + beta);
                a2 = (2*i - 1 + alpha + beta)*(alpha*alpha-beta*beta);
                a3 = (2*i + alpha + beta)*(2*i + alpha + beta - 1)*(2*i + alpha + beta - 2);
                a4 = 2* ( i-1+alpha)*(i-1+beta)*(2*i + alpha+beta);

		/* compute jacobi's Polinom of the ith degree using the reccurence formula explained above*/
		        jacobi[i] = jacobi[i] + (jacobi[i-1] * (a2));   /*jacobi[i]+=jacobi[i-1]*a2 */
		        jacobi[i] = jacobi[i] + (jacobi[i-1] * (multiplyBy * (a3)));   /*jacobi[i]+=jacobi[i-1]*multiplyBy*a3 */
                jacobi[i] = jacobi[i] + (jacobi[i-2] * (-1*a4));   /* jacobi[i]+=jacobi[i-2]*(-a4) */
		        jacobi[i] = jacobi[i] * (1.0/a1);   /* jacobi[i]/=a1 */
	    } 
	//}	


	switch (PolinomType) 
		{
		
			case LEGENDRE:
			    tekuci = legendre[degree];
				PolinomValue = tekuci(isolatedPoint);
				break;
			case HERMITE:
			    tekuci = hermite[degree];
				PolinomValue = tekuci(isolatedPoint);
				break;
			case CHEBYSHEV:
			    tekuci = chebyshev[degree];
				PolinomValue = tekuci(isolatedPoint);
				break;
			case JACOBI:
			    tekuci = jacobi[degree];
				PolinomValue = tekuci(isolatedPoint);
				break;
		} 
			

}
~OrthogonalPolinom(void)
{
delete legendre; //[MAXORDER]
delete hermite;
delete chebyshev;
delete jacobi;
}
};