Category: C Applications

Blog Entry © Wednesday, May 27, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Grade School Arithmetic

Technical Report Grade School ArithmeticDownload
#pragma once
#include <stdint.h>

/* Algorithm due to Microsft's Coilot 
function udiv_restoring(N, D, n) :
    R = 0
    Q = 0

    negD = (~D + 1)

    for i from n - 1 down to 0
    {
        R = (R << 1) | ((N >> i) & 1)

    T = R + negD

    if MSB(T) == 0:
        R = T
        Q = Q | (1 << i)

    return (Q, R) 
 */

class Arithmetic
{
public:
    static bool udiv_restoring(
        uint32_t numer,
        uint32_t denom,
        uint32_t& quo,
        uint32_t& rem,
        int n);
    static bool umul_shift_add(
        uint32_t a,
        uint32_t b,
        uint64_t& product,
        int n);
};


#include <cstdint>
#include "Arithmetic.h"

static inline uint32_t mask_n(int bits) {
    return (bits >= 32) ? 0xFFFFFFFFu : ((1u << bits) - 1u);
}

static inline uint32_t msb(uint32_t x, int bits) {
    // returns top bit of a 'bits'-wide value
    return (x >> (bits - 1)) & 1u;
}

bool Arithmetic::udiv_restoring(
    uint32_t numer,
    uint32_t denom,
    uint32_t& quo,
    uint32_t& rem,
    int n)
{
    if (denom == 0 || n <= 0 || n > 32) return false;
    
    if (numer == 0)
    {
        quo = rem = 0;
        return true;
    }

    quo = 0;
    rem = 0;

    if (n == 32) {
        uint64_t R = 0;
        uint64_t D = (uint64_t)denom;
        uint64_t maskW = (1ull << 33) - 1ull;           // 33-bit mask
        uint64_t negD = ((~D) + 1ull) & maskW;          // 33-bit two's complement

        for (int i = n - 1; i >= 0; --i) {
            R = ((R << 1) | ((numer >> i) & 1u)) & maskW;

            uint64_t T = (R + negD) & maskW;         // R - D

            // Sign bit is bit 32 (the 33rd bit)
            if (((T >> 32) & 1ull) == 0ull) {
                R = T;
                quo |= (1u << i);
            }
        }
        rem = (uint32_t)(R & 0xFFFFFFFFu);
        return true;
    }

    // n < 32 case: we can keep everything in uint32_t using (n+1) bits
    uint32_t maskN = mask_n(n);
    uint32_t maskW = mask_n(n + 1);

    uint32_t N = numer & maskN;
    uint32_t D = denom & maskN;

    // Two's complement of D in (n+1) bits
    uint32_t Dw = D;                  // placed in low bits of (n+1)-wide register
    uint32_t negD = ((~Dw) + 1u) & maskW;

    uint32_t R = 0;

    for (int i = n - 1; i >= 0; --i) {
        R = ((R << 1) | ((N >> i) & 1u)) & maskW;

        uint32_t T = (R + negD) & maskW;            // trial subtract: R - D (in w bits)

        if (msb(T, n + 1) == 0) {                   // non-negative in (n+1) bits
            R = T;
            quo |= (1u << i);
        }
    }

    rem = R & maskN;                                // remainder fits in n bits
    return true;
}

bool Arithmetic::umul_shift_add(
    uint32_t a,
    uint32_t b,
    uint64_t& product,
    int n)
{
    if (n <= 0 || n > 32) return false;

    uint64_t A = a;    // promote to avoid overflow
    uint32_t B = b;

    product = 0;

    for (int i = 0; i < n; ++i) {
        if (B & 1u) {
            product += A;
        }

        A <<= 1;
        B >>= 1;
    }

    return true;
}


#include <chrono>
#include <cstdint>
#include <iostream>
#include <limits>
#include <random>
#include <string>
#include "Arithmetic.h"

namespace {

    constexpr int TESTS_PER_N = 200000;

    uint32_t make_mask(int n) {
        return (n == 32) ? 0xFFFFFFFFu : ((1u << n) - 1u);
    }

    void clear_bad_input() {
        std::cin.clear();
        std::cin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
    }

