/******************************************************** * 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 #include #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= 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