Tag: coding

Blog Entry © Thursday, January 23, 2025, by James Pate Williams, Jr. Merge Sort Verus Quick Sort

== Menu ==
1 Side-by-Side Tests
2 Timing Comparisons
3 Exit
Enter an option: 1
Enter a PRNG Seed >= 1: 1
0.001251        0.001251        0.001251
0.563585        0.014985        0.014985
0.193304        0.174108        0.174108
0.808740        0.193304        0.193304
0.585009        0.303995        0.303995
0.479873        0.350291        0.350291
0.350291        0.479873        0.479873
0.895962        0.513535        0.513535
0.822840        0.563585        0.563585
0.746605        0.585009        0.585009
0.174108        0.710501        0.710501
0.858943        0.746605        0.746605
0.710501        0.808740        0.808740
0.513535        0.822840        0.822840
0.303995        0.858943        0.858943
0.014985        0.895962        0.895962
== Menu ==
1 Side-by-Side Tests
2 Timing Comparisons
3 Exit
Enter an option:

== Menu ==
1 Side-by-Side Tests
2 Timing Comparisons
3 Exit
Enter an option: 2
Enter a PRNG Seed >= 1: 1
merge sort mean runtime (ms) = 523
quick sort mean runtime (ms) = 435
merge sort std dev (s) = 0.033457
quick sort std dev (s) = 0.027816
== Menu ==
1 Side-by-Side Tests
2 Timing Comparisons
3 Exit
Enter an option:
// MergeVersusQuick.c : This file contains the 'main' function.
// Program execution begins and ends there.
// mergesort is from "A Numerical Library in C for Scientists
// and Engineers" by H. T. Lau Translated from ALGOL NUMAL
// QuickSort is from "Introduction to Algorithms" by Thomas H.
// Cormen, Charles E. Leiserson, and Ronald L. Rivest

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define LENGTH1			17		// static side-by-side test
#define LENGTHM1		16		// upper index
#define LENGTH2		  4097		// time test array length
#define LENGTHM2      4096		// upper inndex
#define NUMBER_TESTS  4096		// number of timing tests

void mergesort(float a[], int p[], int low, int up)
{
	int* allocate_integer_vector(int, int);
	void free_integer_vector(int*, int);
	void merge(int, int, int, int[], float[], int[]);
	int i, lo, step, stap, umlp1, umsp1, rest, restv, * hp;

	hp = allocate_integer_vector(low, up);
	for (i = low; i <= up; i++) p[i] = i;
	restv = 0;
	umlp1 = up - low + 1;
	step = 1;
	do {
		stap = 2 * step;
		umsp1 = up - stap + 1;
		for (lo = low; lo <= umsp1; lo += stap)
			merge(lo, step, step, p, a, hp);
		rest = up - lo + 1;
		if (rest > restv && restv > 0)
			merge(lo, rest - restv, restv, p, a, hp);
		restv = rest;
		step *= 2;
	} while (step < umlp1);
	free_integer_vector(hp, low);
}

int* allocate_integer_vector(int l, int u)
{
	/*  Allocates an integer vector of range [l..u].  */

	void system_error(char*);
	int* p;

	p = (int*)malloc((unsigned)(u - l + 1) * sizeof(int));
	if (!p) system_error("Failure in allocate_integer_vector().");
	return p - l;
}

void free_integer_vector(int* v, int l)
{
	/*  Frees an integer vector of range [l..u].  */

	free((char*)(v + l));
}

void system_error(char error_message[])
{
	void exit(int);

	printf("%s", error_message);
	exit(1);
}

void merge(int lo, int ls, int rs, int p[], float a[], int hp[])
{
	/* this procedure is used internally by MERGESORT */

	int l, r, lout, rout, i, pl, pr;

	l = lo;
	r = lo + ls;
	lout = rout = 0;
	i = lo;
	do {
		pl = p[l];
		pr = p[r];
		if (a[pl] > a[pr]) {
			hp[i] = pr;
			r++;
			rout = (r == lo + ls + rs);
		}
		else {
			hp[i] = pl;
			l++;
			lout = (l == lo + ls);
		}
		i++;
	} while (!(lout || rout));
	if (rout) {
		for (i = lo + ls - 1; i >= l; i--) p[i + rs] = p[i];
		r = l + rs;
	}
	for (i = r - 1; i >= lo; i--) p[i] = hp[i];
}

