Category: C++ Computer Applications

Blog Entry © Monday, March 23, 2026, by James Pate Williams, Jr. Gauss-Legendre Quadrature

Blog Entry Monday March 23 2026Download
#pragma once

/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Numerical Analysis: An Algorithmic
* Approach" by S. D. Conte and Carl de Boor
* Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1
* Translated to better C++ on Saturday, March 21,
* 2026
*/

#include <string>
#include <vector>

class Mueller
{
private:
	static double epsilon1, epsilon2;
public:
	static double f(std::string name,
		double x, int n, int m,
		std::vector<double>& zeros);
	static double g(
		std::string name,
		std::string graph,
		double x, int n, int m,
		std::vector<double>& zeros);
	static double OrthogonalPolynomialRoot(
		std::string name, 
		std::string graph, int n,
		double epsilon1, double epsilon2,
		int maxIt, int m, double x0,
		double x1, double x2, std::string& msg,
		std::vector<double>& zeros);
	static void OrthogonalPolynomialRoots(
		std::string name,
		std::string graph, int n, int nr,
		unsigned int seed, std::string& msg,
		std::vector<double>& zeros);
};


#include "Mueller.h"
#include <cmath>
#include <iostream>
#include <string>
#include <vector>

double Mueller::epsilon1 = 1.0e-12;
double Mueller::epsilon2 = 1.0e-15;

double Mueller::f(std::string name,
    double x, int n, int m,
    std::vector<double>& zeros)
{
    double a1n = 1, a2n = 1, a3n = 1, a4n = 1;
    double f0 = 1, f1 = x, f2 = 1;

    if (name == "Laguerre")
        f1 = -x + 1;

    for (int i = 2; i <= n; i++)
    {
        if (name == "Chebyshev")
        {
            a1n = 1.0;
            a2n = 0.0;
            a3n = 2.0;
            a4n = 1.0;
        }

        else if (name == "Laguerre")
        {
            a1n = i - 1 + 1;
            a2n = 2 * (i - 1) + 1;
            a3n = -1.0;
            a4n = i - 1;
        }

        else if (name == "Legendre")
        {
            a1n = i + 1 - 1;
            a2n = 0.0;
            a3n = 2.0 * (i - 1) + 1;
            a4n = i - 1;
        }

        f2 = ((a2n + a3n * x) * f1 - a4n * f0) / a1n;
        f0 = f1;
        f1 = f2;

        // deflate the function

        for (int i = 0; i < m; i++)
        {
            double z = zeros[i];
            double denom = x - z;

            if (fabs(denom) < epsilon2)
                zeros[i] = 1.001 * z;

            else
                f2 /= denom;
        }
    }

    return f2;
}

double Mueller::g(std::string name, 
    std::string graph, double x,
    int n, int m, std::vector<double>& zeros)
{
    double f0 = 1, f1 = x, f2 = 1;

    if (graph != "Extrema")
        return f(name, x, n, m, zeros);

    else
    {
        if (name == "Laguerre") {
            if (x == 1.0)
                x -= 0.0000001;

            else if (x == -1.0)
                x += 0.0000001;
        }

        double g0x = 1.0, g1x = 1.0, g2x = 1.0;

        if (name == "Chebyshev") {
            g0x = n;
            g1x = -n * x;
            g2x = 1.0 - x * x;
        }

        else if (name == "Laguerre") {
            g0x = -n;
            g1x = n;
            g2x = x;
        }

        else if (name == "Legendre")
        {
            g0x = n;
            g1x = -n * x;
            g2x = 1.0 - x * x;
        }

        if (fabs(g2x) < epsilon2)
        {
            std::cout << "g2x is too small!";
            return 0.0;
        }

        f0 = f(name, x, n - 1, 0, zeros);
        f1 = f(name, x, n, 0, zeros);
        f2 = (g1x * f1 + g0x * f0) / g2x;

        // deflate the function

        for (int i = 0; i < m; i++)
        {
            double z = zeros[i];
            double denom = x - z;

            if (fabs(denom) == 0.0)
                zeros[i] = 1.001 * z;

            else
                f2 /= denom;
        }

        return f2;
    }
}

double Mueller::OrthogonalPolynomialRoot(
    std::string name, 
    std::string graph, int n,
    double epsilon1, double epsilon2,
    int maxIt, int m, double x0,
    double x1, double x2, std::string& msg,
    std::vector<double>& zeros)
{
    double x3 = 0.0;
    int i = 2;

    while (true)
    {
        double h1 = x1 - x0;

        if (fabs(h1) < epsilon2)
        {
            msg = "h1 is too small!";
            return 0.0;
        }

        double h2 = x2 - x1;

        if (fabs(h2) < epsilon2)
        {
            msg = "h2 is too small!";
            return 0.0;
        }

        double f0 = g(name, graph, x0, n, m, zeros);
        double f1 = g(name, graph, x1, n, m, zeros);
        double f2 = g(name, graph, x2, n, m, zeros);
        double g1 = (f1 - f0) / h1;
        double g2 = (f2 - f1) / h2;

        if (fabs(f2) < epsilon2)
        {
            x3 = x2;
            break;
        }

        if (fabs(h1 + h2) < epsilon2)
        {
            msg = "h1 + h2 is too small!";
            return 0.0;
        }

        double k2 = (g2 - g1) / (h1 + h2);
        double c1 = g2 + h2 * k2;

        double d = c1 * c1 - 4.0 * f2 * k2;

        if (d < 0.0)
            d = 0.0;

        double s1 = sqrt(d);
        double m1 = c1 - s1;
        double p1 = c1 + s1;
        double denom;

        if (fabs(m1) > fabs(p1))
            denom = m1;

        else
            denom = p1;

        if (fabs(denom) < epsilon2)
        {
            msg = "denom is too small!";
            return 0.0;
        }

        double h3 = -2.0 * f2 / denom;

        x3 = x2 + h3;

        double f3 = g(name, graph, x3, n, m, zeros);

        if (fabs(f3) < epsilon2)
            break;

        if (fabs(h3) < epsilon2)
        {
            msg = "h3 is too small!";
            return 0.0;
        }

        double g3 = (f3 - f2) / h3;

        double delta = fabs(x3 - x2);

        if (delta < epsilon1)
            break;

        i++;

        if (i > maxIt)
            break;

        x0 = x1;
        x1 = x2;
        x2 = x3;

        if (fabs(x1 - x0) < epsilon1)
            break;

        if (fabs(x2 - x1) < epsilon1)
            break;
    }

    return x3;
}

void Mueller::OrthogonalPolynomialRoots(
    std::string name, 
    std::string graph, int n, int nr,
    unsigned int seed, std::string& msg,
    std::vector<double>& zeros) {
    double epsilon1 = 1.0e-10;
    epsilon2 = 1.0e-20;

    double x0 = -0.5, x1 = 0.0, x2 = 0.5;
    double z;
    int maxIt = 128;

    zeros.resize(nr);
    srand(seed);

    for (int m = 0; m < nr; m++) {
        x0 = (double)rand() / RAND_MAX;
        x1 = (double)rand() / RAND_MAX;
        x2 = (double)rand() / RAND_MAX;

        z = zeros[m] = OrthogonalPolynomialRoot(
            name, graph, n, epsilon1, epsilon2,
            maxIt, m, x0, x1, x2, msg, zeros);
    }
}


// Quadrature1d.cpp (c) Saturday, March 21, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Reference: https://dlmf.nist.gov/3.5#x

#include <iomanip>
#include <iostream>
#include <vector>
#include "Mueller.h"

static void SelectionSort(
    std::vector<double>& x,
    std::vector<double>& w,
    int n) {
    for (int i = 0; i < n - 1; i++) {
        for (int j = i + 1; j < n; j++) {
            if (x[i] > x[j]) {
                double t = x[i];
                x[i] = x[j];
                x[j] = t;
                t = w[i];
                w[i] = w[j];
                w[j] = t;
            }
        }
    }
}

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

static double GaussLegendre(
    std::vector<double> x, std::vector<double> w,
    double (*fx)(double), int n) {
    double integral = 0.0;

    for (int i = 0; i < n; i++) {
        integral += fx(x[i]) * w[i];
    }

    return integral;
}

static double Integrate(double (*fx)(double), int n,
    unsigned int seed) {
    std::string msg;
    std::vector<double> x(n);
    std::vector<double> weights(n);

    Mueller::OrthogonalPolynomialRoots(
        "Legendre", "Zeros", n, n, seed, msg, x);

    for (int i = 0; i < n; i++) {
        double xi = x[i];
        double dpn = Mueller::g("Legendre", "Extrema",
            xi, n, 0, x);
        weights[i] = 2.0 / ((1.0 - xi * xi) * dpn * dpn);
    }

    SelectionSort(x, weights, n);
    std::cout << "Abscissas (x) and Weights (w)" << std::endl;
    std::cout << "for Gauss-Legendre Quadrature";
    std::cout << std::endl;
    std::cout << "x\t\t\tw" << std::endl;

    for (int i = 0; i < n; i++) {
        std::cout << std::fixed << std::setprecision(15);
        std::cout << std::setw(18);
        std::cout << x[i] << '\t';
        std::cout << std::fixed << std::setprecision(15);
        std::cout << std::setw(18);
        std::cout << weights[i] << std::endl;
    }

    return GaussLegendre(x, weights, fx, n);
}

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

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

    x = a + h2;

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

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

int main() {
    unsigned int seed = 1;
    
    for (int n = 10; n <= 40; n += 10)
    {
        double integral = Integrate(fx, n, seed++);

        std::cout << std::fixed << std::setprecision(15);
        std::cout << std::setw(18);
        std::cout << n << '\t' << integral << std::endl;
    }

    for (int n = 50000; n <= 100000; n += 10000)
    {
        double integral = SimpsonsRule(n, -1, 1, fx);

        std::cout << "Simpson ";
        std::cout << n << "\t\t";
        std::cout << std::fixed << std::setprecision(15);
        std::cout << std::setw(18); 
        std::cout << integral << std::endl;
    }

    return 0;
}

Blog Input © Saturday, March 21, 2026, by James Pate Williams, Jr., Three Model Universes

Blog Entry Saturday March 21 2026Download
#pragma once

class UniversalModels
{
public:
	static void deSitterUniverse(
		double epsilon, double lambda, double ct,
		double A, double B, double& K);
	static void RadiationUniverse(
		double A, double epsilon,
		double deltaT, double& K);
	static void FriedmannUniverse(
		int epsilon, double M,
		double cT, double& K);
};


#include "UniversalModels.h"
#include <math.h>

void UniversalModels::deSitterUniverse(
	double epsilon, double lambda,
	double ct, double A, double B, double& K)
{
	if (lambda > 0)
	{
		if (epsilon == +1)
		{
			K = cosh(B * ct) / B;
		}

		else if (epsilon == -1)
		{
			K = sinh(B * ct) / B;
		}

		else
		{
			K = A * exp(B * ct);
		}
	}

	else if (lambda < 0)
	{
		if (epsilon == -1)
		{
			lambda = -3.0 * B * B;
			K = cos(B * ct) / B;
		}
	}
}

void UniversalModels::RadiationUniverse(
	double A, double epsilon,
	double deltaT, double& K)
{
	double c = 1.0, kappa = 1.0;

	if (epsilon == 0)
	{
		K = sqrt(2.0 * c * sqrt(kappa * A / 3.0) * deltaT);
	}

	else if (epsilon == -1)
	{
		K = sqrt(c * c * pow(deltaT, 2.0) + 2 * c *
			sqrt(kappa * A / 3.0) * deltaT);
	}

	else if (epsilon == +1)
	{
		K = sqrt(-c * c * pow(deltaT, 2.0) + 2 * c *
			sqrt(kappa * A / 3.0) * deltaT);
	}
}

void UniversalModels::FriedmannUniverse(
	int epsilon, double M, double cT, double& K)
{
	if (epsilon == 0)
	{
		K = M * cT * cT / 12.0;
	}

	else if (epsilon == -1)
	{
		K = M * (cosh(cT) - 1.0) / 6.0;
	}

	else if (epsilon == +1)
	{
		K = M * (1.0 - cos(cT)) / 6.0;
	}
}


// ModelUniverses.cpp (c) Thursday, March 19, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Reference: "General Relativity An Introduction
// to the theory of the gravitational field "
// (c) 1982 by Hans Stephani translated by Martin
// Pollack and John Stewart

#include "framework.h"
#include "ModelUniverses.h"
#include "UniversalModels.h"
#include <math.h>
#include <string>
#include <vector>

#define MAX_LOADSTRING 100

typedef struct tagPoint2dDeSitter
{
    double ct, K;
} Point2dDeSitter, * PPoint2dDeSitter;

typedef struct tagPoint2dRadiation
{
    double deltaT, K;
} Point2dRadiation, * PPoint2dRadiation;

typedef struct tagPoint2dFriedmann
{
    double cT, K;
} Point2dFriedmann, * PPoint2dFriedmann;


// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
bool deSitter, radiation, friedmann;            // universal model type
std::wstring fTitle, xTitle, yTitle;            // graph titles
std::vector<Point2dDeSitter> pointsDeSitter;    // graph points
std::vector<Point2dRadiation> pointsRadiation;  // graph points
std::vector<Point2dFriedmann> pointsFriedmann;  // graph points

