PowerSeriesCombinA.cin

From TORI
Jump to navigation Jump to search

PowerSeriesCombinA.cin is the C++ routine that combines two asymptotic series, analogy of routine PowerSeriesCombine.cin

Specification: void PowerSeriesCombinA( double *a, double *b, double *c, int N)

a,b,c are arrays of length N+1. The 0th elements of these arrays in not used in this computation. The first element of array c is supposed to be unity.
Array a should be declared to store N+1 elements.
Array b should store at least N coefficients of power series \( B(z)= \sum_{n=1}^N b_n z^n + o(z^N) \)
Array c should store at least N coefficients of power series \( C(z)= \sum_{n=1}^N c_n z^n + o(z^N) \)

Then, after to call   PowerSeriesCombinA(a,b,c,N);
array a is supposed to provide N coefficients a[1] .. a[N] of the power series \(\ A(z) = B(C(z)) = \sum_{n=1}^N a_n z^n + o(z^N) \).

WARNING: Call "Arbogast(b,c,n)" assumes that the inner series is tangent to identity, i.e., c[1]==1.
Therefore, call "PowerSeriesCombinA(a,b,c,n)" also requires c[1]==1.
If this condition is not satisfied, use PowerSeriesCombine.cin instead.

PowerSeriesCombinA.cin

//Arbogast from «Arbogast.cin» should be also loaded.

//
//The main module should #include "Arbogast.cin" before to #include "PowerSeriesCombinA.cin"
void PowerSeriesCombinA(double *a,
                        double *b,
                        double *c,
                        int N)
{
    a[0] = b[0];

    a[1] = b[1] * c[1];

    for(int n=2; n<=N; n++)
    {
        a[n] = b[1]*c[n]
               + Arbogast(b, c, n);
    }
}
//

Example 1

//

#include <math.h>
#include  <stdio.h>
//#include  "PowerSeriesInversion.cin"
//#include  "PowerSeriesCombine.cin"
#include  "Arbogast.cin"
#include  "PowerSeriesCombinA.cin"

int main(){ int n; int N=5; int N1=N+1;
double a[N1]; 
double b[N1]; 
double c[N1]; 
             for(n=1;n<N1;n++) {      a[n]=0.;}
double f=1.; for(n=1;n<N1;n++) {f/=n; b[n]=f;}
double k=-1; for(n=1;n<N1;n++) {k*=-1;c[n]=(0.+k)/n;}
printf("before:\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",a[n]); printf("\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",b[n]); printf("\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",c[n]); printf("\n");
//printf("\nPowerSeriesInversion(a,b,N):\n");
//	    PowerSeriesInversion(a,b,N);
printf("\nPowerSeriescombinA(a,b,c,N):\n");
	    PowerSeriesCombinA(a,b,c,N);
for(n=0;n<N1;n++)  printf("%21.18lf ",a[n]); printf("\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",b[n]); printf("\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",c[n]); printf("\n");
printf("\nPowerSeriescombinA(a,c,b,N):\n");
	    PowerSeriesCombinA(a,c,b,N);
for(n=0;n<N1;n++)  printf("%21.18lf ",a[n]); printf("\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",b[n]); printf("\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",c[n]); printf("\n");
printf("\nPowerSeriescombinA(a,b,b,N):\n");
	    PowerSeriesCombinA(a,b,b,N);
for(n=0;n<N1;n++)  printf("%21.18lf ",a[n]); printf("\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",b[n]); printf("\n");
for(n=0;n<N1;n++)  printf("%21.18lf ",c[n]); printf("\n");
}
//<pre>
does
<pre>
before:
 0.000000000000000000  0.000000000000000000  0.000000000000000000  0.000000000000000000  0.000000000000000000  0.000000000000000000 
 0.000000000000000000  1.000000000000000000  0.500000000000000000  0.166666666666666657  0.041666666666666664  0.008333333333333333 
 0.000000000000000000  1.000000000000000000 -0.500000000000000000  0.333333333333333315 -0.250000000000000000  0.200000000000000011 

