PowerSeriesCombine.cin

From TORI
Jump to navigation Jump to search

PowerSeriesCombine.cin is the C++ routine that combines two asymptotic series.

Specification: void PowerSeriesCombine( 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.

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   PowerSeriesCombine(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) \).

PowerSeriesCombine.cin

//
void PowerSeriesCombine(double *a,
                        double *b,
                        double *c,
                        int N)
{
    // a(z) = b(c(z))
    // b[0..N], c[0..N]
    // assume c[0]=0

    a[0] = b[0];

    for(int n=1; n<=N; n++)
    {
        double sum = 0.0;

        for(int k=1; k<=n; k++)
        {
            // compute coefficient of z^n in c(z)^k

            double power[N+1];	// = {0};
            double next[N+1];	//  = {0};

            // initialize
            for(int i=0;i<=n;i++)
                power[i] = c[i];

            // raise to k-th power
            for(int p=2; p<=k; p++)
            {
                for(int i=0;i<=n;i++)
                    next[i] = 0.0;

                for(int i=0;i<=n;i++)
                    for(int j=0;j<=n;j++)
                        if(i+j<=n)
                            next[i+j] += power[i]*c[j];

                for(int i=0;i<=n;i++)
                    power[i] = next[i];
            }

            sum += b[k] * power[n];
        }

        a[n] = sum;
    }
}
//

Example 1

//
#include <math.h>
#include  <stdio.h>
//#include  "PowerSeriesInversion.cin"
#include  "PowerSeriesCombine.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("\nPowerSeriesCombine(a,b,c,N):\n");
	    PowerSeriesCombine(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("\nPowerSeriesCombine(a,b,b,N):\n");
	    PowerSeriesCombine(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");
}

//

does

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 

PowerSeriesCombine(a,b,c,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 

PowerSeriesCombine(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 

Notes, attribution, modifications

The routine PowerSeriesCombine is generated by ChatGPT.
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 should allow 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