/****************************************************
* This program finds the simple elements of a poly- *
* nomial fraction the denominator of which is given *
* as a product of irreducible polynomials.          *
* ------------------------------------------------- *
* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp *
* and A. Reverchon (vol 2), Eyrolles, Paris, 1988"  *
* [BIBLI 05].                                       *
* ------------------------------------------------- *
* SAMPLE RUNS:                                      *
*                                                   *
*                      2x9 +1                       *
* EXAMPLE:  F(x) = ---------------                  *
*                  x3 (x2 +x +1)^2                  *
*                                                   *
*                                                   *
* FIND THE SIMPLE ELEMENTS OF A POLYNOMIAL FRACTION:*
*                                                   *
* Enter numerator: 2x9 +1                           *
*                                                   *
* Enter factor 1                                    *
* Power (0 to exit): 3                              *
* X2 coefficient: 0                                 *
* X  coefficient: 1                                 *
* Constant      : 0                                 *
* Enter factor 2                                    *
* Power (0 to exit): 2                              *
* X2 coefficient: 1                                 *
* X  coefficient: 1                                 *
* Constant      : 1                                 *
* Enter factor 3                                    *
* Power (0 to exit): 0                              *
*                                                   *
* Integer part:                                     *
*                                                   *
*       2                                           *
* 2.00 X  - 4.00 X + 2.00                           *
*                                                   *
* Factor 1                                          *
* A(1, 1) =          1.0000000000                   *
* A(1, 2) =         -2.0000000000                   *
* A(1, 3) =          1.0000000000                   *
*                                                   *
* Factor 2                                          *
* A(2, 1) =         -3.0000000000                   *
* B(2, 1) =          3.0000000000                   *
* A(2, 2) =          3.0000000000                   *
* A(2, 2) =          0.0000000000                   *
*                                                   *
*                                                   *
* So the result is:                                 *
*                                                   *
*                        1   -2   1      -3 +3x     *
* F(x) = (2x2 -4x +2) + (- + -- + --) + (------ +   *
*                        x   x2   x3     x2+x+1     *
*       <integer part>   <first factor>  <second    *
*                                                   *
*   3 + 0x                                          *
* ----------)                                       *
* (x2+x+1)^2                                        *
* factor>                                           *
*                                                   *
*                       C++ version by J-P Moreau.  *
*                           (www.jpmoreau.fr)       *
*****************************************************
Explanations:
------------
1. Theory
   ------
Let us consider a polynomial fraction F(x)=P(x)/Q(x) the denominator
of which is given as a product of irreducible polynomials. Such a
polynomial has a degree=1, or =2 but with no real roots.

So the polynomial Q(x) must be given in the form:

   Q(x) = H1(x)^a1.H2(x)^a2...Hn(x)^an   (^ stands for power)

where the Hi(x) polynomials are irreducible.

We can then prove (not done here) that there exists a unique decomposition
of fraction F(x) in the form:

   F(x) = E(x) + B1(x) + B2(x) + ... + Bn(x)

where E(x) is a polynomial called "integer part" of F(x) and the Bi(x)
take the form:

           j=aj
   Bi(x) = Sum (Lij(x)/Hj(x))  
           j=1

where degree of Lij(x) polynomial strictly < degree of Hj(x).

So, for "1st type" elements, of degree=1, Lij(x) is a scalar.
For "2nd type" elements, Lij(x) has a degree <= 1 and has the form:

   aij + bij x.

Let us take an example:

             2x9 + 1
   F(x) = --------------
          x^3 (x2+x+1)^2

The irreducible polynomials are: H1(x)=x with a1=3
and H2(x)=x2+x+1 with a2=2.

So, according to the theory, F(x) can be written in a unique way as:
 
                 a11   a12   a13   a21 + b21x   a22 + b22x
   F(x) = E(x) + --- + --- + --- + ---------- + ----------             
                  x    x2    x3      x2+x+1     (x2+x+1)^2

                 [    B1(x)    ]   [        B2(x)         ]

Finding the simple elements of F(x) is nothing more than to find the
integer part E(x) and all the coefficients aij and bij corresponding
to the different factors.

2. The program
   -----------

The fraction F(x) is no longer given as two polynomials P(x)/Q(x) as
in the other programs concerning polynomial fractions, but as 3 parameters:

   - numer, numerator polynomial of F(x),
   - nbfactor, number of denominator factors,
   - factor, table of ten variables of C++ type ONEFACTOR.

The type ar_onefactor is defined as follows:   typedef struct {
                                                         double a,b,c;
                                                         int power;
			                                           } ar_onefactor;
where a, b, c are the real coefficients of Hi(x)=ax2+bx+c
(a=0 if degree=1) and power is ai.

The output of Function SimpleElements() is:

    - intpol, polynomial giving the "integer part" E(x),
    - U, a vector of size mlig (calculated by the program)
	  containing all the coefficients aij and bij, in the very order 
	  in which they appear in the above decomposition.

The result of the Function SimpleElements() is TRUE if the calculation
is successful, else FALSE.

The program first verifies if input data are valid, chiefly the
irreducibility of denominator polynomial factors.

If data are correct, the program identifies initial fraction to the
sum of simple elements that must be obtained, according to theory,
for several x values with an arbitrary step of 0.25. So we obtain a
linear system to be solved by a classical Gauss-Jordan method in
double precision.

As output, the program gives the integer part and, for each factor,
the corresponding coefficients a(i,j) and b(i,j).

For accuracy problems, the size MX of the linear system is limited
to 20. 
               (adapted from above reference by J-P Moreau, Paris).
-------------------------------------------------------------------*/
#include <stdio.h>
#include <malloc.h>
#include <math.h>

