Category: C++ Computer Applications

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

Blog Entry Saturday May 30 2026Download
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

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

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

    negD = (~D + 1)

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

    T = R + negD

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

    return (Q, R) 
 */

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


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

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

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

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

    quo = 0;
    rem = 0;

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

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

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

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

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

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

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

    uint32_t R = 0;

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

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

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

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

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

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

    product = 0;

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

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

    return true;
}


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

namespace {

    constexpr int TESTS_PER_N = 200000;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

} // namespace

int main() {
    bool verbose = true;

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

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

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

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

Blog Entry Saturday May 16 2026Download

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

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

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

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

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

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

        sum += f(x);
    }

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

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

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

    return factorial;
}

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

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

        sign *= -1;
    }

    sign = 1;

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

        sign *= -1;
    }

    return sumB - sumA;
}

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

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

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

static double SimpsonsRule(
    int n, double a, double b, double(*fx)(double)) {
    double h = (b - a) / n;
    double h2 = 2.0 * h;
    double s = 0.0;
    double t = 0.0;
    double x = a + h;

    for (int i = 1; i < n; i += 2) {
        s += fx(x);
        x += h2;
    }

    x = a + h2;

    for (int i = 2; i < n; i += 2) {
        t += fx(x);
        x += h2;
    }

    return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}

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

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

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

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

    if (nRow < 2)
        return;

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

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

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

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

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

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

    if (nRow < 3) {
        return;
    }

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

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

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

            T[k][j] = ratio;
        }

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

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

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

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

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

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

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

Blog Entry Wednesday May 13 2026Download

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

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

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

Blog Entry Tuesday May 5 2026Download

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

Butterfly curve (transcendental) – Wikipedia

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

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

#define MAX_LOADSTRING 100

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

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
std::vector<Point3d> points;

// Forward declarations of functions included in this code module:
ATOM                MyRegisterClass(HINSTANCE hInstance);
BOOL                InitInstance(HINSTANCE, int);
LRESULT CALLBACK    WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    About(HWND, UINT, WPARAM, LPARAM);

int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    UNREFERENCED_PARAMETER(hPrevInstance);
    UNREFERENCED_PARAMETER(lpCmdLine);

    // TODO: Place code here.

    // Initialize global strings
    LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
    LoadStringW(hInstance, IDC_BUTTERFLYEQUATION, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, nullptr, 0, 0))
    {
        if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
        {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }

    return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
    WNDCLASSEXW wcex = { 0 };

    wcex.cbSize = sizeof(WNDCLASSEX);

    wcex.style          = CS_HREDRAW | CS_VREDRAW;
    wcex.lpfnWndProc    = WndProc;
    wcex.cbClsExtra     = 0;
    wcex.cbWndExtra     = 0;
    wcex.hInstance      = hInstance;
    wcex.hIcon          = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_BUTTERFLYEQUATION));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_BUTTERFLYEQUATION);
    wcex.lpszClassName  = szWindowClass;
    wcex.hIconSm        = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

    return RegisterClassExW(&wcex);
}

//
//   FUNCTION: InitInstance(HINSTANCE, int)
//
//   PURPOSE: Saves instance handle and creates main window
//
//   COMMENTS:
//
//        In this function, we save the instance handle in a global variable and
//        create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   hInst = hInstance; // Store instance handle in our global variable

   HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

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

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

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

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

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

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

