/* mpf_deriv.c  multiple precision derivative coefficients */
/*              has mpz_deriv and  mpf_deriv does rderiv   */

#include <stdio.h>
#include "mpf_deriv.h" /* includes gmp and digits */

static mpf_t atemp;
static mpf_t btemp;
static int debug = 0;
void mpz_deriv(int order, int npoints, int point, mpz_t a[], mpz_t bb);

void mpf_rderiv(int order, int npts, int point, mpf_t h, mpf_t cx[])
{
  int i, j;
  mpz_t a[50], b;

  mpz_init(b);
  for(i=0; i<50; i++) mpz_init(a[i]);

  mpf_set_default_prec(digits*3.32); /* in bits */
  mpf_init(atemp);
  mpf_init(btemp);

  mpz_deriv(order, npts, point, a, b);

  mpf_set_z(btemp, b);
  for(i=0; i<npts; i++)
  {
    mpf_set_z(atemp, a[i]);
    mpf_div(cx[i], atemp, btemp);
    for(j=0; j<order; j++) mpf_div(cx[i], cx[i], h);
  }
} /* end mpf_rderiv   rderiv */

void mpz_deriv(int order, int npoints, int point, mpz_t a[], mpz_t bb)
{
  /* compute the exact coefficients to numerically compute a derivative */
  /* of order 'order' using 'npoints' at ordinate point,                */
  /* where order>=1, npoints>=order+1, 0 <= point < npoints,            */
  /* 'a' array returned with numerator of coefficients,                 */
  /* 'bb' returned with denominator of h^order                          */
  
  mpz_t h[50];          /* coefficients of h */
  mpz_t x[50];          /* coefficients of x, variable for differentiating */
  mpz_t numer[50][50];  /* numerator of a term numer[term][pos] */
  mpz_t denom[50];      /* denominator of coefficient */
  int i, j, k, term, jorder, iat, jterm;
  mpz_t ipower, isum, r, tmp, tmp2;

  mpz_init(ipower);
  mpz_init(isum);
  mpz_init(r);
  mpz_init(tmp);
  mpz_init(tmp2);

  if(npoints<=order || npoints>=49)
  {
    printf("ERROR in call to deriv, npoints=%d >=49 || < order=%d+1=%d \n\n",
           npoints, order, order+1);
    return;
  }

  for(i=0; i<=npoints; i++)
  {
    mpz_init(a[i]);
    mpz_init(h[i]);
    mpz_init(x[i]);
    mpz_init(denom[i]);
    for(j=0; j<=npoints; j++)
      mpz_init(numer[i][j]);
  }

  for(term=0; term<npoints; term++)
  {
    mpz_set_si(denom[term], 1);
    j = 0;
    for(i=0; i<npoints-1; i++)
    {
      if(j==term) j++; /* no index j in jth term */
      mpz_set_si(numer[term][i], -j); /* from (x-x[j]) */
      mpz_set_si(tmp, term-j);
      mpz_mul(denom[term], denom[term], tmp);
      j++;
    }
  } /* end term setting up numerator and denominator */
  if(debug>0)
    for(i=0; i<npoints; i++)
    {
      gmp_printf("denom[%d]=%Zd \n", i, denom[i]);
      for(j=0; j<npoints-1; j++)
        gmp_printf("numer[%d][%d]=%Zd \n", i, j, numer[i][j]);
    }
        
  /* substitute xj = x0 + j*h, just keep j value in numer[term][i] */
  /* don't need to count x0 because same as x's */
  /* done by above setup */
      
  /* multiply out polynomial in x with coefficients of h */
  /* starting from (x+j1*h)*(x+j2*h) * ... */
  /* then differentiate d/dx */
  for(term=0; term<npoints; term++)
  {
    mpz_set_si(x[0], 1); /* initialize */
    mpz_set_si(x[1], 0);
    mpz_set_si(h[1], 0);
    k = 1;
    for(i=0; i<npoints-1; i++)
    {
      for(j=0; j<k; j++)            /* multiply by (x + j h) */
      {
        mpz_mul(h[j], x[j], numer[term][i]); /* multiply by constant x[j] */
      }
      for(j=0; j<k; j++)         /* multiply by x */
      {
        mpz_add(h[j+1], h[j+1], x[j]); /* multiply by x and add */
      }
      k++;
      for(j=0; j<k; j++) mpz_set(x[j], h[j]);
      mpz_set_si(x[k], 0);
      mpz_set_si(h[k], 0);
    }
    for(j=0; j<k; j++)
    {
      mpz_set(numer[term][j], x[j]);
    } /* have p(x) for this 'term' */

    /* differentiate 'order' number of times */
    for(jorder=0; jorder<order; jorder++)
    {
      for(i=1; i<k; i++) /* differentiate p'(x) */
      {
        mpz_mul_si(numer[term][i-1], numer[term][i], i);
      }
      k--;
    } /* have p^(order) for this term */
    if(debug>0)
      for(i=0; i<k; i++)
        gmp_printf("p^(%d)(x)= %Zd x^%d at term=%d \n",
               order, numer[term][i], i, term);
  } /* end 'term' loop  for one 'order', for one 'npoints' */

  /* substitute each point for x, evaluate, get coefficients of f(x[j]) */
  iat = point; /* evaluate at x[iat], substitute */
  for(jterm=0; jterm<npoints; jterm++) /* get each term */
  {
    mpz_set_si(ipower, iat);
    mpz_set(isum, numer[jterm][0]);
    for(j=1; j<k; j++)
    {
      mpz_mul(tmp, ipower, numer[jterm][j]);
      mpz_add(isum, isum, tmp);
      mpz_mul_si(ipower, ipower, iat);
    }
    mpz_set(a[jterm], isum);
    if(debug>0)
      gmp_printf("f^(%d)(x[%d]) = (1/h^%d) (%Zd/%Zd f(x[%d]) + \n",
             order, iat, order, a[jterm], denom[jterm], jterm);
  }
  if(debug>0) printf("\n");
  mpz_set_si(bb, 0);
  for(jterm=0; jterm<npoints; jterm++) /* clean up each term */
  {
    if(mpz_cmp_si(a[jterm], 0)==0) continue;
    mpz_gcd(r, a[jterm], denom[jterm]);
    mpz_div(a[jterm], a[jterm], r);
    mpz_div(denom[jterm], denom[jterm], r);
    if(mpz_cmp_si(denom[jterm],0)<0)
    {
      mpz_neg(denom[jterm], denom[jterm]);
      mpz_neg(a[jterm], a[jterm]);
    }
    if(mpz_cmp(denom[jterm],bb)>0) mpz_set(bb, denom[jterm]);
                                                  /* largest denominator */
  }
  for(jterm=0; jterm<npoints; jterm++) /* check for divisor */
  {
    mpz_mul(tmp, a[jterm], bb);
    mpz_div(r, tmp, denom[jterm]);
    mpz_mul(tmp2, r, denom[jterm]);
    if(mpz_cmp(tmp2, tmp) != 0)
    {
      mpz_gcd(r, bb, denom[jterm]);
      mpz_mul(tmp, bb, denom[jterm]);
      mpz_div(bb, tmp, r);
    }
  }    
  for(jterm=0; jterm<npoints; jterm++) /* adjust for divisor */
  {
    mpz_mul(tmp, a[jterm], bb);
    mpz_div(a[jterm], tmp, denom[jterm]);
    if(debug>0)
      gmp_printf("f^(%d)(x[%d])=(1/%Ddh^%d)(%Dd f(x[%d]) + \n",
             order, iat, bb, order, a[jterm], jterm);
  } /* end computing terms of coefficients */
} /* end mpz_deriv */