    template <typename TrialFn>
    double run_suite(const char* label, TrialFn trial, bool verbose) {
        std::mt19937 rng(12345); // deterministic

        auto t0 = std::chrono::high_resolution_clock::now();

        for (int n = 1; n <= 32; ++n) {
            const uint32_t mask = make_mask(n);

            for (int i = 0; i < TESTS_PER_N; ++i) {
                if (!trial(rng, mask, n)) {
                    std::cout << label << ": FAILED (n=" << n << ", i=" << i << ")\n";
                    return -1.0;
                }
            }

            if (verbose) {
                std::cout << "n=" << n << " passed\n";
            }
        }

        auto t1 = std::chrono::high_resolution_clock::now();
        double secs = std::chrono::duration<double>(t1 - t0).count();

        std::cout << label << " runtime = " << secs << " sec\n";
        return secs;
    }

    bool trial_division(std::mt19937& rng, uint32_t mask, int n) {
        const uint32_t numer = rng() & mask;
        const uint32_t denom = (rng() & mask) | 1u; // ensure non-zero

        uint32_t q = 0, r = 0;
        if (!Arithmetic::udiv_restoring(numer, denom, q, r, n)) {
            std::cout << "Failure numer=" << numer << " denom=" << denom << "\n";
            return false;
        }

        const uint32_t q2 = numer / denom;
        const uint32_t r2 = numer % denom;

        if (q != q2 || r != r2) {
            std::cout << "Mismatch n=" << n
                << " numer=" << numer
                << " denom=" << denom
                << " got q=" << q << " r=" << r
                << " expected q=" << q2 << " r=" << r2 << "\n";
            return false;
        }
        return true;
    }

    bool trial_multiplication(std::mt19937& rng, uint32_t mask, int n) {
        const uint32_t a = rng() & mask;
        const uint32_t b = rng() & mask;

        uint64_t prod = 0;
        if (!Arithmetic::umul_shift_add(a, b, prod, n)) {
            std::cout << "Failure a=" << a << " b=" << b << "\n";
            return false;
        }

        const uint64_t expected = static_cast<uint64_t>(a) * static_cast<uint64_t>(b);
        if (prod != expected) {
            std::cout << "Mismatch n=" << n
                << " a=" << a
                << " b=" << b
                << " got=" << prod
                << " expected=" << expected << "\n";
            return false;
        }
        return true;
    }

} // namespace

int main() {
    bool verbose = true;

    while (true) {
        std::cout << "\nArithmetic Lab\n";
        std::cout << "1. Test Division (restoring)\n";
        std::cout << "2. Test Multiplication (shift-add)\n";
        std::cout << "3. Run ALL tests\n";
        std::cout << "4. Toggle verbose (currently " << (verbose ? "ON" : "OFF") << ")\n";
        std::cout << "5. Exit\n";
        std::cout << "Choice: ";

        int choice = 0;
        if (!(std::cin >> choice)) {
            clear_bad_input();
            std::cout << "Invalid input. Please enter a number.\n";
            continue;
        }

        if (choice == 1) {
            run_suite("Division test", trial_division, verbose);
        }
        else if (choice == 2) {
            run_suite("Multiplication test", trial_multiplication, verbose);
        }
        else if (choice == 3) {
            const double d = run_suite("Division test", trial_division, verbose);
            if (d >= 0.0) run_suite("Multiplication test", trial_multiplication, verbose);
        }
        else if (choice == 4) {
            verbose = !verbose;
        }
        else if (choice == 5) {
            return 0;
        }
        else {
            std::cout << "Invalid choice.\n";
        }
    }
}

Blog Entry © Saturday, May 16, 2026, by James Pate Williams, Jr. Some More Linear Algebra Examples

Blog Entry Saturday May 16 2026Download

Blog Entry © Thursday, May 14, 2026, by James Pate Williams, Jr. More Numerical Integration Results

Blog Entry Thursday May 14 2026Download
// NumericalIntegrals.cpp (c) Thursday, May 14, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD

#include <iomanip>
#include <iostream>
#include <vector>
#include <stdlib.h>

static double f(double x) {
    return sin(x);
}