// Forward declarations of functions included in this code module:
ATOM                MyRegisterClass(HINSTANCE hInstance);
BOOL                InitInstance(HINSTANCE, int);
LRESULT CALLBACK    WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    About(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    DataInputDialogDeSitter(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    DrawGraphDialogDeSitter(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    DataInputDialogRadiation(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    DrawGraphDialogRadiation(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    DataInputDialogFriedmann(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    DrawGraphDialogFriedmann(HWND, UINT, WPARAM, LPARAM);

int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    UNREFERENCED_PARAMETER(hPrevInstance);
    UNREFERENCED_PARAMETER(lpCmdLine);

    // TODO: Place code here.

    // Initialize global strings
    LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
    LoadStringW(hInstance, IDC_MODELUNIVERSES, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_MODELUNIVERSES));

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, nullptr, 0, 0))
    {
        if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
        {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }

    return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
    WNDCLASSEXW wcex = { 0 };

    wcex.cbSize = sizeof(WNDCLASSEX);

    wcex.style          = CS_HREDRAW | CS_VREDRAW;
    wcex.lpfnWndProc    = WndProc;
    wcex.cbClsExtra     = 0;
    wcex.cbWndExtra     = 0;
    wcex.hInstance      = hInstance;
    wcex.hIcon          = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_MODELUNIVERSES));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_MODELUNIVERSES);
    wcex.lpszClassName  = szWindowClass;
    wcex.hIconSm        = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

    return RegisterClassExW(&wcex);
}

//
//   FUNCTION: InitInstance(HINSTANCE, int)
//
//   PURPOSE: Saves instance handle and creates main window
//
//   COMMENTS:
//
//        In this function, we save the instance handle in a global variable and
//        create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   hInst = hInstance; // Store instance handle in our global variable

   HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    switch (message)
    {
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);
            // Parse the menu selections:
            switch (wmId)
            {
            case IDM_ABOUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
                break;
            case IDM_DE_SITTER:
                deSitter = true;
                fTitle = L"de Sitter Universe";
                xTitle = L"Bct";
                yTitle = L"K";
                DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG_DE_SITTER), 
                    hWnd, DataInputDialogDeSitter);
                break;
            case IDM_RADIATION:
                radiation = true;
                fTitle = L"Radiation Universe";
                xTitle = L"t - t0";
                yTitle = L"K";
                DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG_RADIATION), hWnd,
                    DataInputDialogRadiation);
                break;
            case IDM_FRIEDMANN:
                friedmann = true;
                fTitle = L"Friedmann Universe";
                xTitle = L"cT";
                yTitle = L"K";
                DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG_FRIEDMANN), hWnd,
                    DataInputDialogFriedmann);
                break;
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
        {
            PAINTSTRUCT ps;
            HDC hdc = BeginPaint(hWnd, &ps);
            // TODO: Add any drawing code that uses hdc here...
            EndPaint(hWnd, &ps);
        }
        break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

INT_PTR CALLBACK DataInputDialogDeSitter(
    HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    static WCHAR data[128] = L"";
    static HWND hCombo_A = GetDlgItem(hDlg, IDC_COMBO_A_1);
    static HWND hCombo_B = GetDlgItem(hDlg, IDC_COMBO_B);
    static HWND hCombo_ct = GetDlgItem(hDlg, IDC_COMBO_CT);
    static HWND hCombo_epsilon = GetDlgItem(hDlg, IDC_COMBO_EPSILON_1);
    static HWND hCombo_lambda = GetDlgItem(hDlg, IDC_COMBO_LAMBDA);

    switch (message)
    {
    case WM_INITDIALOG:
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"1");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"2");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"3");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"4");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"5");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"6");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"7");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"8");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"9");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"10");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_SETCURSEL, 0, 0);
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"1");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"2");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"3");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"4");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"5");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"6");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"7");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"8");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"9");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"10");
        SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_SETCURSEL, 0, 0);
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"1");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"2");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"3");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"4");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"5");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"6");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"7");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"8");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"9");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"10");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_SETCURSEL, 0, 0);
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_1, CB_ADDSTRING, 0, (LPARAM)L"-1");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_1, CB_ADDSTRING, 0, (LPARAM)L"0");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_1, CB_ADDSTRING, 0, (LPARAM)L"+1");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_1, CB_SETCURSEL, 0, 0);
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-5");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-4");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-3");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-2");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-1");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"1");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"2");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"3");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"4");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"5");
        SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_SETCURSEL, 0, 0);
        return (INT_PTR)TRUE;
    case WM_COMMAND:
    {
        if (LOWORD(wParam) == IDC_BUTTON_DRAW_DE_SITTER)
        {
            char line[128] = "";
            GetDlgItemTextA(hDlg, IDC_COMBO_A_1, line, 128);
            double A = atof(line);
            GetDlgItemTextA(hDlg, IDC_COMBO_B, line, 128);
            double B = atof(line);
            GetDlgItemTextA(hDlg, IDC_COMBO_CT, line, 128);
            double CT = atof(line);
            int epsilon = GetDlgItemInt(hDlg, IDC_COMBO_EPSILON_1, FALSE, TRUE);
            GetDlgItemTextA(hDlg, IDC_COMBO_LAMBDA, line, 128);
            double lambda = atof(line);
			double ct = 0.0, K = 0.0;
			double delta_ct = CT / 256.0;
			int index = 0;

            pointsDeSitter.clear();

			while (ct <= CT)
			{
                UniversalModels::deSitterUniverse(
                    epsilon, lambda, ct, A, B, K);
				Point2dDeSitter pt = { 0 };
				pt.ct = ct;
				pt.K = K;
				ct += delta_ct;
                pointsDeSitter.push_back(pt);
			}

			DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG_DE_SITTER),
				hDlg, DrawGraphDialogDeSitter);
            return (INT_PTR)TRUE;
        }

        else if (LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
    }
    }

    return (INT_PTR)FALSE;
}

static void FindMinMaxDeSitter(
    double& xmin,
    double& xmax,
    double& ymin,
    double& ymax)
{
    xmin = pointsDeSitter[0].ct;
    xmax = pointsDeSitter[0].ct;
    ymin = pointsDeSitter[0].K;
    ymax = pointsDeSitter[0].K;

    for (int i = 1; i < pointsDeSitter.size(); i++)
    {
        Point2dDeSitter pt = pointsDeSitter[i];

        if (pt.ct < xmin)
            xmin = pt.ct;
        if (pt.ct > xmax)
            xmax = pt.ct;
        if (pt.K < ymin)
            ymin = pt.K;
        if (pt.K > ymax)
            ymax = pt.K;
    }
}

static void DrawTitles(
    HDC hdc, RECT clientRect,
    const std::wstring& fTitle,
    const std::wstring& xTitle,
    const std::wstring& yTitle,
    int sx0, int sx1, int sy0, int sy1)
{
    HFONT hCustomFont = CreateFont(
        -MulDiv(10, GetDeviceCaps(hdc, LOGPIXELSY), 72), // Height in logical units
        0,                      // Width (0 = default)
        0,                      // Escapement
        0,                      // Orientation
        FW_BOLD,                // Weight (700 = bold)
        FALSE,                  // Italic
        FALSE,                  // Underline
        FALSE,                  // StrikeOut
        DEFAULT_CHARSET,        // Charset
        OUT_DEFAULT_PRECIS,     // Output precision
        CLIP_DEFAULT_PRECIS,    // Clipping precision
        DEFAULT_QUALITY,        // Quality
        FIXED_PITCH | FF_MODERN,// Pitch and family
        L"Courier New"          // Typeface name
    );

    HFONT hOldFont = (HFONT)SelectObject(hdc, hCustomFont);

    SIZE sz;
    int width = clientRect.right - clientRect.left;

    // Title: fTitle
    GetTextExtentPoint32W(hdc, fTitle.c_str(), (int)fTitle.length(), &sz);
    int w = sz.cx;
    int h = sz.cy;
    TextOutW(hdc, (width - w) / 2, h, fTitle.c_str(), (int)fTitle.length());
    // draw x-axis title
    GetTextExtentPoint32W(hdc, xTitle.c_str(), (int)xTitle.length(), &sz);
    w = sz.cx;
    TextOutW(hdc, sx0 + (sx1 - sx0 - w) / 2, sy1 + 2 * h, 
        xTitle.c_str(), (int)xTitle.length());
    // draw y-axis title
    GetTextExtentPoint32W(hdc, yTitle.c_str(), (int)yTitle.length(), &sz);
    w = sz.cx;
    TextOutW(hdc, sx1 + w / 5, sy0 + (sy1 - sy0) / 2 - h / 2, 
        yTitle.c_str(), (int)yTitle.length());
    SelectObject(hdc, hOldFont); // Restore original font
}

static void DrawFormattedText(HDC hdc, char loctext[], RECT rect)
{
    // Draw the text with formatting options
    DrawTextA(hdc, loctext, (int)strlen(loctext), &rect, DT_SINGLELINE | DT_NOCLIP);
}

static INT_PTR CALLBACK DrawGraphDialogDeSitter(
    HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    char line[256] = "";
    switch (message)
    {
    case WM_INITDIALOG:
        strcpy_s(line, 256, "Drawing Dialog deSitter Universe");
        SetWindowTextA(hDlg, line);
        return (INT_PTR)TRUE;
    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    case WM_PAINT:
        double xmax = 0, xmin = 0;
        double ymax = 0, ymin = 0;
        FindMinMaxDeSitter(xmin, xmax, ymin, ymax);
        double h = 0, pi = 0, plm = 0, theta = 0;
        float xSpan = (float)(xmax - xmin);
        float ySpan = (float)(ymax - ymin);
        RECT rect = { };
        GetClientRect(hDlg, &rect);
        float width = (float)(rect.right - rect.left + 1);
        float height = (float)(rect.bottom - rect.top - 32 + 1);
        float sx0 = 2.0f * width / 16.0f;
        float sx1 = 14.0f * width / 16.0f;
        float sy0 = 2.0f * height / 16.0f;
        float sy1 = 14.0f * height / 16.0f;
        float deltaX = xSpan / 8.0f;
        float deltaY = ySpan / 8.0f;
        float xSlope = (sx1 - sx0) / xSpan;
        float xInter = (float)(sx0 - xSlope * xmin);
        float ySlope = (sy0 - sy1) / ySpan;
        float yInter = (float)(sy0 - ySlope * ymax);
        float px = 0, py = 0, sx = 0, sy = 0;
        PAINTSTRUCT ps;
        POINT wPt = { };
        HDC hdc = BeginPaint(hDlg, &ps);
        int i = 0;
        float x = (float)xmin;
        float y = (float)ymax;
        DrawTitles(hdc, rect, fTitle, xTitle, yTitle,
            (int)sx0, (int)sx1, (int)sy0, (int)sy1);
        px = (float)pointsDeSitter[0].ct;
        py = (float)pointsDeSitter[0].K;
        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        MoveToEx(hdc, (int)sx, (int)sy0, &wPt);

        while (i <= 8)
        {
            sx = xSlope * x + xInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
            LineTo(hdc, (int)sx, (int)sy1);
            char numberStr[128] = "";
            sprintf_s(numberStr, 128, "%5.4lf", x);
            SIZE size = { };
            GetTextExtentPoint32A(
                hdc,
                numberStr,
                (int)strlen(numberStr),
                &size);
            RECT textRect = { };
            textRect.left = (long)(sx - size.cx / 2.0f);
            textRect.right = (long)(sx + size.cx / 2.0f);
            textRect.top = (long)sy1;
            textRect.bottom = (long)(sy1 + size.cy / 2.0f);
            DrawFormattedText(hdc, numberStr, textRect);
            x += deltaX;
            i++;
        }

        i = 0;
        y = (float)ymin;

        while (i <= 8)
        {
            sy = ySlope * y + yInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
            LineTo(hdc, (int)sx, (int)sy);

            if (i != 0)
            {
                char numberStr[128] = "";
                sprintf_s(numberStr, 128, "%+5.3lf", y);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    numberStr,
                    (int)strlen(numberStr),
                    &size);
                RECT textRect = { 0 };
                textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
                textRect.right = (long)(sx0 - size.cx / 2.0f);
                textRect.top = (long)(sy - size.cy / 2.0f);
                textRect.bottom = (long)(sy + size.cy / 2.0f);
                DrawFormattedText(hdc, numberStr, textRect);
            }

            y += deltaY;
            i++;
        }

        HGDIOBJ bPenNew = NULL;
        HGDIOBJ hPenOld = NULL;

        bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
        hPenOld = SelectObject(hdc, bPenNew);

        HRGN clipRegion = CreateRectRgn(
            (int)sx0, (int)sy0,              // Left, Top
            (int)(sx1), (int)(sy1)           // Right, Bottom
        );

        // Apply clipping region

        SelectClipRgn(hdc, clipRegion);

        px = (float)pointsDeSitter[0].ct;
        py = (float)pointsDeSitter[0].K;
        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        wPt.x = wPt.y = 0;
        MoveToEx(hdc, (int)sx, (int)sy, &wPt);

        for (size_t j = 1; j < pointsDeSitter.size(); j++)
        {
            px = (float)pointsDeSitter[j].ct;
            py = (float)pointsDeSitter[j].K;
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            LineTo(hdc, (int)sx, (int)sy);
        }

        DeleteObject(bPenNew);
        bPenNew = NULL;

        SelectObject(hdc, hPenOld);
        DeleteObject(clipRegion);
        EndPaint(hDlg, &ps);
        return (INT_PTR)TRUE;
    }

    return (INT_PTR)FALSE;
}