#include "polynoms.h"
#include "polfract.h"


#define  MAXC     40       //maximum number of lines of square matrix M
#define  MACH_EPS 1.2e-16  //machine dependant smallest real number

#define  MX       20       //maximum size of system to solve with acceptable accuracy
#define  STEP     0.25     //increment step for x


typedef     struct {
              double a,b,c;
              int power;
			} ar_onefactor;


typedef ar_onefactor TENFACTORS[11]; //)     
                                     //)
typedef double MATRIX[MAXC][MAXC];   //)  index 0 not used here
                                     //)
typedef double VECTOR[MAXC];         //)


//global variables
ar_polynom   numer,integer;
TENFACTORS   factor;
int          nbfactor;
VECTOR       U;
int          i,j,k;


//*** operations on polynomials (see polynoms.cpp) ***

// P(X) * Q(X) = R(X) - R can be either P or Q
bool MultPolynom(ar_polynom *P, ar_polynom *Q, ar_polynom *R) {
  int i,j, n;
  ar_number u;
  ar_polynom *nr;
  nr = (ar_polynom *) calloc(1,sizeof(ar_polynom));
  nr->degree=P->degree+Q->degree;
  if (nr->degree>AR_MAXPOL) return FALSE;  // R degree is too big
  for (n=0; n<=nr->degree; n++)  {
    if (!SetNumber(&nr->coeff[n],"0")) return FALSE;
    for (i=0; i<=P->degree; i++) {
      j=n-i;
      if (j>=0 && j<=Q->degree) {
        if (!MultNumber(P->coeff[i],Q->coeff[j],&u)) return FALSE;
        if (!AddNumber(nr->coeff[n],u,&nr->coeff[n])) return FALSE;
      }
    }
  }
  //copy nr in R
  R->degree=nr->degree;
  for (i=0; i<=nr->degree; i++)  R->coeff[i]=nr->coeff[i];
  free(nr);
  return TRUE;
}


//Euclidian division of two polynomials
bool DivPolynom(ar_polynom *P, ar_polynom *Q, ar_polynom *H, ar_polynom *R) {
  int i,j;
  ar_number u;
  ar_polynom *nh, *nr;
  nh = (ar_polynom *) calloc(1,sizeof(ar_polynom));
  nr = (ar_polynom *) calloc(1,sizeof(ar_polynom));
  //The Q polynomial must be <> zero
  if (Q->degree==0 && Q->coeff[0].value==0) return FALSE;;
  nr->degree=P->degree;
  for (i=0; i<=P->degree; i++)  nr->coeff[i]=P->coeff[i];
  nh->degree = P->degree - Q->degree;
  if (nh->degree<0) {
    nh->degree=0; if (!SetNumber(&nh->coeff[0],"0")) return FALSE;
  }
  else {
    for (i=nh->degree; i>=0; i--)  {
      if (!DivNumber(nr->coeff[nr->degree],Q->coeff[Q->degree],
                       &nh->coeff[i])) return FALSE;
	  for (j=i; j<=nr->degree; j++) {
        if (!MultNumber(nh->coeff[i],Q->coeff[j-i], &u)) return FALSE;
        u.p=-u.p; u.value=-u.value;
        if (!AddNumber(nr->coeff[j],u, &nr->coeff[j])) return FALSE;
      }
      if (nr->degree > 0)  nr->degree--;
    }
    while (fabs(nr->coeff[nr->degree].value) < AR_SMALL && 
              nr->degree>0)  nr->degree--;
  }
  H->degree=nh->degree;
  for (i=0; i<=nh->degree; i++)  H->coeff[i]=nh->coeff[i];
  R->degree=nr->degree;
  for (i=0; i<=nr->degree; i++)  R->coeff[i]=nr->coeff[i];
  free(nh); free(nr);
  return TRUE;
} //DivPolynom()  