int Partition(float a[], int lo, int hi)
{
	int pivotIndex = lo + (hi - lo) / 2;
	float x = a[pivotIndex];
	float t = x;

	a[pivotIndex] = a[hi];
	a[hi] = t;

	int storeIndex = lo;

	for (int i = lo; i < hi; i++)
	{
		if (a[i] < x)
		{
			t = a[i];
			a[i] = a[storeIndex];
			a[storeIndex++] = t;
		}
	}

	t = a[storeIndex];
	a[storeIndex] = a[hi];
	a[hi] = t;

	return storeIndex;
}

void DoQuickSort(float a[], int n, int p, int r)
{
	if (p < r)
	{
		int q = Partition(a, p, r);

		DoQuickSort(a, n, p, q - 1);
		DoQuickSort(a, n, q + 1, r);
	}
}

void QuickSort(float a[], int n)
{
	DoQuickSort(a, n, 1, n);
}

double runtime[2][NUMBER_TESTS];
float A[LENGTH1] = { 0 }, a[LENGTH1] = { 0 };
float AA[LENGTH2] = { 0 }, aa[LENGTH2] = { 0 };
float b[LENGTH1] = { 0 };
int n = NUMBER_TESTS;
int pp[LENGTH2];

int main()
{
	while (1)
	{
		float x;
		int i, j, option = 0, p[LENGTH1], seed = 1;

		printf("== Menu ==\n");
		printf("1 Side-by-Side Tests\n");
		printf("2 Timing Comparisons\n");
		printf("3 Exit\n");
		printf("Enter an option: ");
		scanf_s("%d", &option);

		if (option == 3)
			break;

		if (option == 1 || option == 2)
		{
			printf("Enter a PRNG Seed >= 1: ");
			scanf_s("%u", &seed);
			srand(seed);

			if (option == 1)
			{
				for (i = 1; i <= LENGTHM1; i++)
				{
					x = (float)rand() / RAND_MAX;
					A[i] = x;
					a[i] = x;
					b[i] = x;
				}

				mergesort(A, p, 1, LENGTHM1);
				QuickSort(a, LENGTHM1);

				for (i = 1; i <= LENGTHM1; i++)
					printf("%f\t%f\t%f\n", b[i], A[p[i]], a[i]);
			}

			else if (option == 2)
			{
				double mean[2] = { 0 }, median[2] = { 0 }, stdDev[2] = { 0 };
				double Sx[2] = { 0 };
				
				for (j = 0; j < n; j++)
				{
					for (i = 1; i <= LENGTHM2; i++)
					{
						x = (float)rand() / RAND_MAX;
						AA[i] = x;
						aa[i] = x;
					}
					
					clock_t start1 = clock();
					mergesort(AA, pp, 1, LENGTHM2);
					clock_t finis1 = clock();
					clock_t start2 = clock();
					QuickSort(aa, LENGTHM2);
					clock_t finis2 = clock();
					runtime[0][j] = ((double)finis1 - start1) / 
						CLOCKS_PER_SEC;
					runtime[1][j] = ((double)finis2 - start2) / 
						CLOCKS_PER_SEC;
					mean[0] += runtime[0][j];
					mean[1] += runtime[1][j];
				}

				for (j = 0; j < n; j++)
				{
					Sx[0] = pow((double)runtime[0][j] - mean[0], 2.0) / ((double)n - 1);
					Sx[1] = pow((double)runtime[1][j] - mean[1], 2.0) / ((double)n - 1);
				}

				stdDev[0] = (float)sqrt(Sx[0]);
				stdDev[1] = (float)sqrt(Sx[1]);

				printf("merge sort mean runtime (ms) = %3.0lf\n", 1.0e6 * mean[0] / n);
				printf("quick sort mean runtime (ms) = %3.0lf\n", 1.0e6 * mean[1] / n);
				printf("merge sort std dev (s) = %f\n", stdDev[0]);
				printf("quick sort std dev (s) = %f\n", stdDev[1]);
			}
		}
	}

	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 Friday, July 19, 2024, Easy Internet Math “Puzzle” (c) James Pate Williams, Jr.

Internet Math PuzzleDownload
#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);
}