INT_PTR CALLBACK DataInputDialogRadiation(
    HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    static WCHAR data[128] = L"";
    static HWND hCombo_B = GetDlgItem(hDlg, IDC_COMBO_B);
    static HWND hCombo_ct = GetDlgItem(hDlg, IDC_COMBO_DELTA_T);
    static HWND hCombo_epsilon = GetDlgItem(hDlg, IDC_COMBO_EPSILON_2);
    static HWND hCombo_lambda = GetDlgItem(hDlg, IDC_COMBO_LAMBDA);

    switch (message)
    {
    case WM_INITDIALOG:
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e1");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e2");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e3");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e4");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e5");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"6");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"7");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"8");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"9");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"10");
        SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_SETCURSEL, 0, 0);
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-5");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-4");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-3");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-2");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-1");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+0");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+1");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+2");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+3");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+4");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+5");
        SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_SETCURSEL, 0, 0);
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_2, CB_ADDSTRING, 0, (LPARAM)L"-1");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_2, CB_ADDSTRING, 0, (LPARAM)L"0");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_2, CB_ADDSTRING, 0, (LPARAM)L"+1");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_2, CB_SETCURSEL, 0, 0);
        return (INT_PTR)TRUE;
    case WM_COMMAND:
    {
        if (LOWORD(wParam) == IDC_BUTTON_DRAW_RADIATION)
        {
            char line[128] = "";
            GetDlgItemTextA(hDlg, IDC_COMBO_A_2, line, 128);
            double A = atof(line);
            GetDlgItemTextA(hDlg, IDC_COMBO_DELTA_T, line, 128);
            double deltaT = atof(line);
            int epsilon = GetDlgItemInt(hDlg, IDC_COMBO_EPSILON_2, FALSE, TRUE);
            double step = deltaT / 256.0, t = 0.0, K = 0.0;
            int index = 0;

            pointsRadiation.clear();

            while (t <= deltaT)
            {
                UniversalModels::RadiationUniverse(
                    A, epsilon, t, K);
                Point2dRadiation pt = { 0 };
                pt.deltaT = t;
                pt.K = K;
                t += step;
                pointsRadiation.push_back(pt);
            }

            DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG_RADIATION),
                hDlg, DrawGraphDialogRadiation);
            return (INT_PTR)TRUE;
        }

        else if (LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
    }
    }

    return (INT_PTR)FALSE;
}

static void FindMinMaxRadiation(
    double& xmin,
    double& xmax,
    double& ymin,
    double& ymax)
{
    xmin = pointsRadiation[0].deltaT;
    xmax = pointsRadiation[0].deltaT;
    ymin = pointsRadiation[0].K;
    ymax = pointsRadiation[0].K;

    for (int i = 1; i < pointsRadiation.size(); i++)
    {
        Point2dRadiation pt = pointsRadiation[i];

        if (pt.deltaT < xmin)
            xmin = pt.deltaT;
        if (pt.deltaT > xmax)
            xmax = pt.deltaT;
        if (pt.K < ymin)
            ymin = pt.K;
        if (pt.K > ymax)
            ymax = pt.K;
    }
}

static INT_PTR CALLBACK DrawGraphDialogRadiation(
    HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    char line[256] = "";
    switch (message)
    {
    case WM_INITDIALOG:
        strcpy_s(line, 256, "Drawing Dialog Radiation Universe");
        SetWindowTextA(hDlg, line);
        return (INT_PTR)TRUE;
    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    case WM_PAINT:
        double xmax = 0, xmin = 0;
        double ymax = 0, ymin = 0;
        FindMinMaxRadiation(xmin, xmax, ymin, ymax);
        double h = 0, pi = 0, plm = 0, theta = 0;
        float xSpan = (float)(xmax - xmin);
        float ySpan = (float)(ymax - ymin);
        RECT rect = { };
        GetClientRect(hDlg, &rect);
        float width = (float)(rect.right - rect.left + 1);
        float height = (float)(rect.bottom - rect.top - 32 + 1);
        float sx0 = 2.0f * width / 16.0f;
        float sx1 = 14.0f * width / 16.0f;
        float sy0 = 2.0f * height / 16.0f;
        float sy1 = 14.0f * height / 16.0f;
        float deltaX = xSpan / 8.0f;
        float deltaY = ySpan / 8.0f;
        float xSlope = (sx1 - sx0) / xSpan;
        float xInter = (float)(sx0 - xSlope * xmin);
        float ySlope = (sy0 - sy1) / ySpan;
        float yInter = (float)(sy0 - ySlope * ymax);
        float px = 0, py = 0, sx = 0, sy = 0;
        PAINTSTRUCT ps;
        POINT wPt = { };
        HDC hdc = BeginPaint(hDlg, &ps);
        int i = 0;
        float x = (float)xmin;
        float y = (float)ymax;
        DrawTitles(hdc, rect, fTitle, xTitle, yTitle,
            (int)sx0, (int)sx1, (int)sy0, (int)sy1);
        px = (float)pointsRadiation[0].deltaT;
        py = (float)pointsRadiation[0].K;
        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        MoveToEx(hdc, (int)sx, (int)sy0, &wPt);

        while (i <= 8)
        {
            sx = xSlope * x + xInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
            LineTo(hdc, (int)sx, (int)sy1);
            char numberStr[128] = "";
            sprintf_s(numberStr, 128, "%5.4lf", x);
            SIZE size = { };
            GetTextExtentPoint32A(
                hdc,
                numberStr,
                (int)strlen(numberStr),
                &size);
            RECT textRect = { };
            textRect.left = (long)(sx - size.cx / 2.0f);
            textRect.right = (long)(sx + size.cx / 2.0f);
            textRect.top = (long)sy1;
            textRect.bottom = (long)(sy1 + size.cy / 2.0f);
            DrawFormattedText(hdc, numberStr, textRect);
            x += deltaX;
            i++;
        }

        i = 0;
        y = (float)ymin;

        while (i <= 8)
        {
            sy = ySlope * y + yInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
            LineTo(hdc, (int)sx, (int)sy);

            if (i != 0)
            {
                char numberStr[128] = "";
                sprintf_s(numberStr, 128, "%+5.3lf", y);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    numberStr,
                    (int)strlen(numberStr),
                    &size);
                RECT textRect = { 0 };
                textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
                textRect.right = (long)(sx0 - size.cx / 2.0f);
                textRect.top = (long)(sy - size.cy / 2.0f);
                textRect.bottom = (long)(sy + size.cy / 2.0f);
                DrawFormattedText(hdc, numberStr, textRect);
            }

            y += deltaY;
            i++;
        }

        HGDIOBJ bPenNew = NULL;
        HGDIOBJ hPenOld = NULL;

        bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
        hPenOld = SelectObject(hdc, bPenNew);

        HRGN clipRegion = CreateRectRgn(
            (int)sx0, (int)sy0,              // Left, Top
            (int)(sx1), (int)(sy1)           // Right, Bottom
        );

        // Apply clipping region

        SelectClipRgn(hdc, clipRegion);

        px = (float)pointsRadiation[0].deltaT;
        py = (float)pointsRadiation[0].K;
        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        wPt.x = wPt.y = 0;
        MoveToEx(hdc, (int)sx, (int)sy, &wPt);

        for (size_t j = 1; j < pointsRadiation.size(); j++)
        {
            px = (float)pointsRadiation[j].deltaT;
            py = (float)pointsRadiation[j].K;
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            LineTo(hdc, (int)sx, (int)sy);
        }

        DeleteObject(bPenNew);
        bPenNew = NULL;

        SelectObject(hdc, hPenOld);
        DeleteObject(clipRegion);
        EndPaint(hDlg, &ps);
        return (INT_PTR)TRUE;
    }

    return (INT_PTR)FALSE;
}

INT_PTR CALLBACK DataInputDialogFriedmann(
    HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    static WCHAR data[128] = L"";
    static HWND hCombo_M = GetDlgItem(hDlg, IDC_COMBO_M);
    static HWND hCombo_CT1 = GetDlgItem(hDlg, IDC_COMBO_CT1);
    static HWND hCombo_epsilon = GetDlgItem(hDlg, IDC_COMBO_EPSILON_3);

    switch (message)
    {
    case WM_INITDIALOG:
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"1");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"2");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"3");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"4");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"5");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"6");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"7");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"8");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"9");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"10");
        SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_SETCURSEL, 0, 0);
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"1");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"2");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"3");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"4");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"5");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"6");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"7");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"8");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"9");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"10");
        SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_SETCURSEL, 0, 0);
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_3, CB_ADDSTRING, 0, (LPARAM)L"-1");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_3, CB_ADDSTRING, 0, (LPARAM)L"0");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_3, CB_ADDSTRING, 0, (LPARAM)L"+1");
        SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_3, CB_SETCURSEL, 0, 0);
        return (INT_PTR)TRUE;
    case WM_COMMAND:
    {
        if (LOWORD(wParam) == IDC_BUTTON_DRAW_FRIEDMANN)
        {
            char line[128] = "";
            GetDlgItemTextA(hDlg, IDC_COMBO_M, line, 128);
            double M = atof(line);
            GetDlgItemTextA(hDlg, IDC_COMBO_CT1, line, 128);
            double cT = atof(line);
            int epsilon = GetDlgItemInt(hDlg, IDC_COMBO_EPSILON_3, FALSE, TRUE);
            double step = cT / 256.0, t = 0.0, K = 0.0;
            int index = 0;

            pointsFriedmann.clear();

            while (t <= cT)
            {
                UniversalModels::FriedmannUniverse(
                    epsilon, M, t, K);
                Point2dFriedmann pt = { 0 };
                pt.cT = t;
                pt.K = K;
                t += step;
                pointsFriedmann.push_back(pt);
            }

            DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG_FRIEDMANN),
                hDlg, DrawGraphDialogFriedmann);
            return (INT_PTR)TRUE;
        }

        else if (LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
    }
    }

    return (INT_PTR)FALSE;
}

static void FindMinMaxFriedmann(
    double& xmin,
    double& xmax,
    double& ymin,
    double& ymax)
{
    xmin = pointsFriedmann[0].cT;
    xmax = pointsFriedmann[0].cT;
    ymin = pointsFriedmann[0].K;
    ymax = pointsFriedmann[0].K;

    for (int i = 1; i < pointsFriedmann.size(); i++)
    {
        Point2dFriedmann pt = pointsFriedmann[i];

        if (pt.cT < xmin)
            xmin = pt.cT;
        if (pt.cT > xmax)
            xmax = pt.cT;
        if (pt.K < ymin)
            ymin = pt.K;
        if (pt.K > ymax)
            ymax = pt.K;
    }
}