/******************************************                                     
*  SOLVING A LINEAR MATRIX SYSTEM AX = B  *
* with  Gauss-Jordan  method  using full  *
* pivoting at each step. During the pro-  *
* cess, original AA matrix  and BB vector *
* are destroyed to spare storage location.*
* --------------------------------------- *
* INPUTS:    AA   MATRIX N*N              *                                     
*            BB   VECTOR N                *                                     
* --------------------------------------- *                                     
* OUTPUTS:   AA   INVERSE OF AA N*N       *                                     
*            DET  DETERMINANT OF AA       *                                     
*            BB   SOLUTION VECTOR N       *
* --------------------------------------- *
* Note: Index zero not used here.         *                                     
******************************************/
void MATINV(int N,MATRIX AA,VECTOR BB,double *DET)   { 	
      VECTOR PC,PL,CS;                                                  
      double PV,PAV,temp,TT;   
      int I,IK,J,JK,K;
                                                         
// Initializations :                                                             
      *DET= 1.0; 
      for (I=1; I<=N; I++)  {                                                                
        PC[I]= 0.0;                                                                
        PL[I]= 0.0;                                                                
        CS[I]= 0.0;                                                                
      }              
// Main loop                                                                         
      for (K=1; K<=N; K++) {                                                                  
// Searching greatest pivot :                                               
        PV=AA[K][K];                                                              
        IK=K;                                                                    
        JK=K;                                                                    
        PAV= fabs(PV);                                                            
        for (I=K; I<=N; I++)                                                                
          for (J=K; J<=N; J++)  {     
            temp = fabs(AA[I][J]);                                                        
            if (temp > PAV) {                                      
              PV=AA[I][J];                                                        
              PAV= fabs(PV);                                                      
              IK=I;                                                              
              JK=J;                                                              
            }                                                               
          }                                                                 
                                                                               
// Search terminated, the pivot is in location I=IK, J=JK.                     
// Memorizing pivot location:  		  

		PC[K]=JK;                                                                
        PL[K]=IK;                                                                
                                                                               
// DETERMINANT DET is actualised                                              
// If DET=0, ERROR MESSAGE and STOP                                           
// Machine dependant EPSMACH equals here 1e-20
                                                                               
        if (IK!=K) *DET=-*DET;                                                   
        if (JK!=K) *DET=-*DET;                                                   
        *DET=*DET*PV;  
        temp= fabs(*DET);                                                            
        if (temp < MACH_EPS) {                                          
          printf("\n The determinant equals ZERO !!!\n");                                                              
          return;                                                                  
        }                                                                   
                                                                               
// POSITIONNING PIVOT IN K,K:		
        if(IK!=K)                                                         
          for (I=1; I<=N; I++) {                                                              
// EXCHANGE LINES IK and K of matrix AA:			
            TT=AA[IK][I];                                                         
            AA[IK][I]=AA[K][I];                                                    
            AA[K][I]=TT;                                                          
          }                                                                 
          TT=BB[IK];                                                           
          BB[IK]=BB[K];                                                      
          BB[K]=TT;                                                            
                                                                           
// PIVOT is at correct line		  
        if(JK!=K)                                                         
          for (I=1; I<=N; I++) {                                                              
// EXCHANGE COLUMNS JK and K of matrix AA:
            TT=AA[I][JK];                                                         
            AA[I][JK]=AA[I][K];                                                    
            AA[I][K]=TT;                                                          
          }                                                                 
                                                                           
// The PIVOT is at correct column.                                              
// and is located in K,K.                                                   
                                                                               
// Column K of matrix AA is stored in CS vector                             
// then column K is set to zero. 		  
        for (I=1; I<=N; I++) {                                                                
          CS[I]=AA[I][K];                                                         
          AA[I][K]= 0.0;                                                          
        }                                                                   
                                                                               
        CS[K]= 0.0;                                                                
        AA[K][K]= 1.0;                                                              
// Line K of matrix AA is modified:		
        temp= fabs(PV);                                          
        if(temp < MACH_EPS) {                                            
          printf("\n PIVOT IS TOO SMALL - STOP\n");                               
          return;                                                                  
        }                                                                   
        for (I=1; I<=N; I++)                                                                
          AA[K][I]=AA[K][I]/PV;                                                    
                                                                           
        BB[K]=BB[K]/PV;                                                  
                                                                           
// Other lines of matrix AA are modified:		  
        for (J=1; J<=N; J++) {                                                                
          if (J==K) J++;                                                  
          for (I=1; I<=N; I++)                                                              
// Line J of matrix AA is modified:
            AA[J][I]=AA[J][I]-CS[J]*AA[K][I];
          BB[J]=BB[J]-CS[J]*BB[K];                                     
        }                                                                   
// Line K is ready
      } //End of K loop

// MATRIX AA INVERSION IS DONE - REARRANGEMENT OF MATRIX AA
                                                                               
// EXCHANGE LINES                                                                        
      for (I=N; I>0; I--) {                                                               
        IK=(int) PC[I];                                                                
        if (IK==I) goto fin1;                                                   
// EXCHANGE LINES I and PC(I) of matrix AA:
        for (J=1; J<=N; J++) {                                                                
          TT=AA[I][J];                                                            
          AA[I][J]=AA[IK][J];                                                      
          AA[IK][J]=TT;                                                           
        }
        TT=BB[I];                                                          
        BB[I]=BB[IK];                                                    
        BB[IK]=TT;                                                         
// NO MORE EXCHANGE is NECESSARY                                                      
// Go to next line                                                 
        fin1: ;     			  
      } // for i                                                                     
                                                                               
// EXCHANGE COLUMNS:  	  
      for (J=N; J>0; J--) {                                                                         
        JK=(int) PL[J];                                                                
        if (JK==J) goto fin2;                                                   
// EXCHANGE COLUMNS J ET PL(J) of matrix AA:
        for (I=1; I<=N; I++) {                                                                
          TT=AA[I][J];                                                            
          AA[I][J]=AA[I][JK];                                                      
          AA[I][JK]=TT;                                                           
        } 
        fin2: ;                                                                  
// NO MORE EXCHANGE is NECESSARY                                                      
// Go to next column.
      }                                                                     
// REARRANGEMENT TERMINATED.                                                        
      return;  
  } //Matinv()       


