Rexx Emulation by James Pate Williams, Jr.

Rexx was a popular computer programming language on IBM main frame computers. I found an interesting website: https://en.wikipedia.org/wiki/Rexx#NUMERIC\

I use a very slowly converging infinite series and Fourier series to compute Pi.

// https://en.wikipedia.org/wiki/Rexx#NUMERIC\
// Translated from Rexx code
// by James Pate Williams, Jr.
// Added a lot of functionality
// Copyrighted February 7 - 10, 2023

#include "FourierSeries.h" 
#include <stdlib.h>
#include <iomanip>
#include <iostream>
#include <chrono>
using namespace std;
typedef chrono::high_resolution_clock Clock;

void Option3(int digits, int N)
{
    long double i = 0;
    long double sgn = 1;
    long double sum = 0;
    long double pi0 = Pi;
    long double pi1 = 0.0;

    while (i <= N)
    {
        sum += sgn / (2.0 * i + 1);
        sgn *= -1;
        i = i + 1.0;
    }

    pi1 = 4.0 * sum;
    cout << setw(static_cast<std::streamsize>((int)digits) + 2);
    cout << setprecision(digits) << pi1 << '\t';
    cout << setprecision(digits) << (fabsl(pi0 - pi1)) << endl;
}

void Option4(int digits, int N)
{
    long double* a = new long double[N + 1];
    long double* b = new long double[N + 1];
    memset(a, 0, (N + 1) * sizeof(long double));
    memset(b, 0, (N + 1) * sizeof(long double));
    FourierSeries::CreateCosSeries(N, a);
    FourierSeries::CreateSinSeries(N, b);
    long double pi = 4.0 * FourierSeries::Series(N, 1.0, a, b);
    cout << setw(static_cast<std::streamsize>((int)digits) + 2);
    cout << setprecision(digits) << pi << '\t';
    cout << setprecision(digits) << (fabsl(pi - Pi)) << endl;
    delete[] a;
    delete[] b;
}

int main()
{
    int choice, digits, N = 0, N1 = 0;

    while (1)
    {
        cout << "==Menu==" << endl;
        cout << "1 Compute Sqrt(2)" << endl;
        cout << "2 Compute Exp(1)" << endl;
        cout << "3 Compute Pi" << endl;
        cout << "4 Compute Pi (Fourier Series)" << endl;
        cout << "5 Generate Option 3 Table" << endl;
        cout << "6 Generate Option 4 Table" << endl;
        cout << "7 Exit" << endl;
        cin >> choice;
        if (choice == 7)
            break;
        cout << "  Digits = ";
        cin >> digits;
        if (choice >= 3)
        {
            cout << "  Terms  = ";
            cin >> N;
            N1 = N + 1;
        }
        auto start_time = Clock::now();
        if (choice == 1)
        {
            long double n = 2;
            long double r = 1;
            while (1)
            {
                long double rr = (n / r + r) / 2;
                if (rr == r)
                    break;
                r = rr;
            }

            cout << setprecision(digits) << r << endl;
        }
        else if (choice == 2)
        {
            long double e = 2.5;
            long double f = 0.5;
            long double n = 3;
            do
            {
                f /= n;
                long double ee = e + f;
                if (ee == e)
                    break;
                e = ee;
                n += 1;
            } while (1);

            cout << setprecision(digits) << e << endl;
        }
        else if (choice == 3)
            Option3(digits, N);
        else if (choice == 4)
            Option4(digits, N);
        else if (choice == 5)
        {
            for (int n = 8; n <= N; n *= 2)
                Option3(digits, n);
        }
        else if (choice == 6)
        {
            for (int n = 8; n <= N; n *= 2)
                Option4(digits, n);
        }
        auto end_time = Clock::now();
        cout << "runtime in milliseconds = ";
        cout << std::chrono::duration_cast<std::chrono::milliseconds>
            (end_time - start_time).count();
        cout << endl;
        cout << "runtime in nanoseconds  = ";
        cout << std::chrono::duration_cast<std::chrono::nanoseconds>
            (end_time - start_time).count();
        cout << endl;
    }

    return 0;
}
#pragma once

const long double c = 2.0e9;
const long double Pi = 3.1415926535897932384626433832795;

