/**************************************************** * GCD and SCM of two polynomials * * * * ------------------------------------------------- * * Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp * * and A. Reverchon (vol 2), Eyrolles, Paris, 1988" * * [BIBLI 05]. * * ------------------------------------------------- * * SAMPLE RUN: * * * * GCD AND SCM OF TWO POLYNOMIALS: * * * * P(x) = x5 +5x4 +7x3 +5x2 +x -1 * * Q(x) = x4 +4x3 -7x +2 * * * * * * Greatest common divisor: * * * * + X2 +3X -1 * * * * * * Smallest common multiple: * * * * + X7 +6 X6 +10 X5 +2 X4 -8 X3 -10 X2 -3 X +2 * * * * * * ------------------------------------------------- * * Functions used (of module Polynoms): * * * * AddNumber(), EnterPolynom(), DivNumber(), * * DisplayPolynom(), MultNumber() and SetNumber(). * * * * C++ version by J-P Moreau. * * (To be linked with polynoms.cpp).* * (www.jpmoreau.fr) * ***************************************************** Explanations: The Greatest common divisor (GCD) and the Smallest Common multiple (SCM) of two polynomials are defined but for a coefficient. That is to say, if a polynomial R(x) is the GCD of P(x) and Q(x), any polynomial µ*R(x) is also a GCD. To find the GCD, after permutating P(x) and Q(x), if the degree of P is lower than the degree of Q, we first perform the Euclidian division P/Q: P=Q.H1+R1 (see program Divpol.pas). if R1<>0, we then divide Q by R1: Q=R1.H2+R2. if R2<>0, we continue until Rn=0. Rn-1 is then the GCD of P(x) and Q(x). To find the SCM of P(x) and Q(x), we use the relation: P(x).Q(x)=GCD(P,Q).SCM(P,Q). ---------------------------------------------------------------*/ #include #include #include #include "polynoms.h" ar_polynom *P, *Q, *R; // 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)); //verify that P and Q are not void if (P->degree==0 && P->coeff[0].value==0) return FALSE; if (Q->degree==0 && Q->coeff[0].value==0) return FALSE; 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; } 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; } // SCM of two polynomials bool SCMPolynom(ar_polynom *P,ar_polynom *Q,ar_polynom *R) { ar_polynom *u, *v; u = (ar_polynom *) calloc(1,sizeof(ar_polynom)); v = (ar_polynom *) calloc(1,sizeof(ar_polynom)); if (!GCDPolynom(P,Q,u)) return FALSE; if (!MultPolynom(P,Q,v)) return FALSE; if (!DivPolynom(v,u,R,v)) return FALSE; free(u); free(v); return TRUE; } void main() { P = (ar_polynom *) calloc(1,sizeof(ar_polynom)); Q = (ar_polynom *) calloc(1,sizeof(ar_polynom)); R = (ar_polynom *) calloc(1,sizeof(ar_polynom)); printf("\n GCD AND SCM OF TWO POLYNOMIALS:\n\n"); if (!EnterPolynom(" P(X) = ", P)) return; if (!EnterPolynom(" Q(X) = ", Q)) return; printf("\n\n"); if (GCDPolynom(P,Q,R)) { printf(" Greatest common divisor:\n"); DisplayPolynom(R); } else printf(" Error in looking for GCD."); printf("\n\n"); if (SCMPolynom(P,Q,R)) { printf(" Smallest common multiple:\n"); DisplayPolynom(R); } else printf(" Error in looking for SCM."); printf("\n\n"); free(P); free(Q); free(R); } // end of file gcdpol.cpp