//find the simple elements - See explanations above
bool SimpleElements(ar_polynom *numer,int nbfactor,TENFACTORS factor,
					ar_polynom *intpol, VECTOR U)   {

  int    i,j,k,h;
  double deter,x,s,t,w;
  double ta,tb,tc;
  ar_number y,z;
  bool   ok;
  MATRIX M;
  int    mcol,mlig,ulig,vlig;
  ar_polynom  inter,denom,R;
  ulig=0;
  if (nbfactor<1 || nbfactor>10) return FALSE;
  //verify that each factor is not reducible
  for (i=1; i<=nbfactor; i++) {
	ta=factor[i].a; tb=factor[i].b; tc=factor[i].c;
    if (ta!=0 && ta*ta >= 4*tb*tc) {
      printf("\n One of the factors is reducible!\n");
      return FALSE;
    }
  }
  //verify that each factor is different from others
  for (i=1; i<=nbfactor; i++)
    for (j=1; j<i; j++) {
      if (factor[i].a==factor[j].a && factor[i].b==factor[j].b
		    && factor[i].c==factor[j].c)  {
        printf("\n At least two factors are identical!\n");
        return FALSE;
      }
    }
  //verify that system size < MX
  ulig=0;
  for (i=1; i<=nbfactor; i++) {
    if (factor[i].a!=0)
      ulig=ulig+2*factor[i].power;
    else
      ulig=ulig+factor[i].power;
  }
  if (ulig>MX) {
    printf("\n Denominator too big!\n");
    return FALSE;
  }
  //verify that fraction is not reducible
  for (i=1; i<=nbfactor; i++) {
    if (factor[i].a!=0)  denom.degree=2; else denom.degree=1;
    denom.coeff[0].is_real=TRUE; denom.coeff[0].value=factor[i].c;
    denom.coeff[1].is_real=TRUE; denom.coeff[1].value=factor[i].b;
    denom.coeff[2].is_real=TRUE; denom.coeff[2].value=factor[i].a;
    if (!DivPolynom(numer,&denom,&inter,&R)) return FALSE;
    if (R.degree==0 && fabs(R.coeff[0].value)<AR_SMALL) {
      printf("\n Fraction is reducible!\n");
      return FALSE;
    }
  }
  //calculate denominator polynomial
  denom.degree=0;
  denom.coeff[0].is_real=TRUE; denom.coeff[0].value=1;
  for (i=1; i<=nbfactor; i++) {
    ta=factor[i].a; tb=factor[i].b; tc=factor[i].c;
    if (ta!=0) inter.degree=2; else inter.degree=1;
    inter.coeff[0].is_real=TRUE; inter.coeff[0].value=tc;
    inter.coeff[1].is_real=TRUE; inter.coeff[1].value=tb;
    inter.coeff[2].is_real=TRUE; inter.coeff[2].value=ta;
    for (j=1; j<=factor[i].power; j++)
      if (!MultPolynom(&denom,&inter,&denom)) return FALSE;
  }

  //calculate integer part
  if (!DivPolynom(numer,&denom,intpol,&R)) return FALSE;

  //create linear system MU = V to be solved
  mlig=ulig; mcol=mlig; vlig=ulig;
  x=-vlig*STEP/2 -STEP;
  for (i=1; i<=vlig; i++) {
    //eliminate an x for which denominator =0
    j=0;
    do {
      x+=STEP; ok=TRUE;
      for (k=1; k<=nbfactor; k++) {
        ta=factor[k].a; tb=factor[k].b; tc=factor[k].c;
        if (fabs(tc+x*(tb+x*ta)) < AR_SMALL)  ok=FALSE;
	  }
    } while (!ok);
    //fill first line of matrix M
    for (k=1; k<=nbfactor; k++) {
      ta=factor[k].a; tb=factor[k].b; tc=factor[k].c;
      y.is_real=TRUE; y.value=x;
      if (!EvaluatePolynom(&denom,&y,&z)) return FALSE;
      w=z.value; t=tc+x*(tb+ta*x); s=1;
      for (h=1; h<=factor[k].power; h++) {
        s*=t; j++; M[i][j]=w/s;
        if (ta!=0) {
          j++; M[i][j]=w*x/s;
        }
      }
    }
    y.is_real=TRUE; y.value=x;
    if (!EvaluatePolynom(&R,&y,&z)) return FALSE;
    U[i]=z.value;
  } //i loop
  //solve linear system MU=V by Gauss-Jordan method
  MATINV(mlig,M,U,&deter);
  //Process succesful 
  return TRUE;
}