class FourierSeries
{
private:
	static long double f(long double x);
	static long double cosTermFunction(int n, long double x);
	static long double sinTermFunction(int n, long double x);
public:
	static void CreateCosSeries(int N, long double a[]);
	static void CreateSinSeries(int N, long double b[]);
	static long double Series(int N, long double x,
		long double a[], long double b[]);
};
#include <math.h>
#include "FourierSeries.h"
#include "Integral.h"

long double FourierSeries::f(long double x)
{
    return atanl(x);
}

long double FourierSeries::cosTermFunction(int n, long double x)
{
    return cosl(n * x * Pi / c) * f(x);
}

long double FourierSeries::sinTermFunction(int n, long double x)
{
    return sinl(n * x * Pi / c) * f(x);
}

void FourierSeries::CreateCosSeries(int N, long double a[])
{
    long double e[7] = { 0 };

    e[1] = e[2] = 1.0e-12;

    for (int n = 0; n <= N; n++)
        a[n] = Integral::integral(
            n, -c, c, cosTermFunction, e, 1, 1) / c;
}

void FourierSeries::CreateSinSeries(int N, long double b[])
{
    long double e[7] = { 0 };

    e[1] = e[2] = 1.0e-12;

    for (int n = 1; n <= N; n++)
        b[n] = Integral::integral(
            n, -c, c, sinTermFunction, e, 1, 1) / c;
}

long double FourierSeries::Series(int N, long double x,
    long double a[], long double b[])
{
    long double sum = 0.0;

    for (int n = 1; n <= N; n++)
        sum += a[n] * cosl(2.0 * n * x) + b[n] * sinl(2.0 * n * x);

    return a[0] / 2.0 + sum / 2.0;
}
#pragma once

// Translated from C source code found in the tome
// "A Numerical Library in C for Scientists and
// Engineers" by H.T. Lau, PhD

class Integral
{
private:
	static long double integralqad(
		int n,
		int transf, long double (*fx)(int, long double), long double e[],
		long double* x0, long double* x1, long double* x2, long double* f0, long double* f1,
		long double* f2, long double re, long double ae, long double b1);
	static void integralint(
		int n,
		int transf, long double (*fx)(int, long double), long double e[],
		long double* x0, long double* x1, long double* x2, long double* f0, long double* f1,
		long double* f2, long double* sum, long double re, long double ae, long double b1,
		long double hmin);
public:
	static long double integral(int n, long double a, long double b,
		long double (*fx)(int, long double), long double e[],
		int ua, int ub);
};
#include <math.h>
#include "Integral.h"

long double Integral::integral(
	int n,
	long double a, long double b,
	long double (*fx)(int, long double), long double e[],
	int ua, int ub)
{
	long double x0, x1, x2, f0, f1, f2, re, ae, b1 = 0, x;

	re = e[1];
	if (ub)
		ae = e[2] * 180.0 / fabsl(b - a);
	else
		ae = e[2] * 90.0 / fabsl(b - a);
	if (ua) {
		e[3] = e[4] = 0.0;
		x = x0 = a;
		f0 = (*fx)(n, x);
	}
	else {
		x = x0 = a = e[5];
		f0 = e[6];
	}
	e[5] = x = x2 = b;
	e[6] = f2 = (*fx)(n, x);
	e[4] += integralqad(n, 0, fx, e, &x0, &x1, &x2, &f0, &f1, &f2, re, ae, b1);
	if (!ub) {
		if (a < b) {
			b1 = b - 1.0;
			x0 = 1.0;
		}
		else {
			b1 = b + 1.0;
			x0 = -1.0;
		}
		f0 = e[6];
		e[5] = x2 = 0.0;
		e[6] = f2 = 0.0;
		ae = e[2] * 90.0;
		e[4] -= integralqad(n, 1, fx, e, &x0, &x1, &x2, &f0, &f1, &f2, re, ae, b1);
	}
	return e[4];
}

long double Integral::integralqad(
	int n,
	int transf, long double (*fx)(int, long double), long double e[],
	long double* x0, long double* x1, long double* x2, long double* f0, long double* f1,
	long double* f2, long double re, long double ae, long double b1)
{
	/* this function is internally used by INTEGRAL */

	long double sum, hmin, x, z;

	hmin = fabs((*x0) - (*x2)) * re;
	x = (*x1) = ((*x0) + (*x2)) * 0.5;
	if (transf) {
		z = 1.0 / x;
		x = z + b1;
		(*f1) = (*fx)(n, x) * z * z;
	}
	else
		(*f1) = (*fx)(n, x);
	sum = 0.0;
	integralint(n, transf, fx, e, x0, x1, x2, f0, f1, f2, &sum, re, ae, b1, hmin);
	return sum / 180.0;
}

