/****************************************************
*        Multiply two polynomial fractions          *
*                                                   *
* ------------------------------------------------- *
* Ref.: "Mathématiques en Turbo-Pascal By M. Ducamp *
* and A. Reverchon (vol 2), Eyrolles, Paris, 1988"  *
* [BIBLI 05].                                       *
* ------------------------------------------------- *
* SAMPLE RUNS:                                      *
*                                                   *
* MULTIPLY TWO POLYNOMIAL FRACTIONS:                *
*                                                   *
* Enter polynomial fraction F1(x):                  *
*                                                   *
* P(x) = x -1                                       *
* Q(x) = x2 +5x -7                                  *
*                                                   *
* Enter polynomial fraction F2(x):                  *
*                                                   *
* P(x) = x +3                                       *
* Q(x) = x2 -1                                      *
*                                                   *
*                                                   *
* + X +3                                            *
* ------------                                      *
* + X3 +6 X2 -2 X -7                                *
*                                                   *
* ------------------------------------------------- *
* Functions used (of unit Polynoms):                *
*                                                   *
*  AddNumber(), EnterPolynom(), DivNumber(),        *
*  DisplayPolynom(), MultNumber() and SetNumber().  *
*                                                   *
*                       C++ version by J-P Moreau.  *
*                           (www.jpmoreau.fr)       *
****************************************************/
#include <stdio.h>
#include <malloc.h>
#include <math.h>

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


ar_fraction F,F1,F2;


// 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;
}


//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)<AR_SMALL)
	  R->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; i<nq->degree; 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;
}

//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;
}

//Multiply two polynomial fractions
bool MultPolFract(ar_fraction *F1,ar_fraction *F2,ar_fraction *F) {
  if (!SimpPolFract(F1,FALSE)) return FALSE;
  if (!SimpPolFract(F2,FALSE)) return FALSE;
  F->numer=F1->numer;
  F1->numer=F2->numer;
  F2->numer=F->numer;
  if (!SimpPolFract(F1,FALSE)) return FALSE;
  if (!SimpPolFract(F2,FALSE)) return FALSE;
  if (!MultPolynom(&F1->numer,&F2->numer,&F->numer)) return FALSE;
  if (!MultPolynom(&F1->denom,&F2->denom,&F->denom)) return FALSE;
  return TRUE;
}


void main()  { 

  printf("\n MULTIPLICATION OF TWO POLYNOMIAL FRACTIONS:\n\n");
  EnterPolFract(" Enter polynomial fraction F1(x)):", &F1);
  printf("\n");
  EnterPolFract(" Enter polynomial fraction F2(x):", &F2);
  printf("\n");

  if (MultPolFract(&F1,&F2,&F))
    //display result F
    DisplayPolFract(&F);
  else
    printf(" Error in multiplication.\n");
  printf("\n\n");

}

//end of file multfrac.cpp