static double MonteCarlo(double a, double b,
    double (*f)(double), int n){ 
    double sum = 0;

    for (int i = 0; i < n; i++) {
        double x = (b - a) * (double)rand() / RAND_MAX;

        sum += f(x);
    }

    return (b - a) * sum / n;
}

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

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

    return factorial;
}

static double Series(double a, double b, int n)
{
    double sumA = 0.0, sumB = 0.0;
    int sign = 1;

    for (int i = 0; i <= n; i++) {
        sumA += sign * pow(a, 2 * i + 2) /
            Factorial(2 * i + 2);

        sign *= -1;
    }

    sign = 1;

    for (int i = 0; i <= n; i++) {      
        sumB += sign * pow(b, 2 * i + 2) /
            Factorial(2 * i + 2);

        sign *= -1;
    }

    return sumB - sumA;
}

static double CompositeTrapezoidalRule(
    double a, double b, int n) {
    double pi = 4.0 * atan(1.0);
    double endPts = 0.5 * (f(a) + f(b));
    double sum = 0, xk = 0.0;
    double h = (b - a) / n;

    for (int k = 1; k <= n - 1; k++) {
        xk = a + k * h;
        sum += f(xk);
    }

    return h * (0.5 * endPts + sum);
}

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;
}

static void Romberg(double a, double b,
    double (*f)(double), int mStart, int nRow,
    std::vector<std::vector<double>>& T) {
    int m = mStart;
    double h = (b - a) / m;
    double sum = 0.5 * (f(a) + f(b));

    if (m > 1) {
        for (int i = 1; i <= m - 1; i++) {
            sum += f(a + i * h);
        }
    }

    T[0][0] = sum * h;

    std::cout << "romberg t-table" << std::endl;
    std::cout << std::fixed;
    std::cout << std::setprecision(5) << T[0][0];
    std::cout << std::endl;

    if (nRow < 2)
        return;

    for (int k = 2; k <= nRow; k++) {
        h /= 2.0;
        m *= 2;
        sum = 0.0;

        for (int i = 1; i <= m; i += 2) {
            sum += f(a + i * h);
        }

        T[k][1] = 0.5 * T[k - 1LL][1] + sum * h;

        for (int j = 1; j <= k - 1; j++) {
            T[k - 1LL][j] = T[k][j] - T[k - 1LL][j];
            T[k][j + 1LL] = T[k][j] + T[k - 1LL][j] /
                (pow(4.0, j) - 1.0);
        }

        for (int j = 1; j <= k; j++) {
            std::cout << std::fixed;
            std::cout << std::setprecision(5);
            std::cout << T[k][j] << '\t';
        }

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

    if (nRow < 3) {
        return;
    }

    std::cout << "table of ratios" << std::endl;
    double ratio = 0.0;

    for (int k = 1; k <= nRow - 2; k++) {
        for (int j = 1; j <= k; j++) {
            if (T[k + 1LL][j] == 0.0) {
                ratio = 0.0;
            }

            else {
                ratio = T[k][j] / T[k + 1LL][j];
            }

            T[k][j] = ratio;
        }

        for (int j = 1; j <= k; j++) {
            std::cout << std::fixed;
            std::cout << std::setprecision(5);
            std::cout << T[k][j] << '\t';
        }

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

double MonteCarloVolume(double R, int n)
{
    double pi = 4.0 * atan(1.0), pi2 = 2.0 * pi;
    double R2 = R * R, sum = 0;

    for (int i = 0; i < n; i++)
    {
        double r = R2 * (double)rand() / RAND_MAX;
        double t = pi * (double)rand() / RAND_MAX;
        double p = pi2 * (double)rand() / RAND_MAX;
        sum += r * r * sin(t);
    }

    return R * pi * pi2 * sum / n;
}

int main()
{
    srand(1);
    std::vector<std::vector<double>> T;
    T.resize(35);
    for (int i = 0; i < 35; i++) {
        T[i].resize(35);
    }
    Romberg(0.0, 2.0, f, 2, 7, T);
    std::cout << std::setprecision(11);
    std::cout << "analytic integral of sine = " << -cos(2.0) + cos(0.0);
    std::cout << std::endl;
    std::cout << "simpson's rule integral   = " << SimpsonsRule(500, 0, 2.0, f);
    std::cout << std::endl;
    std::cout << "monte carlo integral      = " << MonteCarlo(0.0, 2.0, f, 2130);
    std::cout << std::endl;
    std::cout << "infinite series integral  = " << Series(0.0, 2.0, 16);
    std::cout << std::endl;
    double integral = CompositeTrapezoidalRule(0.0, 2.0, 175000000);
    std::cout << "romberg integral          = " << integral << std::endl;
    std::cout << "actual spherical volume   = " << 4.0 * 4.0 * atan(1.0) / 3.0;
    std::cout << std::endl;
    double volume = MonteCarloVolume(1.0, 1000000);
    std::cout << "approx spherical volume   = " << volume;
    std::cout << std::endl;
}

Blog Entry © Wednesday, May 13, 2026, by James Pate Williams, Jr. Adaptive n-Quadrature Versus Monte Carlo Integration

Blog Entry Wednesday May 13 2026Download

Blog Entry © Monday, May 11, 2026, by James Pate Williams, Jr., Laplace Equation in a Solid Cylinder

Blog Entry Monday May 11 2026Download
Laplace Equation in a Solid CylinderDownload

Blog Entry © Tuesday, May 5, 2026, by James Pate Williams, Jr., Solution of the Laplace Equation on a Rectangle

Blog Entry Tuesday May 5 2026Download

Blog Entry © Sunday, April 26, 2026, by James Pate Williams, Jr. Butterfly Curve

Butterfly curve (transcendental) – Wikipedia

// ButterflyEquation.cpp : Defines the entry point for the application.
//

#include "framework.h"
#include "ButterflyEquation.h"
#include <vector>

#define MAX_LOADSTRING 100

typedef struct tagPoint3d
{
    double t, x, y;
} Point3d, * PPoint3d;

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
std::vector<Point3d> 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 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_BUTTERFLYEQUATION, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    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_BUTTERFLYEQUATION));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_BUTTERFLYEQUATION);
    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;
}