static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
    // Draw the text with formatting options
    DrawTextA(hdc, text, -1, &rect, DT_SINGLELINE | DT_NOCLIP);
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    switch (message)
    {
    case WM_CREATE:
        CreateButterflyGraphPoints();
        break;
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);
            // Parse the menu selections:
            switch (wmId)
            {
            case IDM_ABOUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
                break;
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
        {
            PAINTSTRUCT ps;
            HDC hdc = BeginPaint(hWnd, &ps);
            double h = 0, pi = 0, plm = 0, theta = 0;
            double xMax = 0, xMin = 0, yMax = 0, yMin = 0;
            FindMinMax(xMin, xMax, yMin, yMax);
            float xSpan = (float)(xMax - xMin);
            float ySpan = (float)(yMax - yMin);
            RECT rect = { };
            GetClientRect(hWnd, &rect);
            float width = (float)(rect.right - rect.left + 1);
            float height = (float)(rect.bottom - rect.top - 32 + 1);
            float sx0 = 2.0f * width / 16.0f;
            float sx1 = 14.0f * width / 16.0f;
            float sy0 = 2.0f * height / 16.0f;
            float sy1 = 14.0f * height / 16.0f;
            float deltaX = xSpan / 8.0f;
            float deltaY = ySpan / 8.0f;
            float xSlope = (sx1 - sx0) / xSpan;
            float xInter = (float)(sx0 - xSlope * xMin);
            float ySlope = (sy0 - sy1) / ySpan;
            float yInter = (float)(sy0 - ySlope * yMax);
            float px = 0, py = 0, sx = 0, sy = 0;
            POINT wPt = { };
            int i = 0;
            float x = (float)xMin;
            float y = (float)yMax;
            px = x;
            py = y;
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
            char buffer[128] = { };

            while (i <= 8)
            {
                sx = xSlope * x + xInter;
                wPt.x = wPt.y = 0;
                MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
                LineTo(hdc, (int)sx, (int)sy1);

                sprintf_s(buffer, "%5.4lf", x);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    buffer,
                    (int)strlen(buffer),
                    &size);
                RECT textRect = { };
                textRect.left = (long)(sx - size.cx / 2.0f);
                textRect.right = (long)(sx + size.cx / 2.0f);
                textRect.top = (long)sy1;
                textRect.bottom = (long)(sy1 + size.cy / 2.0f);
                DrawFormattedText(hdc, buffer, textRect);
                x += deltaX;
                i++;
            }

            i = 0;
            y = (float)yMin;

            while (i <= 8)
            {
                sy = ySlope * y + yInter;
                wPt.x = wPt.y = 0;
                MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
                LineTo(hdc, (int)sx, (int)sy);

                if (i != 0)
                {
                    sprintf_s(buffer, "%+5.3lf", y);
                    SIZE size = { };
                    GetTextExtentPoint32A(
                        hdc,
                        buffer,
                        (int)strlen(buffer),
                        &size);
                    RECT textRect = { };
                    textRect.left = (long)(sx0 - size.cx - size.cx / 5.0f);
                    textRect.right = (long)(sx0 - size.cx / 2.0f);
                    textRect.top = (long)(sy - size.cy / 2.0f);
                    textRect.bottom = (long)(sy + size.cy / 2.0f);
                    DrawFormattedText(hdc, buffer, textRect);
                }

                y += deltaY;
                i++;
            }

            HGDIOBJ bPenNew = NULL;
            HGDIOBJ hPenOld = NULL;

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

            px = (float)points[0].x;
            py = (float)points[0].y;
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx, (int)sy, &wPt);

            for (size_t j = 1; j < points.size(); j++)
            {
                px = (float)points[j].x;
                py = (float)points[j].y;
                sx = xSlope * px + xInter;
                sy = ySlope * py + yInter;
                LineTo(hdc, (int)sx, (int)sy);
            }

            SelectObject(hdc, hPenOld);
            DeleteObject(bPenNew);
            EndPaint(hWnd, &ps);
        }
        break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

Blog Entry © Sunday, April 19, 2026, by James Pate Williams, Jr., Scattering from a Spherically Symmetric Potential

Blog Entry Sunday April 19 2026Download

Vector Analysis by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition© 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Pages 48 and 50

Vector Analysis Chapter 1 by James Pate WilliamsDownload
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider 

#include <iostream>
#include <vector>

static double InnerProduct(
	std::vector<double> A,
	std::vector<double> B,
	int n)
{
	double sum = 0.0;

	for (int i = 0; i < n; i++)
		sum += A[i] * B[i];

	return sum;
}

static void VectorProduct(
	std::vector<double> A,
	std::vector<double> B,
	std::vector<double>&C)
{
	C.resize(3);
	C[0] = A[1] * B[2] - A[2] * B[1];
	C[1] = A[2] * B[0] - A[0] * B[2];
	C[2] = A[0] * B[1] - A[1] * B[0];
}

static double TripleProduct(
	std::vector<double> A,
	std::vector<double> B,
	std::vector<double> C)
{
	double sum0 = 0.0, sum1 = 0.0;

	sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
	sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
	return sum0 - sum1;
}