void Integral::integralint(
	int n,
	int transf, long double (*fx)(int, long double), long double e[],
	long double* x0, long double* x1, long double* x2, long double* f0, long double* f1,
	long double* f2, long double* sum, long double re, long double ae, long double b1,
	long double hmin)
{
	/* this function is internally used by INTEGRALQAD of INTEGRAL */

	int anew;
	long double x3, x4, f3, f4, h, x, z, v, t;

	x4 = (*x2);
	(*x2) = (*x1);
	f4 = (*f2);
	(*f2) = (*f1);
	anew = 1;
	while (anew) {
		anew = 0;
		x = (*x1) = ((*x0) + (*x2)) * 0.5;
		if (transf) {
			z = 1.0 / x;
			x = z + b1;
			(*f1) = (*fx)(n, x) * z * z;
		}
		else
			(*f1) = (*fx)(n, x);
		x = x3 = ((*x2) + x4) * 0.5;
		if (transf) {
			z = 1.0 / x;
			x = z + b1;
			f3 = (*fx)(n, x) * z * z;
		}
		else
			f3 = (*fx)(n, x);
		h = x4 - (*x0);
		v = (4.0 * ((*f1) + f3) + 2.0 * (*f2) + (*f0) + f4) * 15.0;
		t = 6.0 * (*f2) - 4.0 * ((*f1) + f3) + (*f0) + f4;
		if (fabsl(t) < fabsl(v) * re + ae)
			(*sum) += (v - t) * h;
		else if (fabsl(h) < hmin)
			e[3] += 1.0;
		else {
			integralint(n, transf, fx, e, x0, x1, x2, f0, f1, f2, sum,
				re, ae, b1, hmin);
			*x2 = x3;
			*f2 = f3;
			anew = 1;
		}
		if (!anew) {
			*x0 = x4;
			*f0 = f4;
		}
	}
}

Backpropagation Artificial Neural Networks for Approximating Simple Functions by James Pate Williams, Jr.

The functions to be calculated are the exclusive or (XOR) function, the eight input and output identity function, and four non-linear three input and one output functions. I used the logistic function that represents 0 as 0.1 and 1 as 0.9. The three input functions use 8 hidden units, a tolerance of 1.0e-9, and maximum number of epochs as 1,000,000. The backpropagation feed-forward neural network algorithm is from Tom M. Mitchell’s classic textbook “Machine Learning” which has a copyright of 1997.

Georgia Climate Data by James Pate Williams, Jr.

Abstract

I became interested in attempting to predict temperatures in the state of Georgia way back in the day, i.e., the date was November 1 – 4, 2009. I found a neat website for the temperatures from January 1895 to December 2001. Unfortunately, the website no longer exists online, however, I saved the data to an extinct PC of mine and a USB solid state drive. I used two methods to attempt predictions of the annual temperatures from January 2002 to December 2025. The first algorithm is polynomial least squares curve fitting [1]. The second method is a radial basis function neural network [2] [3] that is trained by an evolutionary hill-climber of my design and implementation [4] [5]. The applications were written in one of my favorite computer programming languages, namely, C# [6].

Methodologies

I created a polynomial least squares dynamic-link library (DLL) using Gaussian elimination with pivoting and an inverse matrix calculation utilizing an upper-triangular matrix [7]. At first the driver application was capable of fitting a 1-degree to 100-degrees polynomial. I found that my algorithm was only valid for a 1-degree to 76-degrees fitting polynomial. I used degrees: 1, 4, 8, 16, 32, 64, and 76.