void CreateButterflyGraphPoints()
{
    double p = 4.0 * atan(1.0);
    double h = 12.0 * p / 1024.0;
    double t = 0.0;

    for (int i = 0; i <= 1024; i++)
    {
        Point3d pt = { 0 };
        double c = 2.0 * cos(4.0 * t);
        double d = pow(sin(t / 12.0), 5.0);
        double x = sin(t) * (exp(cos(t)) - c - d);
        double y = cos(t) * (exp(cos(t)) - c - d);
        pt.t = t;
        pt.x = x;
        pt.y = y;
        points.push_back(pt);
        t += h;
    }
}

static void FindMinMax(
    double& xMin, double& xMax,
    double& yMin, double& yMax)
{
    // uses global 2D double points structure

    xMin = yMin = DBL_MAX;
    xMax = yMax = DBL_MIN;

    for (size_t i = 0; i < points.size(); i++)
    {
        Point3d pt = points[i];
        double x = pt.x;
        double y = pt.y;

        if (x < xMin)
            xMin = x;
        if (x > xMax)
            xMax = x;
        if (y < yMin)
            yMin = y;
        if (y > yMax)
            yMax = y;
    }
}

static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
    // Draw the text with formatting options
    DrawTextA(hdc, text, -1, &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:
        CreateButterflyGraphPoints();
        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);
            double h = 0, pi = 0, plm = 0, theta = 0;
            double xMax = 0, xMin = 0, yMax = 0, yMin = 0;
            FindMinMax(xMin, xMax, yMin, yMax);
            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 = { };
            int i = 0;
            float x = (float)xMin;
            float y = (float)yMax;
            px = x;
            py = y;
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
            char buffer[128] = { };

            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);

                sprintf_s(buffer, "%5.4lf", x);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    buffer,
                    (int)strlen(buffer),
                    &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, buffer, 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)
                {
                    sprintf_s(buffer, "%+5.3lf", y);
                    SIZE size = { };
                    GetTextExtentPoint32A(
                        hdc,
                        buffer,
                        (int)strlen(buffer),
                        &size);
                    RECT textRect = { };
                    textRect.left = (long)(sx0 - size.cx - size.cx / 5.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, buffer, textRect);
                }

                y += deltaY;
                i++;
            }

            HGDIOBJ bPenNew = NULL;
            HGDIOBJ hPenOld = NULL;

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

            px = (float)points[0].x;
            py = (float)points[0].y;
            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 < points.size(); j++)
            {
                px = (float)points[j].x;
                py = (float)points[j].y;
                sx = xSlope * px + xInter;
                sy = ySlope * py + yInter;
                LineTo(hdc, (int)sx, (int)sy);
            }

            SelectObject(hdc, hPenOld);
            DeleteObject(bPenNew);
            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 © Sunday, April 19, 2026, by James Pate Williams, Jr., Scattering from a Spherically Symmetric Potential

