/********************************************************
*  Calculate a limited development of a real function   *
*  f(x)og(x) at x=0 up at order 25, knowing the limited *
*  developments of f(x) and g(x).                       *
* ----------------------------------------------------- *
* SAMPLE RUN:                                           *
*                                                       *
* Limited development of a real function f(x)og(x):     *
*                                                       *
* Function to develop: exp(sin(x))                      *
*                                                       *
* Limited development of f(x)=exp(x) is:                *
* 1 + x + x^2/2! + x^3/3! + ... + x^n/n! + ...          *
*                                                       *
* Limited development of g(x)=sin(x) is:                *
* x - x^3/3! +x^5/5! - ... + (-1)^n x^2n+1/(2n+1)! + ...*
*                                                       *
* Order of development (max=25): 10                     *
*                                                       *
* Input coefficients of limited dev. of f(x):           *
*  a0 = 1                                               * 
*  a1 = 1                                               *
*  a2 = 0.5                                             *
*  a3 = 0.16666667                                      *
*  a4 = 0.04166667                                      *
*  a5 = 0.00833333                                      *
*  a6 = 0.00138889                                      *
*  a7 = 0.00019841                                      *
*  a8 = 0.00002480                                      *
*  a9 = 0.00000276                                      *
*  a10 = 0.000000276                                    *
*                                                       *
* Input coefficients of limited dev. of g(x):           *
*  b0 = 0                                               * 
*  b1 = 1                                               *
*  b2 = 0                                               *
*  b3 = -0.16666667                                     *
*  b4 = 0                                               *
*  b5 = 0.00833333                                      *
*  b6 = 0                                               *
*  b7 = -0.00019841                                     *
*  b8 = 0                                               *
*  b9 = 0.00000276                                      *
*  b10 = 0                                              *                          
*                                                       *
* The coefficients of limited dev. of f(x)og(x) are:    *
*  c0 =  1.0000000                                      *
*  c1 =  1.0000000                                      *
*  c2 =  0.5000000                                      *
*  c3 =  0.0000000                                      *
*  c4 = -0.1250000                                      *
*  c5 = -0.0666667                                      *
*  c6 = -0.0041667                                      *
*  c7 =  0.0111111                                      *
*  c8 =  0.0053819                                      *
*  c9 =  0.0001764                                      *
*  c10 = -0.0008132                                     *
*                                                       *
* ----------------------------------------------------- *
* Ref.: "Mathematiques en Turbo Pascal By Alain         *
*        Reverchon and Marc Ducamp, Editions Eyrolles,  *
*        Paris, 1991" [BIBLI 03].                       *
*                                                       *
*                            C++ version by J-P Moreau. *
*                                (www.jpmoreau.fr)      *
********************************************************/
#include <stdio.h>
#include <math.h>

#define  SIZE  25

typedef double Tbinome[SIZE+1][SIZE+1];

int i,m;
double T1[SIZE+1],T2[SIZE+1],R[SIZE+1];


    double ProductBin(int i,double *t2,int *id, Tbinome binome)   {
      int k,p,q; double u;
      u=t2[id[1]];
      p=i; q=1;
      for (k=2; k<=i; k++) {
        u *= t2[id[k]];
        if (id[k]==id[k-1])
          q++;
        else {
          u = u*binome[p][q];
          p = p-q;
          q = 1;
        }
      }
      return u;
	}


void ar_dlcomp(int n, double *t1, double *t2, double *res)  {
  int i,j,m,pos; double x,s;
  Tbinome binome;
  int id[SIZE+1];
  if (n > SIZE) return;
  if (t2[0] != 0) return;
  res[0]=t1[0];
  //calculate coefficients of binome
  binome[1][0]=1;
  binome[1][1]=1; 
  for (i=2; i<=n; i++) {
    binome[i][0]=1;
    binome[i][i]=1;
    for (j=1; j<i; j++)
      binome[i][j]=binome[i-1][j-1]+binome[i-1][j];
  }
  //calculate coefficients of gof
  for (m=1; m<=n; m++) {
    x=t1[1]*t2[m];
    for (i=2; i<=m; i++) {
      //calculate s=gamma(i,m)
      s=0;
      for (j=1; j<i; j++) id[j]=1;
      pos=1;
      do {
        for (j=i-pos+1; j<i; j++) id[j]=id[i-pos];
        id[i]=m;
        for (j=1; j<i; j++) id[i]=id[i]-id[j];
        while (id[i] >= id[i-1]) {
          pos=1;
          s += ProductBin(t2,id,binome);
          id[i]=id[i]-1;
          id[i-1]=id[i-1]+1;
        }
        pos++;
        if (pos < i)  id[i-pos]=id[i-pos]+1;
      } while (pos<i); 
      x += s*t1[i];
    }
    res[m]=x;
  }
}


void main()  {

  printf("\n Limited development of a real function f(x)og(x):\n\n");
  printf(" Function to develop: exp(sin(x))\n\n");
  printf(" Limited development of f(x)=exp(x) is:\n");
  printf("  1 + x + x^2/2! + x^3/3! + ... + x^n/n! + ...\n\n");
  printf(" Limited development of g(x)=sin(x) is:\n");
  printf("  x - x^3/3! +x^5/5! - ... + (-1)^n x^2n+1/(2n+1)! + ...\n\n");
  printf(" Order of development (max=25): "); scanf("%d", &m);
  printf("\n Input coefficients of limited dev. of f(x):\n");
  for (i=0; i<=m; i++) {
    printf("  a%d = ",i); scanf("%lf", &T1[i]);
  }
  printf("\n Input coefficients of limited dev. of g(x):\n");
  for (i=0; i<=m; i++) {
    printf("  b%d = ",i); scanf("%lf", &T2[i]);
  }
  printf("\n");

  ar_dlcomp(m,T1,T2,R);

  printf(" The coefficients of limited dev. of f(x)og(x) are:\n");
  for (i=0; i<=m; i++) 
    printf("  c%d = %10.7f\n", i, R[i]);
  printf("\n\n");

}

// end of file ltddev3.cpp