Predicted Statewide Georgia Annual Average Temperatures in Degrees Fahrenheit
Poly Degree14816326476Model
YearT (F)T (F)T (F)T (F)T (F)T (F)T (F)Average
200264.263.163.863.662.463.062.963.3
200364.263.263.663.464.163.763.763.7
200464.263.363.463.464.563.964.063.8
200564.263.363.463.564.263.964.063.8
200664.263.463.463.663.863.763.863.7
200764.263.563.463.763.363.663.663.6
200864.263.663.463.763.063.463.463.5
200964.263.663.463.762.963.363.363.5
201064.163.763.563.762.963.363.263.5
201164.163.863.663.663.163.363.263.5
201264.163.863.763.663.363.463.363.6
201364.163.963.763.663.663.563.563.7
201464.163.963.863.663.963.663.663.8
201564.164.063.963.664.163.863.863.9
201664.164.063.963.664.363.963.964.0
201764.164.164.063.764.464.064.164.1
201864.164.164.163.764.464.164.264.1
201964.164.264.163.964.464.264.264.2
202064.164.264.264.064.364.364.364.2
202164.164.364.264.164.264.364.364.2
202264.164.364.264.264.264.364.364.2
202364.164.364.364.364.164.364.364.2
202464.064.464.364.564.064.364.364.3
202564.064.464.364.664.064.364.364.3

The second method was a radial basis function neural network that was trained by an evolutionary hill-climber of my design and implementation. I used 8, 16, 32, 64 basis functions. The population of the hill-climber was 16 and generations 262,144.

Predicted Statewide Georgia Annual Average Temperatures in Degrees Fahrenheit
Basis8163264Model
YearT (F)T (F)T (F)T (F)Average
200265.464.265.465.265.1
200365.564.265.465.265.1
200465.564.265.565.365.1
200565.564.265.565.365.1
200665.664.265.565.365.2
200765.664.265.665.365.2
200865.664.265.665.465.2
200965.764.265.665.465.2
201065.764.265.665.465.2
201165.764.265.765.465.3
201265.764.265.765.565.3
201365.864.265.765.565.3
201465.864.265.865.565.3
201565.864.365.865.565.4
201665.964.365.865.665.4
201765.964.365.865.665.4
201865.964.365.965.665.4
201966.064.365.965.665.5
202066.064.365.965.765.5
202166.064.366.065.765.5
202266.064.366.065.765.5
202366.164.366.065.765.5
202466.164.366.165.865.6
202566.164.366.165.865.6

I trust the polynomial fitting results more than the radial basis function neural network values. The polynomial fitting had mean square errors < 1 whereas the radial basis function neural network had mean square errors between 1 and 9. The general trends were increasing temperatures which are in line with the theorized global warming.

References

[1]H. T. Lau, A Numerical Library in C for Scientists and Engineers, Boca Raton: CRC Press, 1995.
[2]T. M. Mitchell, Machine Learning, Boston: WCB McGraw-Hill, 1997.
[3]A. P. Engelbrecht, Computational Intelligence An Introduction, Hoboken: John Wiley and Sons, 2002.
[4]Z. Michalewicz, Genetic Algorithms + Data Structures = Evolutionary Programs 3rd Edition, Berlin: Springer, 1999.
[5]D. B. Fogel, Evolutionary Computation Toward a New Philosophy of Machine Intelligence, Piscataway: IEEE Press, 2000.
[6]C. Petzold, Programming Windows with C#, Redmond: Microsoft Press, 2002.
[7]S. D. Conte and C. d. Boor, Elementary Numerical An Algorithmic Approach, New York: McGraw-Hill Book Company, 1980.

Yet Another Revisitation of Reproducing Ordnance Pamphlet 770 by James Pate Williams, Jr.

This is another attempt to reproduce the United States Navy’s Ordnance Pamphlet 770: https://eugeneleeslover.com/USN-GUNS-AND-RANGE-TABLES/OP-770-1.html which contains ballistic tables for the battleship USS Iowa (BB-61) artillery (16-inch/50 caliber) and is dated October 1941. My C# Windows desktop application is capable of calculating the elevation from range table which has the columns range in yards, angle of elevation in degrees and minutes, positive angle of fall in degrees and minutes, time of flight in seconds, apogee also called summit in feet, striking velocity in feet per second, and energy in foot pound force. Three corrections can be applied to the trajectory: trunnion height in feet, acceleration of gravity correction, and the curvature of the Earth correction (Vincenty calculation). The first image below is the ballistic settings interface. The second image is the uncorrected table. The third image is the application of a trunnion height of 32 feet. The fourth image is the curvature of the Earth correction. The fifth image is the trunnion height of 32 feet and Vincenty corrections. It is to be noticed that the striking velocity and kinetic energy are the only non-monotonically increasing or decreasing data fields.