Blog Entry Sunday April 19 2026Download

Blog Entry © April 1, 2026, by James Pate Williams, Jr. More Romberg Integration Results

Blog Entry Wednesday April 1 2026Download
/*
  Author:  Pate Willimas c 1995

  The following program is a translation of the FORTRAN77
  subroutine found in "Elementary Numerical Analysis 
  An Algorithmic Approach Third Edition" (c) 1980 by S. D.
  Conte and Carl de Boor pages 343-344. The program
  uses Romberg extrapolation to find the integral of a
  function. Also see: https://dlmf.nist.gov/3.5#E10
*/

#include <math.h>
#include <stdio.h>
#include <vector>
#include "Bessel.h"

typedef double real;

long fe;
real t[10][10];

real f(real x)
{
	fe++;
	return(expl(-x * x));
}

real g(real x)
{
	fe++;
	std::vector<real> j(4);
	Bessel::bessj(x, 0, j);
	return exp(-x) * j[0];
}

real romberg(
	real a, real b,
	int start, int row,
	real (*fx)(real))
{
	int i, /*j, */k, m;
	real h, ratio, sum;

	m = start;
	h = (b - a) / m;
	sum = 0.5 * (fx(a) + fx(b));
	if (m > 1)
		for (i = 1; i < m; i++)
			sum += fx(a + i * h);
	t[1][1] = h * sum;
	//printf_s("Romberg T-Table\n%9.7lf\n", t[1][1]);
	if (row < 2) return(t[1][1]);
	for (k = 2; k <= row; k++)
	{
		h = 0.5 * h;
		m *= 2;
		sum = 0.0;
		for (i = 1; i <= m; i += 2)
			sum += fx(a + h * i);
		t[k][1] = 0.5 * t[k - 1][1] + sum * h;
		for (i = 1; i < k; i++)
		{
			t[k - 1][i] = t[k][i] - t[k - 1][i];
			t[k][i + 1] = t[k][i] - t[k - 1][i] /
				(powl(4.0, (real)i) - 1.0);
			//for (j = 1; j < k; j++)
				//printf_s("%9.7lf ", t[k][j]);
		}
		//printf_s("\n");
	}
	if (row < 3) return(t[2][2]);
	//printf_s("Table of ratios\n");
	for (k = 1; k <= row - 2; k++) {
		for (i = 1; i <= k; i++)
		{
			if (t[k + 1][i] == 0.0)
				ratio = 0.0;
			else
				ratio = t[k][i] / t[k + 1][i];
			t[k][i] = ratio;
			//printf_s("%5.2lf ", ratio);
		}
		//printf_s("\n");
	}
	return(t[row][row - 1]);
}

real CompositeTrapezoidalRule(
	real a, real b, int n,
	real(*fx)(real)) {
	real endPts = 0.5 * (fx(a) + fx(b));
	real sum = 0, xk = 0.0;
	real h = (b - a) / n;
	for (int k = 1; k <= n - 1; k++)
	{
		xk = a + k * h;
		sum += fx(xk);
	}
	return h * (endPts + sum);
}

double CompositeSimpsonsRule(
	double a, double b, int n,
	real(*fx)(real))
{
	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;
}

