Category: Computer Science
Comparison of Two Modern Primality Tests: AKS and APR by James Pate Williams, Jr.
Numerical Differentiation in C# by James Pate Williams, Jr. Translated from FORTRAN 77 Formulas
using System;
using System.Windows.Forms;
namespace NumericalDifferentiation
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
comboBox1.SelectedIndex = 0;
}
private double F1(double x)
{
return x * (x * (x - 2) + 4) - 1;
}
private double F2(double x)
{
return x * Math.Exp(-x) + 1;
}
private double F3(double x)
{
return x * x * Math.Sin(x) + 1;
}
private double DF1(double x)
{
return 3 * x * x - 4 * x + 4;
}
private double DF2(double x)
{
return Math.Exp(-x) - x * Math.Exp(-x);
}
private double DF3(double x)
{
return 2 * x * Math.Sin(x) + x * x * Math.Cos(x);
}
private double DDF1(double x)
{
return 6 * x - 4;
}
private double DDF2(double x)
{
return -2 * Math.Exp(-x) + x * Math.Exp(-x);
}
private double DDF3(double x)
{
return
2 * Math.Sin(x) + 2 * x * Math.Cos(x) +
2 * x * Math.Cos(x) - x * x * Math.Sin(x);
}
private void button1_Click(object sender, EventArgs e)
{
double a = (double)numericUpDown1.Value;
double b = (double)numericUpDown2.Value;
int n = (int)numericUpDown3.Value;
double deriv1c, deriv1f;
double deriv2c, deriv2f;
double x = a, y, u, v;
double error1c, error1f, error2c, error2f;
double h = (b - a) / n;
if (comboBox1.SelectedIndex == 0)
textBox1.Text += "x\t f'(x) e f'(x) 1 f'(x) 2 " +
"error1 error2 f''(x)e f''(x)1 f''(x)2 " +
"error1 error2\r\n";
else if (comboBox1.SelectedIndex == 1)
textBox1.Text += "x\t g'(x) e g'(x) 1 g'(x) 2 " +
"error1 error2 g''(x)e g''(x)1 g''(x)2 " +
"error1 error2\r\n";
else if (comboBox1.SelectedIndex == 2)
textBox1.Text += "x\t h'(x) e h'(x) 1 h'(x) 2 " +
"error1 error2 h''(x)e h''(x)1 h''(x)2 " +
"error1 error2\r\n";
while (x <= b)
{
switch (comboBox1.SelectedIndex)
{
case 0:
y = F1(x);
u = DF1(x);
v = DDF1(x);
deriv1c = Differentiation.Derivative1sc(x, h, F1);
deriv1f = Differentiation.Derivative1sf(x, h, F1);
deriv2c = Differentiation.Derivative2sc(x, h, F1);
deriv2f = Differentiation.Derivative2sf(x, h, F1);
error1c = Math.Abs(deriv1c - u);
error1f = Math.Abs(deriv1f - u);
error2c = Math.Abs(deriv2c - v);
error2f = Math.Abs(deriv2f - v);
textBox1.Text +=
x.ToString("F5").PadLeft(8) + " " +
u.ToString("F5").PadLeft(8) + " " +
deriv1c.ToString("F5").PadLeft(8) + " " +
deriv1f.ToString("F5").PadLeft(8) + " " +
error1c.ToString("F5").PadLeft(8) + " " +
error1f.ToString("F5").PadLeft(8) + " " +
v.ToString("F5").PadLeft(8) + " " +
deriv2c.ToString("F5").PadLeft(8) + " " +
deriv2f.ToString("F5").PadLeft(8) + " " +
error2c.ToString("F5").PadLeft(8) + " " +
error2f.ToString("F5").PadLeft(8) + "\r\n";
break;
case 1:
y = F2(x);
u = DF2(x);
v = DDF2(x);
deriv1c = Differentiation.Derivative1sc(x, h, F2);
deriv1f = Differentiation.Derivative1sf(x, h, F2);
deriv2c = Differentiation.Derivative2sc(x, h, F2);
deriv2f = Differentiation.Derivative2sf(x, h, F2);
error1c = Math.Abs(deriv1c - u);
error1f = Math.Abs(deriv1f - u);
error2c = Math.Abs(deriv2c - v);
error2f = Math.Abs(deriv2f - v);
textBox1.Text +=
x.ToString("F5").PadLeft(8) + " " +
u.ToString("F5").PadLeft(8) + " " +
deriv1c.ToString("F5").PadLeft(8) + " " +
deriv1f.ToString("F5").PadLeft(8) + " " +
error1c.ToString("F5").PadLeft(8) + " " +
error1f.ToString("F5").PadLeft(8) + " " +
v.ToString("F5").PadLeft(8) + " " +
deriv2c.ToString("F5").PadLeft(8) + " " +
deriv2f.ToString("F5").PadLeft(8) + " " +
error2c.ToString("F5").PadLeft(8) + " " +
error2f.ToString("F5").PadLeft(8) + "\r\n";
break;
case 2:
y = F3(x);
u = DF3(x);
v = DDF3(x);
deriv1c = Differentiation.Derivative1sc(x, h, F3);
deriv1f = Differentiation.Derivative1sf(x, h, F3);
deriv2c = Differentiation.Derivative2sc(x, h, F3);
deriv2f = Differentiation.Derivative2sf(x, h, F3);
error1c = Math.Abs(deriv1c - u);
error1f = Math.Abs(deriv1f - u);
error2c = Math.Abs(deriv2c - v);
error2f = Math.Abs(deriv2f - v);
textBox1.Text +=
x.ToString("F5").PadLeft(8) + " " +
u.ToString("F5").PadLeft(8) + " " +
deriv1c.ToString("F5").PadLeft(8) + " " +
deriv1f.ToString("F5").PadLeft(8) + " " +
error1c.ToString("F5").PadLeft(8) + " " +
error1f.ToString("F5").PadLeft(8) + " " +
v.ToString("F5").PadLeft(8) + " " +
deriv2c.ToString("F5").PadLeft(8) + " " +
deriv2f.ToString("F5").PadLeft(8) + " " +
error2c.ToString("F5").PadLeft(8) + " " +
error2f.ToString("F5").PadLeft(8) + "\r\n";
break;
default:
break;
}
x += h;
}
textBox1.Text += comboBox1.SelectedItem;
}
private void button2_Click(object sender, EventArgs e)
{
textBox1.Text = string.Empty;
}
}
}
// Derivative formulas translated from FORTRAN 77
// source code found in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 18 - 19, 2023
using System;
namespace NumericalDifferentiation
{
// Numerically compute first and second
// ordinary derivatives
class Differentiation
{
static public double Derivative1sc(
double a, double h, Func<double, double> f)
{
return (f(a + h) - f(a - h)) / (h + h);
}
static public double Derivative1sf(
double a, double h, Func<double, double> f)
{
return (-3.0 * f(a) + 4.0 * f(a + h) - f(a + h + h)) / (h + h);
}
static public double Derivative2sc(
double a, double h, Func<double, double> f)
{
return (f(a - h) - 2.0 * f(a) + f(a + h)) / (h + h);
}
static public double Derivative2sf(
double a, double h, Func<double, double> f)
{
return (f(a) - 2.0 * f(a + h) + f(a + h + h)) / (h + h);
}
}
}
Adams-Moulton Predictor Corrector Method to Solve a One-Dimensional Ordinary Differential Equation Translated from FORTRAN 77 Code by James Pate Williams, Jr.
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
namespace AdamsMoultonMethod1
{
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
}
private double Solution(double x)
{
return Math.Exp(x) - 1.0 - x;
}
private void button1_Click(object sender, EventArgs e)
{
int numberSteps = (int)numericUpDown1.Value;
AdamsMoulton.AdamsMoultonMethod(
numberSteps, Solution, textBox1);
}
private void button2_Click(object sender, EventArgs e)
{
textBox1.Text = string.Empty;
}
}
}
// Adams-Moulton Method translated from FORTRAN 77
// source code found in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 18, 2023
using System;
using System.Windows.Forms;
namespace AdamsMoultonMethod1
{
public struct Point
{
public double x, y;
}
class AdamsMoulton
{
static private string FormatNumber(double x)
{
return x.ToString("E11").PadLeft(12);
}
static private string FormatLine(
int n, double x, double y, double delta, double error)
{
return
n.ToString("D2") + " " +
FormatNumber(x) + " " +
FormatNumber(y) + " " +
FormatNumber(delta) + " " +
FormatNumber(error) + "\r\n";
}
static public void AdamsMoultonMethod(
int numberSteps, Func<double, double> solution,
TextBox textBox)
{
double h = 1.0 / numberSteps;
double xn = 0.0, yn = 0.0, yofxn;
double xBegin = 0.0, yBegin = 0.0;
double fnPredict, ynPredict;
double delta = 0.0, error = 0.0;
double[] f = new double[numberSteps];
int n = 0;
textBox.Text += "n\txn\tyn\tdelta\terror\r\n";
textBox.Text += FormatLine(n, xn, yn, delta, error);
f[1] = xBegin + yBegin;
for (n = 1; n <= 3; n++)
{
xn = xBegin + n * h;
yn = solution(xn);
f[n + 1] = xn + yn;
textBox.Text += FormatLine(n, xn, yn, delta, error);
}
for (n = 4; n <= numberSteps; n++)
{
ynPredict =
yn + h / 24.0 * (55.0 * f[4] - 59.0 * f[3] +
37.0 * f[2] - 9.0 * f[1]);
xn = xBegin + n * h;
fnPredict = xn + ynPredict;
yn += h / 24.0 * (9.0 * fnPredict + 19.0 * f[4]
- 5.0 * f[3] + f[2]);
delta = (yn - ynPredict) / 14.0;
f[1] = f[2];
f[2] = f[3];
f[3] = f[4];
f[4] = xn + yn;
yofxn = solution(xn);
error = yn - yofxn;
textBox.Text += FormatLine(n, xn, yn, delta, error);
}
}
}
}
Determinant Calculations: Bareiss Algorithm and Gaussian Elimination by James Pate Williams, Jr.
Below are three screenshots of two methods of calculating the determinant of a matrix, namely the Bareiss Algorithm and Gaussian Elimination:
using System;
using System.Diagnostics;
using System.Windows.Forms;
namespace MatrixInverseComparison
{
public partial class MainForm : Form
{
public MainForm()
{
InitializeComponent();
}
static private string FormatNumber(double x)
{
string result = string.Empty;
if (x > 0)
result += x.ToString("F5").PadLeft(10);
else
result += x.ToString("F5").PadLeft(10);
return result;
}
private void MultiplyPrintMatricies(
double[,] A, double[,] B, int n)
{
double[,] I = new double[n, n];
textBox1.Text += "Matrix Product:\r\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
double sum = 0.0;
for (int k = 0; k < n; k++)
sum += A[i, k] * B[k, j];
I[i, j] = sum;
}
}
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
textBox1.Text += FormatNumber(I[i, j]) + " ";
}
textBox1.Text += "\r\n";
}
}
private void PrintMatrix(string title, double[,] A, int n)
{
textBox1.Text += title + "\r\n";
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
textBox1.Text += FormatNumber(A[i, j]) + " ";
}
textBox1.Text += "\r\n";
}
}
private void Compute(double[,] MI, int n)
{
double determinantGE = 1;
double[] b = new double[n];
double[] x = new double[n];
double[,] MB1 = new double[n, n];
double[,] MB2 = new double[n, n];
double[,] MG = new double[n, n];
double[,] MS = new double[n, n];
double[,] IG = new double[n, n];
double[,] IS = new double[n, n];
int[] pivot = new int[n];
Stopwatch sw = new();
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
MB1[i, j] = MB2[i, j] = MG[i, j] = MS[i, j] = MI[i, j];
}
PrintMatrix("Initial Matrix: ", MI, n);
sw.Start();
int flag = DirectMethods.Factor(MG, n, pivot);
if (flag != 0)
{
for (int i = 0; i < n; i++)
determinantGE *= MG[i, i];
determinantGE *= flag;
}
sw.Stop();
PrintMatrix("Gaussian Elimination Final:", MG, n);
for (int i = 0; i < n; i++)
b[i] = 0;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
b[j] = 1.0;
DirectMethods.Substitute(MG, b, x, n, pivot);
for (int k = 0; k < n; k++)
IG[k, j] = x[k];
b[j] = 0.0;
}
}
PrintMatrix("Gaussian Elimination Inverse:", IG, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinantGE) + "\r\n";
MultiplyPrintMatricies(IG, MI, n);
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
sw.Start();
double determinant1 = BareissAlgorithm.Determinant1(MB1, n);
sw.Stop();
PrintMatrix("Bareiss Algorithm Final 1:", MB1, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinant1) + "\r\n";
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
sw.Start();
double determinant2 = BareissAlgorithm.Determinant2(MB2, n);
sw.Stop();
PrintMatrix("Bareiss Algorithm Final 2:", MB2, n);
textBox1.Text += "Determinant: " +
FormatNumber(determinant2) + "\r\n";
textBox1.Text += "Runtime (MS) = " +
sw.ElapsedMilliseconds + "\r\n";
}
private void button1_Click(object sender, EventArgs e)
{
int n = (int)numericUpDown1.Value;
int seed = (int)numericUpDown2.Value;
double[,] A = new double[n, n];
Random random = new Random(seed);
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
double q = n * random.NextDouble();
while (q == 0.0)
q = n * random.NextDouble();
A[i, j] = q;
}
}
Compute(A, n);
}
private void button2_Click(object sender, EventArgs e)
{
textBox1.Text = string.Empty;
}
}
}
using System;
namespace MatrixInverseComparison
{
public class DirectMethods
{
// Substitute and Factor translated from FORTRAN 77
// source code found in "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte and Carl de
// Boor. Translator: James Pate Williams, Jr. (c)
// August 14 - 17, 2023
static public void Substitute(
double[,] w,
double[] b,
double[] x,
int n,
int[] pivot)
{
double sum;
int i, j, n1 = n - 1;
if (n == 1)
{
x[0] = b[0] / w[0, 0];
return;
}
// forward substitution
x[0] = b[pivot[0]];
for (i = 0; i < n; i++)
{
for (j = 0, sum = 0.0; j < i; j++)
sum += w[i, j] * x[j];
x[i] = b[pivot[i]] - sum;
}
// backward substitution
x[n1] /= w[n1, n1];
for (i = n - 2; i >= 0; i--)
{
for (j = i + 1, sum = 0.0; j < n; j++)
sum += w[i, j] * x[j];
x[i] = (x[i] - sum) / w[i, i];
}
}
// Factor returns +1 if an even number of exchanges
// Factor returns -1 if an odd number of exchanges
// Factor retrurn 0 if matrix is singular
static public int Factor(
double[,] w, int n, int[] pivot)
// returns 0 if matrix is singular
{
double awikod, col_max, ratio, row_max, temp;
double[] d = new double[n];
int flag = 1, i, i_star, j, k;
for (i = 0; i < n; i++)
{
pivot[i] = i;
row_max = 0;
for (j = 0; j < n; j++)
row_max = Math.Max(row_max, Math.Abs(w[i, j]));
if (row_max == 0)
{
flag = 0;
row_max = 1;
}
d[i] = row_max;
}
if (n <= 1)
return flag;
// factorization
for (k = 0; k < n - 1; k++)
{
// determine pivot row the row i_star
col_max = Math.Abs(w[k, k]) / d[k];
i_star = k;
for (i = k + 1; i < n; i++)
{
awikod = Math.Abs(w[i, k]) / d[i];
if (awikod > col_max)
{
col_max = awikod;
i_star = i;
}
}
if (col_max == 0)
flag = 0;
else
{
if (i_star > k)
{
// make k the pivot row by
// interchanging with i_star
flag *= -1;
i = pivot[i_star];
pivot[i_star] = pivot[k];
pivot[k] = i;
temp = d[i_star];
d[i_star] = d[k];
d[k] = temp;
for (j = 0; j < n; j++)
{
temp = w[i_star, j];
w[i_star, j] = w[k, j];
w[k, j] = temp;
}
}
// eliminate x[k]
for (i = k + 1; i < n; i++)
{
w[i, k] /= w[k, k];
ratio = w[i, k];
for (j = k + 1; j < n; j++)
w[i, j] -= ratio * w[k, j];
}
}
}
if (w[n - 1, n - 1] == 0)
flag = 0;
return flag;
}
}
}
namespace MatrixInverseComparison
{
// One implementation is based on https://en.wikipedia.org/wiki/Bareiss_algorithm
// Another perhaps better implementation is found on the following webpage
// https://cs.stackexchange.com/questions/124759/determinant-calculation-bareiss-vs-gauss-algorithm
class BareissAlgorithm
{
static public double Determinant1(double[,] M, int n)
{
double M00;
for (int k = 0; k < n; k++)
{
if (k - 1 == -1)
M00 = 1;
else
M00 = M[k - 1, k - 1];
for (int i = k + 1; i < n; i++)
{
for (int j = k + 1; j < n; j++)
{
M[i, j] = (
M[i, j] * M[k, k] -
M[i, k] * M[k, j]) / M00;
}
}
}
return M[n - 1, n - 1];
}
static public double Determinant2(double[,] A, int dim)
{
if (dim <= 0)
{
return 0;
}
double sign = 1;
for (int k = 0; k < dim - 1; k++)
{
//Pivot - row swap needed
if (A[k, k] == 0)
{
int m;
for (m = k + 1; m < dim; m++)
{
if (A[m, k] != 0)
{
double[] tempRow = new double[dim];
for (int i = 0; i < dim; i++)
tempRow[i] = A[m, i];
for (int i = 0; i < dim; i++)
A[m, i] = A[k, i];
for (int i = 0; i < dim; i++)
A[k, i] = tempRow[i];
sign = -sign;
break;
}
}
//No entries != 0 found in column k -> det = 0
if (m == dim)
{
return 0;
}
}
//Apply formula
for (int i = k + 1; i < dim; i++)
{
for (int j = k + 1; j < dim; j++)
{
A[i, j] = A[k, k] * A[i, j] - A[i, k] * A[k, j];
if (k >= 1)
{
A[i, j] /= A[k - 1, k - 1];
}
}
}
}
return sign * A[dim - 1, dim - 1];
}
}
}
Microsoft C# Constraint Satisfaction Dynamic Link Library (DLL) and Test Application
Constraint Satisfaction DLL Source Code Files Source Code File Lines of Code Arc.cs 42 Common.cs 214 Label.cs 62 NQPArcConsisitencies.cs 148 NQPBacktrack.cs 51 NQPSosicGu.cs 234 Vertex.cs 18 YakooGCPWCSA.cs 160 Total 929 Source Code File Lines of Code GCPGraphForm.cs 173 MainForm.cs 42 NQPGraphForm.cs 137 NQPTableForm.cs 223 Total 575 Source 1. Foundations of Constraint Satisfaction by Edward Tsang
Sosic and Gu Fast N-Queens Problem Solver Implemented by James Pate Williams, Jr.
I implemented the Sosic and Gu fast N-Queens solver first in C++ and later in C#. The five screen shots below were created by latest C# application. None of my other algorithms can compete with this method.
More N-Queens Puzzle Results by James Pate Williams, Jr.
N-Queens Tests 8 Queens Microseconds 50 Samples Algorithm Minimum Average Maximum Std Dev Arc-Consistency 3 1846 4148 14220 2294 Arc-Consistency 4 3850 8987 27903 5022 Backjump 175 1057 3873 883 Backmark 198 737 3323 596 Backtracking 171 789 3132 652 Hill-Climber 1991 2097 2913 157 9 Queens Microseconds 50 Samples Algorithm Minimum Average Maximum Std Dev Arc-Consistency 3 2628 5785 15748 3202 Arc-Consistency 4 7492 18057 56587 11114 Backjump 208 1410 6184 1325 Backmark 214 866 4468 880 Backtracking 210 1269 6156 1186 Hill-Climber 1309 1398 1953 102 10 Queens Microseconds 50 Samples Algorithm Minimum Average Maximum Std Dev Arc-Consistency 3 3820 10680 30708 6963 Arc-Consistency 4 12624 46645 148221 31212 Backjump 248 3970 16736 3805 Backmark 283 1504 7613 1433 Backtracking 227 2857 12882 3022 Hill-Climber 24575 25608 33563 2037
My hill-climber results are pretty disappointing. I have gotten some decent runs of back-marking for n = 50 and n = 60 queens. I show the n = 60 queens experiment below.
==Menu==
**Instance Submenu**
1 Single Instance Test
2 Multiple Instances Test
3 Exit
Enter an option '1' to '3': 2
**Algorithm Submenu**
1 Arc-Consistency 3
2 Arc-Consistency 4
3 Backjump
4 Backmark
5 Chronological Bactracking
6 Evolutionary Hill-Climber
7 Exit
Enter an algorithm ('1' to '7'): 4
Enter the number of queens for option '2': 60
Enter the number of samples to be analyzed: 10
Minimum Runtime in Microseconds: 4819
Average Runtime in Microseconds: 10419330
Maximum Runtime in Micorseconds: 45439444
Standard Sample Deviation: 16665688
Notice the misspelling of microseconds.
Six C++ Algorithms to Solve the N-Queens Puzzle (Problem) by James Pate Williams, Jr.
Most of the six algorithms were developed back in the years 1999 to 2001 when I was a Master of Software Engineering and Doctoral candidate at Auburn University. Professor Gerry V. Dozier was my Master research advisor. The algorithms developed were in alphabetical order:
- Arc-Consistency 3
- Arc-Consistency 4
- Backjump
- Backmark
- Chronological Backtracking
- Evolutionary Hill-Climber (My algorithm)
The first five algorithms were stated in pseudocode in the excellent treatise Foundations of Constraint Satisfaction by Edward Tsang. Below are six runs of my test C++ application (Win32 Console App) using 8 queens and 50 samples.
The Graph Coloring Problem by James Pate Williams, Jr.
The Graph Coloring Problem and N-Queens Puzzle are both NP-Complete constraint satisfaction problems. There does not exist polynomial time algorithms to solve large instances of either problem. The Graph Coloring Problem is given a graph of n-vertices and m-edges find the minimum number of colors such that no two connected vertices are of the same color. The N-Queens Puzzle is given a n-by-n chessboard and n-queens place the queens such that no two queens are attacking one another. A queen in chess can move any number of squares diagonally, horizontally, or vertically.
As a Master of Software Engineering graduate student in the era 1998 to 2000, I majored in the subdiscipline of artificial intelligence known as constraint satisfaction. I developed evolutionary hill-climbing algorithms to solve constraint satisfaction problems. I specifically solved instances of the Graph Coloring Problem and the N-Queens Puzzle. I compared my algorithms with other researchers’ algorithms. One of the algorithms tested was Makato Yokoo’s Weak-Commitment Search Algorithm (WCSA). Yakoo wrote a very nice book named “Distributed Constraint Satisfaction” which has excellent flow-charts of several algorithms including the Min-Conflicts Search Algorithm (MCSA) and his WCSA. I implemented both algorithms in the C++ and Java computer languages.
I rewrote the preceding two algorithms in the C# computer language sometime in the epoch 2008 – 2009. Below are three instances of the Graph Coloring Problem with four, eight, and twelve vertices solved by the Yakoo’s WCSA.
Numerical Explorations 