/**************************************************** * 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 * * * * * * 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 #include #include #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; jMX) { 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)