void Integrate(int trial, real(*fx)(real))
{
	real a = 0.0, b = 0.0;
	int row, start;

	if (trial == 0)
	{
		a = 0.0;
		b = 1.0;
	}
	else
	{
		a = 0.0;
		b = 30.0;
	}

	if (trial == 0)
		printf("Romberg integration of f(x) = exp(- x * x)\n");
	else
		printf("Romberg integration of g(x) = exp(-x) * J0(x)\n");
	printf("number of trapezoidal intervals = ");
	scanf_s("%d", &start);
	printf("number of rows in table (<= 8) = ");
	scanf_s("%d", &row);
	fe = 0;
	printf_s("integral = %14.11lf\tevals = %d\n",
		romberg(a, b, start, row, fx), fe);
	printf("number of trapezoidal intervals = ");
	scanf_s("%d", &row);
	fe = 0;
	printf_s("integral = %14.11lf\tevals = %d\n",
		CompositeTrapezoidalRule(a, b, row, fx), fe);
	fe = 0;
	printf("number of Simpson's   intervals = ");
	scanf_s("%d", &row);
	printf_s("integral = %14.11lf\tevals = %d\n",
		CompositeSimpsonsRule(a, b, row, fx), fe);
	if (trial == 0)
		printf_s("integral = %14.11lf\n", 0.74682413279);
	else
	{
		real pi = 4.0 * atan(1.0);
		real si = sin(pi * 45.0 / 180.0);
		printf_s("integral = %14.11lf\n", si);
	}
}

int main(void)
{
	Integrate(0, f);
	Integrate(1, g);
	return(0);
}

Blog Entry © Sunday, March 22, 2026, by James Pate Williams, Jr. Mueller’s Method

Blog Entry Sunday March 22 2026Download
/*
* 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 "Elementary 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 
*/

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

//#define DEBUG

static _Lcomplex HornersMethod(_Lcomplex coeff[], _Lcomplex z, int degree) {
	int i = 0;
	_Lcomplex c = coeff[degree];

	for (i = degree; i >= 1; i--) {
		_Lcomplex product = _LCmulcc(c, z);
		c._Val[0] = product._Val[0] + coeff[i - 1]._Val[0];
		c._Val[1] = product._Val[1] + coeff[i - 1]._Val[1];
	}
#ifdef DEBUG
	_Lcomplex sum = { 0 };
	for (i = 0; i <= degree; i++) {
		_Lcomplex term =
			_LCmulcc(cpowl(z, _LCbuild(i, 0.0)), coeff[i]);
		sum._Val[0] += term._Val[0];
		sum._Val[1] += term._Val[1];
	}

	long double delta = fabsl(cabsl(c) - cabsl(sum));

	if (delta > 1.0e-12)
		exit(-5);
#endif
	return c;
}

static void comdiv(
	long double xr, long double xi,
	long double yr, long double yi,
	long double* zr, long double* zi)
{
	long double h, d;

	if (fabs(yi) < fabs(yr)) {
		if (yi == 0.0) {
			*zr = xr / yr;
			*zi = xi / yr;
		}
		else {
			h = yi / yr;
			d = h * yi + yr;
			*zr = (xr + h * xi) / d;
			*zi = (xi - h * xr) / d;
		}
	}
	else {
		h = yr / yi;
		d = h * yr + yi;
		*zr = (xr * h + xi) / d;
		*zi = (xi * h - xr) / d;
	}
}

#ifdef DEBUG
static _Lcomplex MyComplexDivide(_Lcomplex numer, _Lcomplex denom) {
	long double norm2 =
		denom._Val[0] * denom._Val[0] +
		denom._Val[1] * denom._Val[1];
	_Lcomplex result = { 0 };

	result._Val[0] = (
		numer._Val[0] * denom._Val[0] +
		numer._Val[1] * denom._Val[1]) / norm2;
	result._Val[1] = (
		numer._Val[1] * denom._Val[0] -
		numer._Val[0] * denom._Val[1]) / norm2;
	return result;
}
#endif

static _Lcomplex ComplexDivide(_Lcomplex numer, _Lcomplex denom) {
	_Lcomplex result = { 0 };

	comdiv(
		numer._Val[0], numer._Val[1],
		denom._Val[0], denom._Val[1],
		&result._Val[0], &result._Val[1]);
#ifdef DEBUG
	_Lcomplex myResult = MyComplexDivide(numer, denom);
	long double delta = fabsl(cabsl(result) - cabsl(myResult));

	if (delta > 1.0e-12)
		exit(-6);
#endif
	return result;
}

