Category: Computer Science
Blog Entry © Thursday, July 31, 2025, Numerically Solving Four Second Order Linear Ordinary Differential Equation Boundary Value Problems by James Pate Williams, Jr.
Blog Entry (c) Wednesday, July 30, 2025, by James Pate Williams, Jr. Four Root Finding Algorithms and C Source
function exp(-x) - sin(0.5 * pi * x)
The interval found by bisection is:
a = +0.4435735341 b = +0.4435735341
bisection # iterations = 50
The root by regula falsi:
x = +0.4435735341 f(x) = -0.0000000000
regula falsi # iterations = 8
The root by the secant method:
x = +0.4435735341 f(x) = +0.0000000000
secant method # iterations = 9
The root by Newton's Method:
x = +0.4435735341 f(x) = -0.0000000000
Newton's Method # iterations = 6
D:\roots\x64\Release\roots.exe (process 28152) exited with code 0 (0x0).
Press any key to close this window . . .
#include <math.h>
#include <stdio.h>
typedef double real;
static real f(real x)
{
double pi2 = 2.0 * atan(1.0);
return(exp(-x) - sin(pi2 * x));
}
static real g(real x)
{
double pi2 = 2.0 * atan(1.0);
return(-exp(-x) - pi2 * cos(pi2 * x));
}
static int bisection(
real(*f)(real), real* a, real* b, real xtol, int* flag)
{
real error, fa, fm, xm;
int n = 0;
fa = (*f)(*a);
if (fa * (*f)(*b) > 0.0)
{
*flag = -1;
return(n);
}
error = fabsl(*b - *a);
while (error > xtol)
{
n++;
error *= 0.5;
if (error <= xtol)
{
*flag = 0;
return(n);
}
xm = 0.5 * (*a + *b);
if (xm + error == xm)
{
*flag = 1;
return(n);
}
fm = (*f)(xm);
if (fa * fm > 0.0)
{
*a = xm;
fa = fm;
}
else
*b = xm;
}
*flag = 2;
return(n);
}
static int regula(
real (*f)(real), real a, real b, real xtol,
real ftol, int ntol, real* w, int* flag)
{
int n = 0;
real fa, fb, fw, signfa, prvsfw;
fa = f(a);
if (fa >= 0.0) signfa = +1.0; else signfa = -1.0;
fb = f(b);
if (signfa * fb >= 0.0)
{
*flag = -1;
return n;
}
*w = a;
fw = fa;
for (n = 0; n <= ntol; n++)
{
if (fabs(a - b) <= xtol)
{
*flag = 0;
return n;
}
if (fabs(fw) <= ftol)
{
*flag = 1;
return n;
}
*w = (fa * b - fb * a) / (fa - fb);
if (fw >= 0.0) prvsfw = +1.0; else prvsfw = -1.0;
fw = f(*w);
if (signfa * fw > 0.0)
{
a = *w;
fa = fw;
if (fw * prvsfw > 0.0) fb = 0.5 * fb;
}
else
{
b = *w;
fb = fw;
if (fw * prvsfw > 0.0) fa = 0.5 * fa;
}
}
*flag = 2;
return n;
}
static int secantMethod(
real(*f)(real), real xtol, real ftol,
int ntol, real xm1, real x0, real* w)
{
real fm1 = f(xm1), f0 = f(x0);
real df = fabs(fm1 - f0), f1 = 0.0;
real dx = fabs(xm1 - x0), x1 = 0.0;
int n = 0;
while (n < ntol && df > ftol && dx > xtol)
{
x1 = (f0 * xm1 - fm1 * x0) / (f0 - fm1);
f1 = f(x1);
df = fabs(f(x1) - f0);
dx = fabs(x1 - x0);
fm1 = f0;
f0 = f1;
xm1 = x0;
x0 = x1;
n++;
}
*w = x1;
return n;
}
static int NewtonsMethod(
real(*f)(real), real(*g)(real),
real xtol, real ftol, int ntol,
real* w)
{
// f is the function
// g is the function's derivative
// xtol is root's tolerance
// ftol is the function's tolerance
// ntol is the maximum # of iterations
real f1 = 0.0, g1 = 0.0, x0 = *w, x1 = 0.0;
real f0 = f(x0), g0 = g(x0);
real deltaX = DBL_MAX, deltaF = DBL_MAX;
int n = 0;
while (n < ntol && deltaX > xtol && deltaF > ftol)
{
x1 = x0 - f0 / g0;
f1 = f(x1);
g1 = g(x1);
deltaX = fabs(x1 - x0);
deltaF = fabs(f1 - f0);
f0 = f1;
g0 = g1;
x0 = x1;
n++;
}
*w = x1;
return n;
}
int main(void)
{
int flag = 0, ntol = 0;
real a0 = 0, b0 = 0, ftol = 0, w1 = 0, w2 = 1.0;
real a1 = 0, b1 = 0, w3 = 0, xtol = 0;
a0 = 0.0;
b0 = 1.0;
ntol = 128;
ftol = 1.0e-15;
xtol = 1.0e-15;
int its1 = bisection(f, &a0, &b0, xtol, &flag);
a1 = 0.0;
b1 = 1.0;
int its2 = regula(f, a1, b1, xtol, ftol, ntol, &w1, &flag);
int its3 = secantMethod(f, ftol, xtol, ntol, 0.0, 1.0, &w2);
int its4 = NewtonsMethod(f, g, xtol, ftol, ntol, &w3);
printf("function exp(-x) - sin(0.5 * pi * x)\n");
printf("The interval found by bisection is:\n");
printf("a = %+13.10lf b = %+13.10lf\n", a0, b0);
printf("bisection # iterations = %ld\n", its1);
printf("The root by regula falsi:\n");
printf("x = %+13.10lf f(x) = %+13.10lf\n", w1, f(w1));
printf("regula falsi # iterations = %ld\n", its2);
printf("The root by the secant method:\n");
printf("x = %+13.10lf f(x) = %+13.10lf\n", w2, f(w2));
printf("secant method # iterations = %ld\n", its3);
printf("The root by Newton's Method:\n");
printf("x = %+13.10lf f(x) = %+13.10lf\n", w2, f(w3));
printf("Newton's Method # iterations = %ld\n", its4);
return(0);
}
Blog Entry © Tuesday, July 29, 2025, Double and Triple Monte Carlo Integration by James Pate Williams, Jr.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
static double randomRange(double lo, double hi)
{
return (hi - lo) * (double)rand() / RAND_MAX + lo;
}
static double integrand(double r, double w)
{
return pow(r, 4.0) * (2.0 - r) * w * w * exp(-r);
}
static double StarkEffectIntegral(double E, int N)
{
double sum = 0.0;
for (int i = 0; i <= N; i++)
{
double r = randomRange(0.0, 100.0);
double w = randomRange(-1.0, 1.0);
sum += integrand(r, w);
}
return 100.0 * 2.0 * E * sum / (16.0 * (N - 1));
}
static void firstOrderStarkEffect(double E)
{
double exact = -3.0 * E;
int N[9] = {
1000000, 2000000, 3000000, 4000000,
5000000, 6000000, 7000000, 8000000,
9000000 };
for (int n = 0; n < 9; n++)
{
int iN = N[n];
double integ = StarkEffectIntegral(E, iN);
double error = 100.0 * fabs(integ - exact) / fabs(exact);
printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
iN, integ, error);
}
printf("exact value = %13.10lf\n", exact);
}
static double ee1(int N, double R, double Z)
{
double pi = 4.0 * atan(1.0);
double sum = 0.0;
for (int i = 0; i <= N; i++)
{
double r1 = randomRange(1.0e-25, R);
double r2 = randomRange(0.0, r1);
sum += R * r1 * r1 * exp(-2.0 * Z * (r1 + r2)) * r2 * r2;
}
return 16.0 * pi * pi * sum / (N - 1);
}
static double ee2(int N, double R, double Z)
{
double pi = 4.0 * atan(1.0);
double sum = 0.0;
for (int i = 0; i <= N; i++)
{
double r1 = randomRange(1.0e-25, R);
double r2 = randomRange(r1, R);
sum += R * (R - r2) * r2 * exp(-2.0 * Z * (r1 + r2)) * r1 * r1;
}
return 16.0 * pi * pi * sum / (N - 1);
}
static void firstOrderHelium(double Z)
{
double pi = 4.0 * atan(1.0), R = 25.0;
double exact = 5.0 * pi * pi / (8.0 * pow(Z, 5.0));
int N[9] = {
1000000, 2000000, 3000000, 4000000,
5000000, 6000000, 7000000, 8000000,
9000000 };
for (int n = 0; n < 9; n++)
{
int iN = N[n];
double integ = ee1(iN, R, Z) + ee2(iN, R, Z);
double error = 100.0 * fabs(integ - exact) / fabs(exact);
printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
iN, integ, error);
}
printf("exact value = %13.10lf\n", exact);
}
int main(void)
{
firstOrderStarkEffect(2.0);
firstOrderHelium(2.0);
return 0;
}
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
static double randomRange(double lo, double hi)
{
return (hi - lo) * (double)rand() / RAND_MAX + lo;
}
static double f(double x, double y, double z)
{
return pow(sin(x), 2.0) + y * sin(z);
}
static double g(double x, double y, double z)
{
return x + y * z * z;
}
static double integral(
double x0, double x1,
double y0, double y1,
double z0, double z1,
double (*f)(double, double, double),
int N)
{
double sum = 0.0;
for (int n = 0; n <= N; n++)
{
double x = randomRange(x0, x1);
double y = randomRange(y0, y1);
double z = randomRange(z0, z1);
sum += f(x, y, z);
}
return (x1 - x0) * (y1 - y0) * (z1 - z0) *
sum / (N - 1);
}
int main(void)
{
double pi = 4.0 * atan(1.0);
double x0 = 0.0, x1 = pi;
double y0 = 0.0, y1 = 1.0;
double z0 = 0.0, z1 = pi;
double exact = 0.5 * pi * (2.0 + pi);
int N[9] = {
1000000, 2000000, 3000000, 4000000,
5000000, 6000000, 7000000, 8000000,
9000000 };
printf("integrand pow(sin(x), 2.0) + y * sin(z)\n");
printf("x = 0 to pi, y = 0 to 1, z = 0 to pi\n");
for (int n = 0; n < 9; n++)
{
int iN = N[n];
double integ = integral(
x0, x1, y0, y1, z0, z1, f, iN);
double error = 100.0 * fabs(integ - exact) / fabs(exact);
printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
iN, integ, error);
}
printf("exact value = %13.10lf\n", exact);
x0 = -1.0;
x1 = 5.0;
y0 = 2.0;
y1 = 4.0;
z0 = 0.0;
z1 = 1.0;
exact = 36.0;
printf("integrand x + y * z * z\n");
printf("x = -1 to 5, y = 2 to 4, z = 0 to 1\n");
for (int n = 0; n < 9; n++)
{
int iN = N[n];
double integ = integral(
x0, x1, y0, y1, z0, z1, g, iN);
double error = 100.0 * fabs(integ - exact) / fabs(exact);
printf("N = %4ld\tintegral = %13.10lf\t%% error = %13.10lf\n",
iN, integ, error);
}
printf("exact value = %13.10lf\n", exact);
return 0.0;
}
Blog Entry © Sunday, July 27, 2025, A Bit of Programming Nostalgia Prime Number Related Programs by James Williams, Jr.
/*
Author: Pate Williams c 1995
The following program is a solution to problem 18.15
in Pascalgorithms by Edwin D. Reilly and Francis D.
Federighi page 627. The program uses Simpson's rule
to calculate the number of primes less than or equal
a given number.
*/
#include <math.h>
#include <stdio.h>
typedef double real;
static real f(real x)
{
return(1.0 / log(x));
}
static real simpson(int n, real a, real b)
{
int i;
real evensum, h, oddsum, twoh, x;
if (n % 2 == 1) n = n - 1;
h = (b - a) / n;
twoh = h + h;
x = h + a;
oddsum = 0.0;
for (i = 1; i <= n / 2; i++)
{
oddsum += f(x);
x = twoh + x;
}
x = twoh + a;
evensum = 0.0;
for (i = 1; i <= n / 2 - 1; i++)
{
evensum += f(x);
x = twoh + x;
}
return(h / 3.0 * (f(a) + f(b) + 4.0 * oddsum + 2.0 * evensum));
}
int main(void)
{
int i, n, Nmaximum = 0, Nminimum = 0, Nstep;
printf("n = "); scanf_s("%d", &n);
printf("N minimum = "); scanf_s("%d", &Nminimum);
printf("N maximum = "); scanf_s("%d", &Nmaximum);
printf("N step = "); scanf_s("%d", &Nstep);
printf("\n");
printf("----------------------------------------\n");
printf("Min\t\tMax\t\tprimes\n");
printf("----------------------------------------\n");
for (i = Nminimum; i <= Nmaximum; i += Nstep)
{
printf("%8d\t%8d\t%8.0lf\n", Nminimum, i + Nstep,
simpson(n, Nminimum, i + Nstep));
}
printf("----------------------------------------\n");
return(0);
}
n = 1024
N minimum = 0
N maximum = 10000000
N step = 1000000
----------------------------------------
Min Max primes
----------------------------------------
0 1000000 78551
0 2000000 148923
0 3000000 216788
0 4000000 283122
0 5000000 348361
0 6000000 412754
0 7000000 476461
0 8000000 539590
0 9000000 602224
0 10000000 664424
0 11000000 726239
----------------------------------------
D:\PrimeCounter\x64\Release\PrimeCounter.exe (process 51884) exited with code 0 (0x0).
Press any key to close this window . . .
/*
Author: Pate Williams c 1995
The following is a translation of the Pascal program
sieve found in Pascalgorithms by Edwin D. Reilly and
Francis D. Federighi page 652. This program uses sets
to represent the sieve (see C Programming Language An
Applied Perspective by Lawrence Miller and Alec Qui-
lici pages 160 - 162).
*/
#include <math.h>
#include <stdio.h>
#define _WORD_SIZE 32
#define _VECT_SIZE 524288
#define SET_MIN 0
#define SET_MAX 16777215
typedef unsigned long SET[_VECT_SIZE];
typedef long ELEMENT;
typedef unsigned long LONG;
SET set;
static int get_bit_pos(int* long_ptr, int* bit_ptr,
ELEMENT element)
{
*long_ptr = element / _WORD_SIZE;
*bit_ptr = element % _WORD_SIZE;
return(element >= SET_MIN && element <= SET_MAX);
}
static void set_bit(ELEMENT element, int inset)
{
int bit, word;
if (get_bit_pos(&word, &bit, element))
{
if (inset > 0)
set[word] |= (01 << bit);
else
set[word] &= ~(01 << bit);
}
}
static int get_bit(ELEMENT element)
{
int bit, word;
return(get_bit_pos(&word, &bit, element) ?
(set[word] >> bit) & 01 : 0);
}
static void set_Add(ELEMENT element)
{
set_bit(element, 1);
}
static void set_Del(ELEMENT element)
{
set_bit(element, 0);
}
static int set_Mem(ELEMENT element)
{
return get_bit(element);
}
static void primes(long n)
{
long c, i, inc, k;
double x;
set_Add(2);
for (i = 3; i <= n; i++)
if ((i + 1) % 2 == 0)
set_Add(i);
else
set_Del(i);
c = 3;
do
{
i = c * c;
inc = c + c;
while (i <= n)
{
set_Del(i);
i = i + inc;
}
c += 2;
while (set_Mem(c) == 0) c += 1;
} while (c * c <= n);
k = 0;
for (i = 2; i <= n; i++)
if (set_Mem(i) == 1) k++;
x = n / log(n) - 5.0;
x = x + exp(1.0 + 0.15 * log(n) * sqrt(log(n)));
printf("%8ld\t%8ld\t%8.0lf\n", n, k, x);
}
int main(void)
{
long n = 100L;
printf("----------------------------------------\n");
printf("n\t\tprimes\t\ttheory\n");
printf("----------------------------------------\n");
do
{
primes((int)n);
n = 10L * n;
} while (n < (long)SET_MAX);
printf("----------------------------------------\n");
return(0);
}
----------------------------------------
n primes theory
----------------------------------------
100 25 29
1000 168 181
10000 1229 1261
100000 9592 9634
1000000 78498 78396
10000000 664579 665060
----------------------------------------
D:\Sieve\x64\Release\Sieve.exe (process 60092) exited with code 0 (0x0).
Press any key to close this window . . .
Blog Entry © Sunday, July 27, 2025, A Bit of Programming Nostalgia by James Williams, Jr.
Blog Entry Thursday, July 24, 2025, © James Pate Williams, Jr. Comparison of Data Between Two Sources
Blog Entry © Sunday, July 20, 2025, by James Pate Williams, Jr. Example and Two Exercises from a Numerical Analysis Textbook
Happy 59th Anniversary of the historic moon landing by Neil Armstrong, Buzz Aldrin, and Michael Collins in Apolla 11.
Newtonian Gravity and Cosmology © Friday, July 18, 2025, by James Pate Williams, Jr. in Conjunction with the Microsoft Copilot AI
Some General Relativity Mathematics by James Pate Williams, Jr. (c) July 13, 2025
// SomeRelativityMath.cpp : Defines the entry point for the application.
// https://arxiv.org/pdf/2506.09946
#include "pch.h"
#include "framework.h"
#include "SomeRelativityMath.h"
#include <stdio.h>
#include <vector>
#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
char text[8192];
std::vector<POINT2D> points1, points2;
std::vector<POINT2D> points3, points4;
// 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 DrawGraph(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_SOMERELATIVITYMATH, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_SOMERELATIVITYMATH));
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_SOMERELATIVITYMATH));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_SOMERELATIVITYMATH);
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 void FindMinMax(
double& xMin, double& xMax,
double& yMin, double& yMax)
{
// uses global 2D double points structure
xMin = yMin = DBL_MAX;
xMax = yMax = DBL_MIN;
for (size_t i = 0; i < points1.size(); i++)
{
POINT2D pt = points1[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;
}
for (size_t i = 0; i < points2.size(); i++)
{
POINT2D pt = points2[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;
}
for (size_t i = 0; i < points3.size(); i++)
{
POINT2D pt = points3[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;
}
for (size_t i = 0; i < points4.size(); i++)
{
POINT2D pt = points4[i];
double x = pt.x;
double y = pt.y;
if (x < xMin)
xMin = x;
if (x > xMax)
xMax = x;
if (y < yMin)
yMin = y;
if (y > yMax)
yMax = y;
}
}
static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
// Draw the text with formatting options
DrawTextA(hdc, text, -1, &rect, DT_SINGLELINE | DT_NOCLIP);
}
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
switch (message)
{
case WM_CREATE:
{
double P1 = 0, P2 = 0, P3 = 0, P4 = 0;
int numberPts = 1024;
double x0 = 1, h = 0.5 / numberPts, x = x0;
points1.resize(numberPts);
points2.resize(numberPts);
points3.resize(numberPts);
points4.resize(numberPts);
for (int i = 0; i < numberPts; i++)
{
P1 = x + 1.0;
P2 = 3.0 * x * x + 4.0 * x + 2.0;
P3 = 15.0 * x * x * x + 16.0 * x * x + 6.0 * x;
P4 = 105.0 * x * x * x * x + 156.0 * x * x * x +
94.0 * x * x + 24.0 * x;
POINT2D pt = { 0 };
pt.x = x;
pt.y = log(P1);
points1[i] = pt;
pt.y = log(P2);
points2[i] = pt;
pt.y = log(P3);
points3[i] = pt;
pt.y = log(P4);
points4[i] = pt;
x += h;
}
break;
}
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
double h = 0, pi = 0, plm = 0, theta = 0;
double xMax = 0, xMin = 0, yMax = 0, yMin = 0;
FindMinMax(xMin, xMax, yMin, yMax);
float xSpan = (float)(xMax - xMin);
float ySpan = (float)(yMax - yMin);
RECT rect = { };
GetClientRect(hWnd, &rect);
float width = (float)(rect.right - rect.left + 1);
float height = (float)(rect.bottom - rect.top - 32 + 1);
float sx0 = 2.0f * width / 16.0f;
float sx1 = 14.0f * width / 16.0f;
float sy0 = 2.0f * height / 16.0f;
float sy1 = 14.0f * height / 16.0f;
float deltaX = xSpan / 8.0f;
float deltaY = ySpan / 8.0f;
float xSlope = (sx1 - sx0) / xSpan;
float xInter = (float)(sx0 - xSlope * xMin);
float ySlope = (sy0 - sy1) / ySpan;
float yInter = (float)(sy0 - ySlope * yMax);
float px = 0, py = 0, sx = 0, sy = 0;
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hWnd, &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, "%3.2f", x);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
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, "%3.1f", y);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
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 rPenNew = NULL;
HGDIOBJ hPenOld = NULL;
rPenNew = CreatePen(PS_SOLID, 2, RGB(255, 0, 0));
hPenOld = SelectObject(hdc, rPenNew);
px = (float)points1[0].x;
py = (float)points1[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 < points1.size(); j++)
{
px = (float)points1[j].x;
py = (float)points1[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(rPenNew);
HGDIOBJ gPenNew = CreatePen(PS_SOLID, 2, RGB(0, 255, 0));
hPenOld = SelectObject(hdc, gPenNew);
px = (float)points2[0].x;
py = (float)points2[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 < points2.size(); j++)
{
px = (float)points2[j].x;
py = (float)points2[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(gPenNew);
HGDIOBJ bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
px = (float)points3[0].x;
py = (float)points3[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 < points3.size(); j++)
{
px = (float)points3[j].x;
py = (float)points3[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(bPenNew);
HGDIOBJ pPenNew = CreatePen(PS_SOLID, 2, RGB(128, 0, 128));
hPenOld = SelectObject(hdc, pPenNew);
px = (float)points4[0].x;
py = (float)points4[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 < points4.size(); j++)
{
px = (float)points4[j].x;
py = (float)points4[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(pPenNew);
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 INT_PTR CALLBACK DrawGraph(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
char line[256] = { };
switch (message)
{
case WM_INITDIALOG:
GetWindowTextA(hDlg, line, 256);
strcat_s(text, 8192, " ");
strcat_s(text, 8192, line);
SetWindowTextA(hDlg, text);
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(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 = (sx1 - sx0) / xSpan;
float xInter = (float)(sx0 - xSlope * xMin);
float ySlope = (sy0 - sy1) / ySpan;
float yInter = (float)(sy0 - ySlope * yMax);
float px = 0, py = 0, sx = 0, sy = 0;
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hDlg, &ps);
int i = 0;
float x = (float)xMin;
float y = (float)yMax;
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, "%3.0f", x);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
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, "%3.0f", y);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
buffer,
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 rPenNew = NULL;
HGDIOBJ hPenOld = NULL;
rPenNew = CreatePen(PS_SOLID, 2, RGB(255, 0, 0));
hPenOld = SelectObject(hdc, rPenNew);
px = (float)points1[1].x;
py = (float)points1[1].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 < points1.size(); j++)
{
px = (float)points1[j].x;
py = (float)points1[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(rPenNew);
px = (float)points2[1].x;
py = (float)points2[1].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
HGDIOBJ gPenNew = CreatePen(PS_SOLID, 2, RGB(0, 255, 0));
hPenOld = SelectObject(hdc, gPenNew);
px = (float)points2[1].x;
py = (float)points2[1].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 < points2.size(); j++)
{
px = (float)points2[j].x;
py = (float)points2[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(gPenNew);
px = (float)points3[1].x;
py = (float)points3[1].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
HGDIOBJ bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
px = (float)points3[1].x;
py = (float)points3[1].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 < points3.size(); j++)
{
px = (float)points3[j].x;
py = (float)points3[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(bPenNew);
px = (float)points4[1].x;
py = (float)points4[1].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
HGDIOBJ tPenNew = CreatePen(PS_SOLID, 2, RGB(128, 0, 128));
hPenOld = SelectObject(hdc, tPenNew);
px = (float)points4[1].x;
py = (float)points4[1].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 < points4.size(); j++)
{
px = (float)points4[j].x;
py = (float)points4[j].y;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
SelectObject(hdc, hPenOld);
DeleteObject(tPenNew);
EndPaint(hDlg, &ps);
return (INT_PTR)FALSE;
break;
}
return (INT_PTR)FALSE;
}