/* Arithmetic operations on polynomials with rational coefficients
 *
 * In the following descriptions a, b, c are polynomials of degree
 * na, nb, nc respectively.  The degree of a polynomial cannot
 * exceed a run-time value FMAXPOL.  An operation that attempts
 * to use or generate a polynomial of higher degree may produce a
 * result that suffers truncation at degree FMAXPOL.  The value of
 * FMAXPOL is set by calling the function
 *
 *     polini( maxpol );
 *
 * where maxpol is the desired maximum degree.  This must be
 * done prior to calling any of the other functions in this module.
 * Memory for internal temporary polynomial storage is allocated
 * by polini().
 *
 * Each polynomial is represented by an array containing its
 * coefficients, together with a separately declared integer equal
 * to the degree of the polynomial.  The coefficients appear in
 * ascending order; that is,
 *
 *                                        2                      na
 * a(x)  =  a[0]  +  a[1] * x  +  a[2] * x   +  ...  +  a[na] * x  .
 *
 *   wrapper functions to the following:
 *
 * `a', `b', `c' are arrays of fracts.
 * fpoleva( a, na, &x, &sum );	Evaluate polynomial a(t) at t = x.
 * fpoladd( a, na, b, nb, c );	c = b + a, nc = max(na,nb)
 * fpolsub( a, na, b, nb, c );	c = b - a, nc = max(na,nb)
 * fpolmul( a, na, b, nb, c );	c = b * a, nc = na+nb
 *
 *
 * Division:
 *
 * i = fpoldiv( a, na, b, nb, c );	c = b / a, nc = FMAXPOL
 *
 * returns i = the degree of the first nonzero coefficient of a.
 * The computed quotient c must be divided by x^i.  An error message
 * is printed if a is identically zero.
 *
 *
 * Change of variables:
 * If a and b are polynomials, and t = a(x), then
 *     c(t) = b(a(x))
 * is a polynomial found by substituting a(x) for t.  The
 * subroutine call for this is
 *
 * fpolsbt( a, na, b, nb, c );
 *
 *
 * Notes:
 * fpoldiv() is an integer routine; fpoleva() is double.
 * Any of the arguments a, b, c may refer to the same array.
 *
 */

#include <stdio.h>
#include "mconf.h"
#ifndef NULL
#define NULL 0
#endif
typedef struct{
	double n;
	double d;
	}fract;

#ifdef ANSIPROT
extern void * malloc ( long );
extern void free ( void * );
#else
void * malloc();
void free ();
#endif

int FMAXPOL = 0;
extern int FMAXPOL;


void fpoladd_wrap( an, ad, na, bn, bd, nb, cn, cd, nc)
double an[], ad[], bn[], bd[], cn[], cd[];
int na, nb, nc;
{
  extern void fpoladd(  fract a[], int na, fract b[], int nb, fract c[]);
  fract *a, *b, *c;
  int j;
  a = (fract *) malloc( (na+1) * sizeof (fract) ); 
  b = (fract *) malloc( (nb+1) * sizeof (fract) ); 
  c = (fract *) malloc( (nc+1) * sizeof (fract) ); 
  for (j=0; j<=na; j++) {
    a[j].n = an[j];
    a[j].d = ad[j];
  }
  for (j=0; j<=nb; j++) {
    b[j].n = bn[j];
    b[j].d = bd[j];
  }
  for (j=0; j<=nc; j++) {
    c[j].n = 0;
    c[j].d = 1;
  }
  fpoladd(a, na, b, nb, c);
  for (j=0; j<=nc; j++) {
    cn[j] = c[j].n;
    cd[j] = c[j].d;
  }
  free(a);
  free(b);
  free(c);
}

void fpolsub_wrap( an, ad, na, bn, bd, nb, cn, cd, nc)
double an[], ad[], bn[], bd[], cn[], cd[];
int na, nb, nc;
{
  extern void fpolsub(  fract a[], int na, fract b[], int nb, fract c[]);
  fract *a, *b, *c;
  int j;
  a = (fract *) malloc( (na+1) * sizeof (fract) ); 
  b = (fract *) malloc( (nb+1) * sizeof (fract) ); 
  c = (fract *) malloc( (nc+1) * sizeof (fract) ); 
  for (j=0; j<=na; j++) {
    a[j].n = an[j];
    a[j].d = ad[j];
  }
  for (j=0; j<=nb; j++) {
    b[j].n = bn[j];
    b[j].d = bd[j];
  }
  for (j=0; j<=nc; j++) {
    c[j].n = 0;
    c[j].d = 1;
  }
  fpolsub(a, na, b, nb, c);
  for (j=0; j<=nc; j++) {
    cn[j] = c[j].n;
    cd[j] = c[j].d;
  }
  free(a);
  free(b);
  free(c);
}