static int Deflate(
	_Lcomplex coeff[], _Lcomplex zero,
	_Lcomplex* fzero, _Lcomplex* fzrdfl,
	_Lcomplex zeros[], int i, int* count,
	int degree) {
	_Lcomplex denom = { 0 };

	(*count)++;

	*fzero = HornersMethod(coeff, zero, degree);
	*fzrdfl = *fzero;

	if (i < 1)
		return 0;

	for (int j = 1; j <= i; j++) {
		denom._Val[0] = zero._Val[0] - zeros[j - 1]._Val[0];
		denom._Val[1] = zero._Val[1] - zeros[j - 1]._Val[1];

		if (cabsl(denom) == 0.0) {
			zeros[i] = _LCmulcr(zero, 1.001);
			return 1;
		}

		else
			*fzrdfl = ComplexDivide(*fzrdfl, denom);
	}

	return 0;
}

static void Mueller(
	_Lcomplex coeff[], _Lcomplex zeros[],
	double epsilon1, double epsilon2,
	int degree, int fnreal, int maxIts, int n, int nPrev) {
	double eps1 = max(epsilon1, 1.0e-12);
	double eps2 = max(epsilon2, 1.0e-20);
	int count = 0, i = 0;
	_Lcomplex c = { 0 };
	_Lcomplex denom = { 0 };
	_Lcomplex divdf1 = { 0 };
	_Lcomplex divdf2 = { 0 };
	_Lcomplex divdf1p = { 0 };
	_Lcomplex fzero = { 0 };
	_Lcomplex fzr = { 0 };
	_Lcomplex fzdfl = { 0 };
	_Lcomplex fzrdfl = { 0 };
	_Lcomplex fzrprv = { 0 };
	_Lcomplex four = _LCbuild(4.0, 0.0);
	_Lcomplex h = { 0 };
	_Lcomplex hprev = { 0 };
	_Lcomplex sqr = { 0 };
	_Lcomplex zero = { 0 };
	_Lcomplex p5 = _LCbuild(0.5, 0.0);
	_Lcomplex zeropp5 = { 0 };
	_Lcomplex zeromp5 = { 0 };
	_Lcomplex diff = { 0 };
	_Lcomplex tadd = { 0 };
	_Lcomplex tmul = { 0 };
	_Lcomplex umul = { 0 };
	_Lcomplex vmul = { 0 };

	for (i = nPrev; i < n; i++) {
		count = 0;

	Label1:

		zero = zeros[i];
		h = p5;

		zeropp5._Val[0] = zero._Val[0] + p5._Val[0];
		zeropp5._Val[1] = zero._Val[1] + p5._Val[1];

		if (Deflate(
			coeff, zeropp5, &fzr, &divdf1p,
			zeros, i, &count, degree))
			goto Label1;

		zeromp5._Val[0] = zero._Val[0] - p5._Val[0];
		zeromp5._Val[1] = zero._Val[1] - p5._Val[1];

		if (Deflate(
			coeff, zeromp5, &fzr, &fzrprv,
			zeros, i, &count, degree))
			goto Label1;

		hprev._Val[0] = -1.0;
		hprev._Val[1] = 0.0;
		diff._Val[0] = fzrprv._Val[0] - divdf1p._Val[0];
		diff._Val[1] = fzrprv._Val[1] - divdf1p._Val[1];
		if (cabsl(hprev) == 0)
			exit(-2);
		divdf1p = ComplexDivide(diff, hprev);

		if (Deflate(
			coeff, zero, &fzr, &fzrdfl,
			zeros, i, &count, degree))
			goto Label1;

	Label2:

		diff._Val[0] = fzrdfl._Val[0] - fzrprv._Val[0];
		diff._Val[1] = fzrdfl._Val[1] - fzrprv._Val[1];
		if (cabsl(h) == 0)
			exit(-3);
		divdf1 = ComplexDivide(diff, h);
		diff._Val[0] = divdf1._Val[0] - divdf1p._Val[0];
		diff._Val[1] = divdf1._Val[1] - divdf1p._Val[1];
		tadd._Val[0] = h._Val[0] + hprev._Val[0];
		tadd._Val[1] = h._Val[1] + hprev._Val[1];
		if (cabsl(tadd) == 0)
			exit(-3);
		divdf2 = ComplexDivide(diff, tadd);
		hprev = h;
		divdf1p = divdf1;
		tmul = _LCmulcc(h, divdf2);
		c._Val[0] = divdf1._Val[0] + tmul._Val[0];
		c._Val[1] = divdf1._Val[1] + tmul._Val[1];
		tmul = _LCmulcc(c, c);
		umul = _LCmulcc(four, fzrdfl);
		vmul = _LCmulcc(umul, divdf2);
		sqr._Val[0] = tmul._Val[0] - vmul._Val[0];
		sqr._Val[1] = tmul._Val[1] - vmul._Val[1];

		if (fnreal && sqr._Val[0] < 0.0)
		{
			sqr._Val[0] = 0.0;
			sqr._Val[1] = 0.0;
		}

		sqr = csqrtl(sqr);

		if ((c._Val[0] * sqr._Val[0] + c._Val[1] * sqr._Val[1]) < 0.0) {
			denom._Val[0] = c._Val[0] - sqr._Val[0];
			denom._Val[1] = c._Val[1] - sqr._Val[1];
		}
		else {
			denom._Val[0] = c._Val[0] + sqr._Val[0];
			denom._Val[1] = c._Val[1] + sqr._Val[1];
		}
		if (cabsl(denom) <= 0.0)
		{
			denom._Val[0] = 1.0;
			denom._Val[1] = 0.0;
		}
		if (cabsl(denom) == 0)
			exit(-4);
		tmul = _LCmulcr(fzrdfl, -2.0);
		h = ComplexDivide(tmul, denom);
		fzrprv = fzrdfl;
		zero._Val[0] = zero._Val[0] + h._Val[0];
		zero._Val[1] = zero._Val[1] + h._Val[1];

		if (count > maxIts)
			goto Label4;

	Label3:

		if (Deflate(
			coeff, zero, &fzr, &fzrdfl,
			zeros, i, &count, degree))
			goto Label1;

		if (cabsl(h) < eps1 * cabsl(zero))
			goto Label4;

		if (max(cabsl(fzr), cabsl(fzdfl)) < eps2)
			goto Label4;

		if (cabsl(fzrdfl) >= 10.0 * cabsl(fzrprv)) {
			h = _LCmulcr(h, 0.5);
			zero._Val[0] = zero._Val[0] - h._Val[0];
			zero._Val[1] = zero._Val[1] - h._Val[1];
			goto Label3;
		}

		else
			goto Label2;

	Label4:

		zeros[i] = zero;
	}
}

