harmonic.c

#include <stdio.h>

void Series_Sum(double sum[]);

int main()
{
    int i;
    double x, sum[3001];
	
    Series_Sum( sum );

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

    return 0;
}

// Order
#define O 10
// Number of iterations
#define M 80
#include <math.h>
double fact[O-1] = {1};
double cal_sum(double x);

void Series_Sum(double sum[])
{
    // Calculating the integeral x
    sum[0] = 3.14159265358979*3.14159265358979/6;
    for(int x = 1; x <= 300; x++)
    {
        double result = 0.;
        for(int i = x; i > 0; i--)
        {
            result += 1./i;
        }
        result /= x;
        sum[10*x] = result;
    }
    
    // Factorial list
    for(int i = 1; i <= O - 2; i++)
    {
        fact[i] = fact[i-1]*i;
    }

    // Calculating else
    double x = 0.1;
    for(int i = 1; i < 100; i++, x += 0.1)
    {
        if(!(i%10)) continue;
        sum[i] = cal_sum(x);
    }

    x = 10.1;
    for(int i = 101; i < 3000; i++, x += 0.1)
    {
        if(!(i%10)) continue;
        double x2 = x*x;
        sum[i] =(-420 + x2*(
                  231 + x2*(
                 -220 + x2*(
                  462 + x2*(
                -4620 + x *(
                27720 + x *(
                55440*(0.57721566490153286060651209+log(x)))))))))
                /(55440*x2*x2*x2*x2*x2*x);
    }
}

double cal_sum(double x)
{
    double result = 0.;

    // Calculating the core sum
    for(int k = M; k > 0; k--)
    {
        double product = x + k;
        for(int i = 0; i <= O - 2; i++)
        {
            product *= k + i;
        }
        result += 1./product;
    }
    result += pow(M+1, 1-O)/(O-1);
    // Now that the result is the core sum
    for(int i = O - 2; i > 0; i--)
    {
        result *= i - x;
        result += 1./i/fact[i];
    }
    return result;
}

Version 1 of “Numerical Summation of a Series”