Category: Number Theory

// ThreeIntegrals.cpp
// Copyright (c) Sunday, October 11, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Integrals due to Gauss, Cauchy, and Dirichlet

#include <chrono>
#include <iomanip>
#include <iostream>
#include <vector>

static double Gauss(double x)
{
    return exp(-x * x);
}

static double Cauchy(double x)
{
    return 1.0 / (x * x + 1.0);
}

static double Dirichlet(double x)
{
    return sin(x) / x;
}

static double SimpsonsRule(
    double a, double b, double(*fx)(double), int nSteps)
{
    double h = (b - a) / nSteps;
    double h2 = 2.0 * h;
    double s = 0.0;
    double t = 0.0;
    double x = a + h;

    for (int i = 1; i < nSteps; i += 2)
    {
        s += fx(x);
        x += h2;
    }

    x = a + h2;

    for (int i = 2; i < nSteps; i += 2)
    {
        t += fx(x);
        x += h2;
    }

    return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}

static double ComputePi1(double bound)
{
    return SimpsonsRule(-bound, bound, Gauss, 16384);
}

static double ComputePi2(double bound, int nSteps)
{
    return 2.0 * SimpsonsRule(0.0, bound, Cauchy, nSteps);
}

static double ComputePi3(double bound, int nSteps)
{
    double zero = 1.0e-32;
    double integral1 = SimpsonsRule(-bound, -zero, Dirichlet, nSteps);
    double integral2 = SimpsonsRule(+zero, +bound, Dirichlet, nSteps);
    return integral1 + integral2;
}

static double nFactorial(int n)
{
    double factorial = 1.0;

    for (int i = 1; i <= n; i++)
        factorial *= i;

    return factorial;
}

typedef struct tagGene {
    int R; int nSteps;
} Gene, * PGene;

static double HillClimber(
    int fNumber, int generations, int population, Gene& gene)
{
    double pi = 4.0 * atan(1.0), pi2 = 0.0, pi3 = 0.0;
    std::vector<double> fitness(population);
    std::vector<std::vector<Gene>> genome;

    genome.resize(population);

    for (int i = 0; i < population; i++)
    {
        genome[i].resize(1);

        Gene g = { 0 };
        
        if (fNumber == 2)
        {
            g.R = rand() % 100000;
            g.nSteps = rand() % 100000;

            double pi2 = ComputePi2(g.R, g.nSteps);
            fitness[i] = fabs(pi - pi2);

            if (fitness[i] <= 1.0e-15)
            {
                gene = g;
                return pi2;
            }

            genome[i][0] = g;
        }

        else if (fNumber == 3)
        {
            g.R = rand() % 10000000;
            g.nSteps = rand() % 10000000;

            double pi3 = ComputePi3(g.R, g.nSteps);
            fitness[i] = fabs(pi - pi3);

            if (fitness[i] <= 1.0e-15)
            {
                gene = g;
                return pi3;
            }

            genome[i][0] = g;
        }
    }

    for (int g = 0; g < generations; g++)
    {
        double childFitness = DBL_MAX;
        int parent1 = rand() % population;
        int parent2 = rand() % population;
        int child = parent1 > parent2 ? parent1 : parent2;
        Gene childGene = genome[child][0];

        if (fNumber == 2)
        {
            pi2 = ComputePi2(childGene.R, childGene.nSteps);
            childFitness = fabs(pi - pi2);

            if (childFitness <= 1.0e-15)
            {
                gene = childGene;
                return childFitness;
            }
        }

        else if (fNumber == 3)
        {
            pi3 = ComputePi3(childGene.R, childGene.nSteps);
            childFitness = fabs(pi - pi3);

            if (childFitness <= 1.0e-15)
            {
                gene = childGene;
                return childFitness;
            }
        }

        double maxFitness = fitness[0];

        for (int i = 1; i < population; i++)
            if (fitness[i] > maxFitness)
                maxFitness = fitness[i];

        std::vector<int> maxIndices;

        for (int i = 0; i < population; i++)
            if (fitness[i] == maxFitness)
                maxIndices.push_back(i);

        int selection = maxIndices[rand() % maxIndices.size()];

        fitness[selection] = childFitness;
        genome[selection][0] = childGene;
    }

    double minFitness = fitness[0];

    for (int i = 1; i < population; i++)
        if (fitness[i] < minFitness)
            minFitness = fitness[i];

    std::vector<int> minIndices;

    for (int i = 0; i < population; i++)
        if (fitness[i] == minFitness)
            minIndices.push_back(i);

    int selection = minIndices[rand() % minIndices.size()];

    gene = genome[selection][0];

    if (fNumber == 2)
        return ComputePi2(gene.R, gene.nSteps);
    else if (fNumber == 3)
        return ComputePi3(gene.R, gene.nSteps);
    else
        return 0.0;
}

int main()
{
    Gene gene = { 0 };
    std::cout << std::fixed << std::setprecision(12);
    std::cout << "Function\t\tValue of Pi\tRuntime (mu-s)" << std::endl;
    auto start0 = std::chrono::high_resolution_clock::now();
    double pi0 = 4.0 * atan(1.0);
    auto stop0 = std::chrono::high_resolution_clock::now();
    auto duration0 = std::chrono::duration_cast<std::chrono::microseconds>(stop0 - start0);
    std::cout << "Pi [4 * atan(1)]   = " << pi0 << '\t' << duration0.count() << std::endl;
    // Gauss Integral
    auto start1 = std::chrono::high_resolution_clock::now();
    double pi1 = pow(ComputePi1(6.0e2), 2.0);
    auto stop1 = std::chrono::high_resolution_clock::now();
    auto duration1 = std::chrono::duration_cast<std::chrono::microseconds>(stop1 - start1);
    std::cout << "Gauss Integral     = " << pi1 << '\t' << duration1.count() << std::endl;
    // Cauchy Integral
    auto start2 = std::chrono::high_resolution_clock::now();
    double pi2 = HillClimber(2, 200000, 3000, gene);
    auto stop2 = std::chrono::high_resolution_clock::now();
    auto duration2 = std::chrono::duration_cast<std::chrono::microseconds>(stop2 - start2);
    std::cout << "Cauchy Integral    = " << pi2 << '\t' << duration2.count() << std::endl;
    // Dirichlet integral
    auto start3 = std::chrono::high_resolution_clock::now();
    double pi3 = HillClimber(3, 50000, 1000, gene);
    auto stop3 = std::chrono::high_resolution_clock::now();
    auto duration3 = std::chrono::duration_cast<std::chrono::microseconds>(stop3 - start3);
    std::cout << "Dirichlet Integral = " << pi3 << '\t' << duration3.count() << std::endl;
	return 0;
}