Georgia Climate Data – Temperatures in Degrees Fahrenheit by James Pate Williams, Jr.

I became interested in attempting to predict average annual temperatures in the state of Georgia way back in the day. I found a neat website for temperatures for January 1895 to December 2001. Unfortunately, the website no longer exists online, however, I saved the data to an extinct PC of mine and a USB solid state drive. I used several methods to attempt predictions of the annual temperatures from 2002 to 2025. The first algorithm is polynomial least squares utilizing a 32-degree polynomial and a 75-degree polynomial. The application was written in one of my favorite programming languages, namely, C#.

Parabolic Cylinder Functions by James Pate Williams, Jr.

Solution to the differential equation: d2y/dx2 = -0.25 * x * x * y – a * y which is valid for all real and complex numbers. We examine the real solutions to the second order ordinary differential equations. See “Handbook of Mathematical Functions” by Milton Abramowitz and Irene A. Stegun, Chapter 19. Parabolic Cylinder functions. Page 686 Equation 19.2.7 Recurrence formula and 19.2.5 Series solutions.

Win32 C Orthogonal Polynomials

I learned about the Laguerre and Legendre polynomials when I first read “Introduction to Quantum Mechanics” by Pauling and Wilson way back in the early 1970s. I later learned about the Chebyshev and other orthogonal polynomials. Beginning on March 30, 2015, I created yet another application to graph various orthogonal polynomials in C#. A few days ago, I wrote a Win32 C application to graph the Chebyshev, Laguerre, and Legendre orthogonal polynomials. At a later date, I will probably add Hermite and Jacobi polynomials.


Win32 C Romberg Extrapolation by James Pate Williams, Jr.

// Romberg.cpp : Defines the entry point for the application.
// https://learn.microsoft.com/en-us/windows/win32/controls/create-a-simple-combo-box

#include "stdafx.h"
#include <CommCtrl.h>
#include <math.h>
#include <stdio.h>
#include "Romberg.h"

#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

// 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    RombergDialog(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_ROMBERG, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, NULL, 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;

    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_ROMBERG));
    wcex.hCursor        = LoadCursor(NULL, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_ROMBERG);
    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, NULL, NULL, hInstance, NULL);

   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 ID_OPTION_ROMBERGINTEGRATION:
				DialogBox(hInst, MAKEINTRESOURCE(IDD_DIALOG1), hWnd, RombergDialog);
				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;
}

double f0(double x)
{
	return exp(-x * x);
}

double f1(double x)
{
	return x * x;
}

double f2(double x)
{
	double pi = 4.0 * atan(1.0);
	return sin(101.0 * pi * x);
}

double f3(double x)
{
	double pi = 4.0 * atan(1.0);
	return 1.0 + sin(10.0 * pi * x);
}

double f4(double x)
{
	return fabs(x - 1.0 / 3.0);
}

double f5(double x)
{
	return sqrt(x);
}

double f6(double x)
{
	double result = 1.0;

	if (x != 0.0)
		result = sin(x) / x;

	return result;
}

/* https://en.wikipedia.org/wiki/Romberg%27s_method#Implementation */

char table[8192];
double R1[2048], R2[2048];

void print_row(size_t i, double* R)
{
	char buffer[128];
	sprintf_s(buffer, 128, "R[%2zu] = ", i);
	strcat_s(table, 8192, buffer);

	for (size_t j = 0; j <= i; ++j)
	{
		sprintf_s(buffer, 128, "%.*lf ", 16, R[j]);
		strcat_s(table, 8192, buffer);
	}

	strcat_s(table, 8192, "\r\n");
}

double RombergExtrapolation(
	double(*f)(double),
	double a, double b,
	size_t max_steps,
	double acc)
{
	double *Rp = &R1[0], *Rc = &R2[0];
	double h = (b - a);
	Rp[0] = (f(a) + f(b))*h*.5;

	print_row(0, Rp);

	for (size_t i = 1; i < max_steps; ++i)
	{
		h /= 2.;
		double c = 0;
		size_t ep = 1 << (i - 1);
		
		for (size_t j = 1; j <= ep; ++j)
		{
			c += f(a + (2 * j - 1)*h);
		}

		Rc[0] = h*c + .5*Rp[0];

		for (size_t j = 1; j <= i; ++j)
		{
			double n_k = pow(4, j);
			Rc[j] = (n_k*Rc[j - 1] - Rp[j - 1]) / (n_k - 1);
		}

		print_row(i, Rc);

		if (i > 1 && fabs(Rp[i - 1] - Rc[i]) < acc)
		{
			return Rc[i];
		}

		double *rt = Rp;
		Rp = Rc;
		Rc = rt;
	}

	return Rp[max_steps - 1];
}