void fpolmul_wrap( an, ad, na, bn, bd, nb, cn, cd, nc)
double an[], ad[], bn[], bd[], cn[], cd[];
int na, nb, nc;
{
  extern void fpolmul(  fract a[], int na, fract b[], int nb, fract c[]);
  fract *a, *b, *c;
  int j;
  a = (fract *) malloc( (na+1) * sizeof (fract) ); 
  b = (fract *) malloc( (nb+1) * sizeof (fract) ); 
  c = (fract *) malloc( (nc+1) * sizeof (fract) ); 
  for (j=0; j<=na; j++) {
    a[j].n = an[j];
    a[j].d = ad[j];
  }
  for (j=0; j<=nb; j++) {
    b[j].n = bn[j];
    b[j].d = bd[j];
  }
  for (j=0; j<=nc; j++) {
    c[j].n = 0;
    c[j].d = 1;
  }
  fpolmul(a, na, b, nb, c);
  for (j=0; j<=nc; j++) {
    cn[j] = c[j].n;
    cd[j] = c[j].d;
  }
  free(a);
  free(b);
  free(c);
}

int fpoldiv_wrap( an, ad, na, bn, bd, nb, cn, cd, nc)
double an[], ad[], bn[], bd[], cn[], cd[];
int na, nb, nc;
{
  extern int fpoldiv(  fract a[], int na, fract b[], int nb, fract c[]);
  fract *a, *b, *c;
  int j, ret;
  a = (fract *) malloc( (na+1) * sizeof (fract) ); 
  b = (fract *) malloc( (nb+1) * sizeof (fract) ); 
  c = (fract *) malloc( (nc+1) * sizeof (fract) ); 
  for (j=0; j<=na; j++) {
    a[j].n = an[j];
    a[j].d = ad[j];
  }
  for (j=0; j<=nb; j++) {
    b[j].n = bn[j];
    b[j].d = bd[j];
  }
  for (j=0; j<=nc; j++) {
    c[j].n = 0;
    c[j].d = 1;
  }
  ret = fpoldiv(a, na, b, nb, c);
  for (j=0; j<=nc; j++) {
    cn[j] = c[j].n;
    cd[j] = c[j].d;
  }
  free(a);
  free(b);
  free(c);
  return ret;
}

void fpoleva_wrap( an, ad, na, x, s)
double an[], ad[];
int na;
fract *x, *s;
{
  extern void fpoleva(  fract a[], int na, fract *x, fract *s);
  fract *a;
  int j;
  a = (fract *) malloc( (na+1) * sizeof (fract) ); 
  for (j=0; j<=na; j++) {
    a[j].n = an[j];
    a[j].d = ad[j];
  }
  s->n = 0.0;
  s->d = 1.0;
  fpoleva(a, na, x, s);
  free(a);
}

void fpolsbt_wrap( an, ad, na, bn, bd, nb, cn, cd, nc)
double an[], ad[], bn[], bd[], cn[], cd[];
int na, nb, nc;
{
  extern void fpolsbt(  fract a[], int na, fract b[], int nb, fract c[]);
  fract *a, *b, *c;
  int j;
  a = (fract *) malloc( (na+1) * sizeof (fract) ); 
  b = (fract *) malloc( (nb+1) * sizeof (fract) ); 
  c = (fract *) malloc( (nc+1) * sizeof (fract) ); 
  for (j=0; j<=na; j++) {
    a[j].n = an[j];
    a[j].d = ad[j];
  }
  for (j=0; j<=nb; j++) {
    b[j].n = bn[j];
    b[j].d = bd[j];
    }
  for (j=0; j<=nc; j++) {
    c[j].n = 0;
    c[j].d = 1;
  }
  fpolsbt(a, na, b, nb, c);
  for (j=0; j<=nc; j++) {
    cn[j] = c[j].n;
    cd[j] = c[j].d;
  }
  free(a);
  free(b);
  free(c);
}