Tag: C/C++
The programming language I utilized.
Blog Entry © Monday, November 17, 2025, by James Pate Williams, Jr. An Elitist Evolutionary Hill Climber to Solve the 8-Tile Puzzle
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;
}
Rice-Golomb Encoder and Decoder Copyright (c) Thursday, April 3, 2025, to Sunday, April 6, 2025, by James Pate Williams, Jr. BA, BS, Master of Software Engineering, Doctor of Philosophy Computer Science
Online references:
https://en.wikipedia.org/wiki/Golomb_coding
// Rice-Golomb Encoder and Decoder
// Copyright (c) Thursday, April 3, 2025
// by James Pate Williams, Jr.
// BA, BS, Master of Software Engineering
// Doctor of Philosophy Computer Science
// Online references:
// https://en.wikipedia.org/wiki/Golomb_coding
// https://ntrs.nasa.gov/api/citations/19790014634/downloads/19790014634.pdf
#include <iostream>
#include <string>
#include <vector>
//#include <stdlib.h>
bool Encode(const char* NChars, size_t NCharsCount,
long M, long long& N, std::vector<char>& qBits,
std::vector<char>& rBits, unsigned int& qSize, unsigned int& rSize,
long long& q, long long& r, unsigned int& NSize) {
N = NChars[0] - (long long)'0';
for (unsigned int i = 1; i < NCharsCount; i++) {
N = 10 * N + (long long)NChars[i] - (long long)'0';
}
q = N / M;
r = N % M;
qSize = 0;
while (qSize < q) {
qBits.push_back('1');
qSize++;
}
qBits.push_back('0');
qSize++;
rSize = 0;
unsigned int b = (unsigned int)floor(log2(M));
if (b > 62) {
return false;
}
long long p = (long long)pow(2, b + 1);
if (r < p - M) {
long long rr = r;
while (rr > 0) {
long long digit = (rr & 1) == 1 ? 1 : 0;
rBits.push_back((char)digit + '0');
rSize++;
rr >>= 1;
}
rBits.push_back('0');
rSize++;
}
else {
long long rr = r + p - M;
while (rSize < b + 1) {
long long digit = rr & 1 ? 1 : 0;
rBits.push_back((char)digit + '0');
rSize++;
rr >>= 1;
}
}
long long rValue = rBits[0];
for (size_t i = 1; i < rSize; i++) {
rValue = rValue * 2 + rBits[i];
}
long long NBitCount = 0;
while (N > 0) {
N >>= 1;
NBitCount++;
}
std::cout << "q-bits size = " << qSize << std::endl;
std::cout << "r-bits size = " << rSize << std::endl;
std::cout << "N-bits size = " << qSize + rSize << std::endl;
std::cout << "N-Chars * 8-Bits per Char = " << NCharsCount * 8 << std::endl;
std::cout << "% Compression = " << 100.0 * (1.0 - (qSize + rSize) /
(NCharsCount * 8.0)) << std::endl;
return true;
}
void Decode(long long M, long long& N,
std::vector<char> qBits, std::vector<char> rBits,
unsigned int& qSize, unsigned int& rSize,
long long& q, long long& r) {
int count = 0;
while (qBits[count] != '0') {
count++;
}
q = count;
int c = (int)rSize - 1;
unsigned int b = (unsigned int)floor(log2(M));
long long p = (long long)pow(2, b + 1);
long long s = 0;
r = rBits[c--] - (long long)'0';
do {
r = 2 * r + rBits[c] - (long long)'0';
c--;
} while (c >= 0);
if (r < p - M) {
s = r;
}
else {
s = r + p - M;
c = 1;
r = rBits[0] - (long long)'0';
while (c < (int)(b + 1)) {
r = 2 * r + rBits[c] - (long long)'0';
c++;
}
s = r;
}
r = s;
N = q * M + r;
}
int main() {
char line[128] = { };
size_t NSize = 0, qSize = 0, rSize = 0;
long long M = 10, N = 42, q = -1, r = -1;
std::vector<char> qBits, rBits;
std::cout << "M = ";
std::cin.getline(line, 127);
std::string str1(line);
M = std::stoi(str1);
std::cout << "N = ";
std::cin.getline(line, 127);
std::string str2(line);
Encode(str2.c_str(), strlen(str2.c_str()), M, N,
qBits, rBits, qSize, rSize, q, r, NSize);
std::cout << "q = " << q << std::endl;
std::cout << "r = " << r << std::endl;
std::cout << "q-size = " << qSize << std::endl;
std::cout << "r-size = " << rSize << std::endl;
std::cout << "q ";
for (unsigned int i = 0; i < qSize; i++) {
std::cout << qBits[i] << ' ';
}
std::cout << std::endl;
std::cout << "r ";
for (int i = (int)rSize - 1; i >= 0; i--) {
std::cout << rBits[i] << ' ';
}
std::cout << std::endl;
Decode(M, N, qBits, rBits, qSize, rSize, q, r);
std::cout << "q = " << q << std::endl;
std::cout << "r = " << r << std::endl;
std::cout << "q-size = " << qSize << std::endl;
std::cout << "r-size = " << rSize << std::endl;
std::cout << "q ";
for (unsigned int i = 0; i < qSize; i++) {
std::cout << qBits[i] << ' ';
}
std::cout << std::endl;
std::cout << "r ";
for (int i = rSize - 1; i >= 0; i--) {
std::cout << rBits[i] << ' ';
}
std::cout << std::endl;
std::cout << "N = " << N << std::endl;
return 0;
}
M = 64
N = 1027
q-bits size = 17
r-bits size = 3
N-bits size = 20
N-Chars * 8-Bits per Char = 32
% Compression = 37.5
q = 16
r = 3
q-size = 17
r-size = 3
q 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
r 0 1 1
q = 16
r = 3
q-size = 17
r-size = 3
q 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
r 0 1 1
N = 1027
Chapter One Straight-Line Program Interpreter from “Modern Compiler Implementation in Java Second Edition” (c) 2002 by Andrew W. Appel, Translation to C++ by James Pate Williams, Jr. on Thursday, April 3, 2025
Wayback in the Spring Semester of 2006, after I was awarded my Doctor of Philosophy Degree in Computer Science, I partially audited a Compiler Design Course. Due to my mental aberrations, I was unable to complete the course. The instructor was on a Sabbatical from the United States Air Force Academy in Colorado Springs, Colorado. The textbook we used, and I still have a copy, was “Modern Compiler Implementation in Java Second Edition” © 2002 by Andrew W. Appel. Below is a translation from Java to C++ that I just completed.
// Chapter One program translated from Java to C++ by
// James Pate Williams, Jr. (c) Wednesday April 3, 2025
// Reference: "Modern Complier Implementation in Java
// Second Edition" (c) 2002 by Andrew W. Appel
#ifndef _SLPInterpreter_H
#include <iostream>
#include <stack>
#include <string>
#include <vector>
class TableEntry {
public:
std::string symbol, value;
TableEntry(std::string symbol, std::string value) {
this->symbol = symbol;
this->value = value;
}
};
std::stack<std::string> sStack;
std::vector<TableEntry> symbolTable;
class Exp {
public:
Exp() { };
virtual ~Exp() { };
};
std::stack<Exp> eStack;
class ExpList {
public:
ExpList() { };
virtual ~ExpList() { };
};
class Stm {
public:
Stm() { };
virtual ~Stm() { };
};
class CompoundStm : public Stm {
public:
Stm stm1, stm2;
CompoundStm(Stm stm1, Stm stm2) {
this->stm1 = stm1;
this->stm2 = stm2;
};
};
class AssignStm : public Stm {
public:
std::string id;
Exp exp;
AssignStm(std::string id, Exp exp) {
this->id = id;
this->exp = exp;
bool found = false;
for (int i = 0; !found && i < (int)symbolTable.size(); i++) {
if (symbolTable[i].symbol == id) {
found = true;
}
}
if (!found) {
symbolTable.push_back(TableEntry(id, ""));
}
};
void Print() {
std::cout << this->id << ' ';
};
};
class PrintStm : public Stm {
public:
ExpList exps;
PrintStm(ExpList exps) {
this->exps = exps;
};
};
class IdExp : public Exp {
public:
std::string id;
IdExp(std::string id) {
this->id = id;
Print();
TableEntry te(id, "");
};
void Print() {
std::cout << id << ' ';
};
};
class NumExp : public Exp {
public:
int num;
NumExp(int num) {
this->num = num;
Print();
char buffer[128] = { };
_itoa_s(num, buffer, 127, 10);
sStack.push(std::string(buffer));
};
void Print() {
std::cout << num << ' ';
};
};
enum class ArithmeticOp {
Plus, Minus, Times, Div
};
class OpExp : public Exp {
public:
Exp left, right;
ArithmeticOp op;
OpExp(Exp left, ArithmeticOp op, Exp right) {
this->left = left;
this->op = op;
this->right = right;
std::string ops = "";
switch (op) {
case ArithmeticOp::Plus:
ops = "+";
break;
case ArithmeticOp::Minus:
ops = "-";
break;
case ArithmeticOp::Times:
ops = "*";
break;
case ArithmeticOp::Div:
ops = "/";
break;
};
std::cout << ops << std::endl;
eStack.push(left);
eStack.push(right);
sStack.push(ops);
};
};
class EseqExp : public Exp {
public:
Stm stm; Exp exp;
EseqExp(Stm stm, Exp exp) {
this->stm = stm;
this->exp = exp;
};
};
class PairExpList : public ExpList {
public:
Exp head;
ExpList tail;
PairExpList(Exp head, ExpList tail) {
this->head = head;
this->tail = tail;
};
};
class LastExpList : public ExpList {
public:
Exp head;
LastExpList(Exp head) {
this->head = head;
};
};
#endif _SLPInterpreter_H
int main() {
int a = 0, b = 0;
Stm prog(CompoundStm(AssignStm("a",
OpExp(NumExp(5), ArithmeticOp::Plus, NumExp(3))),
CompoundStm(AssignStm("b",
EseqExp(PrintStm(PairExpList(IdExp("a"),
LastExpList(OpExp(IdExp("a"),
ArithmeticOp::Minus, NumExp(1))))),
OpExp(NumExp(10), ArithmeticOp::Times, IdExp("a")))),
PrintStm(LastExpList(IdExp("b"))))));
bool first = true;
int result = 0;
//sStack.push("0");
while (!sStack.empty()) {
std::string lts, ops, rts;
if (first) {
ops = sStack.top();
sStack.pop();
lts = sStack.top();
sStack.pop();
rts = sStack.top();
sStack.pop();
first = false;
}
else {
lts = sStack.top();
sStack.pop();
ops = sStack.top();
sStack.pop();
rts = sStack.top();
sStack.pop();
}
int lvi = std::stoi(lts);
int rvi = std::stoi(rts);
if (ops == "+") {
result = lvi + rvi;
}
else if (ops == "-") {
result = lvi - rvi;
}
else if (ops == "*") {
result = lvi * rvi;
}
else if (ops == "/") {
result = lvi / rvi;
}
char ascii[128] = { };
_itoa_s(result, ascii, 10);
if (sStack.size() != 0) {
sStack.push(std::string(ascii));
}
}
std::cout << "Result = " << result << std::endl;
return 0;
}
Blog Entry © Sunday, March 29, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Slater Determinant Coefficients for Z = 2 to 4
Enter the atomic number Z (2 to 6 or 0 to quit): 2
2 1 1 + a(1)b(2)
1 0 0 - a(2)b(1)
# Even Permutations = 1
Enter the atomic number Z (2 to 6 or 0 to quit): 3
6 3 1 + a(1)b(2)c(3)
5 2 0 - a(1)b(3)c(2)
4 2 0 - a(2)b(1)c(3)
3 1 1 + a(2)b(3)c(1)
2 1 1 + a(3)b(1)c(2)
1 0 0 - a(3)b(2)c(1)
# Even Permutations = 3
Enter the atomic number Z (2 to 6 or 0 to quit): 4
24 12 0 + a(1)b(2)c(3)d(4)
23 11 1 - a(1)b(2)c(4)d(3)
22 11 1 - a(1)b(3)c(2)d(4)
21 10 0 + a(1)b(3)c(4)d(2)
20 10 0 + a(1)b(4)c(2)d(3)
19 9 1 - a(1)b(4)c(3)d(2)
18 9 1 - a(2)b(1)c(3)d(4)
17 8 0 + a(2)b(1)c(4)d(3)
16 8 0 + a(2)b(3)c(1)d(4)
15 7 1 - a(2)b(3)c(4)d(1)
14 7 1 - a(2)b(4)c(1)d(3)
13 6 0 + a(2)b(4)c(3)d(1)
12 6 0 + a(3)b(1)c(2)d(4)
11 5 1 - a(3)b(1)c(4)d(2)
10 5 1 - a(3)b(2)c(1)d(4)
9 4 0 + a(3)b(2)c(4)d(1)
8 4 0 + a(3)b(4)c(1)d(2)
7 3 1 - a(3)b(4)c(2)d(1)
6 3 1 - a(4)b(1)c(2)d(3)
5 2 0 + a(4)b(1)c(3)d(2)
4 2 0 + a(4)b(2)c(1)d(3)
3 1 1 - a(4)b(2)c(3)d(1)
2 1 1 - a(4)b(3)c(1)d(2)
1 0 0 + a(4)b(3)c(2)d(1)
# Even Permutations = 12
Enter the atomic number Z (2 to 6 or 0 to quit):
// AOPermutations.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Saturday, March 29, 2025
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Signs of the atomic orbitals in a Slater Determinant
#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
int main()
{
char alpha[] = { 'a', 'b', 'c', 'd', 'e', 'f' }, line[128] = {};
int factorial[7] = { 1, 1, 2, 6, 24, 120, 720 };
while (true)
{
int col = 0, counter = 0, row = 0, sign = 1, t = 0, Z = 0, zfact = 0;
int numberEven = 0;
std::cout << "Enter the atomic number Z (2 to 6 or 0 to quit): ";
std::cin.getline(line, 127);
std::string str(line);
Z = std::stoi(str);
if (Z == 0)
{
break;
}
if (Z < 2 || Z > 6)
{
std::cout << "Illegal Z, please try again" << std::endl;
continue;
}
zfact = factorial[Z];
std::vector<char> orb(Z);
std::vector<int> tmp(Z), vec(Z);
for (int i = 0; i < Z; i++)
{
orb[i] = alpha[i];
vec[i] = i + 1;
}
do
{
for (int i = 0; i < (int)vec.size(); i++)
{
tmp[i] = vec[i];
}
t = 0;
do
{
t++;
} while (std::next_permutation(tmp.begin(), tmp.end()));
std::cout << t << '\t' << t / 2 << '\t';
std::cout << (t / 2 & 1) << '\t';
if (Z == 2 || Z == 3)
{
if ((t / 2 & 1) == 0)
{
std::cout << "-\t";
}
else
{
std::cout << "+\t";
numberEven++;
}
}
else
{
if ((t / 2 & 1) == 1)
{
std::cout << "-\t";
}
else
{
std::cout << "+\t";
numberEven++;
}
}
for (int i = 0; i < Z; i++)
{
std::cout << orb[i] << '(' << vec[i] << ')';
}
row++;
std::cout << std::endl;
if (zfact != 2 && row == zfact)
{
std::cout << std::endl;
break;
}
row %= Z;
} while (std::next_permutation(vec.begin(), vec.end()));
std::cout << "# Even Permutations = ";
std::cout << numberEven << std::endl;
}
return 0;
}
Blog Entry (c) Wednesday, November 6, 2024, by James Pate Williams, Jr. Small Angular Momentum Quantum Numbers Gaunt Coefficients
// GauntCoefficients.cpp (c) Monday, November 4, 2024
// by James Pate Williams, Jr., BA, BS, MSWE, PhD
// Computes the Gaunt angular momentum coefficients
// Also the Wigner-3j symbols are calculated
// https://en.wikipedia.org/wiki/3-j_symbol
// https://doc.sagemath.org/html/en/reference/functions/sage/functions/wigner.html#
// https://www.geeksforgeeks.org/factorial-large-number/
#include <iostream>
using namespace std;
typedef long double real;
real pi;
// iterative n-factorial function
real Factorial(int n)
{
real factorial = 1;
for (int i = 2; i <= n; i++)
factorial *= i;
if (n < 0)
factorial = 0;
return factorial;
}
real Delta(int lt, int rt)
{
return lt == rt ? 1.0 : 0.0;
}
real Wigner3j(
int j1, int j2, int j3,
int m1, int m2, int m3)
{
real delta = Delta(m1 + m2 + m3, 0) *
powl(-1.0, j1 - j2 - m3);
real fact1 = Factorial(j1 + j2 - j3);
real fact2 = Factorial(j1 - j2 + j3);
real fact3 = Factorial(-j1 + j2 + j3);
real denom = Factorial(j1 + j2 + j3 + 1);
real numer = delta * sqrt(
fact1 * fact2 * fact3 / denom);
real fact4 = Factorial(j1 - m1);
real fact5 = Factorial(j1 + m1);
real fact6 = Factorial(j2 - m2);
real fact7 = Factorial(j2 + m2);
real fact8 = Factorial(j3 - m3);
real fact9 = Factorial(j3 + m3);
real sqrt1 = sqrtl(
fact4 * fact5 * fact6 * fact7 * fact8 * fact9);
real sumK = 0;
int K = (int)fmaxl(0, fmaxl((real)j2 - j3 - m1,
(real)j1 - j3 + m2));
int N = (int)fminl((real)j1 + j2 - j3,
fminl((real)j1 - m1, (real)j2 + m2));
for (int k = K; k <= N; k++)
{
real f0 = Factorial(k);
real f1 = Factorial(j1 + j2 - j3 - k);
real f2 = Factorial(j1 - m1 - k);
real f3 = Factorial(j2 + m2 - k);
real f4 = Factorial(j3 - j2 + m1 + k);
real f5 = Factorial(j3 - j1 - m2 + k);
sumK += powl(-1.0, k) / (f0 * f1 * f2 * f3 * f4 * f5);
}
return numer * sqrt1 * sumK;
}
real GauntCoefficient(
int l1, int l2, int l3, int m1, int m2, int m3)
{
real factor = sqrtl(
(2.0 * l1 + 1.0) *
(2.0 * l2 + 1.0) *
(2.0 * l3 + 1.0) /
(4.0 * pi));
real wigner1 = Wigner3j(l1, l2, l3, 0, 0, 0);
real wigner2 = Wigner3j(l1, l2, l3, m1, m2, m3);
return factor * wigner1 * wigner2;
}
int main()
{
pi = 4.0 * atanl(1.0);
cout << "Gaunt(1, 0, 1, 1, 0, 0) = ";
cout << GauntCoefficient(1, 0, 1, 1, 0, 0);
cout << endl;
cout << "Gaunt(1, 0, 1, 1, 0, -1) = ";
cout << GauntCoefficient(1, 0, 1, 1, 0, -1);
cout << endl;
real number = -1.0 / 2.0 / sqrtl(pi);
cout << "-1.0 / 2.0 / sqrt(pi) = ";
cout << number << endl;
return 0;
}
Blog Entry (c) Tuesday, July 23, 2024, by James Pate Williams, Jr. Mueller’s Method for Finding the Complex and/or Real Roots of a Complex and/or Real Polynomial
I originally implemented this algorithm in FORTRAN IV in the Summer Quarter of 1982 at the Georgia Institute of Technology. I was taking a course named “Scientific Computing I” taught by Professor Gunter Meyer. I made a B in the class. Later in 2015 I re-implemented the recipe in C# using Visual Studio 2008 Professional. VS 2015 did not have support for complex numbers nor large integers. In December of 2015 I upgraded to Visual Studio 2015 Professional which has support for big integers and complex numbers. I used Visual Studio 2019 Community version for this project. Root below should be function.
Degree (0 to quit) = 2
coefficient[2].real = 1
coefficient[2].imag = 0
coefficient[1].real = 1
coefficient[1].imag = 0
coefficient[0].real = 1
coefficient[0].imag = 0
zero[0].real = -5.0000000000e-01 zero[0].imag = 8.6602540378e-01
zero[1].real = -5.0000000000e-01 zero[1].imag = -8.6602540378e-01
root[0].real = 0.0000000000e+00 root[0].imag = -2.2204460493e-16
root[1].real = 3.3306690739e-16 root[1].imag = -7.7715611724e-16
Degree (0 to quit) = 3
coefficient[3].real = 1
coefficient[3].imag = 0
coefficient[2].real = 0
coefficient[2].imag = 0
coefficient[1].real = -18.1
coefficient[1].imag = 0
coefficient[0].real = -34.8
coefficient[0].imag = 0
zero[0].real = -2.5026325486e+00 zero[0].imag = -8.3036679880e-01
zero[1].real = -2.5026325486e+00 zero[1].imag = 8.3036679880e-01
zero[2].real = 5.0052650973e+00 zero[2].imag = 2.7417672687e-15
root[0].real = 0.0000000000e+00 root[0].imag = 1.7763568394e-15
root[1].real = 3.5527136788e-14 root[1].imag = -1.7763568394e-14
root[2].real = 2.8421709430e-14 root[2].imag = 1.5643985575e-13
Degree (0 to quit) = 5
coefficient[5].real = 1
coefficient[5].imag = 0
coefficient[4].real = 2
coefficient[4].imag = 0
coefficient[3].real = 3
coefficient[3].imag = 0
coefficient[2].real = 4
coefficient[2].imag = 0
coefficient[1].real = 5
coefficient[1].imag = 0
coefficient[0].real = 6
coefficient[0].imag = 0
zero[0].real = -8.0578646939e-01 zero[0].imag = 1.2229047134e+00
zero[1].real = -8.0578646939e-01 zero[1].imag = -1.2229047134e+00
zero[2].real = 5.5168546346e-01 zero[2].imag = 1.2533488603e+00
zero[3].real = 5.5168546346e-01 zero[3].imag = -1.2533488603e+00
zero[4].real = -1.4917979881e+00 zero[4].imag = 1.8329656063e-15
root[0].real = 8.8817841970e-16 root[0].imag = 4.4408920985e-16
root[1].real = -2.6645352591e-15 root[1].imag = -4.4408920985e-16
root[2].real = 8.8817841970e-16 root[2].imag = 1.7763568394e-15
root[3].real = 3.4638958368e-14 root[3].imag = -1.4210854715e-14
root[4].real = 8.8817841970e-16 root[4].imag = 2.0710031449e-14
Blog Entry Friday, July 19, 2024, Easy Internet Math “Puzzle” (c) James Pate Williams, Jr.
#include <math.h>
#include <iostream>
using namespace std;
long double f(long double x)
{
return powl(8.0, x) - powl(2.0, x) -
2.0 * (powl(6.0, x) - powl(3.0, x));
}
long double g(long double x)
{
return powl(8.0, x) * logl(8.0) - powl(2.0, x) * logl(2.0) -
2.0 * (powl(6.0, x) * logl(6.0) - powl(3.0, x) * logl(3.0));
}
long double Newton(long double x, int maxIts, int& iterations)
{
long double x0 = x;
long double x1 = 0.0;
iterations = 0;
while (true) {
long double dx = 0.0;
long double fx = f(x0);
long double gx = g(x0);
x1 = x0 - fx / gx;
dx = fabsl(x1 - x0);
iterations++;
if (dx < 1.0e-15)
break;
if (fabsl(fx) < 1.0e-15)
break;
if (iterations == maxIts)
break;
x0 = x1;
}
return x1;
}
int main() {
int iterations = 0, maxIts;
long double x0 = 0.0, x1 = 0.0;
while (true) {
cout << "x0 = ";
cin >> x0;
if (x0 == 0)
break;
cout << "maximum iterations = ";
cin >> maxIts;
x1 = Newton(x0, maxIts, iterations);
cout << "x1 = " << x1 << endl;
cout << "iterations = ";
cout << iterations << endl;
}
return 0;
}
Blog Entry Monday, June 24, 2024 (c) James Pate Williams, Jr. Computing Binomial Coefficients and Pascal’s Triangle in the C Language
Enter n (<= 18) below:
5
Enter k (<= 18) below:
0
1 1
Enter n (<= 18) below:
5
Enter k (<= 18) below:
1
5 5
Enter n (<= 18) below:
5
Enter k (<= 18) below:
2
10 10
Enter n (<= 18) below:
0
Enter n (<= 18) below:
0
Pascal's Triangle:
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
1 7 21 35 35 21 7 1
1 8 28 56 70 56 28 8 1
1 9 36 84 126 126 84 36 9 1
1 10 45 120 210 252 210 120 45 10 1
1 11 55 165 330 462 462 330 165 55 11 1
1 12 66 220 495 792 924 792 495 220 66 12 1
1 13 78 286 715 1287 1716 1716 1287 715 286 78 13 1
1 14 91 364 1001 2002 3003 3432 3003 2002 1001 364 91 14 1
1 15 105 455 1365 3003 5005 6435 6435 5005 3003 1365 455 105 15 1
1 16 120 560 1820 4368 8008 11440 12870 11440 8008 4368 1820 560 120 16 1
1 17 136 680 2380 6188 12376 19448 24310 24310 19448 12376 6188 2380 680 136 17 1
1 18 153 816 3060 8568 18564 31824 43758 48620 43758 31824 18564 8568 3060 816 153 18 1
C:\Users\james\source\repos\BinomialCoefficeint\Debug\BinomialCoefficeint.exe (process 40028) exited with code 0.
Press any key to close this window . . .
// BinomialCoefficient.c (c) Monday, June 24, 2024
// by James Pate Williams, Jr. BA, BS, MSwE, PhD
#include <stdio.h>
#include <stdlib.h>
typedef long long ll;
ll** Binomial(ll n)
{
ll** C = (ll**)calloc(n + 1, sizeof(ll*));
if (C == NULL)
exit(-1);
for (int i = 0; i < n + 1; i++)
{
C[i] = (ll*)calloc(n + 1, sizeof(ll));
if (C[i] == NULL)
exit(-1);
}
if (n >= 0)
{
C[0][0] = 1;
}
if (n >= 1)
{
C[1][0] = 1;
C[1][1] = 1;
}
if (n >= 2)
{
for (int i = 2; i <= n; i++)
{
for (int j = 2; j <= n; j++)
{
C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
}
}
}
return C;
}
ll Factorial(ll n)
{
ll fact = 1;
if (n > 1)
{
for (int i = 2; i <= n; i++)
fact = i * fact;
}
return fact;
}
ll BC(ll n, ll k)
{
return Factorial(n) / (Factorial(n - k) * Factorial(k));
}
int main()
{
int i = 0, j = 0;
ll** C = Binomial(20);
while (1)
{
char buffer[256] = { '\0' };
printf_s("Enter n (<= 18) below:\n");
scanf_s("%s", buffer, sizeof(buffer));
printf_s("\n");
ll n = atoll(buffer);
if (n == 0)
break;
printf_s("Enter k (<= 18) below:\n");
scanf_s("%s", buffer, sizeof(buffer));
printf_s("\n");
ll k = atoll(buffer);
printf_s("%lld\t%lld\n\n", C[n + 2][k + 2], BC(n, k));
}
printf_s("Pascal's Triangle:\n\n");
for (i = 2; i <= 20; i++)
{
for (j = 2; j <= 20; j++)
if (C[i][j] != 0)
printf_s("%5lld ", C[i][j]);
printf_s("\n");
}
for (i = 0; i <= 20; i++)
free(C[i]);
free(C);
}
Numerical Explorations 