static void Exercises_Section_1_13_1_Triple_Products()
{
	// triple products (a)
	std::vector<double> A1 = { 2, 0, 0 };
	std::vector<double> B1 = { 0, 3, 0 };
	std::vector<double> C1 = { 0, 0, 5 };
	double tp_1 = TripleProduct(A1, B1, C1);
	std::cout << "Exercise 1 (a)" << '\t';
	std::cout << tp_1 << std::endl;
	// triple products (b)
	std::vector<double> A2 = { 1, 1, 1 };
	std::vector<double> B2 = { 3, 1, 0 };
	std::vector<double> C2 = { 0, -1, 5 };
	std::cout << "Exercise 1 (b)" << '\t';
	double tp_2 = TripleProduct(A2, B2, C2);
	std::cout << tp_2 << std::endl;
	// triple products (c)
	std::vector<double> A3 = { 2, -1, 1 };
	std::vector<double> B3 = { 1, 1, 1 };
	std::vector<double> C3 = { 2, 0, 3 };
	double tp_3 = TripleProduct(A3, B3, C3);
	std::cout << "Exercise 1 (c)" << '\t';
	std::cout << tp_3 << std::endl;
	// triple products (d)
	std::vector<double> A4 = { 0, 0, 1 };
	std::vector<double> B4 = { 1, 0, 0 };
	std::vector<double> C4 = { 0, 1, 0 };
	double tp_4 = TripleProduct(A4, B4, C4);
	std::cout << "Exercise 1 (d)" << '\t';
	std::cout << tp_4 << std::endl;
	// volume of a parallelpipped
	std::vector<double> A5 = { 3, 4, 1 };
	std::vector<double> B5 = { 2, 3, 4 };
	std::vector<double> C5 = { 0, 0, 5 }; 
	double tp_5 = TripleProduct(A5, B5, C5);
	std::cout << "Exercise 2" << '\t';
	std::cout << tp_5 << std::endl;
	// volume of a parallelpipped
	std::vector<double> A6 = { 3, 2, 1 };
	std::vector<double> B6 = { 4, 2, 1 };
	std::vector<double> C6 = { 0, 1, 4 };
	std::vector<double> D6 = { 0, 0, 7 };
	std::vector<double> AB = { -1, 0, 0 };
	std::vector<double> AC = { 3, 1, -3 };
	std::vector<double> AD = { 3, 2, -6 };
	std::cout << "Exercise 3" << '\t';
	double tp_6 = TripleProduct(AB, AC, AD);
	std::cout << tp_6 << std::endl;
	// volume of a tetrahedron
	std::vector<double> AB1 = { 1, 1, 0 };
	std::vector<double> AC1 = { 1, -1, 0 };
	std::vector<double> AD1 = { 0, 0, 2 };
	std::cout << "Exercise 4" << '\t';
	double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
	std::cout << fabs(tp_7) << std::endl;
	std::vector<double> P1 = { 0, 0, 0 };
	std::vector<double> P2 = { 1, 1, 0 };
	std::vector<double> P3 = { 3, 4, 0 };
	std::vector<double> P4 = { 4, 5, 0 };
	std::vector<double> P5 = { 0, 0, 1 };
	std::vector<double> Q1 = { 1, 1, 0 };
	std::vector<double> Q2 = { 2, 3, 0 };
	std::vector<double> Q3 = { 1, 1, 1 };
	std::cout << "Exercise 5" << '\t';
	double tp_8 = TripleProduct(Q1, Q2, Q3);
	std::cout << fabs(tp_8) << std::endl;
	std::vector<double> A10 = { 1, 1, 1 };
	std::vector<double> B10 = { 2, 4, -1 };
	std::vector<double> C10 = { 1, 1, 3 };
	double tp_10 = TripleProduct(A10, B10, C10);
	std::vector<double> D10 = { 0 };
	VectorProduct(A10, B10, D10);
	double magnitude = sqrt(InnerProduct(D10, D10, 3));
	std::cout << "Exercise 10" << '\t';
	std::cout << tp_10 / magnitude << '\t';
	std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
	std::vector<double> A11 = { 1, 1, 1 };
	std::vector<double> B11 = { 2, 4, -1 };
	std::vector<double> C11 = { 1, 1, 3 };
	std::vector<double> D11 = { 3, 2, 1 };
	std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
	VectorProduct(A11, B11, AB11);
	VectorProduct(B11, C11, BC11);
	VectorProduct(C11, A11, CA11);
	VectorProduct(BC11, CA11, BCCA11);
	double Q11 = InnerProduct(AB11, BCCA11, 3);
	double A_x = A11[0], A_y = A11[1], A_z = A11[2];
	double B_x = B11[0], B_y = B11[1], B_z = B11[2];
	double C_x = C11[0], C_y = C11[1], C_z = C11[2];
	double term1 = +(A_y * B_z - A_z * B_y) * (B_z * C_x - B_x * C_z) * (C_x * A_y - C_y * A_x);
	double term2 = -(A_y * B_z - A_z * B_y) * (B_x * C_y - B_y * C_x) * (C_z * A_x - C_x * A_z);
	double term3 = +(A_z * B_x - A_x * B_z) * (B_x * C_y - B_y * C_x) * (C_y * A_z - C_z * A_y);
	double term4 = -(A_z * B_x - A_x * B_z) * (B_y * C_z - B_z * C_y) * (C_x * A_y - C_y * A_x);
	double term5 = +(A_x * B_y - A_y * B_x) * (B_y * C_z - B_z * C_y) * (C_z * A_x - C_x * A_z);
	double term6 = -(A_x * B_y - A_y * B_x) * (B_z * C_x - B_x * C_z) * (C_y * A_z - C_z * A_y);
	double P11 = term1 + term2 + term3 + term4 + term5 + term6;
	std::cout << "Q = (A x B) . (B x C) x (C x A) = " << Q11 << std::endl;
	std::cout << "P = (A x B) . (B x C) x (C x A) = " << P11 << std::endl;
}