static INT_PTR CALLBACK DrawGraphDialogFriedmann(
    HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    char line[256] = "";
    switch (message)
    {
    case WM_INITDIALOG:
        strcpy_s(line, 256, "Drawing Dialog Friedmann Universe");
        SetWindowTextA(hDlg, line);
        return (INT_PTR)TRUE;
    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    case WM_PAINT:
        double xmax = 0, xmin = 0;
        double ymax = 0, ymin = 0;
        FindMinMaxFriedmann(xmin, xmax, ymin, ymax);
        double h = 0, pi = 0, plm = 0, theta = 0;
        float xSpan = (float)(xmax - xmin);
        float ySpan = (float)(ymax - ymin);
        RECT rect = { };
        GetClientRect(hDlg, &rect);
        float width = (float)(rect.right - rect.left + 1);
        float height = (float)(rect.bottom - rect.top - 32 + 1);
        float sx0 = 2.0f * width / 16.0f;
        float sx1 = 14.0f * width / 16.0f;
        float sy0 = 2.0f * height / 16.0f;
        float sy1 = 14.0f * height / 16.0f;
        float deltaX = xSpan / 8.0f;
        float deltaY = ySpan / 8.0f;
        float xSlope = (sx1 - sx0) / xSpan;
        float xInter = (float)(sx0 - xSlope * xmin);
        float ySlope = (sy0 - sy1) / ySpan;
        float yInter = (float)(sy0 - ySlope * ymax);
        float px = 0, py = 0, sx = 0, sy = 0;
        PAINTSTRUCT ps;
        POINT wPt = { };
        HDC hdc = BeginPaint(hDlg, &ps);
        int i = 0;
        float x = (float)xmin;
        float y = (float)ymax;
        DrawTitles(hdc, rect, fTitle, xTitle, yTitle,
            (int)sx0, (int)sx1, (int)sy0, (int)sy1);
        px = (float)pointsFriedmann[0].cT;
        py = (float)pointsFriedmann[0].K;
        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        MoveToEx(hdc, (int)sx, (int)sy0, &wPt);

        while (i <= 8)
        {
            sx = xSlope * x + xInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
            LineTo(hdc, (int)sx, (int)sy1);
            char numberStr[128] = "";
            sprintf_s(numberStr, 128, "%5.4lf", x);
            SIZE size = { };
            GetTextExtentPoint32A(
                hdc,
                numberStr,
                (int)strlen(numberStr),
                &size);
            RECT textRect = { };
            textRect.left = (long)(sx - size.cx / 2.0f);
            textRect.right = (long)(sx + size.cx / 2.0f);
            textRect.top = (long)sy1;
            textRect.bottom = (long)(sy1 + size.cy / 2.0f);
            DrawFormattedText(hdc, numberStr, textRect);
            x += deltaX;
            i++;
        }

        i = 0;
        y = (float)ymin;

        while (i <= 8)
        {
            sy = ySlope * y + yInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
            LineTo(hdc, (int)sx, (int)sy);

            if (i != 0)
            {
                char numberStr[128] = "";
                sprintf_s(numberStr, 128, "%+5.3lf", y);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    numberStr,
                    (int)strlen(numberStr),
                    &size);
                RECT textRect = { 0 };
                textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
                textRect.right = (long)(sx0 - size.cx / 2.0f);
                textRect.top = (long)(sy - size.cy / 2.0f);
                textRect.bottom = (long)(sy + size.cy / 2.0f);
                DrawFormattedText(hdc, numberStr, textRect);
            }

            y += deltaY;
            i++;
        }

        HGDIOBJ bPenNew = NULL;
        HGDIOBJ hPenOld = NULL;

        bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
        hPenOld = SelectObject(hdc, bPenNew);

        HRGN clipRegion = CreateRectRgn(
            (int)sx0, (int)sy0,              // Left, Top
            (int)(sx1), (int)(sy1)           // Right, Bottom
        );

        // Apply clipping region

        SelectClipRgn(hdc, clipRegion);

        px = (float)pointsFriedmann[0].cT;
        py = (float)pointsFriedmann[0].K;
        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        wPt.x = wPt.y = 0;
        MoveToEx(hdc, (int)sx, (int)sy, &wPt);

        for (size_t j = 1; j < pointsFriedmann.size(); j++)
        {
            px = (float)pointsFriedmann[j].cT;
            py = (float)pointsFriedmann[j].K;
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            LineTo(hdc, (int)sx, (int)sy);
        }

        DeleteObject(bPenNew);
        bPenNew = NULL;

        SelectObject(hdc, hPenOld);
        DeleteObject(clipRegion);
        EndPaint(hDlg, &ps);
        return (INT_PTR)TRUE;
    }

    return (INT_PTR)FALSE;
}


//Microsoft Visual C++ generated resource script.
//
#include "resource.h"

#define APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
#ifndef APSTUDIO_INVOKED
#include "targetver.h"
#endif
#define APSTUDIO_HIDDEN_SYMBOLS
#include "windows.h"
#undef APSTUDIO_HIDDEN_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
#undef APSTUDIO_READONLY_SYMBOLS

#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU)
LANGUAGE 9, 1

/////////////////////////////////////////////////////////////////////////////
//
// Icon
//

// Icon with lowest ID value placed first to ensure application icon
// remains consistent on all systems.

IDI_MODELUNIVERSES       ICON         "ModelUniverses.ico"
IDI_SMALL               ICON         "small.ico"

/////////////////////////////////////////////////////////////////////////////
//
// Menu
//

IDC_MODELUNIVERSES MENU
BEGIN
    POPUP "&Begin"
    BEGIN
        MENUITEM "&de Sitter Universe", IDM_DE_SITTER
        MENUITEM "&Radiation Universe", IDM_RADIATION
        MENUITEM "&Friedmann Universe", IDM_FRIEDMANN
        MENUITEM SEPARATOR
        MENUITEM "E&xit",               IDM_EXIT
    END
    POPUP "&Help"
    BEGIN
        MENUITEM "&About ...",          IDM_ABOUT
    END
END


/////////////////////////////////////////////////////////////////////////////
//
// Accelerator
//

IDC_MODELUNIVERSES ACCELERATORS
BEGIN
    "?",            IDM_ABOUT,              ASCII,  ALT
    "/",            IDM_ABOUT,              ASCII,  ALT
END


/////////////////////////////////////////////////////////////////////////////
//
// Dialog
//

IDD_ABOUTBOX DIALOGEX 0, 0, 170, 62
STYLE DS_SETFONT | DS_MODALFRAME | DS_FIXEDSYS | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "About ModelUniverses"
FONT 8, "MS Shell Dlg"
BEGIN
    ICON            IDI_MODELUNIVERSES,IDC_STATIC,14,14,21,20
    LTEXT           "ModelUniverses, Version 1.0",IDC_STATIC,42,14,114,8,SS_NOPREFIX
    LTEXT           "Copyright (c) 2026",IDC_STATIC,42,26,114,8
    DEFPUSHBUTTON   "OK",IDOK,113,41,50,14,WS_GROUP
END

IDD_DATA_INPUT_DIALOG_DE_SITTER DIALOGEX 0, 0, 280, 200
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog de Sitter Universe"
FONT 10, "Courier New", 700
BEGIN
LTEXT       "A:", IDC_STATIC, 10, 10, 100, 10
COMBOBOX    IDC_COMBO_A_1, 120, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT       "B:", IDC_STATIC, 10, 40, 100, 10
COMBOBOX    IDC_COMBO_B, 120, 40, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT       "ct:", IDC_STATIC, 10, 70, 100, 10
COMBOBOX    IDC_COMBO_CT, 120, 70, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT       "epsilon:",
IDC_STATIC, 10, 100, 50, 60
COMBOBOX    IDC_COMBO_EPSILON_1, 120, 100, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT       "Lambda:",
IDC_STATIC, 10, 130, 50, 60
COMBOBOX    IDC_COMBO_LAMBDA, 120, 130, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
PUSHBUTTON  "Draw", IDC_BUTTON_DRAW_DE_SITTER, 10, 160, 50, 16
PUSHBUTTON  "Cancel", IDCANCEL, 210, 160, 50, 16
END

IDD_DRAW_GRAPH_DIALOG_DE_SITTER DIALOGEX 0, 0, 600, 310
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Draw Graph Dialog de Sitter Universe"
FONT 10, "Courier New", 700
BEGIN
END

IDD_DATA_INPUT_DIALOG_RADIATION DIALOGEX 0, 0, 280, 170
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog Radiation Universe"
FONT 10, "Courier New", 700
BEGIN
LTEXT       "A:", IDC_STATIC, 10, 10, 100, 10
COMBOBOX    IDC_COMBO_A_2, 120, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT       "t - t0:", IDC_STATIC, 10, 40, 100, 10
COMBOBOX    IDC_COMBO_DELTA_T, 120, 40, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT       "epsilon:",
IDC_STATIC, 10, 70, 50, 60
COMBOBOX    IDC_COMBO_EPSILON_2, 120, 70, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
PUSHBUTTON  "Draw", IDC_BUTTON_DRAW_RADIATION, 10, 100, 50, 16
PUSHBUTTON  "Cancel", IDCANCEL, 210, 100, 50, 16
END

IDD_DRAW_GRAPH_DIALOG_RADIATION DIALOGEX 0, 0, 600, 310
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Draw Graph Dialog Radiation Universe"
FONT 10, "Courier New", 700
BEGIN
END

IDD_DATA_INPUT_DIALOG_FRIEDMANN DIALOGEX 0, 0, 280, 170
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog Friedmann Universe"
FONT 10, "Courier New", 700
BEGIN
LTEXT       "M:", IDC_STATIC, 10, 10, 100, 10
COMBOBOX    IDC_COMBO_M, 120, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT       "cT:", IDC_STATIC, 10, 40, 100, 10
COMBOBOX    IDC_COMBO_CT1, 120, 40, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT       "epsilon:",
IDC_STATIC, 10, 70, 50, 60
COMBOBOX    IDC_COMBO_EPSILON_3, 120, 70, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
PUSHBUTTON  "Draw", IDC_BUTTON_DRAW_FRIEDMANN, 10, 100, 50, 16
PUSHBUTTON  "Cancel", IDCANCEL, 210, 100, 50, 16
END

IDD_DRAW_GRAPH_DIALOG_FRIEDMANN DIALOGEX 0, 0, 600, 310
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Draw Graph Dialog Friedmann Universe"
FONT 10, "Courier New", 700
BEGIN
END

/////////////////////////////////////////////////////////////////////////////
//
// DESIGNINFO
//

#ifdef APSTUDIO_INVOKED
GUIDELINES DESIGNINFO
BEGIN
    IDD_ABOUTBOX, DIALOG
    BEGIN
        LEFTMARGIN, 7
        RIGHTMARGIN, 163
        TOPMARGIN, 7
        BOTTOMMARGIN, 55
    END
END
#endif    // APSTUDIO_INVOKED

#ifdef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// TEXTINCLUDE
//
1 TEXTINCLUDE
BEGIN
    "resource.h\0"
END

2 TEXTINCLUDE
BEGIN
    "#ifndef APSTUDIO_INVOKED\r\n"
    "#include ""targetver.h""\r\n"
    "#endif\r\n"
    "#define APSTUDIO_HIDDEN_SYMBOLS\r\n"
    "#include ""windows.h""\r\n"
    "#undef APSTUDIO_HIDDEN_SYMBOLS\r\n"
    "\0"
END

3 TEXTINCLUDE
BEGIN
    "\r\n"
    "\0"
END

#endif    // APSTUDIO_INVOKED

/////////////////////////////////////////////////////////////////////////////
//
// String Table
//

STRINGTABLE
BEGIN
   IDC_MODELUNIVERSES   "MODELUNIVERSES"
   IDS_APP_TITLE       "ModelUniverses"
END

#endif
/////////////////////////////////////////////////////////////////////////////



#ifndef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//

/////////////////////////////////////////////////////////////////////////////
#endif    // not APSTUDIO_INVOKED


//{{NO_DEPENDENCIES}}
// Microsoft Visual C++ generated include file.
// Used by ModelUniverses.rc

#define IDS_APP_TITLE			103

#define IDR_MAINFRAME			128
#define IDD_MODELUNIVERSES_DIALOG	102
#define IDD_ABOUTBOX			103
#define IDM_ABOUT				104
#define IDM_EXIT				105
#define IDI_MODELUNIVERSES		107
#define IDI_SMALL				108
#define IDC_MODELUNIVERSES		109
#define IDC_MYICON				2
#ifndef IDC_STATIC
#define IDC_STATIC				-1
#endif
// Next default values for new objects
//
#ifdef APSTUDIO_INVOKED
#ifndef APSTUDIO_READONLY_SYMBOLS

#define _APS_NO_MFC					130
#define _APS_NEXT_RESOURCE_VALUE	129
#define _APS_NEXT_COMMAND_VALUE		32771
#define _APS_NEXT_CONTROL_VALUE		1000
#define _APS_NEXT_SYMED_VALUE		110
#endif
#endif

#define IDD_DATA_INPUT_DIALOG_DE_SITTER	1000
#define IDC_COMBO_A_1				1010
#define IDC_COMBO_B					1020
#define IDC_COMBO_CT				1030
#define IDC_COMBO_EPSILON_1			1040
#define IDC_COMBO_LAMBDA			1050
#define IDC_BUTTON_DRAW_DE_SITTER	1060
#define IDD_DRAW_GRAPH_DIALOG_DE_SITTER	2000
#define IDD_DATA_INPUT_DIALOG_RADIATION	3000
#define IDC_COMBO_A_2				3010
#define IDC_COMBO_DELTA_T			3020
#define IDC_COMBO_EPSILON_2			3030
#define IDC_BUTTON_DRAW_RADIATION	3040
#define IDD_DRAW_GRAPH_DIALOG_RADIATION	4000
#define IDD_DATA_INPUT_DIALOG_FRIEDMANN 5000
#define IDC_COMBO_M					5010
#define IDC_COMBO_CT1				5020
#define IDC_COMBO_EPSILON_3			5030
#define IDC_BUTTON_DRAW_FRIEDMANN	5040
#define IDD_DRAW_GRAPH_DIALOG_FRIEDMANN	6000
#define IDM_DE_SITTER				32771
#define IDM_RADIATION				32772
#define IDM_FRIEDMANN				32773

Blog Entry © Thursday, March 19, 2026, by James Pate Williams, Jr. Solution of Laplace’s Two-Dimensional Equation on a Square Boundary

LaplaceFD2dDownload
#pragma once

class Laplace
{
private:
    static int N;

public:
    Laplace(int N);
    double** Solve(
        double L,
        double epsilon,
        int& iterations,
        double (*f)(double));
    void PrintSolution(
        double L,
        double epsilon,
        double** u,
        double** w,
        int& iterations);
};


#include "Laplace.h"
#include <float.h>
#include <math.h>
#include <iomanip>
#include <iostream>

int Laplace::N = 0;

static double fxy(double x, double y)
{
    return x * y * exp(-0.5 * (x * x + y * y));
}

static double gxy(double x, double y)
{
    return x * sin(y);
}

Laplace::Laplace(int N)
{
    this->N = N;
}