PowerSeriescombinA(a,b,c,N):
 0.000000000000000000  1.000000000000000000  0.000000000000000000 -0.000000000000000056  0.000000000000000000  0.000000000000000000 
 0.000000000000000000  1.000000000000000000  0.500000000000000000  0.166666666666666657  0.041666666666666664  0.008333333333333333 
 0.000000000000000000  1.000000000000000000 -0.500000000000000000  0.333333333333333315 -0.250000000000000000  0.200000000000000011 

PowerSeriescombinA(a,c,b,N):
 0.000000000000000000  1.000000000000000000  0.000000000000000000 -0.000000000000000028 -0.000000000000000076  0.000000000000000016 
 0.000000000000000000  1.000000000000000000  0.500000000000000000  0.166666666666666657  0.041666666666666664  0.008333333333333333 
 0.000000000000000000  1.000000000000000000 -0.500000000000000000  0.333333333333333315 -0.250000000000000000  0.200000000000000011 

PowerSeriescombinA(a,b,b,N):
 0.000000000000000000  1.000000000000000000  1.000000000000000000  0.833333333333333259  0.624999999999999889  0.433333333333333348 
 0.000000000000000000  1.000000000000000000  0.500000000000000000  0.166666666666666657  0.041666666666666664  0.008333333333333333 
 0.000000000000000000  1.000000000000000000 -0.500000000000000000  0.333333333333333315 -0.250000000000000000  0.200000000000000011 

Example 2

//
#include <cmath>
//#include <cstring>  // for memset
#include <stdio.h>
#include "Arbogast.cin"
#include "PowerSeriesInversion.cin"
//#include "PowerSeriesCombine.cin"
#include "PowerSeriesCombinA.cin"
double B=sqrt(2.); // 1.4142135623730951 
double Q=4.;
double S=log(B);
int main() { int n;
int N=36; int N1=N+1;
double t[N1], a[N1], b[N1],c[N1], d[N1], e[N1];
t[0]=4.;
for(n=1;n<N1;n++) t[n]=t[n-1]*S/n;

a[0]=0.;
a[1]=1.;;
for(n=2;n<N1;n++) a[n] = - Arbogast(t, a, n)/( pow(t[1],n)-t[1] );

PowerSeriesInversion(b,a,N);

PowerSeriesInversion(c,b,N);

PowerSeriesCombinA(d,b,a,N);

PowerSeriesCombinA(e,a,b,N);

for(n=0;n<N1;n++) printf("%02d %19.16lf %19.16lf %19.16lf %19.16lf %19.16lf\n",n,b[n],a[n],c[n],d[n],e[n]);
}
//

does:

00  0.0000000000000000  0.0000000000000000  0.0000000000000000  0.0000000000000000  0.0000000000000000
01  1.0000000000000000  1.0000000000000000  1.0000000000000000  1.0000000000000000  1.0000000000000000
02  0.4485874311952611 -0.4485874311952611 -0.4485874311952611  0.0000000000000000  0.0000000000000000
03  0.1903722467978066  0.2120891200549195  0.2120891200549195  0.0000000000000000  0.0000000000000000
04  0.0778295765369682 -0.1021843675069716 -0.1021843675069716  0.0000000000000000  0.0000000000000000
05  0.0309358603057080  0.0496986830373718  0.0496986830373717  0.0000000000000000  0.0000000000000000
06  0.0120221257690659 -0.0243075903261195 -0.0243075903261195  0.0000000000000000  0.0000000000000000
07  0.0045849888965618  0.0119330883965108  0.0119330883965108 -0.0000000000000000  0.0000000000000000
08  0.0017207423310578 -0.0058736976420089 -0.0058736976420089 -0.0000000000000000 -0.0000000000000000
09  0.0006368109038799  0.0028968672871058  0.0028968672871058  0.0000000000000000 -0.0000000000000000
10  0.0002327696003031 -0.0014309081060793 -0.0014309081060793 -0.0000000000000000 -0.0000000000000000
11  0.0000841455118381  0.0007076637148566  0.0007076637148567  0.0000000000000000 -0.0000000000000001
12  0.0000301156464937 -0.0003503296158730 -0.0003503296158732 -0.0000000000000000  0.0000000000000002
13  0.0000106807458130  0.0001735756004664  0.0001735756004668 -0.0000000000000000 -0.0000000000000003
14  0.0000037565713616 -0.0000860610119292 -0.0000860610119299  0.0000000000000000  0.0000000000000004
15  0.0000013111367785  0.0000426959089013  0.0000426959089023 -0.0000000000000000 -0.0000000000000004
16  0.0000004543791625 -0.0000211930290684 -0.0000211930290696 -0.0000000000000000  0.0000000000000004
17  0.0000001564298463  0.0000105244225997  0.0000105244226009  0.0000000000000000 -0.0000000000000003
18  0.0000000535232764 -0.0000052285174362 -0.0000052285174370  0.0000000000000000 -0.0000000000000001
19  0.0000000182077863  0.0000025984499950  0.0000025984499946  0.0000000000000000  0.0000000000000009
20  0.0000000061604765 -0.0000012917821009 -0.0000012917820980 -0.0000000000000000 -0.0000000000000024
21  0.0000000020737193  0.0000006423759966  0.0000006423759892 -0.0000000000000000  0.0000000000000049
22  0.0000000006946828 -0.0000003195233018 -0.0000003195232871  0.0000000000000000 -0.0000000000000085
23  0.0000000002316515  0.0000001589712671  0.0000001589712419 -0.0000000000000000  0.0000000000000133
24  0.0000000000769124 -0.0000000791094995 -0.0000000791094604 -0.0000000000000000 -0.0000000000000190
25  0.0000000000254309  0.0000000393753741  0.0000000393753185 -0.0000000000000000  0.0000000000000250
26  0.0000000000083756 -0.0000000196019782 -0.0000000196019049 -0.0000000000000000 -0.0000000000000303
27  0.0000000000027482  0.0000000097599628  0.0000000097598731  0.0000000000000000  0.0000000000000337
28  0.0000000000008985 -0.0000000048603109 -0.0000000048602094 -0.0000000000000000 -0.0000000000000340
29  0.0000000000002927  0.0000000024207097  0.0000000024206050  0.0000000000000000  0.0000000000000297
30  0.0000000000000951 -0.0000000012058126 -0.0000000012057171 -0.0000000000000000 -0.0000000000000197
31  0.0000000000000308  0.0000000006007192  0.0000000006006495  0.0000000000000000  0.0000000000000028
32  0.0000000000000099 -0.0000000002993051 -0.0000000002992820 -0.0000000000000000  0.0000000000000220
33  0.0000000000000032  0.0000000001491436  0.0000000001491913 -0.0000000000000000 -0.0000000000000556
34  0.0000000000000010 -0.0000000000743259 -0.0000000000744721 -0.0000000000000000  0.0000000000000986
35  0.0000000000000003  0.0000000000370440  0.0000000000373185 -0.0000000000000000 -0.0000000000001511
36  0.0000000000000001 -0.0000000000184644 -0.0000000000188983 -0.0000000000000000  0.0000000000002127

Notes, attribution, modifications

The routine PowerSeriesCombinA is generated by ChatGPT as «PowerSeriesCombine» and renamed by Editor to PowerSeriesCombinA at PowerSeriesCombinA.cin in order to avoid confusion with PowerSeriesCombine.cin, for the case, if some code, in order to compare the performance of the two implementations, will have to call both PowerSeriesCombinA and PowerSeriesCombine.
The honest use is assumed; at the reuse, attribute the source.
The attribution helps to trace bugs, if any.

The 0th elements of these arrays in not used in this computation.
Such a style of numeration allows the straight-forward translation of the code from C++ to some Fortran.

In the current version, the variables are declared to be double.
Editor expects, in this routine, all the specifications «double» can be replaced to something more suitable (for example, «float» or «complex double»).
At the moment of the loading, such a modification is not yet tested.

References