Tag: C/C++
The programming language I utilized.
Blog Entry © Saturday – Tuesday, May 30 – June 2, 2026, by James Pate Williams, Jr., theMicrosoft Bing Copilot, the M365 Copilot Partial Reproduction of Figures 8 and 9 fromChapter 2 of Quantum Mechanics Third Edition by Leonard I. Schiff
class Figure
{
public:
static void ComputeFigure(
bool eight,
double& xi2, double& eta2,
std::vector<double>& radius,
std::vector<double>& xi,
std::vector<double>& eta,
std::vector<double>& vix,
std::vector<double>& viy,
std::vector<double>& energy1,
double (*f)(double),
double (*g)(double));
};
#include "Figure.h"
#include <algorithm>
#include <cmath>
#include <iterator>
#include <set>
#include <vector>
// Function to compute intersection of two sorted containers of doubles with tolerance
template <typename InputIt1, typename InputIt2, typename OutputIt>
void set_intersection_with_tolerance(InputIt1 first1, InputIt1 last1,
InputIt2 first2, InputIt2 last2,
OutputIt d_first, double tolerance) {
while (first1 != last1 && first2 != last2) {
double a = *first1;
double b = *first2;
if (std::fabs(a - b) <= tolerance) {
// Values are considered equal within tolerance
*d_first++ = /*a; // or */ (a + b) / 2.0; //if you want averaged value
++first1;
++first2;
}
else if (a < b - tolerance) {
++first1;
}
else {
++first2;
}
}
}
static double f(double xi)
{
return xi * tan(xi);
}
static double g(double xi)
{
return -xi * 1.0 / tan(xi);
}
static void Intersection(
bool eight, std::vector<double>& intersection)
{
double del = 0.0, radius2 = 0.0;
double tolerance = 0.0;
double (*h)(double);
int ilimit = 0;
std::set<double> eta, rad;
if (eight)
{
ilimit = 10000;
radius2 = 1.0;
tolerance = 1.0e-1;
h = f;
}
else
{
ilimit = 10000;
radius2 = 4.0;
tolerance = 1.0e-1;
h = g;
}
del = radius2 / ilimit;
intersection.clear();
for (size_t i = 0; i < ilimit; i++)
{
double xi0 = i * del;
if (eight && xi0 >= 0.65)
{
double et = h(xi0);
double r2 = xi0 * xi0 + et * et;
if (fabs(r2 - radius2) < tolerance)
{
eta.insert(et);
rad.insert(r2);
}
}
else if (!eight && xi0 >= 1.8125)
{
double et = h(xi0);
double r2 = xi0 * xi0 + et * et;
if (fabs(r2 - radius2) < tolerance)
{
eta.insert(et);
rad.insert(r2);
}
}
}
if (eight)
{
size_t count = 0, index = 0;
for (double val : eta)
{
if (index == eta.size() - 1)
{
intersection.push_back(val);
break;
}
count++;
index++;
}
}
else
{
int count = 24, index = 0;
for (double val : eta)
{
if (index == eta.size() - count)
{
intersection.push_back(val);
break;
}
index++;
}
}
/*set_intersection_with_tolerance(
eta.begin(), eta.end(),
rad.begin(), rad.end(),
std::back_inserter(intersection),
tolerance);*/
}
void Figure::ComputeFigure(
bool eight,
double& xi2, double& eta2,
std::vector<double>& radius,
std::vector<double>& xi,
std::vector<double>& eta,
std::vector<double>& vix,
std::vector<double>& viy,
std::vector<double>& intersection,
double (*f)(double),
double (*g)(double))
{
double xi0 = 0.0;
double xi1 = 3.5;
double eta0 = 0.0, eta1 = 0.0;
double (*h)(double);
xi.clear();
eta.clear();
if (eight)
h = f;
else
{
xi0 = 1.5;
h = g;
}
eta0 = h(xi0);
eta1 = h(xi1);
double deltaXi = (xi1 - xi0) / 10000.0;
int count = 0;
vix.clear();
viy.clear();
for (int j = 0; j <= 10000; j++)
{
double x = j * deltaXi;
double hx = h(x);
double vx = x * x + hx * hx;
if (eight && x >= xi0 && x <= xi1)
{
if (count == 0 && x >= xi0)
{
xi.push_back(x);
eta.push_back(hx);
count = 1;
}
else
{
xi.push_back(x);
eta.push_back(hx);
}
for (int k = 0; k < 2; k++)
{
double r = radius[k];
vix.push_back(xi[xi.size() - 1]);
viy.push_back(r * r);
}
}
else if (!eight && x >= xi0 && x <= xi1)
{
if (count == 0 && x >= xi0)
{
xi.push_back(x);
eta.push_back(hx);
count = 1;
}
else if (count == 1)
{
xi.push_back(x);
eta.push_back(hx);
}
if (xi.size() > 0)
{
for (int k = 0; k < 1; k++)
{
double r = radius[k];
vix.push_back(xi[xi.size() - 1]);
viy.push_back(r * r);
}
}
}
if (eight && xi[j] >= 0.0 && eta[j] >= 3.5)
break;
else if (!eight &&
xi.size() > 0 &&
xi[xi.size() - 1] > 1.5 &&
eta.size() > 0 &&
eta[eta.size() - 1] >= 3.5)
break;
}
Intersection(eight, intersection);
}
// SchiffChapter2Fig8and9.cpp : Defines the entry point for the application.
//
#include "framework.h"
#include "SchiffChapter2Fig8and9.h"
#include <time.h>
#include <algorithm>
#include <vector>
#include "Figure.h"
#define MAX_LOADSTRING 100
typedef struct tagPoint2d
{
double x, y;
} Point2d, * PPoint2d;
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
WCHAR line[8192], text[8192]; // wide character buffers
WCHAR title[65536]; // window title
bool eight; // true = plot figure 8
double xi2, eta2; // energy cordinates
std::vector<Point2d> points; // plotting 2d points
std::vector<double> radius; // radius vector { 1.0, 2.0 }
std::vector<double> xi; // Greek x-coordinate
std::vector<double> eta; // Greek y-coordinate
std::vector<double> V0; // potential of the well
std::vector<double> vix; // potential x-coordinate
std::vector<double> viy; // potential y-coordinate
std::vector<double> energy1; // energy eigenvalues
// 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 DrawEtaDialog(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_SCHIFFCHAPTER2FIG8AND9, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_SCHIFFCHAPTER2FIG8AND9));
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_SCHIFFCHAPTER2FIG8AND9));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_SCHIFFCHAPTER2FIG8AND9);
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 fx(double x)
{
return x * tan(x);
}
static double gx(double x)
{
return -x * 1.0 / tan(x);
}
//
// 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_FIGURE8:
{
eight = true;
text[0] = L'\0';
clock_t clock0 = clock();
radius.push_back(1.0);
radius.push_back(2.0);
Figure::ComputeFigure(
eight, xi2, eta2, radius, xi, eta, vix, viy, energy1, fx, gx);
clock_t clock1 = clock() - clock0;
double runtime = (double)clock1 / CLOCKS_PER_SEC;
swprintf_s(line, 8192, L"Runtime in Seconds = %lf\r\n", runtime);
wcscat_s(text, 8192, line);
points.clear();
for (size_t j = 0; j < xi.size(); j++)
{
Point2d pt = { 0 };
pt.x = xi[j];
pt.y = eta[j];
points.push_back(pt);
}
for (size_t i = 0; i < energy1.size(); i++)
{
swprintf_s(line, 8192,
L"E[%zu] = %lf\r\n", i, energy1[i]);
wcscat_s(text, 8192, line);
}
MessageBox(hWnd, text, L"Energy Information",
MB_OK | MB_ICONINFORMATION);
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_ETA_DIALOG), hWnd, DrawEtaDialog);
break;
}
case IDM_FIGURE9:
{
eight = false;
text[0] = L'\0';
clock_t clock0 = clock();
radius.push_back(2.0);
Figure::ComputeFigure(
eight, xi2, eta2, radius, xi, eta, vix, viy, energy1, fx, gx);
clock_t clock1 = clock() - clock0;
double runtime = (double)clock1 / CLOCKS_PER_SEC;
swprintf_s(line, 8192, L"Runtime in Seconds = %lf\r\n", runtime);
wcscat_s(text, 8192, line);
points.clear();
for (size_t j = 0; j < xi.size(); j++)
{
Point2d pt = { 0 };
pt.x = xi[j];
pt.y = eta[j];
points.push_back(pt);
}
for (size_t i = 0; i < energy1.size(); i++)
{
swprintf_s(line, 8192,
L"E[%zu] = %lf\r\n", i, energy1[i]);
wcscat_s(text, 8192, line);
}
MessageBox(hWnd, text, L"Energy Information",
MB_OK | MB_ICONINFORMATION);
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_ETA_DIALOG), hWnd, DrawEtaDialog);
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;
}
static void FindMinMax(
std::vector<Point2d>& points,
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++)
{
Point2d 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);
}
static void DrawQuarterCircleArc(
HDC hdc,
float xSlope, float ySlope,
float xInter, float yInter,
float radius, bool topToRight)
{
auto mapX = [&](float x)
{
return (int)lroundf(xSlope * x + xInter);
};
auto mapY = [&](float y)
{
return (int)lroundf(ySlope * y + yInter);
};
int x1 = mapX(-radius);
int y1 = mapY(+radius);
int x2 = mapX(+radius);
int y2 = mapY(-radius);
int left = min(x1, x2);
int right = max(x1, x2);
int top = min(y1, y2);
int bottom = max(y1, y2);
int xTop = mapX(0.0f);
int yTop = mapY(radius);
int xRight = mapX(radius);
int yRight = mapY(0.0f);
if (topToRight)
{
Arc(hdc, left, top, right, bottom,
xTop, yTop, xRight, yRight);
}
else
{
Arc(hdc, left, top, right, bottom,
xRight, yRight, xTop, yTop);
}
}
INT_PTR CALLBACK DrawEtaDialog(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
SetWindowText(hDlg, title);
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 h = 0, pi = 0, plm = 0, theta = 0;
double xMax = 0, xMin = 0, yMax = 0, yMin = 0;
FindMinMax(points, xMin, xMax, yMin, yMax);
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, xInter, ySlope, yInter;
xSlope = (sx1 - sx0) / xSpan;
xInter = (float)(sx0 - xSlope * xMin);
ySlope = (sy0 - sy1) / ySpan;
yInter = (float)(sy0 - ySlope * yMax);
float px = 0, py = 0, sx = 0, sy = 0;
float vx = 0, vy = 0;
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hDlg, &ps);
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);
HRGN hrgn = CreateRectRgn((int)sx0, (int)sy0, (int)sx1, (int)sy1);
// Select the clipping region into the DC
if (SelectClipRgn(hdc, hrgn) == ERROR) {
MessageBox(hDlg, L"Failed to select clip region",
L"Error", MB_ICONERROR);
return (INT_PTR)FALSE;
}
SelectClipRgn(hdc, hrgn);
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);
}
float radius = 0.0f;
if (eight)
radius = 1.0f;
else
radius = 2.0f;
DrawQuarterCircleArc(
hdc, xSlope, ySlope,
xInter, yInter, radius, false);
if (eight)
{
radius = 2.0f;
DrawQuarterCircleArc(
hdc, xSlope, ySlope,
xInter, yInter, radius, false);
}
SelectObject(hdc, hPenOld);
DeleteObject(bPenNew);
return (INT_PTR)FALSE;
}
return (INT_PTR)FALSE;
}
Blog Entry © Wednesday, May 27, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Grade School Arithmetic
#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 (c) Tuesday, May 26, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Hydrogen-Like Atomic Radial Wave Functions
Included a downloadable PDF and Microsoft Excel Workbook.
Blog Entry (c) Monday, May 25, 2026, by James Pate Williams, Jr. Quantum Mechanical Linear Harmonic Oscillator
Blog Entry © Saturday, May 16, 2026, by James Pate Williams, Jr. Some More Linear Algebra Examples
Blog Entry © Thursday, May 14, 2026, by James Pate Williams, Jr. More Numerical Integration Results
// 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;
}