// demonstration program of function SimpleElements()
void main()  {
  
  printf("\n FIND THE SIMPLE ELEMENTS OF A POLYNOMIAL FRACTION:\n\n");
  EnterPolynom(" Enter numerator: ",&numer);
  printf("\n"); i=0;
  do {
    i++;
    printf(" Enter factor %d\n", i);
    printf(" Power (0 to exit): "); scanf("%d",&factor[i].power);
    if (factor[i].power!=0) {
      printf(" X2 coefficient: "); scanf("%lf",&factor[i].a);
      printf(" X  coefficient: "); scanf("%lf",&factor[i].b);
      printf(" Constant      : "); scanf("%lf",&factor[i].c);
    }
  } while(i!=0 && factor[i].power!=0);
  nbfactor=i-1; printf("\n");

  if (SimpleElements(&numer,nbfactor,factor,&integer,U)) {
    printf("\n Integer part:\n");
    DisplayPolynom(&integer); j=0;
    printf("\n");
    for (i=1; i<=nbfactor; i++) {
      printf("\n Factor %d:\n", i);
      for (k=1; k<=factor[i].power; k++) {
        j++;
        printf(" A(%d,%d) = %22.10f\n", i,k, U[j]);
        if (factor[i].a!=0) {
          j++;
          printf(" B(%d,%d) = %22.10f\n", i,k, U[j]);
        }
      }
	}
  }
  else
    printf(" Error in finding simple elements.\n");
  printf("\n\n");
}

//end of file simpelem.cpp