double** Laplace::Solve(
    double L, double epsilon, int& iterations,
    double (*f)(double))
{
    double h = L / (N + 1.0), uv = 0, max = 0;
    double** u = new double*[N + 2];

    if (u == nullptr)
        exit(-1);

    for (int i = 0; i < N + 2; i++)
    {
        u[i] = new double[N + 2];

        if (u[i] == nullptr)
            exit(-2);

        for (int j = 0; j < N + 2; j++)
            u[i][j] = 0.0;
    }

    double** v = new double* [N + 2];

    if (v == nullptr)
        exit(-3);

    for (int i = 0; i < N + 2; i++)
    {
        v[i] = new double[N + 2];

        if (v[i] == nullptr)
            exit(-4);

        for (int j = 0; j < N + 2; j++)
            v[i][j] = 0.0;
    }

    iterations = 0;

    for (int j = 0; j <= N + 1; j++)
        u[0][j] = 0.0;

    for (int j = 0; j <= N + 1; j++)
        u[N + 1][j] = 0.0;

    for (int i = 0; i <= N + 1; i++)
        u[i][0] = f(i * h);

    for (int i = 0; i <= N + 1; i++)
        u[i][N + 1] = 0.0;

    do
    {
        iterations++;

        for (int i = 1; i <= N; i++)
            for (int j = 1; j <= N; j++)
                u[i][j] = 0.25 * (u[i - 1][j] + u[i + 1][j] + u[i][j - 1] +
                    u[i][j + 1]);

        for (int i = 0; i <= N + 1; i++)
            for (int j = 0; j <= N + 1; j++)
                v[i][j] = u[i][j];

        for (int i = 1; i <= N; i++)
            for (int j = 1; j <= N; j++)
                v[i][j] = 0.25 * (v[i - 1][j] + v[i + 1][j] + v[i][j - 1] +
                    v[i][j + 1]);

        max = DBL_MIN;

        for (int i = 1; i <= N; i++)
        {
            for (int j = 1; j <= N; j++)
            {
                uv = fabs(v[i][j] - u[i][j]);
                max = fmax(uv, max);
            }
        }

        for (int i = 0; i <= N + 1; i++)
            for (int j = 0; j <= N + 1; j++)
                u[i][j] = v[i][j];

    } while (max > epsilon);

    for (int i = 0; i < N + 2; i++)
        delete v[i];

    return u;
}

void Laplace::PrintSolution(
    double L,
    double epsilon,
    double** u,
    double** w,
    int& iterations)
{
    double h = L / (N + 1.0);
    
    std::cout << std::endl;
    std::cout << "x\t\ty\t\tu(x, y)\t\tw(x, y)\t\tpercent error\r\n\r\n";

    for (int i = 0; i <= N + 1; i++)
    {
        double x = i * h;

        for (int j = 0; j <= N + 1; j++)
        {
            double y = j * h;
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << x << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << y << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << u[i][j] << '\t';
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << w[i][j] << '\t';
            double error = 0.0;
            if (w[i][j] != 0.0)
                error = 100.0 * fabs(u[i][j] - w[i][j]) /
                    fabs(w[i][j]);
            std::cout << std::fixed << std::setprecision(5);
            std::cout << std::setw(8);
            std::cout << error << std::endl;
        }
    }

    std::cout << std::endl;
    std::cout << "iterations = " << iterations;
    std::cout << std::endl;
}


// LaplaceEquationFD2d.cpp (c) Tuesday, March 17, 2026 
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// by James Pate Williams, Jr., BA, ,BS, MSwE, PhD
// Reference: Numerical Computation 2 Methods,
// Software And Analysis (c) 1995 by Chrisoph E. Ueberhuber
// Pages 392 - 393 Poisson Matrices
// https://math.libretexts.org/Bookshelves/Differential_Equations/Introduction_to_Partial_Differential_Equations_(Herman)/06%3A_Problems_in_Higher_Dimensions/6.03%3A_Laplaces_Equation_in_2D

#include <math.h>
#include <iostream>
#include "Laplace.h"

double** w;

double SimpsonsRule(double L, double a, double b,
    int Nsr, int n, double (*f)(double, double, int))
{
    double h = (b - a) / Nsr;
    double h2 = 2.0 * h;
    double s = 0.0;
    double t = 0.0;
    double x = a + h;

    for (int i = 1; i < Nsr; i += 2)
    {
        s += f(L, x, n);
        x += h2;
    }

    x = a + h2;

    for (int i = 2; i < Nsr; i += 2)
    {
        t += f(L, x, n);
        x += h2;
    }

    return
        h * (f(L, a, n) + 4 * s + 2 * t + f(L, b, n)) / 3.0;
}

double f(double x)
{
    return x * exp(-x) * cos(x);
}

double g(double L, double x, int n)
{
    double pi = 4.0 * atan(1.0);
    double fx = f(x);
    return fx * sin(n * pi * x / L);
}

double bn(double L, int Nsr, int n)
{
    double cs = 2.0 / L;
    double sr = SimpsonsRule(L, 0.0, L, Nsr, n, g);
    return cs * sr;
}

void ComputeW(double L, int N)
{
    double h = L / (N + 1.0);

    w = new double* [N + 2];

    if (w == nullptr)
        exit(-5);

    for (int i = 0; i < N + 2; i++)
    {
        w[i] = new double[N + 2];

        if (w[i] == nullptr)
            exit(-6);

        for (int j = 0; j < N + 2; j++)
            w[i][j] = 0.0;
    }

    for (int j = 0; j <= N; j++)
        w[0][j] = 0.0;

    for (int j = 0; j <= N; j++)
        w[N + 1][j] = 0.0;

    for (int i = 0; i <= N; i++)
        w[i][0] = f(i * h);

    for (int i = 0; i <= N; i++)
        w[i][N + 1] = 0.0;

    double pi = 4.0 * atan(1.0);
    double x = h, y = 0.0;

    for (int i = 1; i <= N; i++)
    {
        y = h;

        for (int j = 1; j <= N; j++)
        {
            double s = 0.0;

            for (int n = 1; n <= 64; n++)
            {
                double npi = n * pi;
                double an = bn(L, 512, n) / sinh(npi);
                s += an * sin(npi * x / L) * sinh(npi * (L - y) / L);
            }

            w[i][j] = s;
            y += h;
        }

        x += h;
    }
}

int main()
{
    while (true)
    {
        char line[128] = "";
        std::cout << "N = ";
        std::cin.getline(line, 128);
        int N = atoi(line);

        if (N < 4 || N > 128)
        {
            std::cout << "N < 4 || N > 128, try again";
            std::cout << std::endl;
        }

        std::cout << "epsilon = ";
        std::cin.getline(line, 128);

        double epsilon = atof(line);
        Laplace l(N);

        double L = 1.0;
        int iterations = 0;
        ComputeW(L, N);
        double** u = l.Solve(L, epsilon, iterations, f);
        l.PrintSolution(L, epsilon, u, w, iterations);

        for (int i = 0; i < N + 2; i++)
        {
            if (u[i] != nullptr)
                delete u[i];

            if (w[i] != nullptr)
                delete w[i];
        }

        break;
    }

    return 0;
}

Blog Entry © Tuesday, March 17, 2026, by James Pate Williams, Jr., Comparison of Power Series Arctangent and Arcsine Functions with the C++ Built-In Functions

Blog Entry Tuesday March 17 2026Download
#pragma once

//https://scipp-legacy.pbsci.ucsc.edu/~haber/ph116A/taylor11.pdf

class Functions
{
public:
	static void initialize();
	static double factorial(int n);
	static double arccosine(double x, int n);
	static double arccosecant(double x);
	static double arccotangent(double x);
	static double arcsecant(double x);
	static double arcsine(double x, int n);
	static double arctangent(double x, int n);
};


#include "Functions.h"
#include <math.h> 

double pi2 = 0.0;

double Functions::factorial(int n)
{
	double nfactorial = 1.0;

	for (int i = 2; i <= n; i++)
		nfactorial *= i;

	return nfactorial;
}

void Functions::initialize()
{
	pi2 = 2.0 * atan(1.0);
}

double Functions::arccosine(double x, int n)
{
	double sum = pi2 - arcsine(x, n);
	return sum;
}

double Functions::arccosecant(double x)
{
	double sum = 0;
	return sum;
}

double Functions::arccotangent(double x)
{
	double sum = 0;
	return sum;
}

double Functions::arcsecant(double x)
{
	double sum = 0;
	return sum;
}

double Functions::arcsine(double x, int n)
{
	double sum = 0.0;

	if (fabs(x) <= 1.0)
	{
		for (int i = n; i >= 0; i--)
		{
			double ifact = factorial(i);
			double i2 = 2.0 * i, i21 = 2 * i + 1;
			double coeff = factorial(i2) /
				(pow(2, i2) * ifact * ifact * (i21));
			sum += coeff * pow(x, i21);
		}
	}

	return sum;
}

double Functions::arctangent(double x, int n)
{
	double sum = 0.0;

	if (fabs(x) <= 1.0)
	{
		double one = 0.0;

		if (n % 2 == 0)
			one = 1.0;
		else
			one = -1.0;

		for (int i = 2 * n + 1; i >= 0; i--)
		{
			one = -one;
			double i21 = 2 * i + 1.0;
			sum += one * pow(x, i21) / i21;
		}
	}

	return sum;
}


// InvTrigConsoleCPP.cpp (c) Monday, March 16, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD

#include "Functions.h"
#include <math.h>
#include <iomanip>
#include <iostream>

static void CreateTable(
	char title[],
	double a, double b, int n, int nPts,
	double(*fx)(double, int), double(*gx)(double))
{
	double x = a;
	double deltaX = (b - a) / nPts;
	std::cout << title << std::endl;
	std::cout << "# Terms = " << (2 * n + 2) << std::endl;
	std::cout << "x" << '\t' << "fx" << "\t\t";
	std::cout << "sx" << "\t\t" << "error" << std::endl;

	for (int i = 0; i < nPts; i++)
	{
		double f = fx(x, n);
		double s = gx(x);
		double e = 0.0;

		if (fabs(s) != 0)
			e = 100.0 * fabs(f - s) / fabs(s);

		std::cout << std::fixed << std::setprecision(4);
		std::cout << std::setw(4);
		std::cout << x << '\t';
		std::cout << std::fixed << std::setprecision(11);
		std::cout << std::setw(12);
		std::cout << f << '\t';
		std::cout << std::fixed << std::setprecision(11);
		std::cout << std::setw(12);
		std::cout << s << '\t';
		std::cout << e << std::endl;
		x += deltaX;
	}
}

int main()
{
	Functions::initialize();

	while (true)
	{
		char line[128] = "";
		std::cout << "== Menu == " << std::endl;
		std::cout << "1 arcsine" << std::endl;
		std::cout << "2 arctangent" << std::endl;
		std::cout << "3 exit" << std::endl;
		std::cout << "option # : ";
		std::cin.getline(line, 128);
		char option = line[0];

		if (option == '3')
			break;

		if (option < '1' || option > '2')
		{
			std::cout << "invalid option" << std::endl;
			continue;
		}

		std::cout << "# points: ";
		std::cin.getline(line, 128);
		int nPts = atoi(line);

		std::cout << "# terms: ";
		std::cin.getline(line, 128);
		int n = atoi(line);

		if (option == '1')
		{
			char title[] = "Series arcsin Versus C++ asin";
			CreateTable(title, 0.0, 1.0, n, nPts,
				Functions::arcsine, asin);
		}
		else if (option == '2')
		{
			char title[] = "Series arctan Versus C++ atan";
			CreateTable(title, 0.0, 1.0, n, nPts,
				Functions::arctangent, atan);
		}
	}

	return 0;
}

Blog Entry © Sunday March 15, 2026, by James Pate Williams, Jr. Graphing Integer Order Bessel Functions of the First Kind and the real Gamma Function

Blog Entry Sunday March 15 2026Download

Blog Entry © Saturday, March 14, 2026, by James Pate Williams, Jr. Power Series Solution of a Second Order Linear Ordinary Differential Equation

Blog Entry Saturday March 14 2026Download
// SecondOrderLinearODE.cpp (c) Friday March 13, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Solves the Bessel Equation of the First Kind for Order 0
// Using the power series found in the reference below - 
// References: https://math.libretexts.org/Bookshelves/Calculus/Calculus_(OpenStax)/17%3A_Second-Order_Differential_Equations/17.04%3A_Series_Solutions_of_Differential_Equations
// https://mathworld.wolfram.com/BesselDifferentialEquation.html

#include "framework.h"
#include <Windows.h>
#include "SecondOrderLinearODE.h"
#include <string>
#include <vector>

#define MAX_LOADSTRING 100

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
WCHAR fTitle[128];
double chi[1025], eta[1025], coeff[65];

// Forward declarations of functions included in this code module:
ATOM                MyRegisterClass(HINSTANCE hInstance);
BOOL                InitInstance(HINSTANCE, int);
LRESULT CALLBACK    WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    About(HWND, UINT, WPARAM, LPARAM);

