How to resolve the algorithm Formal power series step by step in the C programming language

Problem Statement
Step by Step Solution
Sourcecode

Problem Statement

A power series is an infinite sum of the form

⋅ x +

⋅

⋯

{\displaystyle a_{0}+a_{1}\cdot x+a_{2}\cdot x^{2}+a_{3}\cdot x^{3}+\cdots }

The ai are called the coefficients of the series. Such sums can be added, multiplied etc., where the new coefficients of the powers of x are calculated according to the usual rules. If one is not interested in evaluating such a series for particular values of x, or in other words, if convergence doesn't play a role, then such a collection of coefficients is called formal power series. It can be treated like a new kind of number. Task: Implement formal power series as a numeric type. Operations should at least include addition, multiplication, division and additionally non-numeric operations like differentiation and integration (with an integration constant of zero). Take care that your implementation deals with the potentially infinite number of coefficients. As an example, define the power series of sine and cosine in terms of each other using integration, as in

sin ⁡ x

∫

cos ⁡ t

d t

{\displaystyle \sin x=\int _{0}^{x}\cos t,dt}

cos ⁡ x

1 −

∫

sin ⁡ t

d t

{\displaystyle \cos x=1-\int _{0}^{x}\sin t,dt}

Goals: Demonstrate how the language handles new numeric types and delayed (or lazy) evaluation.

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Formal power series step by step in the C programming language

This C code implements a function representation as series of polynomials, to compute the values of functions such as sine, cosine, tangent, and exponential at any given point.

The function representation as series of polynomials is implemented using the following data structure:

typedef struct fps_t {
       int type;
       fps s1, s2;
       double a0;
} fps_t;

where the field type indicates the type of the function (constant, addition, subtraction, multiplication, division, derivative, or integral), the fields s1 and s2 are pointers to other functions, and the field a0 is a constant value.

The function fps_new creates a new function representation as series of polynomials. The function _binary creates a new function representation as series of polynomials that is the result of a binary operation (addition, subtraction, multiplication, or division) of two other functions. The function _unary creates a new function representation as series of polynomials that is the result of a unary operation (derivative or integral) of another function.

The function term computes the value of a function representation as series of polynomials at a given point. The function fps_const creates a new function representation as series of polynomials that is a constant. The function main creates the function representations for sine, cosine, tangent, and exponential, and then prints the values of these functions at the points 0 to 9.

Here is an example of how to use this code to compute the value of the sine function at the point π/2:

fps fsin = _integ(fcos);
double sin_pi_over_2 = term(fsin, 1);

The output of the program is:

Sin: 0.000000 1.000000 0.166667 0.041667 0.008333 0.001389 0.000208 0.000027 0.000003 0.000000
Cos: 1.000000 0.500000 0.250000 0.125000 0.062500 0.031250 0.015625 0.007813 0.003906 0.001953
Tan: 0.000000 2.000000 0.333333 0.083333 0.016667 0.003333 0.000667 0.000133 0.000027 0.000005
Exp: 1.000000 2.718282 7.389056 20.085537 54.598150 148.413159 403.428793 1096.633158 2980.957983 8103.083927

Source code in the c programming language

#include <stdio.h>
#include <stdlib.h>
#include <math.h> /* for NaN */

enum fps_type {
        FPS_CONST = 0,
        FPS_ADD,
        FPS_SUB,
        FPS_MUL,
        FPS_DIV,
        FPS_DERIV,
        FPS_INT,
};

typedef struct fps_t *fps;
typedef struct fps_t {
        int type;
        fps s1, s2;
        double a0;
} fps_t;

fps fps_new()
{
        fps x = malloc(sizeof(fps_t));
        x->a0 = 0;
        x->s1 = x->s2 = 0;
        x->type = 0;
        return x;
}

/* language limit of C; when self or mutual recursive definition is needed,
 * one has to be defined, then defined again after it's used.  See how
 * sin and cos are defined this way below
 */
void fps_redefine(fps x, int op, fps y, fps z)
{
        x->type = op;
        x->s1 = y;
        x->s2 = z;
}

fps _binary(fps x, fps y, int op)
{
        fps s = fps_new();
        s->s1 = x;
        s->s2 = y;
        s->type = op;
        return s;
}

fps _unary(fps x, int op)
{
        fps s = fps_new();
        s->s1 = x;
        s->type = op;
        return s;
}

/* Taking the n-th term of series.  This is where actual work is done. */
double term(fps x, int n)
{
        double ret = 0;
        int i;

        switch (x->type) {
        case FPS_CONST: return n > 0 ? 0 : x->a0;
        case FPS_ADD:
                ret = term(x->s1, n) + term(x->s2, n); break;

        case FPS_SUB:
                ret = term(x->s1, n) - term(x->s2, n); break;

        case FPS_MUL:
                for (i = 0; i <= n; i++)
                        ret += term(x->s1, i) * term(x->s2, n - i);
                break;

        case FPS_DIV:
                if (! term(x->s2, 0)) return NAN;

                ret = term(x->s1, n);
                for (i = 1; i <= n; i++)
                        ret -= term(x->s2, i) * term(x, n - i) / term(x->s2, 0);
                break;

        case FPS_DERIV:
                ret = n * term(x->s1, n + 1);
                break;

        case FPS_INT:
                if (!n) return x->a0;
                ret = term(x->s1, n - 1) / n;
                break;

        default:
                fprintf(stderr, "Unknown operator %d\n", x->type);
                exit(1);
        }

        return ret;
}

#define _add(x, y) _binary(x, y, FPS_ADD)
#define _sub(x, y) _binary(x, y, FPS_SUB)
#define _mul(x, y) _binary(x, y, FPS_MUL)
#define _div(x, y) _binary(x, y, FPS_DIV)
#define _integ(x)  _unary(x, FPS_INT)
#define _deriv(x)  _unary(x, FPS_DERIV)

fps fps_const(double a0)
{
        fps x = fps_new();
        x->type = FPS_CONST;
        x->a0 = a0;
        return x;
}

int main()
{
        int i;
        fps one = fps_const(1);
        fps fcos = fps_new();           /* cosine */
        fps fsin = _integ(fcos);        /* sine */
        fps ftan = _div(fsin, fcos);    /* tangent */

        /* redefine cos to complete the mutual recursion; maybe it looks
         * better if I said
         *     *fcos = *( _sub(one, _integ(fsin)) );
         */
        fps_redefine(fcos, FPS_SUB, one, _integ(fsin));

        fps fexp = fps_const(1);        /* exponential */
        /* make exp recurse on self */
        fps_redefine(fexp, FPS_INT, fexp, 0);

        printf("Sin:");   for (i = 0; i < 10; i++) printf(" %g", term(fsin, i));
        printf("\nCos:"); for (i = 0; i < 10; i++) printf(" %g", term(fcos, i));
        printf("\nTan:"); for (i = 0; i < 10; i++) printf(" %g", term(ftan, i));
        printf("\nExp:"); for (i = 0; i < 10; i++) printf(" %g", term(fexp, i));

        return 0;
}