/**************************************************** * Nth derivation of a polynomial fraction P(x)/Q(x) * * * * ------------------------------------------------- * * ReF.: "Mathématiques en Turbo-Pascal By M. Ducamp * * and A. Reverchon (vol 2), Eyrolles, Paris, 1988" * * [BIBLI 05]. * * ------------------------------------------------- * * SAMPLE RUNS: * * * * Nth DERIVATION OF A POLYNOMIAL FRACTION: * * * * Enter polynomial fraction F1(x): * * P(x) = 3x2 +2x +1 * * Q(x) = 4x2 +5x +56 * * * * Order: 1 * * * * * * +7 X2 +328 X +107 * * ----------------- * * +16 X4 +40 X3 +473 X2 +560 X +3136 * * * * * * Nth DERIVATION OF A POLYNOMIAL FRACTION: * * * * Enter polynomial fraction F1(x): * * P(x) = x -1 * * Q(x) = x +1 * * * * Order: 5 * * * * * * +240 * * ----------------- * * + X6 +6 X5 +15 X4 +20 X3 +15 X2 +6 X + 1 * * * * ------------------------------------------------- * * Functions used (of unit Polynoms): * * * * AddNumber(), EnterPolynom(), DivNumber(), * * DisplayPolynom(), MultNumber() and SetNumber(). * * * * C++ version by J-P Moreau. * * (www.jpmoreau.fr) * ***************************************************** Explanations: ------------ This program calculates the nth derivative of the polynomial fraction F1=P0/Q0. The first derivative is: F1' = (P0'Q0-P0Q0')/Q0^2 = P1/Q1 The second derivative is: F1" = (F1')' = (P1'Q1-P1Q1')/Q1^2 = P2/Q2 If F1(n) = Pn/Qn, then we can obtain F1(n+1) by the formula: F1(n+1) = (Pn'Qn-PnQn')/Qn^2 = Pn+1/Qn+1 Note: for high values of derivation order n, errors may occur: for instance, let us seek the nth derivative of x-2/x-3. Up to n=6, the results are correct. For n>=7, the results are not correct (the simplification is not finished due to numerical errors caused by successive calls of routine SimpPolFract(). --------------------------------------------------------------*/ #include #include #include #include "polynoms.h" #include "polfract.h" ar_fraction F1,F; int n; //*** 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; } bool DerivPolynom(ar_polynom *P,int n,ar_polynom *R) { int i,j; ar_number u; R->degree=P->degree-n; if (R->degree<0) { R->degree=0; SetNumber(&R->coeff[0], "0"); } else for (i=0; i<=R->degree; i++) { R->coeff[i]=P->coeff[i+n]; for (j=1; j<=n; j++) { //put i+j in u u.is_real=FALSE; u.p=(long) i+j; u.q=1; u.value=(double) i+j; if (!MultNumber(R->coeff[i], u, &R->coeff[i])) return FALSE; } } 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; } //Nth derivative of a polynomial fraction bool DerivPolFract(ar_fraction *F,int order,ar_fraction *G) { ar_polynom P,Q; int i; ar_number a,b; if (!SetNumber(&a,"1")) return FALSE; //a=1 if (!SetNumber(&b,"-1")) return FALSE; //b=-1 for (i=1; i<=order; i++) { if (!MultPolynom(&F->denom,&F->denom,&G->denom)) return FALSE; if (!DerivPolynom(&F->numer,1,&P)) return FALSE; if (!DerivPolynom(&F->denom,1,&Q)) return FALSE; if (!MultPolynom(&F->denom,&P,&P)) return FALSE; if (!MultPolynom(&F->numer,&Q,&Q)) return FALSE; if (!CombiPolynom(&P,&Q,a,b,&G->numer)) return FALSE;; G->denom=F->denom; if (!SimpPolFract(G,FALSE)) return FALSE; if (!MultPolynom(&G->denom,&F->denom,&G->denom)) return FALSE; if (!SimpPolFract(G,FALSE)) return FALSE; *F=*G; } return TRUE; } void main() { printf("\n Nth DERIVATION OF A POLYNOMIAL FRACTION:\n\n"); EnterPolFract(" Enter polynomial fraction F1(x):", &F1); printf("\n"); printf(" Order: "); scanf("%d", &n); if (DerivPolFract(&F1,n,&F)) //display result F DisplayPolFract(&F); else printf("\n Error in derivation of order %d.", n); printf("\n\n"); } //end of file derivfra.cpp