int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    UNREFERENCED_PARAMETER(hPrevInstance);
    UNREFERENCED_PARAMETER(lpCmdLine);

    // TODO: Place code here.

    // Initialize global strings
    LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
    LoadStringW(hInstance, IDC_SECONDORDERLINEARODE, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_SECONDORDERLINEARODE));

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, nullptr, 0, 0))
    {
        if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
        {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }

    return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
    WNDCLASSEXW wcex = { 0 };

    wcex.cbSize = sizeof(WNDCLASSEX);

    wcex.style          = CS_HREDRAW | CS_VREDRAW;
    wcex.lpfnWndProc    = WndProc;
    wcex.cbClsExtra     = 0;
    wcex.cbWndExtra     = 0;
    wcex.hInstance      = hInstance;
    wcex.hIcon          = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_SECONDORDERLINEARODE));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_SECONDORDERLINEARODE);
    wcex.lpszClassName  = szWindowClass;
    wcex.hIconSm        = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

    return RegisterClassExW(&wcex);
}

//
//   FUNCTION: InitInstance(HINSTANCE, int)
//
//   PURPOSE: Saves instance handle and creates main window
//
//   COMMENTS:
//
//        In this function, we save the instance handle in a global variable and
//        create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   hInst = hInstance; // Store instance handle in our global variable

   HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

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

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

    return factorial;
}

static void ComputeCoefficients()
{
    double a0 = 1.0;

    for (int k = 0; k <= 64; k++)
    {
        double kFactorial = Factorial(k);

        coeff[k] = a0 * pow(-1.0, k) /
            (pow(2.0, 2.0 * k) * pow(kFactorial, 2.0));
    }
}

static void ComputePoints()
{
    double chi0 = -15.0;
    double deltaChi = 30.0 / 1024;
 
    for (int i = 0; i <= 1024; i++)
    {
        double sum = 0.0;

        for (int k = 64; k >= 0; k--)
            sum += pow(chi0, 2 * k) * coeff[k];

        chi[i] = chi0;
        eta[i] = sum;

        chi0 += deltaChi;
    }
}

static void FindMinMax(
    double x[],
    double& min,
    double& max)
{
    min = x[0];
    max = x[0];

    for (int i = 1; i <= 1024; i++)
    {
        if (x[i] < min)
            min = x[i];
        if (x[i] > max)
            max = x[i];
    }
}

static void DrawTitles(HDC hdc, RECT clientRect, const std::wstring& fTitle,
    int sx0, int sx1, int sy0, int sy1)
{
    HFONT hCustomFont = CreateFont(
        -MulDiv(10, GetDeviceCaps(hdc, LOGPIXELSY), 72), // Height in logical units
        0,                      // Width (0 = default)
        0,                      // Escapement
        0,                      // Orientation
        FW_BOLD,                // Weight (700 = bold)
        FALSE,                  // Italic
        FALSE,                  // Underline
        FALSE,                  // StrikeOut
        DEFAULT_CHARSET,        // Charset
        OUT_DEFAULT_PRECIS,     // Output precision
        CLIP_DEFAULT_PRECIS,    // Clipping precision
        DEFAULT_QUALITY,        // Quality
        FIXED_PITCH | FF_MODERN,// Pitch and family
        L"Courier New"          // Typeface name
    );

    HFONT hOldFont = (HFONT)SelectObject(hdc, hCustomFont);

    SIZE sz;
    int width = clientRect.right - clientRect.left;

    // Title: fTitle
    GetTextExtentPoint32(hdc, fTitle.c_str(), (int)fTitle.length(), &sz);
    int w = sz.cx;
    int h = sz.cy;
    TextOut(hdc, (width - w) / 2, h, fTitle.c_str(), (int)fTitle.length());

    // x-axis title: "x"
    std::wstring xTitle = L"x";
    GetTextExtentPoint32W(hdc, xTitle.c_str(), (int)xTitle.length(), &sz);
    w = sz.cx;
    TextOutW(hdc, sx0 + (sx1 - sx0 - w) / 2, sy1 + 2 * h, xTitle.c_str(), (int)xTitle.length());

    // y-axis title: "p(x)"
    std::wstring yTitle = L"J0(x)";
    GetTextExtentPoint32W(hdc, yTitle.c_str(), (int)yTitle.length(), &sz);
    w = sz.cx;
    TextOutW(hdc, sx1 + w / 5, sy0 + (sy1 - sy0) / 2 - h / 2, yTitle.c_str(), (int)yTitle.length());

    SelectObject(hdc, hOldFont); // Restore original font
}

static void DrawFormattedText(HDC hdc, char loctext[], RECT rect)
{
    // Draw the text with formatting options
    DrawTextA(hdc, loctext, (int)strlen(loctext), &rect, DT_SINGLELINE | DT_NOCLIP);
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    switch (message)
    {
    case WM_CREATE:
        ComputeCoefficients();
        ComputePoints();
        break;
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);
            // Parse the menu selections:
            switch (wmId)
            {
            case IDM_ABOUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
                break;
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
    {
        PAINTSTRUCT ps;
        HDC hdc = BeginPaint(hWnd, &ps);

        ComputePoints();
        double xMin = 0, yMin = 0;
        double xMax = 0, yMax = 0;
        FindMinMax(chi, xMin, xMax);
        FindMinMax(eta, yMin, yMax);

        double h = 0, pi = 0, plm = 0, theta = 0;
        float xSpan = (float)(xMax - xMin);
        float ySpan = (float)(yMax - yMin);
        RECT rect = { };
        GetClientRect(hWnd, &rect);
        float width = (float)(rect.right - rect.left + 1);
        float height = (float)(rect.bottom - rect.top - 32 + 1);
        float sx0 = 2.0f * width / 16.0f;
        float sx1 = 14.0f * width / 16.0f;
        float sy0 = 2.0f * height / 16.0f;
        float sy1 = 14.0f * height / 16.0f;
        float deltaX = xSpan / 8.0f;
        float deltaY = ySpan / 8.0f;
        float xSlope = (sx1 - sx0) / xSpan;
        float xInter = (float)(sx0 - xSlope * xMin);
        float ySlope = (sy0 - sy1) / ySpan;
        float yInter = (float)(sy0 - ySlope * yMax);
        float px = 0, py = 0, sx = 0, sy = 0;
        POINT wPt = { 0 };
        int i = 0;
        float x = (float)xMin;
        float y = (float)yMax;
        wcscpy_s(fTitle, 128, L"Bessel Function of the First Kind Order 0");
        DrawTitles(hdc, rect, fTitle, (int)sx0, (int)sx1, (int)sy0, (int)sy1);


        px = (float)chi[0];
        py = (float)eta[0];

        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        MoveToEx(hdc, (int)sx, (int)sy0, &wPt);

        while (i <= 8)
        {
            sx = xSlope * x + xInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
            LineTo(hdc, (int)sx, (int)sy1);
            char numberStr[128] = "";
            sprintf_s(numberStr, 128, "%5.4lf", x);
            SIZE size = { };
            GetTextExtentPoint32A(
                hdc,
                numberStr,
                (int)strlen(numberStr),
                &size);
            RECT textRect = { 0 };
            textRect.left = (long)(sx - size.cx / 2.0f);
            textRect.right = (long)(sx + size.cx / 2.0f);
            textRect.top = (long)sy1;
            textRect.bottom = (long)(sy1 + size.cy / 2.0f);
            DrawFormattedText(hdc, numberStr, textRect);
            x += deltaX;
            i++;
        }

        i = 0;
        y = (float)yMin;

        while (i <= 8)
        {
            sy = ySlope * y + yInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
            LineTo(hdc, (int)sx, (int)sy);

            if (i != 0)
            {
                char numberStr[128] = "";
                sprintf_s(numberStr, 128, "%+5.3lf", y);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    numberStr,
                    (int)strlen(numberStr),
                    &size);
                RECT textRect = { 0 };
                textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
                textRect.right = (long)(sx0 - size.cx / 2.0f);
                textRect.top = (long)(sy - size.cy / 2.0f);
                textRect.bottom = (long)(sy + size.cy / 2.0f);
                DrawFormattedText(hdc, numberStr, textRect);
            }

            y += deltaY;
            i++;
        }

        HGDIOBJ nPenNew = NULL;
        HGDIOBJ hPenOld = NULL;

        nPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
        hPenOld = SelectObject(hdc, nPenNew);

        HRGN clipRegion = CreateRectRgn(
            (int)sx0, (int)sy0,              // Left, Top
            (int)(sx1), (int)(sy1)           // Right, Bottom
        );

        // Apply clipping region

        SelectClipRgn(hdc, clipRegion);


        px = (float)chi[0];
        py = (float)eta[0];

        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        wPt.x = wPt.y = 0;
        MoveToEx(hdc, (int)sx, (int)sy, &wPt);

        for (size_t j = 1; j <= 1024; j++)
        {

            px = (float)chi[j];
            py = (float)eta[j];

            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            LineTo(hdc, (int)sx, (int)sy);
        }

        DeleteObject(nPenNew);
        nPenNew = NULL;

        SelectObject(hdc, hPenOld);
        DeleteObject(clipRegion);

        EndPaint(hWnd, &ps);
    }
    break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

Blog Entry © Thursday and Friday, March 12 – 13, 2026, Naive Power Series Computations of the Cosine and Sine Trigonometric Functions and Polynomial Fit to the Trig Functions Using Least Squares

Blog Entry Thursday March 12 2026Download
// TrigLSPolyFit.cpp (c) Thursday, March 12, 2026
// by James Pate Williams, Jr.
// Least Squares Curve Fitting to the Trigonometric
// functions Cosine and Sine

#include <iomanip>
#include <iostream>
#include <vector>
#include "PolyLSFit.h"

static const int HornerN = 50;
double even_factors[HornerN], odd_factors[HornerN + 1];

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

	for (int i = 2; i <= n; i++)
		nfactorial *= i;
	return nfactorial;
}

static void computeEvenFactors()
{
	int one = 1;

	for (int n = 0; n < 50; n += 2)
	{
		even_factors[n] = one / factorial(n);
		one = -one;
	}
}

static void computeOddFactors()
{
	int one = 1;

	for (int n = 1; n < 51; n += 2)
	{
		odd_factors[n] = one / factorial(n);
		one = -one;
	}
}

static double fx_cos(double x, int n)
{
	double sum = even_factors[n];

	for (int i = n - 1; i >= 0; i--)
	{
		sum = sum * x + even_factors[i];
	}

	return sum;
}

static double fx_sin(double x, int n)
{
	double sum = odd_factors[n];

	for (int i = n - 1; i >= 0; i--)
	{
		sum = sum * x + odd_factors[i];
	}

	return sum;
}

static const int M = 33, N = 100, N21 = N + N + 1;
std::vector<double> cosX(N21), cosY(N21), cosCoeff;
std::vector<double> sinX(N21), sinY(N21), sinCoeff;
std::vector<double> pX(M), pY(M);

