Category: C Applications
Blog Entry © Tuesday, June 16, 2026, by James Pate Williams, Jr., More 64-Bit Pseudoprimes Found Using John Pollard’s Factoring with Cubic Integers
Blog Entry © Monday, June 15, 2026, by James Pate Williams, Jr. Filtering a Noisy Signal
#pragma once
#include <complex>
#include <vector>
class Transform
{
public:
static void VandermondeDFT(
int n,
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& y);
static void InverseVandermondeDFT(
int n,
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& y);
static std::vector<std::complex<double>> DFT(
std::vector<double>& x, std::vector<double>& f);
static std::vector<double> InverseDFT(
std::vector<double>& f,
std::vector<std::complex<double>>& X);
/*
* Reference: "Elementary Numerical Analysis:
* An Algorithmic Approach Third Edition" (c)
* 1980 by S. D. Conte and Carl de Boor
* Section 6.5 pages 268 - 277 and Section 6.6
* pages 277 - 283
* Input to FFT
* Z1, Z2 complex n-vectors
* n the length of the vectors
* inzee
* = 1 transform in Z1
* = 2 transform in Z2
* Constructs the discrete Fourier transform in the Cooley-
* Tukey way, but with a twist.
*/
static void FFT(
std::vector<std::complex<double>>& Z1,
int& after, int& now, int& before, int& inzee,
std::vector<std::complex<double>>& Z2);
/*
* This computes an in - place complex - to - complex FFT
* x and y are the real and imaginary arrays of 2^m points.
* dir = 1 gives forward transform
* dir = -1 gives reverse transform
* see http://astronomy.swin.edu.au/~pbourke/analysis/dft/
* Website no longer exists
*/
static void FFT(short dir, int m,
std::vector<double>& x, std::vector<double>& y);
/*
* Reference: "Introduction to Algorithms" by
* Thomas H. Cormen, Charles E. Leiserson, and
* Ronald L. Rivest, pages 794 - 795
*/
static void IterativeFFT(
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& A);
/*
* Reference: "Introduction to Algorithms" by
* Thomas H. Cormen, Charles E. Leiserson, and
* Ronald L. Rivest, page 788
*/
static std::vector<std::complex<double>> RecursiveFFT(
std::vector<std::complex<double>>& a);
};
#include "Transform.h"
void Transform::VandermondeDFT(
int n,
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& y)
{
double pi = 4.0 * atan(1.0);
std::complex<double> z(0.0, 2.0 * pi / n);
std::complex<double> omegaN = exp(z);
std::vector<std::vector<std::complex<double>>> V(n);
for (int k = 0; k < n; k++)
{
V[k].resize(n);
for (int j = 0; j < n; j++)
{
V[k][j] = std::pow(omegaN, k * j);
}
}
for (int k = 0; k < n; k++)
{
std::complex<double> sum = 0.0;
for (int j = 0; j < n; j++)
{
sum += V[k][j] * a[j];
}
y[k] = sum;
}
}
void Transform::InverseVandermondeDFT(
int n,
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& y)
{
double pi = 4.0 * atan(1.0);
std::complex<double> nc = { static_cast<double>(n), 0.0 };
std::complex<double> z(0.0, 2.0 * pi / n);
std::complex<double> omegaN = exp(z);
std::vector<std::vector<std::complex<double>>> invV(n);
for (int k = 0; k < n; k++)
{
invV[k].resize(n);
for (int j = 0; j < n; j++)
{
invV[k][j] = std::pow(omegaN, -k * j);
}
}
for (int k = 0; k < n; k++)
{
std::complex<double> sum = 0.0;
for (int j = 0; j < n; j++)
{
sum += invV[k][j] * y[j];
}
a[k] = sum / nc;
}
}
std::vector<std::complex<double>> Transform::DFT(
std::vector<double>& x, std::vector<double>& f)
{
int length = static_cast<int>(x.size());
double pi = 4.0 * atan(1.0);
double pi2oN = 2.0 * pi / length;
int k, n;
std::vector<double> X(length);
std::vector<double> Y(length);
std::vector<std::complex<double>> Z(length);
f.resize(length);
for (k = 0; k < length; k++)
{
X[k] = Y[k] = 0;
for (n = 0; n < length; n++)
{
X[k] += x[n] * cos(pi2oN * k * n);
Y[k] -= x[n] * sin(pi2oN * k * n);
}
f[k] = pi2oN * k;
X[k] /= length;
Y[k] /= length;
Z[k] = { X[k], Y[k] };
}
return Z;
}
std::vector<double> Transform::InverseDFT(
std::vector<double>& f,
std::vector<std::complex<double>>& X)
{
double imag = 0.0;
int length = static_cast<int>(X.size());
std::vector<double> x(length);
for (int n = 0; n < length; n++)
{
imag = x[n] = 0.0;
for (int k = 0; k < length; k++)
{
x[n] += X[k]._Val[0] * cos(f[k] * n)
- X[k]._Val[1] * sin(f[k] * n);
imag += X[k]._Val[0] * sin(f[k] * n)
+ X[k]._Val[1] * cos(f[k] * n);
}
}
return x;
}
static void FFTStep(
std::vector<std::complex<double>>& Zinp,
int after, int now, int before,
std::vector<std::complex<double>>& Zout)
{
double angle = 0.0, ratio = 0.0;
double twoPi = 2.0 * 4.0 * atan(1.0);
int ia = 0, ib = 0, inp = 0, j = 0;
std::complex<double> arg = 1.0, omega = 0, value = 0;
angle = twoPi / ((now + 1) * (after + 1));
omega = std::complex<double>(cos(angle), -sin(angle));
int address = 1;
for (int i = 1; i <= now; i++)
{
for (int j = 1; j <= after; j++)
{
for (int k = 1; k <= before; k++)
{
address = i * j * k;
if (address < Zout.size())
Zout[address] = { 0.0, 0.0 };
}
}
}
address = 1;
for (int j = 1; j <= now; j++)
{
for (ia = 1; ia <= after; ia++)
{
for (ib = 1; ib <= before; ib++)
{
int address = j * ia * ib;
if (address < Zinp.size())
value = Zinp[address];
for (inp = now - 1; inp >= 1; inp--)
{
address = ia * ib * inp;
if (address < Zinp.size())
value = value * arg + Zinp[address];
}
address = ia * j * ib;
if (address < Zout.size())
Zout[address] = value;
}
arg *= omega;
}
}
}
void Transform::FFT(
std::vector<std::complex<double>>& Z1,
int& after, int& now, int& before, int& inzee,
std::vector<std::complex<double>>& Z2)
{
std::vector<int> prime =
{ 0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
int next = 1, nextmx = 25;
after = 1;
before = (int)Z1.size();
now = 1;
Label10:
if (before / prime[next] * prime[next] < before)
{
next++;
if (next <= nextmx)
goto Label10;
else
{
now = before;
before = 1;
}
}
else
{
now = prime[next];
before /= prime[next];
}
if (inzee == 1)
FFTStep(Z1, after, now, before, Z2);
else
FFTStep(Z2, after, now, before, Z1);
inzee = 3 - inzee;
if (before == 1)
return;
after *= now;
goto Label10;
}
void Transform::FFT(short dir, int m,
std::vector<double>& x, std::vector<double>& y)
{
int n, i, i1, j, k, i2, l, l1, l2;
double c1, c2, tx, ty, t1, t2, u1, u2, z;
// Calculate the number of points
n = 1;
for (i = 0; i < m; i++)
n *= 2;
// Do the bit reversal
i2 = n >> 1;
j = 0;
for (i = 0; i < n - 1; i++)
{
if (i < j)
{
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j)
{
j -= k;
k >>= 1;
}
j += k;
}
// Compute the FFT
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l = 0; l < m; l++)
{
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j = 0; j < l1; j++)
{
for (i = j; i < n; i += l2)
{
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}
// Scaling for forward transform
if (dir == 1)
{
for (i = 0; i < n; i++)
{
x[i] /= n;
y[i] /= n;
}
}
}
static void FFTBase(
std::vector<std::complex<double>> a,
std::vector<std::complex<double>> A)
{
double pi = 4.0 * atan(1.0);
int n = static_cast<int>(a.size());
for (int s = 1; s <= log2(n); s++)
{
int m = static_cast<int>(pow(2, s));
std::complex<double> z(0.0, 2.0 * pi / m);
std::complex<double> omegaM = exp(z);
for (int k = 0; k <= n - 1; k += m)
{
std::complex<double> omega = { 1.0, 0.0 };
for (int j = 0; j <= m / 2 - 1; j++)
{
std::complex<double> t = omega * A[k + j + m / 2];
std::complex<double> u = A[k + j];
std::complex<double> jc = { static_cast<double>(j), 0.0 };
A[k + j] = u + jc;
A[k + j + m / 2] = u - t;
omega *= omegaM;
}
}
}
}
static int Reverse(int k)
{
int digits[32] = { 0 }, i = 0;
while (k > 0)
{
int digit = k & 1;
k >>= 1;
digits[i++] = digit;
}
int result = digits[0];
for (int j = 1; j < i; j++)
result = result * 2 + digits[j];
return result;
}
static void BitReverseCopy(
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& A)
{
int n = static_cast<int>(a.size());
for (int k = 0; k <= n - 1; k++)
A[Reverse(k)] = a[k];
}
void Transform::IterativeFFT(
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& A)
{
BitReverseCopy(a, A);
double pi = 4.0 * atan(1.0);
int n = static_cast<int>(a.size());
for (int s = 1; s <= static_cast<int>(log2(n)); s++)
{
int m = static_cast<int>(pow(2.0, s));
std::complex<double> z(0.0, 2.0 * pi / m);
std::complex<double> omegaM = exp(z);
std::complex<double> omega = { 1.0, 0.0 };
for (int j = 0; j <= m / 2 - 1; j++)
{
for (int k = j; k <= n - 1; k += m)
{
std::complex<double> t = omega * A[k + m / 2];
std::complex<double> u = A[k];
A[k] = u + t;
A[k + m / 2] = u - t;
omega *= omegaM;
}
}
}
}
std::vector<std::complex<double>> Transform::RecursiveFFT(
std::vector<std::complex<double>>& a)
{
int n = static_cast<int>(a.size());
if (n == 1)
return a;
std::vector<std::complex<double>> a0;
std::vector<std::complex<double>> a1;
std::vector<std::complex<double>> y0;
std::vector<std::complex<double>> y1;
for (int i = 0; i <= n - 2; i++)
a0.push_back(a[i]);
for (int i = 1; i <= n - 1; i++)
a1.push_back(a[i]);
y0 = RecursiveFFT(a0);
y1 = RecursiveFFT(a1);
double pi = 4.0 * atan(1.0);
std::complex<double> z(0.0, 2.0 * pi / n);
std::complex<double> omegaN = exp(z);
std::complex<double> omega(1.0, 0.0);
std::vector<std::complex<double>> y(n, 0.0);
for (int k = 0; k <= n / 2 - 1; k++)
{
y[k] = y0[k] + omega * y1[k];
y[(long long)k + n / 2] = y0[k] - omega * y1[k];
omega *= omegaN;
}
return y;
}
// FilteringNoisySignal.cpp : Defines the entry point for the application.
// Copyright (c) Monday, June 15, 2026 by James Pate Williams, Jr.
// Reference: "Numerical Computation 2: Methods, Software and Analysis"
// (c) 1997 by Christoph W. Ueberhuber pages 52-53.
#include "framework.h"
#include "Resource.h"
#include "FilteringNoisySignal.h"
#include "Transform.h"
#include <float.h>
#include <cmath>
#include <vector>
#define MAX_LOADSTRING 100
typedef struct tagPoint2d
{
double t, f;
} 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
char thresholdText[128]; // threshold buffer
char noisePCText[128]; // noise % buffer
double threshold; // noise threshold
double noisePercent; // noise parameter
int Npts = 1024; // number of data points
std::vector<double> originalSignal; // original signal
std::vector<double> perturbedSignalReal; // perturbed signal real part
std::vector<double> perturbedSignalImag; // perturbed signal imag part
std::vector<double> recoveredSignalReal; // recovered signal after filtering
std::vector<double> recoveredSignalImag; // recovered signal after filtering
std::vector<Point2d> points; // plotting data
// 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 InputDialog(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK DrawGraphDialog(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_FILTERINGNOISYSIGNAL, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_FILTERINGNOISYSIGNAL));
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_FILTERINGNOISYSIGNAL));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_FILTERINGNOISYSIGNAL);
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 f(double t)
{
double pi = 4.0 * atan(1.0);
double arg = 2.0 * pi * t;
return 2.0 * sin(arg / 500.0) + cos(arg / 200.0) -
0.5 * sin(arg / 50.0);
}
//
// 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_INPUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_INPUT_DIALOG), hWnd, InputDialog);
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 AddWhiteNoise(
int n, unsigned int seed)
{
double signalMin = DBL_MAX;
double signalMax = DBL_MIN;
std::vector<double> noise(n, 0.0);
perturbedSignalReal.resize(n, 0.0);
srand(seed);
for (int i = 0; i < n; i++)
{
if (originalSignal[i] < signalMin)
signalMin = originalSignal[i];
if (originalSignal[i] > signalMax)
signalMax = originalSignal[i];
}
for (int i = 0; i < n; i++)
noise[i] = (signalMax - signalMin) * rand() /
RAND_MAX + signalMin;
for (int i = 0; i < n; i++)
{
double newSignal = originalSignal[i] +
noisePercent * noise[i];
if (newSignal < signalMin)
newSignal = signalMin;
if (newSignal > signalMax)
newSignal = signalMax;
perturbedSignalReal[i] = newSignal;
}
}
static void Filter(
double threshold, int n, unsigned int seed, HWND hDlg)
{
double pi = 4.0 * atan(1.0);
double step = 2.0 * pi * 200.0 / n, t = 0;
points.clear();
for (int i = 0; i < n; i++)
{
double ft = f(t);
Point2d pt = { t, ft };
points.push_back(pt);
t += step;
}
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG), hDlg,
DrawGraphDialog);
originalSignal.resize(n, 0.0);
for (int i = 0; i < n; i++)
originalSignal[i] = points[i].f;
AddWhiteNoise(n, seed);
points.clear();
t = 0;
for (int i = 0; i < n; i++)
{
double ft = perturbedSignalReal[i];
Point2d pt = { t, ft };
points.push_back(pt);
t += step;
}
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG), hDlg,
DrawGraphDialog);
perturbedSignalImag.resize(n, 0.0);
recoveredSignalReal.resize(n, 0.0);
recoveredSignalImag.resize(n, 0.0);
Transform::FFT(+1, static_cast<int>(log2(n)),
perturbedSignalReal, perturbedSignalImag);
for (int i = 0; i < n; i++)
{
double magnitude2 =
perturbedSignalReal[i] *
perturbedSignalReal[i] +
perturbedSignalImag[i] *
perturbedSignalImag[i];
if (magnitude2 > threshold)
{
recoveredSignalReal[i] = perturbedSignalReal[i];
recoveredSignalImag[i] = perturbedSignalImag[i];
}
}
Transform::FFT(+1, static_cast<int>(log2(n)),
recoveredSignalReal, recoveredSignalImag);
points.clear();
t = 6 * step;
for (int i = 6; i < n - 6; i++)
{
double ft = recoveredSignalReal[i];
Point2d pt = { t, ft };
points.push_back(pt);
t += step;
}
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG), hDlg,
DrawGraphDialog);
}
INT_PTR CALLBACK InputDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
switch (message)
{
case WM_INITDIALOG:
SetDlgItemText(hDlg, IDC_EDIT_THRESHOLD, L"1.0e-2");
SetDlgItemText(hDlg, IDC_EDIT_NOISE, L"0.25");
SetDlgItemText(hDlg, IDC_EDIT_N, L"1024");
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDC_BUTTON_DRAW)
{
GetDlgItemTextA(hDlg, IDC_EDIT_THRESHOLD, thresholdText, 128);
threshold = atof(thresholdText);
GetDlgItemTextA(hDlg, IDC_EDIT_NOISE, noisePCText, 128);
noisePercent = atof(noisePCText);
Npts = GetDlgItemInt(hDlg, IDC_EDIT_N, FALSE, FALSE);
if (LOWORD(wParam) == IDC_BUTTON_DRAW)
{
Filter(threshold, Npts, 1, hDlg);
return (INT_PTR)TRUE;
}
}
if (LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
static void FindMinMax(
double& tMin, double& tMax,
double& fMin, double& fMax)
{
// uses global 2D double points structure
tMin = fMin = DBL_MAX;
tMax = fMax = DBL_MIN;
for (size_t i = 0; i < points.size(); i++)
{
Point2d pt = points[i];
double t = pt.t;
double f = pt.f;
if (t < tMin)
tMin = t;
if (t > tMax)
tMax = t;
if (f < fMin)
fMin = f;
if (f > fMax)
fMax = f;
}
}
static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
// Draw the text with formatting options
DrawTextA(hdc, text, -1, &rect, DT_SINGLELINE | DT_NOCLIP);
}
INT_PTR CALLBACK DrawGraphDialog(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
WCHAR line[256] = { };
switch (message)
{
case WM_INITDIALOG:
SetWindowText(hDlg, L"Graph Dialog");
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;
double tMax = 0, tMin = 0, fMax = 0, fMin = 0;
FindMinMax(tMin, tMax, fMin, fMax);
float tSpan = (float)(tMax - tMin);
float fSpan = (float)(fMax - fMin);
RECT rect = { };
GetClientRect(hDlg, &rect);
float width = (float)(rect.right - rect.left + 1);
float height = (float)(rect.bottom - rect.top - 32 + 1);
float st0 = 2.0f * width / 16.0f;
float st1 = 14.0f * width / 16.0f;
float sf0 = 2.0f * height / 16.0f;
float sf1 = 14.0f * height / 16.0f;
float deltaT = tSpan / 8.0f;
float deltaF = fSpan / 8.0f;
float tSlope = (st1 - st0) / tSpan;
float tInter = (float)(st0 - tSlope * tMin);
float fSlope = (sf0 - sf1) / fSpan;
float fInter = (float)(sf0 - fSlope * fMax);
float pt = 0, pf = 0, st = 0, sf = 0;
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hDlg, &ps);
int i = 0;
float t = (float)tMin;
float f = (float)fMax;
pt = t;
pf = f;
st = tSlope * pt + tInter;
sf = fSlope * pf + fInter;
MoveToEx(hdc, (int)st, (int)sf0, &wPt);
char buffer[128] = { };
while (i <= 8)
{
st = tSlope * t + tInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)st, (int)sf0, &wPt);
LineTo(hdc, (int)st, (int)sf1);
sprintf_s(buffer, "%5.4f", t);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
(int)strlen(buffer),
&size);
RECT textRect = { };
textRect.left = (long)(st - size.cx / 2.0f);
textRect.right = (long)(st + size.cx / 2.0f);
textRect.top = (long)sf1;
textRect.bottom = (long)(sf1 + size.cy / 2.0f);
DrawFormattedText(hdc, buffer, textRect);
t += deltaT;
i++;
}
i = 0;
f = (float)fMin;
while (i <= 8)
{
sf = fSlope * f + fInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)st0, (int)sf, &wPt);
LineTo(hdc, (int)st, (int)sf);
if (i != 0)
{
sprintf_s(buffer, "%+5.3lf", f);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
(int)strlen(buffer),
&size);
RECT textRect = { };
textRect.left = (long)(st0 - size.cx - size.cx / 5.0f);
textRect.right = (long)(st0 - size.cx / 2.0f);
textRect.top = (long)(sf - size.cy / 2.0f);
textRect.bottom = (long)(sf + size.cy / 2.0f);
DrawFormattedText(hdc, buffer, textRect);
}
f += deltaF;
i++;
}
HGDIOBJ bPenNew = NULL;
HGDIOBJ hPenOld = NULL;
bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
pt = (float)points[0].t;
pf = (float)points[0].f;
st = tSlope * pt + tInter;
sf = fSlope * pf + fInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)st, (int)sf, &wPt);
for (size_t j = 1; j < points.size(); j++)
{
pt = (float)points[j].t;
pf = (float)points[j].f;
st = tSlope * pt + tInter;
sf = fSlope * pf + fInter;
LineTo(hdc, (int)st, (int)sf);
}
SelectObject(hdc, hPenOld);
DeleteObject(bPenNew);
return (INT_PTR)FALSE;
}
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_FILTERINGNOISYSIGNAL ICON "FilteringNoisySignal.ico"
IDI_SMALL ICON "small.ico"
/////////////////////////////////////////////////////////////////////////////
//
// Menu
//
IDC_FILTERINGNOISYSIGNAL MENU
BEGIN
POPUP "&Start"
BEGIN
MENUITEM "&Input", IDM_INPUT
MENUITEM SEPARATOR
MENUITEM "E&xit", IDM_EXIT
END
POPUP "&Help"
BEGIN
MENUITEM "&About ...", IDM_ABOUT
END
END
/////////////////////////////////////////////////////////////////////////////
//
// Accelerator
//
IDC_FILTERINGNOISYSIGNAL 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 FilteringNoisySignal"
FONT 8, "MS Shell Dlg"
BEGIN
ICON IDI_FILTERINGNOISYSIGNAL,IDC_STATIC,14,14,21,20
LTEXT "FilteringNoisySignal, 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_INPUT_DIALOG DIALOGEX 0, 0, 310, 120
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Input Dialog"
FONT 10, "Courier New", 700
BEGIN
LTEXT "Threshold", IDC_STATIC, 10, 0, 40, 12
EDITTEXT IDC_EDIT_THRESHOLD, 55, 0, 40, 14, ES_AUTOHSCROLL
LTEXT "Noise %", IDC_STATIC, 10, 15, 40, 12
EDITTEXT IDC_EDIT_NOISE, 55, 15, 40, 14, ES_AUTOHSCROLL
LTEXT "N", IDC_STATIC, 10, 30, 40, 12
EDITTEXT IDC_EDIT_N, 55, 30, 40, 14, ES_AUTOHSCROLL
PUSHBUTTON "Draw", IDC_BUTTON_DRAW, 10, 75, 50, 16
PUSHBUTTON "Cancel", IDCANCEL, 180, 75, 50, 16
END
IDD_DRAW_GRAPH_DIALOG DIALOGEX 0, 0, 410, 310
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Draw Graph Dialog"
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_FILTERINGNOISYSIGNAL "FILTERINGNOISYSIGNAL"
IDS_APP_TITLE "FilteringNoisySignal"
END
#endif
/////////////////////////////////////////////////////////////////////////////
#ifndef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
/////////////////////////////////////////////////////////////////////////////
#endif // not APSTUDIO_INVOKED
//{{NO_DEPENDENCIES}}
// Microsoft Visual C++ generated include file.
// Used by FilteringNoisySignal.rc
#define IDS_APP_TITLE 103
#define IDR_MAINFRAME 128
#define IDD_FILTERINGNOISYSIGNAL_DIALOG 102
#define IDD_ABOUTBOX 103
#define IDM_ABOUT 104
#define IDM_EXIT 105
#define IDI_FILTERINGNOISYSIGNAL 107
#define IDI_SMALL 108
#define IDC_FILTERINGNOISYSIGNAL 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_INPUT_DIALOG 1000
#define IDC_EDIT_THRESHOLD 1010
#define IDC_EDIT_NOISE 1020
#define IDC_EDIT_N 1030
#define IDC_BUTTON_DRAW 1040
#define IDD_DRAW_GRAPH_DIALOG 2000
#define IDM_INPUT 32771
Blog Entry © Thursday – Saturday, June 11 – 13, 2026, by James Pate Williams, Jr. C/C++ Translation of S. D. Conte and Carl de Boor’s Cooley-Tukey Fast Fourier Transform Algorithm and Related Algorithms
#pragma once
#include <complex>
#include <vector>
class Transform
{
public:
static void VandermondeDFT(
int n,
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& y);
static void InverseVandermondeDFT(
int n,
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& y);
static std::vector<std::complex<double>> DFT(
std::vector<double>& x, std::vector<double>& f);
static std::vector<double> InverseDFT(
std::vector<double>& f,
std::vector<std::complex<double>>& X);
/*
* Reference: "Elementary Numerical Analysis:
* An Algorithmic Approach Third Edition" (c)
* 1980 by S. D. Conte and Carl de Boor
* Section 6.5 pages 268 - 277 and Section 6.6
* pages 277 - 283
* Input to FFT
* Z1, Z2 complex n-vectors
* n the length of the vectors
* inzee
* = 1 transform in Z1
* = 2 transform in Z2
* Constructs the discrete Fourier transform in the Cooley-
* Tukey way, but with a twist.
*/
static void FFT(
std::vector<std::complex<double>>& Z1,
int& after, int& now, int& before, int& inzee,
std::vector<std::complex<double>>& Z2);
/*
* This computes an in - place complex - to - complex FFT
* x and y are the real and imaginary arrays of 2^m points.
* dir = 1 gives forward transform
* dir = -1 gives reverse transform
* see http://astronomy.swin.edu.au/~pbourke/analysis/dft/
*/
static void FFT(short dir, int m,
std::vector<double>& x, std::vector<double>& y);
/*
* Reference: "Introduction to Algorithms" by
* Thomas H. Cormen, Charles E. Leiserson, and
* Ronald L. Rivest, pages 794 - 795
*/
static void IterativeFFT(
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& A);
/*
* Reference: "Introduction to Algorithms" by
* Thomas H. Cormen, Charles E. Leiserson, and
* Ronald L. Rivest, page 788
*/
static std::vector<std::complex<double>> RecursiveFFT(
std::vector<std::complex<double>>& a);
};
#include "Transform.h"
void Transform::VandermondeDFT(
int n,
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& y)
{
double pi = 4.0 * atan(1.0);
std::complex<double> z(0.0, 2.0 * pi / n);
std::complex<double> omegaN = exp(z);
std::vector<std::vector<std::complex<double>>> V(n);
for (int k = 0; k < n; k++)
{
V[k].resize(n);
for (int j = 0; j < n; j++)
{
V[k][j] = std::pow(omegaN, k * j);
}
}
for (int k = 0; k < n; k++)
{
std::complex<double> sum = 0.0;
for (int j = 0; j < n; j++)
{
sum += V[k][j] * a[j];
}
y[k] = sum;
}
}
void Transform::InverseVandermondeDFT(
int n,
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& y)
{
double pi = 4.0 * atan(1.0);
std::complex<double> nc = { static_cast<double>(n), 0.0 };
std::complex<double> z(0.0, 2.0 * pi / n);
std::complex<double> omegaN = exp(z);
std::vector<std::vector<std::complex<double>>> invV(n);
for (int k = 0; k < n; k++)
{
invV[k].resize(n);
for (int j = 0; j < n; j++)
{
invV[k][j] = std::pow(omegaN, -k * j);
}
}
for (int k = 0; k < n; k++)
{
std::complex<double> sum = 0.0;
for (int j = 0; j < n; j++)
{
sum += invV[k][j] * y[j];
}
a[k] = sum / nc;
}
}
std::vector<std::complex<double>> Transform::DFT(
std::vector<double>& x, std::vector<double>& f)
{
int length = static_cast<int>(x.size());
double pi = 4.0 * atan(1.0);
double pi2oN = 2.0 * pi / length;
int k, n;
std::vector<double> X(length);
std::vector<double> Y(length);
std::vector<std::complex<double>> Z(length);
f.resize(length);
for (k = 0; k < length; k++)
{
X[k] = Y[k] = 0;
for (n = 0; n < length; n++)
{
X[k] += x[n] * cos(pi2oN * k * n);
Y[k] -= x[n] * sin(pi2oN * k * n);
}
f[k] = pi2oN * k;
X[k] /= length;
Y[k] /= length;
Z[k] = { X[k], Y[k] };
}
return Z;
}
std::vector<double> Transform::InverseDFT(
std::vector<double>& f,
std::vector<std::complex<double>>& X)
{
double imag = 0.0;
int length = static_cast<int>(X.size());
std::vector<double> x(length);
for (int n = 0; n < length; n++)
{
imag = x[n] = 0.0;
for (int k = 0; k < length; k++)
{
x[n] += X[k]._Val[0] * cos(f[k] * n)
- X[k]._Val[1] * sin(f[k] * n);
imag += X[k]._Val[0] * sin(f[k] * n)
+ X[k]._Val[1] * cos(f[k] * n);
}
}
return x;
}
static void FFTStep(
std::vector<std::complex<double>>& Zinp,
int after, int now, int before,
std::vector<std::complex<double>>& Zout)
{
double angle = 0.0, ratio = 0.0;
double twoPi = 2.0 * 4.0 * atan(1.0);
int ia = 0, ib = 0, inp = 0, j = 0;
std::complex<double> arg = 1.0, omega = 0, value = 0;
angle = twoPi / ((now + 1) * (after + 1));
omega = std::complex<double>(cos(angle), -sin(angle));
int address = 1;
for (int i = 1; i <= now; i++)
{
for (int j = 1; j <= after; j++)
{
for (int k = 1; k <= before; k++)
{
address = i * j * k;
if (address < Zout.size())
Zout[address] = { 0.0, 0.0 };
}
}
}
address = 1;
for (int j = 1; j <= now; j++)
{
for (ia = 1; ia <= after; ia++)
{
for (ib = 1; ib <= before; ib++)
{
int address = j * ia * ib;
if (address < Zinp.size())
value = Zinp[address];
for (inp = now - 1; inp >= 1; inp--)
{
address = ia * ib * inp;
if (address < Zinp.size())
value = value * arg + Zinp[address];
}
address = ia * j * ib;
if (address < Zout.size())
Zout[address] = value;
}
arg *= omega;
}
}
}
void Transform::FFT(
std::vector<std::complex<double>>& Z1,
int& after, int& now, int& before, int& inzee,
std::vector<std::complex<double>>& Z2)
{
std::vector<int> prime =
{ 0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 };
int next = 1, nextmx = 25;
after = 1;
before = (int)Z1.size();
now = 1;
Label10:
if (before / prime[next] * prime[next] < before)
{
next++;
if (next <= nextmx)
goto Label10;
else
{
now = before;
before = 1;
}
}
else
{
now = prime[next];
before /= prime[next];
}
if (inzee == 1)
FFTStep(Z1, after, now, before, Z2);
else
FFTStep(Z2, after, now, before, Z1);
inzee = 3 - inzee;
if (before == 1)
return;
after *= now;
goto Label10;
}
void Transform::FFT(short dir, int m,
std::vector<double>& x, std::vector<double>& y)
{
int n, i, i1, j, k, i2, l, l1, l2;
double c1, c2, tx, ty, t1, t2, u1, u2, z;
// Calculate the number of points
n = 1;
for (i = 0; i < m; i++)
n *= 2;
// Do the bit reversal
i2 = n >> 1;
j = 0;
for (i = 0; i < n - 1; i++)
{
if (i < j)
{
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j)
{
j -= k;
k >>= 1;
}
j += k;
}
// Compute the FFT
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l = 0; l < m; l++)
{
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j = 0; j < l1; j++)
{
for (i = j; i < n; i += l2)
{
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}
// Scaling for forward transform
if (dir == 1)
{
for (i = 0; i < n; i++)
{
x[i] /= n;
y[i] /= n;
}
}
}
static void FFTBase(
std::vector<std::complex<double>> a,
std::vector<std::complex<double>> A)
{
double pi = 4.0 * atan(1.0);
int n = static_cast<int>(a.size());
for (int s = 1; s <= log2(n); s++)
{
int m = static_cast<int>(pow(2, s));
std::complex<double> z(0.0, 2.0 * pi / m);
std::complex<double> omegaM = exp(z);
for (int k = 0; k <= n - 1; k += m)
{
std::complex<double> omega = { 1.0, 0.0 };
for (int j = 0; j <= m / 2 - 1; j++)
{
std::complex<double> t = omega * A[k + j + m / 2];
std::complex<double> u = A[k + j];
std::complex<double> jc = { static_cast<double>(j), 0.0 };
A[k + j] = u + jc;
A[k + j + m / 2] = u - t;
omega *= omegaM;
}
}
}
}
static int Reverse(int k)
{
int digits[32] = { 0 }, i = 0;
while (k > 0)
{
int digit = k & 1;
k >>= 1;
digits[i++] = digit;
}
int result = digits[0];
for (int j = 1; j < i; j++)
result = result * 2 + digits[j];
return result;
}
static void BitReverseCopy(
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& A)
{
int n = static_cast<int>(a.size());
for (int k = 0; k <= n - 1; k++)
A[Reverse(k)] = a[k];
}
void Transform::IterativeFFT(
std::vector<std::complex<double>>& a,
std::vector<std::complex<double>>& A)
{
BitReverseCopy(a, A);
double pi = 4.0 * atan(1.0);
int n = static_cast<int>(a.size());
for (int s = 1; s <= static_cast<int>(log2(n)); s++)
{
int m = static_cast<int>(pow(2.0, s));
std::complex<double> z(0.0, 2.0 * pi / m);
std::complex<double> omegaM = exp(z);
std::complex<double> omega = { 1.0, 0.0 };
for (int j = 0; j <= m / 2 - 1; j++)
{
for (int k = j; k <= n - 1; k += m)
{
std::complex<double> t = omega * A[k + m / 2];
std::complex<double> u = A[k];
A[k] = u + t;
A[k + m / 2] = u - t;
omega *= omegaM;
}
}
}
}
std::vector<std::complex<double>> Transform::RecursiveFFT(
std::vector<std::complex<double>>& a)
{
int n = static_cast<int>(a.size());
if (n == 1)
return a;
std::vector<std::complex<double>> a0;
std::vector<std::complex<double>> a1;
std::vector<std::complex<double>> y0;
std::vector<std::complex<double>> y1;
for (int i = 0; i <= n - 2; i++)
a0.push_back(a[i]);
for (int i = 1; i <= n - 1; i++)
a1.push_back(a[i]);
y0 = RecursiveFFT(a0);
y1 = RecursiveFFT(a1);
double pi = 4.0 * atan(1.0);
std::complex<double> z(0.0, 2.0 * pi / n);
std::complex<double> omegaN = exp(z);
std::complex<double> omega(1.0, 0.0);
std::vector<std::complex<double>> y(n, 0.0);
for (int k = 0; k <= n / 2 - 1; k++)
{
y[k] = y0[k] + omega * y1[k];
y[(long long)k + n / 2] = y0[k] - omega * y1[k];
omega *= omegaN;
}
return y;
}
// CooleyTukeyConsole.cpp : This file contains the 'main' function. Program execution begins and ends there.
//
#include <algorithm>
#include <complex>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
#include "Transform.h"
int M = 0, N = 0, Sn = 2048;
static double SimpsonsRule(
double a, double b, int n,
double (*f)(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 += f(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < n; i += 2)
{
t += f(x);
x += h2;
}
return h * (f(a) + 4 * s + 2 * t + f(b)) / 3.0;
}
static double f(double x)
{
return x * x * sin(x);
}
static double AMf(double x)
{
return f(x) * cos(M * x);
}
static double BMf(double x)
{
return f(x) * sin(M * x);
}
static double AM(int m)
{
M = m;
double pi = 4.0 * atan(1.0), twoPi = pi + pi;
double integ = SimpsonsRule(-pi, pi, Sn, AMf);
double value = integ / twoPi;
if (m == 0)
return value;
else
return 2.0 * value;
}
static double BM(int m)
{
M = m;
double pi = 4.0 * atan(1.0), twoPi = pi + pi;
double integ = SimpsonsRule(-pi, pi, Sn, BMf);
return 2.0 * integ / twoPi;
}
static void FormatPrint(double x)
{
if (fabs(x) < 1.0e-12)
x = 0.0;
std::cout << std::setw(13);
std::cout << std::setfill(' ');
std::cout << std::setprecision(10);
std::cout << x << '\t';
}
static void FormatPrint(std::complex<double> z)
{
std::cout << std::setw(13);
std::cout << std::setfill(' ');
std::cout << std::setprecision(10);
std::cout << z << '\t';
}
static void GetCooleyTukeyData(
int n0, int n1, int n2, int n3,
std::vector<double>& a,
std::vector<double>& b,
std::vector<double>& p,
std::vector<double>& x,
std::vector<double>& fx)
{
double pi = 4.0 * atan(1.0), twoPi = pi + pi;
x[0] = fx[0] = 0.0;
for (int i = 1; i <= n0; i++)
{
x[i] = twoPi * i / n0;
fx[i] = f(x[i]);
}
a[0] = AM(0);
for (int m = 1; m <= n0; m++)
{
a[m] = AM(m);
b[m] = BM(m);
}
for (int n = 0; n <= n0; n++)
{
double asum = 0.0;
for (int m = 1; m < n0; m++)
asum += a[m] * cos(m * x[n]);
double bsum = 0.0;
for (int m = 1; m <= n0; m++)
bsum += b[m] * sin(m * x[n]);
p[n] = a[0] / 2.0 + asum + bsum;
}
}
static void GetData(
int n0,
std::vector<double>& a,
std::vector<double>& b,
std::vector<double>& p,
std::vector<double>& x,
std::vector<double>& fx)
{
double pi = 4.0 * atan(1.0), twoPi = pi + pi;
for (int i = 0; i < n0; i++)
{
x[i] = twoPi * i / n0;
fx[i] = f(x[i]);
}
a[0] = AM(0);
for (int m = 1; m < n0; m++)
{
a[m] = AM(m);
b[m] = BM(m);
}
for (int n = 0; n < n0; n++)
{
double asum = 0.0;
for (int m = 1; m < n0; m++)
asum += a[m] * cos(m * x[n]);
double bsum = 0.0;
for (int m = 1; m < n0; m++)
bsum += b[m] * sin(m * x[n]);
p[n] = a[0] / 2.0 + asum + bsum;
}
}
static void TestCooleyTukey(
int n0, int n1, int n2, int n3)
{
int n01 = n0 + 1;
std::vector<double> a(n01), b(n01), p(n01), x(n01), fx(n01);
GetCooleyTukeyData(n0, n1, n2, n3, a, b, p, x, fx);
double pi = 4.0 * atan(1.0), twoPi = pi + pi;
int index1 = 1, inzee = 1;
std::vector<std::complex<double>> pp(n01);
std::vector<std::complex<double>> Z1 = { {0, 0} };
std::vector<std::complex<double>> Z2 = { {0, 0} };
Z1.resize(n01, 0.0);
Z2.resize(n01, 0.0);
for (int i = 1; i <= n0; i++)
Z1[i] = fx[i];
int after = 0, now = 0, before = 0;
Transform::FFT(Z1, after, now, before, inzee, Z2);
std::cout << "Forward Transform" << std::endl;
for (int i = 1; index1 < n01 && i <= after; i++)
{
for (int j = 1; index1 < n01 && j <= now; j++)
{
for (int k = 1; index1 < n01 && k <= before; k++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << index1 << '\t';
FormatPrint(x[index1]);
FormatPrint(fx[index1]);
FormatPrint(Z1[index1]);
FormatPrint(Z2[index1]);
std::cout << std::endl;
index1++;
}
}
}
}
static void TestFFT(int n0)
{
int m = static_cast<int>(log2(n0));
std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0);
GetData(n0, a, b, p, x, fx);
std::vector<double> xx(n0), yy(n0);
for (int i = 0; i < n0; i++)
xx[i] = fx[i];
Transform::FFT(+1, m, xx, yy);
std::cout << "Forward Transform" << std::endl;
for (int i = 0; i < n0; i++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << i << '\t';
FormatPrint(x[i]);
FormatPrint(fx[i]);
FormatPrint(p[i]);
FormatPrint(xx[i]);
FormatPrint(yy[i]);
std::cout << std::endl;
if (i == n0 / 2)
break;
}
Transform::FFT(-1, m, xx, yy);
std::cout << "Backward Transform" << std::endl;
for (int i = 0; i < n0; i++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << i << '\t';
FormatPrint(x[i]);
FormatPrint(fx[i]);
FormatPrint(p[i]);
FormatPrint(xx[i]);
FormatPrint(yy[i]);
std::cout << std::endl;
if (i == n0 / 2)
break;
}
}
static void TestIterativeFFT(int n0)
{
std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0);
GetData(n0, a, b, p, x, fx);
double pi = 4.0 * atan(1.0), twoPi = pi + pi;
std::vector<std::complex<double>> AA(n0);
std::vector<std::complex<double>> aa(n0);
std::vector<std::complex<double>> yy(n0);
for (int i = 0; i < n0; i++)
aa[i] = fx[i];
Transform::IterativeFFT(aa, AA);
std::cout << "Forward Transform" << std::endl;
for (int i = 0; i < n0; i++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << i << '\t';
FormatPrint(x[i]);
FormatPrint(fx[i]);
FormatPrint(AA[i]);
std::cout << std::endl;
}
}
static void TestRecursiveFFT(int n0)
{
std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0);
std::vector<double> xx(n0), yy(n0);
GetData(n0, a, b, p, x, fx);
double pi = 4.0 * atan(1.0), twoPi = pi + pi;
std::vector<std::complex<double>> A = { {0, 0} };
std::vector<std::complex<double>> Y = { {0, 0} };
A.resize(n0);
for (int i = 0; i < n0; i++)
A[i] = fx[i];
Y = Transform::RecursiveFFT(A);
std::cout << "Forward Transform" << std::endl;
for (int i = 0; i < n0; i++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << i << '\t';
FormatPrint(x[i]);
FormatPrint(fx[i]);
FormatPrint(A[i]);
FormatPrint(Y[i]);
std::cout << std::endl;
}
}
static void TestDFT(int n0)
{
std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0), ff(n0);
GetData(n0, a, b, p, x, fx);
std::vector<std::complex<double>> zz = Transform::DFT(fx, ff);
std::cout << "Forward Transform" << std::endl;
for (int i = 0; i < n0; i++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << i << '\t';
FormatPrint(x[i]);
FormatPrint(fx[i]);
FormatPrint(ff[i]);
FormatPrint(zz[i]);
std::cout << std::endl;
}
std::cout << "Inverse Transform" << std::endl;
std::vector<double> inv = Transform::InverseDFT(ff, zz);
for (int i = 0; i < n0; i++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << i << '\t';
FormatPrint(x[i]);
FormatPrint(fx[i]);
FormatPrint(ff[i]);
FormatPrint(inv[i]);
std::cout << std::endl;
}
}
static void TestVandermondeDFT(int n0)
{
std::vector<double> a(n0), b(n0), p(n0), x(n0), fx(n0);
GetData(n0, a, b, p, x, fx);
std::vector<std::complex<double>> aa(n0);
std::vector<std::complex<double>> yy(n0);
for (int i = 0; i < n0; i++)
aa[i] = fx[i];
Transform::VandermondeDFT(n0, aa, yy);
std::cout << "Vandermonde DFT" << std::endl;
for (int i = 0; i < n0; i++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << i << '\t';
FormatPrint(x[i]);
FormatPrint(fx[i]);
FormatPrint(p[i]);
FormatPrint(yy[i]);
std::cout << std::endl;
}
Transform::InverseVandermondeDFT(n0, aa, yy);
std::cout << "Inverse Vandermonde DFT" << std::endl;
for (int i = 0; i < n0; i++)
{
std::cout << std::setw(5);
std::cout << std::setfill(' ');
std::cout << i << '\t';
FormatPrint(x[i]);
FormatPrint(fx[i]);
FormatPrint(p[i]);
FormatPrint(aa[i]);
std::cout << std::endl;
}
}
static int Horner(char line[])
{
int length = static_cast<int>(strlen(line));
int sum = line[0] - '0';
for (int i = 1; i < length; i++)
sum = sum * 10 + line[i] - '0';
return sum;
}
int main()
{
char line[128];
while (true)
{
std::cout << "== Menu ==" << std::endl;
std::cout << "1 Cooley-Tukey" << std::endl;
std::cout << "2 FFT" << std::endl;
std::cout << "3 Iterative FFT" << std::endl;
std::cout << "4 Recursive FFT" << std::endl;
std::cout << "5 DFT" << std::endl;
std::cout << "6 Vandermonde DFT" << std::endl;
std::cout << "7 Exit" << std::endl;
std::cout << "Option 1 - 7 = ";
std::cin.getline(line, 128);
int option = Horner(line);
if (option == 7)
break;
if (option < 1 || option > 7)
{
std::cout << "Unknown Option Number" << std::endl;
continue;
}
if (option == 1)
{
int n0 = 0, n1 = 0, n2 = 0, n3 = 0;
std::cout << "n1 = ";
std::cin.getline(line, 128);
n1 = Horner(line);
std::cout << "n2 = ";
std::cin.getline(line, 128);
n2 = Horner(line);
std::cout << "n3 = ";
std::cin.getline(line, 128);
n3 = Horner(line);
n0 = n1 * n2 * n3;
TestCooleyTukey(n0, n1, n2, n3);
}
else if (option == 2)
{
std::cout << "n0 = ";
std::cin.getline(line, 128);
int n0 = Horner(line);
if (n0 % 2 != 0)
{
std::cout << "n0 must be a power of 2";
std::cout << std::endl;
continue;
}
TestFFT(n0);
}
else if (option == 3)
{
std::cout << "n0 = ";
std::cin.getline(line, 128);
int n0 = Horner(line);
if (n0 % 2 != 0)
{
std::cout << "n0 must be a power of 2";
std::cout << std::endl;
continue;
}
TestIterativeFFT(n0);
}
else if (option == 4)
{
std::cout << "n0 = ";
std::cin.getline(line, 128);
int n0 = Horner(line);
if (n0 % 2 != 0)
{
std::cout << "n0 must be a power of 2";
std::cout << std::endl;
continue;
}
TestRecursiveFFT(n0);
}
else if (option == 5)
{
std::cout << "n0 = ";
std::cin.getline(line, 128);
int n0 = Horner(line);
if (n0 % 2 != 0)
{
std::cout << "n0 must be a power of 2";
std::cout << std::endl;
continue;
}
TestDFT(n0);
}
else if (option == 6)
{
std::cout << "n0 = ";
std::cin.getline(line, 128);
int n0 = Horner(line);
if (n0 % 2 != 0)
{
std::cout << "n0 must be a power of 2";
std::cout << std::endl;
continue;
}
TestVandermondeDFT(n0);
}
}
return 0;
}
Blog Entry © Monday to Tuesday, June 8 – 9, 2026, by James Pate Williams, Jr., DampedNewton’s Method for a System of Equations
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 © 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;
}
Numerical Explorations 