harmonic2.c

#include <stdio.h>
#include <string.h>

void Series_Sum(double sum[]);

int main(void)
{
    int i;
    double x, sum[3001];
    memset((void*)sum, 0xffffffff, sizeof(sum));

    Series_Sum(sum);

    x = 0.0;
    for(i = 0; i < 3001; i++)
    {
        printf("%6.2f %16.20f\n", x + (double)i * 0.10, sum[i]);
    }
    return 0;
}

void Series_Sum(double sum[])
{
    // Calculating the integer x
    sum[0] = 3.14159265358979*3.14159265358979/6;
    for(int x = 1; x <= 300; x++)
    {
        double result = 0.;
        for(int i = 1; i <= x; i++)
        {
            result += 1./i;
        }
        result /= x;
        sum[10*x] = result;
    }
    
    // Calculating else
    double x = 0.1;
    for(int i = 1; i < 3000; i++, x += 0.1)
    {
        if(!(i%10)) continue;
        sum[i] = 0;
        double result = 0.;
        int m = (int)x + 1;
        double f0 = sum[m*10], f1 = sum[m*10-10], f2 = sum[m*10-20];
        int M = 800;

        if(x < 1) 
        {
            int MM = 1500;
            for(int k = MM; k > 0; k--)
            {
                result += 1/(k*(k+m)*(k+x)*(k+m-1));
            }
            result += 1./(3.*MM*MM*MM);
            result = (m-x)*((m-1-x)*result+f1-f0)+f0;
        }
        else
        {
            for(int k = M; k > 0; k--)
            {
                result += 1/((k+x)*k*(k+m)*(k+m-1)*(k+m-2));
            }
            result += 1./(4.*M*M*M*M);
            result = (m-x)*((m-1-x)*((m-2-x)*result*2 + f2 + f0 - 2*f1)/2+f1-f0)+f0;
        }
        sum[i] = result;
    }
}

Version 2 of “Numerical Summation of a Series”