INT_PTR CALLBACK RombergDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
	UNREFERENCED_PARAMETER(lParam);
	WORD word = LOWORD(wParam);
	static TCHAR functions[7][128] =
	{
		TEXT("f(x) = exp(-x * x)"),
		TEXT("f(x) = x * x"),
		TEXT("f(x) = sin(101 * pi * x)"),
		TEXT("f(x) = 1 + sin(10 * pi * x)"),
		TEXT("f(x) = abs(x - 1.0 / 3.0)"),
		TEXT("f(x) = sqrt(x)"),
		TEXT("f(x) = sin(x) / x or 1 if x = 0")
	};
	static TCHAR A[1024];
	static HWND cbWindow;
	static int k = 0;
	static int itemIndex;

	switch (message)
	{
	case WM_INITDIALOG:
		cbWindow = GetDlgItem(hDlg, IDC_COMBO1);
		memset(&A, 0, sizeof(A));
		for (k = 0; k < 6; k++)
		{
			// Add string to combobox.
			SendMessage(cbWindow, CB_ADDSTRING, (WPARAM)k, (LPARAM)functions[k]);
		}
		// Send the CB_SETCURSEL message to display an initial item 
		// in the selection field  
		SendMessage(cbWindow, CB_SETCURSEL, (WPARAM)0, (LPARAM)0);
		return (INT_PTR)TRUE;
	case WM_COMMAND:
		if (word == IDOK || word == IDCANCEL)
		{
			EndDialog(hDlg, word);
			return (INT_PTR)TRUE;
		}
		if (HIWORD(wParam) == CBN_SELCHANGE)
			// If the user makes a selection from the list:
			//   Send CB_GETCURSEL message to get the index of the selected list item.
			//   Send CB_GETLBTEXT message to get the item.
			//   Display the item in a messagebox.
		{
			itemIndex = SendMessage((HWND)lParam, (UINT)CB_GETCURSEL,
				(WPARAM)0, (LPARAM)0);
			TCHAR listItem[256];
			(TCHAR)SendMessage((HWND)lParam, (UINT)CB_GETLBTEXT,
				(WPARAM)itemIndex, (LPARAM)listItem);
		}
		else if (word == IDC_BUTTON1)
		{
			char buffer[128];
			GetDlgItemTextA(hDlg, IDC_EDIT1, buffer, 128);
			double a = atof(buffer);
			GetDlgItemTextA(hDlg, IDC_EDIT2, buffer, 128);
			double b = atof(buffer);
			double integral = 0.0;
			GetDlgItemTextA(hDlg, IDC_EDIT4, buffer, 128);
			double acc = atof(buffer);
			int max_steps = GetDlgItemInt(hDlg, IDC_EDIT5, FALSE, FALSE);

			switch (itemIndex)
			{
			case 0:
				integral = RombergExtrapolation(f0, a, b, max_steps, acc);
				break;
			case 1:
				integral = RombergExtrapolation(f1, a, b, max_steps, acc);
				break;
			case 2:
				integral = RombergExtrapolation(f2, a, b, max_steps, acc);
				break;
			case 3:
				integral = RombergExtrapolation(f3, a, b, max_steps, acc);
				break;
			case 4:
				integral = RombergExtrapolation(f4, a, b, max_steps, acc);
				break;
			case 5:
				integral = RombergExtrapolation(f5, a, b, max_steps, acc);
				break;
			case 6:
				integral = RombergExtrapolation(f6, a, b, max_steps, acc);
				break;
			}

			sprintf_s(buffer, 128, "Integral = %.*lf\r\n", 16, integral);
			strcat_s(table, 8192, buffer);
			SetDlgItemTextA(hDlg, IDC_EDIT7, table);
		}

		else if (word == IDC_BUTTON3)
		{
			table[0] = '\0';
			SetDlgItemTextA(hDlg, IDC_EDIT7, table);
		}
	}
	return (INT_PTR)FALSE;
}