/**************************************************** * Substitution of two polynomial fractions * * F(x) = F1(F2(x)) * * * * ------------------------------------------------- * * Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp * * and A. Reverchon (vol 2), Eyrolles, Paris, 1988" * * [BIBLI 05]. * * ------------------------------------------------- * * SAMPLE RUNS: * * * * SUBSTITUTION OF TWO POLYNOMIAL FRACTIONS * * F(x) = F1(F2(x)) * * * * Enter polynomial fraction F1(x): * * P(x) = x -2 * * Q(x) = 4x +3 * * * * Enter polynomial fraction F2(x): * * P(x) = x -1 * * Q(x) = x +1 * * * * * * - X -3 * * ------------- * * +7 X -1 * * * * * * SUBSTITUTION OF TWO POLYNOMIAL FRACTIONS * * F(x) = F1(F2(x)) * * * * Enter polynomial fraction F1(x): * * P(x) = x3 * * Q(x) = 1 * * * * Enter polynomial fraction F2(x): * * P(x) = 1/2x * * Q(x) = x2 -5 * * * * * * +1/8 X3 * * ------------- * * + X6 -15 X4 +75 X2 -125 * * * * ------------------------------------------------- * * Functions used (of module Polynoms): * * * * AddNumber(), EnterPolynom(), DivNumber(), * * DisplayPolynom(), MultNumber() and SetNumber(). * * * * C++ version by J-P Moreau. * * (www.jpmoreau.fr) * ****************************************************/ #include #include #include #include "polynoms.h" #include "polfract.h" ar_fraction F,F1,F2; ar_number a,b; //*** operations on polynomials *** // 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() //copy polynomial q in polynomial p void polcpy(ar_polynom *p, ar_polynom *q) { int i; p->degree=q->degree; for (i=0; i<=q->degree; i++) p->coeff[i]=q->coeff[i]; } // GCD of two polynomials bool GCDPolynom(ar_polynom *P,ar_polynom *Q, ar_polynom *R) { ar_number z; long pg,pp; int i; ar_polynom *v, *np, *nq; //dynamic memory allocation of polynomials v = (ar_polynom *) calloc(1,sizeof(ar_polynom)); np = (ar_polynom *) calloc(1,sizeof(ar_polynom)); nq = (ar_polynom *) calloc(1,sizeof(ar_polynom)); if (P->degree==0 && P->coeff[P->degree].value==0) return FALSE; if (Q->degree==0 && Q->coeff[Q->degree].value==0) return FALSE; polcpy(np,P); polcpy(nq,Q); polcpy(R,nq); while(R->degree>0) { if (!DivPolynom(np,nq,v,R)) return FALSE; while(R->degree>0 && fabs(R->coeff[R->degree].value)degree--; if (R->degree==0) { if (fabs(R->coeff[0].value)>AR_SMALL) nq->degree=0; } else { polcpy(np,nq); polcpy(nq,R); } for (i=0; idegree; i++) if (!DivNumber(nq->coeff[i],nq->coeff[nq->degree], &nq->coeff[i])) return FALSE; SetNumber(&nq->coeff[nq->degree], "1"); } polcpy(R,nq); if (nq->degree==0) SetNumber(&R->coeff[0], "1"); else { for (pg=0, i=0; i<=R->degree; i++) if (!R->coeff[i].is_real) pg = (pg) ? GCD(pg,R->coeff[i].p) : R->coeff[i].p; for (pp=0, i=0; i<=R->degree; i++) if (!R->coeff[i].is_real) pp = (pp==0) ? R->coeff[i].q : pp * R->coeff[i].q / GCD(pp,R->coeff[i].q); if (pg) { z.is_real = FALSE; z.p=pp; z.q=pg; z.value=(double) pp / (double) pg; for (i=0; i<=R->degree; i++) if (!MultNumber(R->coeff[i], z, &R->coeff[i])) return FALSE; } } free(v); free(np); free(nq); return TRUE; } // a P(X) + b Q(X) = R(X) bool CombiPolynom(ar_polynom *P, ar_polynom *Q, ar_number a, ar_number b, ar_polynom *R) { int i; ar_number u,v; if (Q->degree > P->degree) R->degree=Q->degree; else R->degree=P->degree; if (R->degree > AR_MAXPOL) return FALSE; // degree of R is too big for (i=0; i<=R->degree; i++) { if (!SetNumber(&u,"0")) return FALSE; if (!SetNumber(&v,"0")) return FALSE; if (i<=P->degree) if (!MultNumber(a,P->coeff[i], &u)) return FALSE; if (i<=Q->degree) if (!MultNumber(b,Q->coeff[i], &v)) return FALSE; if (!AddNumber(u,v,&R->coeff[i])) return FALSE; } while (R->degree>0 && fabs(R->coeff[R->degree].value)degree--; return TRUE; } //*** operations on polynomial fractions *** //Simplify a polynomial fraction F(x) = P(x) / Q(x) bool SimpPolFract(ar_fraction *F,bool integer) { ar_polynom P,R; ar_number z; long pg,pp1,pp2; int i; if (!GCDPolynom(&F->numer,&F->denom,&P)) return FALSE; if (!DivPolynom(&F->numer,&P,&F->numer,&R)) return FALSE; if (!DivPolynom(&F->denom,&P,&F->denom,&R)) return FALSE; if (integer) { //integer coefficients asked pg=0; for (i=0; i<=F->numer.degree; i++) if (!F->numer.coeff[i].is_real) if (pg==0) pg=F->numer.coeff[i].p; else pg=GCD(pg,F->numer.coeff[i].p); for (i=0; i<=F->denom.degree; i++) if (!F->denom.coeff[i].is_real) if (pg==0) pg=F->denom.coeff[i].p; else pg=GCD(pg,F->denom.coeff[i].p); pp1=0; for (i=0; i<=F->numer.degree; i++) if (!F->numer.coeff[i].is_real) if (pp1==0) pp1=F->numer.coeff[i].q; else pp1=pp1*F->numer.coeff[i].q/GCD(pp1,F->numer.coeff[i].q); pp2=0; for (i=0; i<=F->denom.degree; i++) if (!F->denom.coeff[i].is_real) if (pp2==0) pp2=F->denom.coeff[i].q; else pp2=pp2*F->denom.coeff[i].q/GCD(pp2,F->denom.coeff[i].q); z.p=int(pp1*pp2/GCD(pp1,pp2)); z.q=pg; if (z.q!=0) { z.value=z.p/z.q; z.is_real=FALSE; for (i=0; i<=F->numer.degree; i++) if (!MultNumber(F->numer.coeff[i],z,&F->numer.coeff[i])) return FALSE; for (i=0; i<=F->denom.degree; i++) if (!MultNumber(F->denom.coeff[i],z,&F->denom.coeff[i])) return FALSE; } } return TRUE; } //Substitution of two polynomial fractions F=F1(F2(x)) bool SubstPolFract(ar_fraction *F1,ar_fraction *F2,ar_fraction *F) { int i; ar_number a,b; ar_polynom P,Q; P.degree=0; for (i=0; inumer.degree; i>=0; i--) { if (!MultPolynom(&P,&F2->numer,&P)) return FALSE; b=F1->numer.coeff[i]; if (!CombiPolynom(&P,&Q,a,b,&P)) return FALSE; if (!MultPolynom(&Q,&F2->denom,&Q)) return FALSE; } F->numer=P; P.degree=0; for (i=0; idenom.degree; i>=0; i--) { if (!MultPolynom(&P,&F2->numer,&P)) return FALSE; b=F1->denom.coeff[i]; if (!CombiPolynom(&P,&Q,a,b,&P)) return FALSE; if (!MultPolynom(&Q,&F2->denom,&Q)) return FALSE; } F->denom=P; for (i=1; i<=(F1->denom.degree-F1->numer.degree); i++) if (!MultPolynom(&F->numer,&F2->denom,&F->numer)) return FALSE; for (i=1; i<= -(F1->denom.degree-F1->numer.degree); i++) if (!MultPolynom(&F->denom,&F2->denom,&F->denom)) return FALSE; if (!SimpPolFract(F,FALSE)) return FALSE; return TRUE; } //*** demonstration program *** void main() { printf("\n SUBSTITUTION OF TWO POLYNOMIAL FRACTIONS\n"); printf(" F = F1(F2(x))\n\n"); EnterPolFract(" Enter polynomial fraction F1(x):", &F1); printf("\n"); EnterPolFract(" Enter polynomial fraction F2(x):", &F2); printf("\n"); if (SubstPolFract(&F1,&F2,&F)) //display result F DisplayPolFract(&F); else //print error message printf(" Error in substitution.\n"); printf("\n\n"); } //end of file substfra.cpp