static void LSPolyFit(int HornerN, int degree)
{
	double pi = 4.0 * atan(1.0);
	double deltaX = 2.0 * pi / N;
	double X = -pi;

	for (int i = 0; i < N21; i++)
	{
		cosX[i] = X;
		cosY[i] = fx_cos(X, HornerN);
		sinX[i] = X;
		sinY[i] = fx_sin(X, HornerN + 1);
		X += deltaX;
	}

	PolyLSFit::leastSquares(cosX, cosY, cosCoeff, degree, N21);
	PolyLSFit::leastSquares(sinX, sinY, sinCoeff, degree, N21);

	std::cout << std::setw(6) << "x";
	std::cout << std::setw(17) << "fit cos";
	std::cout << std::setw(17) << "C++ cos";
	std::cout << std::setw(17) << "% Error";
	std::cout << std::endl;

	deltaX = 2.0 * pi / (M - 1);
	X = -pi;

	for (int i = 0; i < M; i++)
	{
		double Y = cosCoeff[degree];

		for (int j = degree - 1; j >= 0; j--)
		{
			Y = Y * X + cosCoeff[j];
		}

		pX[i] = X;
		pY[i] = Y;

		std::cout << std::fixed << std::setprecision(3);
		std::cout << std::setw(6);
		std::cout << X << ' ';

		if (Y >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);
		
		std::cout << fabs(Y) << ' ';
		double cx = cos(X);

		if (cx >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);
		
		std::cout << fabs(cx) << ' ';
		double pcDifference = 0.0;

		if (fabs(cx) >= 1.0e-12)
			pcDifference = 100.0 * fabs(Y - cx) / fabs(cx);
		else
			pcDifference = -1.0;

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		if (pcDifference != -1.0)
			std::cout << pcDifference << std::endl;
		else
			std::cout << "not a number" << std::endl;

		X += deltaX;
	}

	std::cout << std::endl;

	std::cout << std::setw(6) << "x";
	std::cout << std::setw(17) << "fit sin";
	std::cout << std::setw(17) << "C++ sin";
	std::cout << std::setw(17) << "% Error";
	std::cout << std::endl;

	deltaX = 2.0 * pi / (M - 1);
	X = -pi;

	for (int i = 0; i < M; i++)
	{
		double Y = sinCoeff[degree];

		for (int j = degree - 1; j >= 0; j--)
		{
			Y = Y * X + sinCoeff[j];
		}

		pX[i] = X;
		pY[i] = Y;

		std::cout << std::fixed << std::setprecision(3);
		std::cout << std::setw(6);
		std::cout << X << ' ';

		if (Y >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		std::cout << fabs(Y) << ' ';
		double sx = sin(X);

		if (sx >= 0)
			std::cout << '+';
		else
			std::cout << '-';

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		std::cout << fabs(sx) << ' ';
		double pcDifference = 0.0;

		if (fabs(sx) >= 1.0e-12)
			pcDifference = 100.0 * fabs(Y - sx) / fabs(sx);
		else
			pcDifference = -1.0;

		std::cout << std::fixed << std::setprecision(16);
		std::cout << std::setw(16);

		if (pcDifference != -1.0)
			std::cout << pcDifference << std::endl;
		else
			std::cout << "not a number" << std::endl;

		X += deltaX;
	}

	std::cout << std::endl;
}

int main()
{
	computeEvenFactors();
	computeOddFactors();
	std::cout << "# terms 16" << std::endl;
	LSPolyFit(HornerN, 16);
	std::cout << "# terms 20" << std::endl;
	LSPolyFit(HornerN, 20);
	std::cout << "# terms 24" << std::endl;
	LSPolyFit(HornerN, 24);
	return 0;
}

Blog Entry © Thursday, February 5, 2026, by James Pate Williams, Jr., Three Toom-Cook 3-Way Multiplication Implementations

I translated some of Marco Bodrato ‘s pseudocode found on his website:

Optimal Toom-Cook Polynomial Multiplication, by Marco Bodrato

To Microsoft Win32 C++. I used Microsoft’s Visual Studio 2022 Community Version.

Balanced Toom-Cook 3-Way Multiplication

U = 7

V = 8

W = 56

U = 9

V = 9

W = 81

U = 18

V = 17

W = 306

U = 123

V = 456

W = 56088

U = 1234

V = 5678

W = 7006652

U = 12345

V = 67890

W = 838102050

U = 123456

V = 789012

W = 97408265472

D:\ToomCookBodrato\x64\Debug\ToomCookBodrato.exe (process 30392) exited with code 0 (0x0).

Unbalanced Toom-Cook 3-Way Multiplication

U = 20

V = 5

W = 100

U = 112

V = 20

W = 2240

U = 1234

V = 567

W = 699678

U = 12345

V = 6789

W = 83810205

U = 123456

V = 78901

W = 9740801856

U = 1234567

V = 890123

W = 1098916481741

U = 12345678

V = 9012345

W = 111263509394910

D:\ToomCookBodrato\x64\Debug\ToomCookBodrato.exe (process 13108) exited with code 0 (0x0).

Asymmetrical Squaring, Splitting in Five

U = 8

W = 64

U = 12

W = 144

U = 321

W = 103041

U = 1234

W = 1522756

U = 54321

W = 2950771041

U = 123456

W = 15241383936

D:\ToomCookBodrato\x64\Debug\ToomCookBodrato.exe (process 27592) exited with code 0 (0x0).

Press any key to close this window . . .

The accompanying source code is released under the GNU General Public License, version 3. The full license text is included in the downloadable archive and is also available at:
https://www.gnu.org/licenses/gpl-3.0.html

/*
   Toom–Cook Multiplication (Bodrato-inspired implementation)
   Copyright (C) 2026 James Pate Williams, Jr.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#pragma once

// Algorithm structure inspired by Marco Bodrato’s Toom–Cook notes.
// See: http://www.bodrato.it/toom-cook/
// All implementation code and comments in this file are by
// James Pate Williams, Jr., BA, BS, Master of Software Engineering,
// Doctor of Philosophy in Computer Science

class Multiplication
{

public:

	static long long BalancedToomCook3Way(
		long long b, long long U, long long V);
	static long long UnbalancedToomCook3Way(
		long long b, long long U, long long V);
	static long long AsymeticalSquaring(
		long long b, long long U);
};


/*
   Toom–Cook Multiplication (Bodrato-inspired implementation)
   Copyright (C) 2026 James Pate Williams, Jr.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#include "pch.h"
#include "Multiplication.h"

static void BaseBDigits(
	long long B, long long u,
	std::vector<long long>& digits)
{
	while (u > 0)
	{
		long long d = u % B;
		u = u / B;
		digits.push_back(d);
	}
}

long long Multiplication::BalancedToomCook3Way(
	long long b, long long U, long long V)
{
	double k = 3.0;
	double logU = log(U) / log(b);
	double logV = log(V) / log(b);
	double floorU = floor(floor(logU) / k);
	double floorV = floor(floor(logV) / k);
	double i = fmax(floorU, floorV) + 1.0;
	long long B = static_cast<long long>(pow(b, i));
	long long W = 0;
	std::vector<long long> UDigits;
	std::vector<long long> VDigits;

	for (long long j = 0; j < 5; j++)
	{
		long long x = 0, y = 0;

		if (j == 0)
			x = 1;
		else if (j == 1)
			x = -2;
		else if (j == 2)
			x = 1;
		else if (j == 3)
			x = -1;
		else
			x = 1;

		if (j == 0)
			y = 1;
		else if (j == 1)
			y = 1;
		else if (j == 2)
			y = 1;
		else if (j == 3)
			y = 1;
		else
			y = 0;

		UDigits.clear();
		VDigits.clear();

		BaseBDigits(B, U, UDigits);
		BaseBDigits(B, V, VDigits);

		long long U0 = 0, U1 = 0, U2 = 0;
		long long V0 = 0, V1 = 0, V2 = 0;

		if (UDigits.size() == 3)
		{
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (VDigits.size() == 3)
		{
			V2 = VDigits[2];
			V1 = VDigits[1];
			V0 = VDigits[0];
		}

		if (UDigits.size() == 2)
		{
			U2 = 0;
			U1 = UDigits[1];
			U2 = UDigits[0];
		}

		if (VDigits.size() == 2)
		{
			V2 = 0;
			V1 = VDigits[1];
			V0 = VDigits[0];
		}

		if (UDigits.size() == 1)
		{
			U2 = 0;
			U1 = 0;
			U0 = UDigits[0];
		}

		if (VDigits.size() == 1)
		{
			V2 = 0;
			V1 = 0;
			V0 = VDigits[0];
		}

		if (UDigits.size() == 0 ||
			VDigits.size() == 0)
			break;

		long long x2 = x * x;
		long long y2 = y * y;
		long long xy = x * y;
		
		U = U2 * x2 + U1 * xy + U0 * y2;
		V = V2 * x2 + V1 * xy + V0 * y2;
		
		if (U == 0 || V == 0)
			break;

		long long W3 = U0 + U2;
		long long W2 = V0 + V2;
		long long W0 = W3 - U1;
		long long W4 = W2 - V1;
		
		W3 = W3 + U1;
		W2 = W2 + V1;

		long long W1 = W3 * W2;

		W2 = W0 * W4;
		W0 = ((W0 + U2) << 1) - U0;
		W4 = ((W4 + V2) << 1) - V0;
		W3 = W0 * W4;
		W0 = U0 * V0;
		W4 = U2 * V2;

		W3 = (W3 - W1) / k;
		W1 = (W1 - W2) >> 1;
		W2 = W2 - W0;
		W3 = ((W2 - W3) >> 1) + (W4 << 1);
		W2 = W2 + W1;

		W3 = W4 * x + W3 * y;
		W1 = W2 * x + W1 * y;
		W1 = W1 - W3;

		W = W3 * x * x2 + W1 * x * y2 + W0 * y2 * y2;
	}

	return W;
}

long long Multiplication::UnbalancedToomCook3Way(
	long long b, long long U, long long V)
{
	double k = 3.0;
	double logU = log(U) / log(b);
	double logV = log(V) / log(b);
	double floorU = floor(floor(logU) / k);
	double floorV = floor(floor(logV) / k);
	double i = fmax(floorU, floorV) + 1.0;
	long long B = static_cast<long long>(pow(b, i));
	long long W = 0;
	std::vector<long long> UDigits;
	std::vector<long long> VDigits;

	for (long long j = 0; j < 5; j++)
	{
		long long x = 0, y = 0;

		if (j == 0)
			x = 1;
		else if (j == 1)
			x = -2;
		else if (j == 2)
			x = 1;
		else if (j == 3)
			x = -1;
		else
			x = 1;

		if (j == 0)
			y = 1;
		else if (j == 1)
			y = 1;
		else if (j == 2)
			y = 1;
		else if (j == 3)
			y = 1;
		else
			y = 0;

		long long U0 = 0, U1 = 0, U2 = 0;
		long long V0 = 0, V1 = 0, V2 = 0;
		long long U3 = 0;
		long long x2 = x * x, y2 = y * y;
		long long x3 = x * x2, y3 = y * y2;

		UDigits.clear();
		VDigits.clear();

		BaseBDigits(B, U, UDigits);
		BaseBDigits(B, V, VDigits);

		if (UDigits.size() == 4)
		{
			U3 = UDigits[3];
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 3)
		{
			U3 = 0;
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (VDigits.size() == 3)
		{
			V2 = VDigits[2];
			V1 = VDigits[1];
			V0 = VDigits[0];
		}

		if (UDigits.size() == 2)
		{
			U3 = 0;
			U2 = 0;
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (VDigits.size() == 2)
		{
			V2 = 0;
			V1 = VDigits[1];
			V0 = VDigits[0];
		}

		if (UDigits.size() == 1)
		{
			U3 = 0;
			U2 = 0;
			U1 = 0;
			U0 = UDigits[0];
		}

		if (VDigits.size() == 1)
		{
			V2 = 0;
			V1 = 0;
			V0 = VDigits[0];
		}

		if (UDigits.size() == 0 ||
			VDigits.size() == 0)
			break;

		U = U3 * x3 + U2 * x2 * y + U1 * x * y2 + U0 * y3;
		V = V1 * x + V0 * y;

		if (U == 0 || V == 0)
			break;

		long long W3 = U0 + U2;
		long long W2 = V0 + V1;
		long long W1 = U1 + U3;
		long long W4 = V0 - V1;
		long long W0 = W3 - W1;
		
		W3 = W3 + W1;
		W1 = W3 * W2;
		W2 = W0 * W4;
		W0 = U0 - (((U1 - ((U2 - U3) << 1)) << 1) << 1);
		W4 = W4 - V1;
		W3 = W0 * W4;
		W0 = U0 * V0;
		W4 = U3 * V1;

		W3 = (W3 - W1) / k;
		W1 = (W1 - W2) >> 1;
		W2 = W2 - W0;
		W3 = ((W2 - W3) >> 1) + (W4 << 1);
		W2 = W2 + W1;

		W3 = W4 * x + W3 * y;
		W1 = W2 * x + W1 * y;
		W1 = W1 - W3;

		W = W3 * x3 + W1 * x * y2 + W0 * y2 * y2;
	}

	return W;
}

long long Multiplication::AsymeticalSquaring(
	long long b, long long U)
{
	double k = 5.0;
	double logU = log(U) / log(b);
	double floorU = floor(floor(logU) / k);
	double i = floorU + 1.0;
	long long B = static_cast<long long>(pow(b, i));
	long long W = 0, Wp = 0;
	std::vector<long long> UDigits;

	for (long long j = 0; j < 9; j++)
	{
		long long x = 0;

		if (j == 0)
			x = 0;
		else if (j == 1)
			x = 0;
		else if (j == 2)
			x = -1;
		else if (j == 3)
			x = 1;
		else if (j == 4 || j == 5 || j == 6)
			continue;
		else if (j == 7 || j == 8)
			x = 1;

		UDigits.clear();
		BaseBDigits(B, U, UDigits);
		
		long long x2 = x * x, x3 = x * x2, x4 = x2 * x2;
		long long U0 = 0, U1 = 0, U2 = 0, U3 = 0, U4 = 0;

		if (UDigits.size() == 5)
		{
			U4 = UDigits[4];
			U3 = UDigits[3];
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 4)
		{
			U4 = 0;
			U3 = UDigits[3];
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 3)
		{
			U4 = 0;
			U3 = 0;
			U2 = UDigits[2];
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 2)
		{
			U4 = 0;
			U3 = 0;
			U2 = 0;
			U1 = UDigits[1];
			U0 = UDigits[0];
		}

		if (UDigits.size() == 1)
		{
			U4 = 0;
			U3 = 0;
			U2 = 0;
			U1 = 0;
			U0 = UDigits[0];
		}

		U = U4 * x4 + U3 * x3 + U2 * x2 + U1 * x + U0;

		long long W0 = U0 - U3;
		long long W1 = U3 - U1;
		long long W6 = U1 - U2;
		long long W4 = U1 + U2;
		long long W5 = W6 - U4;
		long long W3 = W5 + (W0 << 1);
		W0 = W0 - W5;
		W6 = W0 + (W6 << 1);
		long long W7 = W6 + W1;
		W5 = W7 + W1;
		long long W8 = W5 + (W4 << 1);
		W4 = W4 - U4;

		long long W2 = (W4) * (W3);
		W4 = (W6) * (W5);
		W3 = (W7) * (W1);
		W1 = U0 * U1 * 2;
		W7 = U3 * U4 * 2;
		W5 = (W8) * (W8); W6 = (W0) * (W0);
		W0 = U0 * U0;
		W8 = U4 * U4;

		W2 = (W4) * (W3);
		W4 = (W6) * (W5);
		W3 = (W7) * (W1);
		W1 = U0 * U1 * 2;
		W7 = U3 * U4 * 2;
		W5 = (W8) * (W8);
		W6 = (W0) * (W0);
		W0 = U0 * U0;
		W8 = U4 * U4;

		W6 = (W6 + W5) >> 1;
		W5 = W5 - W6;
		W4 = (W4 + W6) >> 1;
		W3 = W3 + (W5 >> 1);
		W6 = W6 - W4;
		W5 = W5 - W3 - W1;
		W4 = W4 - W8 - W0;
		W3 = W3 - W7;
		W2 = W2 - W8 - W1 - W7 + W4 + W5;
		W6 = W6 - W2;

		W = W8 * x4 * x4 + W7 * x3 * x4 + W6 * x3 * x3 +
			W5 * x * x4 + W4 * x4 + W3 * x3 + W2 * x2 + W1 * x + W0;
		
		if (W != 0)
			Wp = W;

		if (Wp != 0)
			break;
	}

	return W;
}


/*
   Toom–Cook Multiplication (Bodrato-inspired implementation)
   Copyright (C) 2026 James Pate Williams, Jr.

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/

#include "pch.h"
#include <iostream>
#include "Multiplication.h"

int main()
{
    bool balanced = false, square = true;

    if (balanced && !square)
    {
        long long U = 7;
        long long V = 8;
        long long W = 0;
        // W = U * V = 56
        W = Multiplication::BalancedToomCook3Way(10, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 9;
        V = 9;
        W = Multiplication::BalancedToomCook3Way(10, U, V);
        // W = U * V = 81
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 18;
        V = 17;
        W = Multiplication::BalancedToomCook3Way(100, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 306
        U = 123;
        V = 456;
        W = Multiplication::BalancedToomCook3Way(1000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 56,088
        U = 1234;
        V = 5678;
        W = Multiplication::BalancedToomCook3Way(10000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 7,006,652
        U = 12345;
        V = 67890;
        W = Multiplication::BalancedToomCook3Way(100000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 838,102,050
        U = 123456;
        V = 789012;
        W = Multiplication::BalancedToomCook3Way(1000000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 838,102,050
    }

    else if (!balanced && !square)
    {
        long long U = 20;
        long long V = 5;
        long long W;
        // W = U * V = 100
        W = Multiplication::UnbalancedToomCook3Way(100, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 112;
        V = 20;
        W = Multiplication::UnbalancedToomCook3Way(1000, U, V);
        // W = U * V = 2240
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 1234;
        V = 567;
        W = Multiplication::UnbalancedToomCook3Way(10000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 699,678
        U = 12345;
        V = 6789;
        W = Multiplication::UnbalancedToomCook3Way(100000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 83,810,205
        U = 123456;
        V = 78901;
        W = Multiplication::UnbalancedToomCook3Way(1000000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 97,408,265,472
        U = 1234567;
        V = 890123;
        W = Multiplication::UnbalancedToomCook3Way(10000000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 9,740,885,481,741
        U = 12345678;
        V = 9012345;
        W = Multiplication::UnbalancedToomCook3Way(100000000, U, V);
        std::cout << "U = " << U << std::endl;
        std::cout << "V = " << V << std::endl;
        std::cout << "W = " << W << std::endl;
        // W = U * V = 111,263,509,394,910
    }

    else if (square)
    {
        long long U = 8;
        long long W = Multiplication::AsymeticalSquaring(10, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 12;
        W = Multiplication::AsymeticalSquaring(100, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 321;
        W = Multiplication::AsymeticalSquaring(1000, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 1234;
        W = Multiplication::AsymeticalSquaring(10000, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 54321;
        W = Multiplication::AsymeticalSquaring(100000, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
        U = 123456;
        W = Multiplication::AsymeticalSquaring(1000000, U);
        std::cout << "U = " << U << std::endl;
        std::cout << "W = " << W << std::endl;
    }
    return 0;
}

Blog Entry © Wednesday, February 4, 2026, by James Pate Williams, Jr. and the Microsoft Copilot Two Multiplication Algorithms from Wikipedia

[BEGIN COPILOT TEXT]

## Introduction

This small C++ console program implements two classical multiplication algorithms: 

**long multiplication** (the grade‑school method) and **Karatsuba multiplication**, the first sub‑quadratic multiplication algorithm discovered (Anatoly Karatsuba, 1960).

My goal here is not to optimize or modernize the algorithms, but to preserve their structure in a clear, readable form. The implementations follow the versions presented in the Wikipedia “Multiplication” article as faithfully as possible.

The code is intentionally simple and procedural. It uses 1‑based indexing for digit arrays because that mirrors the mathematical notation more closely than 0‑based indexing.

## Long Multiplication

This function implements the classical long multiplication algorithm in an arbitrary base. 

Digits are stored in reverse order (least significant digit first), which simplifies carry propagation.

The implementation below is intentionally direct and mirrors the textbook algorithm step‑by‑step

## Karatsuba Multiplication (Base 10)

This is a minimal, non‑recursive implementation of the Karatsuba method. 

The inputs are split into high and low parts:

x = x_1 * B ^ m + x_0

y = y_1 * B ^ m + y_0

Karatsuba reduces the number of multiplications from 4 to 3 by computing:

z = x * y = z_2 *B ^ 2m + z_1 * B ^ m + z_0

This implementation follows the Wikipedia pseudocode closely.

[END COPILOT TEXT]

#pragma once

class MyAlgorithms
{

public:

	static void LongMultiplication(
		int base, int p, int q, int& pSize,
		std::vector<int> a,
		std::vector<int> b,
		std::vector<int>& product
	);

	static void KaratsubaBase10(
		int x0, int x1, int y0, int y1,
		int B, int m, long long& z);
};

#include "pch.h"
#include "MyAlgorithms.h"

void MyAlgorithms::LongMultiplication(
	int base, int p, int q, int& pSize,
	std::vector<int> a, std::vector<int> b,
	std::vector<int>& product)
{
	pSize = p + q;
	product.resize(pSize + 1LL);

	for (int b_i = 1; b_i <= q; b_i++)
	{
		int carry = 0;

		for (int a_i = 1; a_i <= p; a_i++)
		{
			product[static_cast<long long>(a_i) + b_i - 1LL] +=
				carry + a[a_i] * b[b_i];
			carry = product[static_cast<long long>(a_i) + b_i - 1] / base;
			product[static_cast<long long>(a_i) + b_i - 1] =
				product[static_cast<long long>(a_i) + b_i - 1] % base;
		}

		product[static_cast<long long>(b_i) + p] = carry;
	}
}

void MyAlgorithms::KaratsubaBase10(
	int B, int m, int x0, int x1,
	int y0, int y1, long long& z)
{
	int pb = static_cast<int>(pow(B, m));
	int x = x1 * pb + x0;
	int y = y1 * pb + y0;
	int z2 = x1 * y1;
	int z1 = x1 * y0 + x0 * y1;
	int z0 = x0 * y0;

	z = z2 * static_cast<int>(pow(B, 2 * m)) + z1 * pb + z0;
}


#include "MyAlgorithms.h"

static void DoLongMultiplication()
{
	int base = 0, p = 0, q = 0, pSize = 0;
	char line[128] = "";
	char inputaStr[128] = "";
	char inputbStr[128] = "";
	char* aReverseStr = nullptr;
	char* bReverseStr = nullptr;
	std::cout << "Enter base = ";
	std::cin.getline(line, 128);
	base = atoi(line);
	std::cout << "a = ";
	std::cin.getline(inputaStr, 128);
	std::cout << "b = ";
	std::cin.getline(inputbStr, 128);
	aReverseStr = _strrev(inputaStr);
	bReverseStr = _strrev(inputbStr);
	p = static_cast<int>(strlen(aReverseStr));
	q = static_cast<int>(strlen(bReverseStr));
	pSize = p + q;
	std::vector<int> a(p + 1);
	std::vector<int> b(q + 1);
	std::vector<int> ab(p + q + 1);
	std::vector<int> product;

	for (int i = 1; i <= p; i++)
		a[i] = aReverseStr[i - 1] - '0';

	for (int i = 1; i <= q; i++)
		b[i] = bReverseStr[i - 1] - '0';

	MyAlgorithms::LongMultiplication(
		base, p, q, pSize, a, b, ab);

	size_t i = ab.size() - 1, j = 1;

	while (i >= 0)
	{
		if (ab[i] == 0)
			i--;
		else
			break;
	}

	product.push_back(0);

	for (j = i; j >= 1; j--)
		product.push_back(ab[j]);

	std::cout << "product = ";

	for (int i = 1; i < product.size(); i++)
		std::cout << product[i];

	std::cout << std::endl;
}

static void DoKaratsuba()
{
	char line[128] = "";
	std::cout << "Enter base = ";
	std::cin.getline(line, 128);
	int B = atoi(line);
	std::cout << "Enter m = ";
	std::cin.getline(line, 128);
	int m = atoi(line);
	std::cout << "x1 = ";
	std::cin.getline(line, 128);
	int x1 = atoi(line);
	std::cout << "x0 = ";
	std::cin.getline(line, 128);
	int x0 = atoi(line);
	std::cout << "y1 = ";
	std::cin.getline(line, 128);
	int y1 = atoi(line);
	std::cout << "y0 = ";
	std::cin.getline(line, 128);
	int y0 = atoi(line);
	long long z = 0;
	MyAlgorithms::KaratsubaBase10(
		B, m, x0, x1, y0, y1, z);
	std::cout << "z = " << z << std::endl;
}

int main()
{
	while (true)
	{
		char line[128] = "";
		std::cout << "== Menu ==" << std::endl;
		std::cout << "1 Long Multiplication" << std::endl;
		std::cout << "2 Karatsuba Multiplication" << std::endl;
		std::cout << "3 Exit" << std::endl;
		std::cout << "Option (1 or 2 or 3) = ";
		std::cin.getline(line, 128);
		char option = line[0];

		if (option == '1')
		{
			DoLongMultiplication();
		}

		else if (option == '2')
		{
			DoKaratsuba();
		}

		else
			break;
	}

	return 0;
}

== Menu ==
1 Long Multiplication
2 Karatsuba Multiplication
3 Exit
Option (1 or 2 or 3) = 1
Enter base = 10
a = 506
b = 208
product = 105248
== Menu ==
1 Long Multiplication
2 Karatsuba Multiplication
3 Exit
Option (1 or 2 or 3) = 2
Enter base = 10
Enter m = 2
x1 = 5
x0 = 6
y1 = 2
y0 = 8
z = 105248
== Menu ==
1 Long Multiplication
2 Karatsuba Multiplication
3 Exit
Option (1 or 2 or 3) = 3

D:\Multiplication\x64\Debug\Multiplication.exe (process 30912) exited with code 0 (0x0).
Press any key to close this window . . .

Blog Entry (c) Tuesday February 3, 2026, by James Pate Williams, Jr. Derivation of the Classical Kinetic Energy Formula

Derivation of the Classical Kinetic Energy Equation from EinsteinDownload
// KineticEnergy.cpp
// Author: James Pate Williams, Jr. and NIST
// Copyright Monday February 2, 2026

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

const double c = 2.99792458E8;
const double mass0 = 1.67492750056E-27;

static double KineticEnergy(
    double v, size_t number,
    double& kilotonsTNT,
    std::vector<double>& terms)
{
    double c2 = c * c, ke = 0;;
    
    if (number >= 0)
        terms[0] = mass0 * c2;
    if (number >= 1)
        terms[1] = 0.5 * mass0 * v * v;
    if (number >= 2)
        terms[2] = 3.0 * mass0 * pow(v, 4.0) / (8.0 * c2);
    if (number >= 3)
        terms[3] = 225.0 * mass0 * pow(v, 6.0) / (720.0 * c2 * c2);

    for (size_t i = 0; i <= number; i++)
        ke += terms[i];

    kilotonsTNT = 2.3900573613767E-13 * ke;
    return ke;
}

int main()
{
    char line[128] = "";

    while (true)
    {
        double ke = 0.0, kilotonsTNT = 0.0, v = 0.0;
        size_t number = 0;

        std::cout << "Enter v = ";
        std::cin.getline(line, 128);
        v = atof(line);

        if (v == 0)
            break;

        std::cout << "# Terms = ";
        std::cin.getline(line, 128);
        number = static_cast<size_t>(atoi(line));
        std::vector<double> terms(number + 1);
        ke = KineticEnergy(v, number, kilotonsTNT, terms);
        std::cout << "KE = " << ke << " joules" << std::endl;
        std::cout << "KE = " << kilotonsTNT << " kilotons TNT";
        std::cout << std::endl;
        std::cout << std::scientific;

        for (size_t i = 0; i <= number; i++)
        {
            std::cout << "KE[" << i << "] = " << terms[i];
            std::cout << " joules" << std::endl;
            kilotonsTNT = 2.3900573613767E-13 * terms[i];
            std::cout << "KE[" << i << "] = " << kilotonsTNT;
            std::cout << " kilotons TNT" << std::endl;
        }
    }
}