static void Exercises_Section_1_14_Vector_Identities()
{
	std::vector<double> A11 = { 1, 1, 1 };
	std::vector<double> B11 = { 2, 4, -1 };
	std::vector<double> C11 = { 1, 1, 3 };
	std::vector<double> D11 = { 3, 2, 1 };
	std::cout << "Section 1.14 page 50 Exercises Exercise 1" << std::endl;
	std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
	std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
	std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
	std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
	std::cout << "TPI1 = (A x B) x (C x D) = [A, B, D]C - [A, B, C]D = " << std::endl;
	double TP1411a = TripleProduct(A11, B11, D11);
	double TP1411b = TripleProduct(A11, B11, C11);
	std::cout << "[A, B, D] = " << TP1411a << std::endl;
	std::cout << "[A, B, C] = " << TP1411b << std::endl;
	std::cout << "TPI1_x = " << TP1411a * C11[0] << std::endl;
	std::cout << "TPI1_y = " << TP1411a * C11[1] << std::endl;
	std::cout << "TPI1_z = " << TP1411a * C11[2] << std::endl;
	std::cout << "TPI2_x = " << TP1411b * D11[0] << std::endl;
	std::cout << "TPI2_y = " << TP1411b * D11[1] << std::endl;
	std::cout << "TPI2_z = " << TP1411b * D11[2] << std::endl;
	std::cout << "RHS1 = [A, B, D]C - [A, B, C]D = " << std::endl;
	std::vector<double> RHS1(3);
	RHS1[0] = TP1411a * C11[0] - TP1411b * D11[0];
	RHS1[1] = TP1411a * C11[1] - TP1411b * D11[1];
	RHS1[2] = TP1411a * C11[2] - TP1411b * D11[2];
	std::cout << "RHS1_x = " << RHS1[0] << std::endl;
	std::cout << "RHS1_y = " << RHS1[1] << std::endl;
	std::cout << "RHS1_z = " << RHS1[2] << std::endl;
	std::vector<double> CD11(3), TD11(3);
	std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
	VectorProduct(A11, B11, AB11);
	VectorProduct(B11, C11, BC11);
	VectorProduct(C11, A11, CA11);
	VectorProduct(BC11, CA11, BCCA11);
	VectorProduct(A11, B11, AB11);
	VectorProduct(C11, D11, CD11);
	VectorProduct(AB11, CD11, TD11);
	std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
	std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
	std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
	std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
	std::cout << "A x B = " << AB11[0] << '\t' << AB11[1] << '\t' << AB11[2] << std::endl;
	std::cout << "C x D = " << CD11[0] << '\t' << CD11[1] << '\t' << CD11[2] << std::endl;
	std::cout << "TD11 = (A x B) x (C x D) = " << std::endl;
	std::cout << "TD11_x = " << TD11[0] << std::endl;
	std::cout << "TD11_y = " << TD11[1] << std::endl;
	std::cout << "TD11_z = " << TD11[2] << std::endl;
	VectorProduct(B11, C11, BC11);
	VectorProduct(C11, A11, CA11);
	VectorProduct(BC11, CA11, D11);
	VectorProduct(A11, B11, AB11);
	double ip12 = InnerProduct(AB11, D11, 3);
	std::cout << "2. Inner Product = " << ip12 << std::endl;
	double tp12 = TripleProduct(A11, B11, C11);
	std::cout << "2. Triple Product ^ 2 = " << tp12 * tp12 << std::endl;
	std::vector<double> ABC11(3), BAC11(3), CAB11(3);
	VectorProduct(A11, BC11, ABC11);
	VectorProduct(B11, CA11, BAC11);
	VectorProduct(C11, AB11, CAB11);
	double zx = ABC11[0] + BAC11[0] + CAB11[0];
	double zy = ABC11[1] + BAC11[1] + CAB11[1];
	double zz = ABC11[2] + BAC11[2] + CAB11[2];
	std::cout << "3. Zero Vector = " << zx << ' ' << zy << ' ' << zz;
	std::cout << std::endl;
}

int main()
{
	Exercises_Section_1_13_1_Triple_Products();
	Exercises_Section_1_14_Vector_Identities();
	return 0;
}