// "Numerical Computation 2: Methods, Software,
// and Analysis" by Christoph W. Ueberhuber
// Chapter 17 Random Numbers
// "The Art of Computer Programming Volume 2"
// "Seminumerical Algorithms Second Edition"
// "Chapter 3 RANDOM NUMBERS" by Donald E. Knuth
// https://en.wikipedia.org/wiki/Mersenne_Twister
using System.Collections.Generic;
namespace SimplePRNGs
{
class PRNGs
{
private long AMxn1, AMxn;
private long AMyn1, AMyn;
private long AMk;
private long AMm = 34359738368;
private long LCG0z0, LCG0z1;
private long LCG1z0, LCG1z1;
private long LCG2z0, LCG2z1, LCG2z2;
private long MFG0z0, MFG0z1;
private readonly long LCG0m = 4294967296;
private readonly long LCG2m = 2147483647;
private readonly List<long> AMV = new();
private long MTindex;
private long[] MT;
private LaggedFibRng fibRng;
public void SetSeedLCG0(long z0)
{
LCG0z0 = z0;
}
public void SetSeedLCG1(long z0)
{
LCG1z0 = z0;
}
public void SetSeedLCG2(long z0, long z1)
{
LCG2z0 = z0;
LCG2z1 = z1;
}
public long LCG0()
{
LCG0z1 = (69069 * LCG0z0) % LCG0m;
LCG0z0 = LCG0z1;
return LCG0z1;
}
public long LCG1()
{
LCG1z1 = (69069 * LCG1z0 + 1) % LCG0m;
LCG1z0 = LCG1z1;
return LCG1z1;
}
public long LCG2()
{
LCG2z2 = (1999 * LCG2z1 + 4444 * LCG2z0) % LCG2m;
LCG2z0 = LCG2z1;
LCG2z1 = LCG2z2;
return LCG2z2;
}
public void SetSeedMFG0(long z0, long z1)
{
MFG0z0 = z0;
MFG0z1 = z1;
}
public long MFG0()
{
long MFG0z2 = (MFG0z1 + MFG0z0) % LCG0m;
MFG0z0 = MFG0z1;
MFG0z1 = MFG0z2;
return MFG0z2;
}
public void ComputeNextXY()
{
AMxn1 = (3141592653 * AMxn + 2718281829) % AMm;
if (AMxn1 < 0)
AMxn1 += AMm;
AMyn1 = (2718281829 * AMyn + 3141592653) % AMm;
if (AMyn1 < 0)
AMyn1 += AMm;
}
public void AMSeed(long k, long x0, long y0)
{
long AMTxn1, AMTxn = x0;
AMxn = x0;
AMyn = y0;
AMk = k;
for (int i = 0; i < k; i++)
{
AMTxn1 = (3141592653 * AMTxn + 2718281829) % AMm;
if (AMTxn1 < 0)
AMTxn1 += AMm;
AMTxn = AMTxn1;
AMV.Add(AMTxn1);
}
}
public long AlgorithmM()
{
ComputeNextXY();
AMxn = AMxn1;
AMyn = AMyn1;
long j = (AMk * AMyn1) / AMm;
long r = AMV[(int)j];
AMV[(int)j] = AMxn1;
if (r < 0)
r += AMm;
return r;
}
public void MTInitialization(long seed)
{
long f = 6364136223846793005;
long n = 312, w = 64;
MTindex = n;
MT = new long[n];
MT[0] = seed;
for (int i = 1; i < n; i++)
MT[i] = f * (MT[i - 1] ^ (MT[i - 1] >> (int)(w - 2))) + i;
}
public long MTExtractNumber()
{
unchecked
{
long n = 312;
long c = (long)0xFFF7EEE000000000;
long b = 0x71D67FFFEDA60000;
long d = 0x5555555555555555;
long u = 29, s = 17, t = 27, l = 43;
if (MTindex == n)
MTTwist();
long y = MT[MTindex];
y ^= ((y >> (int)u) & d);
y ^= ((y << (int)s) & b);
y ^= ((y << (int)t) & c);
y ^= (y >> (int)l);
MTindex++;
return y;
}
}
public void MTTwist()
{
unchecked
{
long n = 312, m = 156, r = 31;
long a = (long)0xB5026F5AA96619E9;
MTindex = n + 1;
long lower_mask = (1 << (int)r) - 1;
long upper_mask = ~lower_mask;
for (int i = 0; i < n; i++)
{
long x = (MT[i] & upper_mask) |
(MT[(i + 1) % n] & lower_mask);
long xA = x >> 1;
if (x % 2 != 0)
xA ^= a;
MT[i] = MT[(i + m) % n] ^ xA;
}
}
MTindex = 0;
}
public void LaggedFibRngSeed(int seed)
{
fibRng = new LaggedFibRng(seed);
}
public long LaggedFibonacci(long modulus)
{
long lo = fibRng.Next();
long hi = fibRng.Next();
long rs = ((hi << 31) | lo) % modulus;
return rs;
}
}
}
//https://learn.microsoft.com/en-us/archive/msdn-magazine/2016/august/test-run-lightweight-random-number-generation
// modified by current author James Pate Williams, Jr. on August 30, 2023
using System.Collections.Generic;
namespace SimplePRNGs
{
public class LaggedFibRng
{
private const int k = 606; // Largest magnitude"-index"
private const int j = 273; // Other "-index"
private const long m = 4294967296; // 2^32
private readonly List<long> vals = null;
private long seed;
public LaggedFibRng(int seed)
{
vals = new List<long>();
for (int i = 0; i < k + 1; ++i)
vals.Add(i);
if (seed % 2 == 0) vals[0] = 11;
// Burn some values away
for (int ct = 0; ct < 1000; ++ct)
{
long dummy = Next();
}
} // ctor
public long Next()
{
// (a + b) mod n = [(a mod n) + (b mod n)] mod n
long left = vals[0] % m; // [x-big]
long right = vals[k - j] % m; // [x-other]
long sum = (left + right) % m; // prevent overflow
if (sum < 0)
seed = sum + m;
else
seed = sum;
vals.Insert(k + 1, seed); // Add new val at end
vals.RemoveAt(0); // Delete now irrelevant [0] val
return seed;
}
}
}
As can be seen the C++ application appears to be much faster than the C# application. Also, in my humble opinion, C++ with header files and C++ source code files is much more elegant than C# with only class source code. I leave it to someone else to add C# interfaces to my class source code files.
A few years ago, I implemented the Shanks-Mestre elliptic curve point counting algorithm in C# using the BigInteger data structure and the algorithm found in Henri Cohen’s textbook A Course in Computational Algebraic Number Theory. I translated the C# application to C++ in August 21-22, 2023. The C++ code uses the long long 64-bit signed integer data type. Below are some results. I also implemented Schoof’s point counting algorithm in C#.
The exclusive or (XOR) function is a very well known simple two input binary function in the world of computer architecture. The truth table has four rows with two input columns 00, 01, 10, and 11. The output column is 0, 1, 1, 0. In other words, the XOR function is false for like inputs and true for unlike inputs. This function is also known in the world of cryptography and is the encryption and decryption function for the stream cipher called the onetime pad. The onetime pad is perfectly secure as long as none of the pad is resused. The original artificial feedforward neural network was unable to correctly generate the outputs of the XOR function. This disadvantage was overcome by the addition of back-propagation to the artificial neural network. The outputs must be adjusted to 0.1 and 0.9 for the logistic also known as the sigmoid activation function to work.
using System;
using System.Windows.Forms;
namespace BPNNTestOne
{
public class BPNeuralNetwork
{
private static double RANGE = 0.1;
private double learningRate, momentum, tolerance;
private double[] deltaHidden, deltaOutput, h, o, p;
private double[,] newV, newW, oldV, oldW, v, w, O, T, X;
private int numberHiddenUnits, numberInputUnits;
private int numberOutputUnits, numberTrainingExamples;
private int maxEpoch;
private Random random;
private TextBox tb;
public double[] P
{
get
{
return p;
}
}
double random_range(double x_0, double x_1)
{
double temp;
if (x_0 > x_1)
{
temp = x_0;
x_0 = x_1;
x_1 = temp;
}
return (x_1 - x_0) * random.NextDouble() + x_0;
}
public BPNeuralNetwork
(
double alpha,
double eta,
double threshold,
int nHidden,
int nInput,
int nOutput,
int nTraining,
int nMaxEpoch,
int seed,
TextBox tbox
)
{
int i, j, k;
random = new Random(seed);
learningRate = eta;
momentum = alpha;
tolerance = threshold;
numberHiddenUnits = nHidden;
numberInputUnits = nInput;
numberOutputUnits = nOutput;
numberTrainingExamples = nTraining;
maxEpoch = nMaxEpoch;
tb = tbox;
h = new double[numberHiddenUnits + 1];
o = new double[numberOutputUnits + 1];
p = new double[numberOutputUnits + 1];
h[0] = 1.0;
deltaHidden = new double[numberHiddenUnits + 1];
deltaOutput = new double[numberOutputUnits + 1];
newV = new double[numberHiddenUnits + 1, numberOutputUnits + 1];
oldV = new double[numberHiddenUnits + 1, numberOutputUnits + 1];
v = new double[numberHiddenUnits + 1, numberOutputUnits + 1];
newW = new double[numberInputUnits + 1, numberHiddenUnits + 1];
oldW = new double[numberInputUnits + 1, numberHiddenUnits + 1];
w = new double[numberInputUnits + 1, numberHiddenUnits + 1];
for (j = 0; j <= numberHiddenUnits; j++)
{
for (k = 1; k <= numberOutputUnits; k++)
{
oldV[j, k] = random_range(-RANGE, +RANGE);
v[j, k] = random_range(-RANGE, +RANGE);
}
}
for (i = 0; i <= numberInputUnits; i++)
{
for (j = 0; j <= numberHiddenUnits; j++)
{
oldW[i, j] = random_range(-RANGE, +RANGE);
w[i, j] = random_range(-RANGE, +RANGE);
}
}
O = new double[numberTrainingExamples + 1, numberOutputUnits + 1];
T = new double[numberTrainingExamples + 1, numberInputUnits + 1];
X = new double[numberTrainingExamples + 1, numberInputUnits + 1];
}
public void SetTrainingExample(double[] t, double[] x, int d)
{
int i, k;
for (i = 0; i <= numberInputUnits; i++)
X[d, i] = x[i];
for (k = 0; k <= numberOutputUnits; k++)
T[d, k] = t[k];
}
private double f(double x)
// squashing function = a sigmoid function
{
return 1.0 / (1.0 + Math.Exp(-x));
}
private double g(double x)
// derivative of the squashing function
{
return f(x) * (1.0 - f(x));
}
public void ForwardPass(double[] x)
{
double sum;
int i, j, k;
for (j = 1; j <= numberHiddenUnits; j++)
{
for (i = 0, sum = 0; i <= numberInputUnits; i++)
sum += x[i] * w[i, j];
h[j] = sum;
}
for (k = 1; k <= numberOutputUnits; k++)
{
for (j = 1, sum = h[0] * v[0, k]; j <= numberHiddenUnits; j++)
sum += f(h[j]) * v[j, k];
o[k] = sum;
p[k] = f(o[k]);
}
}
private void BackwardPass(double[] t, double[] x)
{
double sum;
int i, j, k;
for (k = 1; k <= numberOutputUnits; k++)
deltaOutput[k] = (p[k] * (1 - p[k])) * (t[k] - p[k]);
for (k = 1; k <= numberOutputUnits; k++)
{
newV[0, k] = learningRate * h[0] * deltaOutput[k];
for (j = 1; j <= numberHiddenUnits; j++)
newV[j, k] = learningRate * f(h[j]) * deltaOutput[k];
for (j = 0; j <= numberHiddenUnits; j++)
{
v[j, k] += newV[j, k] + momentum * oldV[j, k];
oldV[j, k] = newV[j, k];
}
}
for (j = 1; j <= numberHiddenUnits; j++)
{
for (k = 1, sum = 0; k <= numberOutputUnits; k++)
sum += v[j, k] * deltaOutput[k];
deltaHidden[j] = g(h[j]) * sum;
for (i = 0; i <= numberInputUnits; i++)
{
newW[i, j] = learningRate * x[i] * deltaHidden[j];
w[i, j] += newW[i, j] + momentum * oldW[i, j];
oldW[i, j] = newW[i, j];
}
}
}
public double[,] Backpropagation(bool intermediate, int number)
{
double error = double.MaxValue;
int d, k, epoch = 0;
while (epoch < maxEpoch && error > tolerance)
{
error = 0.0;
for (d = 1; d <= numberTrainingExamples; d++)
{
double[] x = new double[numberInputUnits + 1];
for (int i = 0; i <= numberInputUnits; i++)
x[i] = X[d, i];
ForwardPass(x);
double[] t = new double[numberOutputUnits + 1];
for (int i = 0; i <= numberOutputUnits; i++)
t[i] = T[d, i];
BackwardPass(t, x);
for (k = 1; k <= numberOutputUnits; k++)
{
O[d, k] = p[k];
error += Math.Pow(T[d, k] - O[d, k], 2);
}
}
epoch++;
error /= (numberTrainingExamples * numberOutputUnits);
if (intermediate && epoch % number == 0)
tb.Text += error.ToString("E10") + "\r\n";
}
tb.Text += "\r\n";
tb.Text += "Mean Square Error = " + error.ToString("E10") + "\r\n";
tb.Text += "Epoch = " + epoch + "\r\n";
return O;
}
}
}
// Learn the XOR Function
// Inputs/Targets
// 0 0 / 0
// 0 1 / 1
// 1 0 / 1
// 1 1 / 0
// The logistic function works
// with 0 -> 0.1 and 1 -> 0.9
//
// BPNNTestOne (c) 2011
// James Pate Williams, Jr.
// All rights reserved.
using System;
using System.Windows.Forms;
namespace BPNNTestOne
{
public partial class MainForm : Form
{
private int seed = 512;
private BPNeuralNetwork bpnn;
public MainForm()
{
InitializeComponent();
int numberInputUnits = 2;
int numberOutputUnits = 1;
int numberHiddenUnits = 2;
int numberTrainingExamples = 4;
double[,] X =
{
{1.0, 0.0, 0.0},
{1.0, 0.0, 0.0},
{1.0, 0.0, 1.0},
{1.0, 1.0, 0.0},
{1.0, 1.0, 1.0}
};
double[,] T =
{
{1.0, 0.0},
{1.0, 0.1},
{1.0, 0.9},
{1.0, 0.9},
{1.0, 0.1}
};
bpnn = new BPNeuralNetwork(0.1, 0.9, 1.0e-12,
numberHiddenUnits, numberInputUnits,
numberOutputUnits, numberTrainingExamples,
5000, seed, textBox1);
for (int i = 1; i <= numberTrainingExamples; i++)
{
double[] x = new double[numberInputUnits + 1];
double[] t = new double[numberOutputUnits + 1];
for (int j = 0; j <= numberInputUnits; j++)
x[j] = X[i, j];
for (int j = 0; j <= numberOutputUnits; j++)
t[j] = T[i, j];
bpnn.SetTrainingExample(t, x, i);
}
double[,] O = bpnn.Backpropagation(true, 500);
textBox1.Text += "\r\n";
textBox1.Text += "x0\tx1\tTarget\tOutput\r\n\r\n";
for (int i = 1; i <= numberTrainingExamples; i++)
{
for (int j = 1; j <= numberInputUnits; j++)
textBox1.Text += X[i, j].ToString("##0.#####") + "\t";
for (int j = 1; j <= numberOutputUnits; j++)
{
textBox1.Text += T[i, j].ToString("####0.#####") + "\t";
textBox1.Text += O[i, j].ToString("####0.#####") + "\t";
}
textBox1.Text += "\r\n";
}
}
}
}