int main(void)
{
	double epsilon1 = 1.0e-15;
	double epsilon2 = 1.0e-15;
	int degree = 0, fnreal = 0, i = 0, maxIts = 1000;
	int n = 0, nPrev = 0;
		
	while (1) {
		_Lcomplex* coeff = NULL;
		_Lcomplex* zeros = NULL;

		printf_s("Degree (0 to quit) = ");
		scanf_s("%d", &degree);

		if (degree == 0)
			break;

		n = degree;
		coeff = calloc(degree + 1, sizeof(_Lcomplex));

		if (coeff == NULL)
			exit(-1);

		zeros = calloc(n, sizeof(_Lcomplex));

		if (zeros == NULL)
			exit(-1);

		for (i = degree; i >= 0; i--) {
			printf_s("coefficient[%d].real = ", i);
			scanf_s("%Lf", &coeff[i]._Val[0]);
			printf_s("coefficient[%d].imag = ", i);
			scanf_s("%Lf", &coeff[i]._Val[1]);
		}

		Mueller(
			coeff, zeros, epsilon1,
			epsilon2, degree, fnreal,
			maxIts, n, nPrev);

		printf_s("\n");

		for (i = 0; i < degree; i++) {
			printf_s("zero[%d].real = %17.10e\t", i, zeros[i]._Val[0]);
			printf_s("zero[%d].imag = %17.10e\n", i, zeros[i]._Val[1]);
		}

		printf_s("\n");

		for (i = 0; i < degree; i++) {
			_Lcomplex func = HornersMethod(coeff, zeros[i], degree);

			printf_s("func[%d].real = %17.10e\t", i, func._Val[0]);
			printf_s("func[%d].imag = %17.10e\n", i, func._Val[1]);
		}

		printf_s("\n");

		free(coeff);
		free(zeros);
	}

	return 0;
}