Simple Pseudo-Random Number Generators in C# Implemented by James Pate Williams, Jr.

// "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

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();

        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;
        }
    }
}

Another Point-Mass Ballistics Application Battleship Iowa 16-Inch Guns in Particular

by James Pate Williams, Jr.

namespace PointMassBallistics
{
    public class TableEntry : IComparable<TableEntry>
    {
        public double range, elevationDegrees, elevationMinutes,
            angleFallDegrees, angleFallMinutes, timeOfFlight,
            strikingVelocity, maximumOrdinate;
        public int CompareTo(TableEntry other)
        {
            if (elevationDegrees < other.elevationDegrees &&
                elevationMinutes < other.elevationMinutes)
                return -1;
            if (elevationDegrees > other.elevationDegrees &&
                elevationMinutes > other.elevationMinutes)
                return +1;
            return 0;
        }
    }
}
// Solves the following system of first order
// ordinary differential equations. Formulas
// are from "Elementary Numerical Analysis:
// An Algorithmic Approach 3rd Edition" by S.
// D. Conte & Carl de Boor (c) 1980 8.12 page 398.
// Extended from two to four equations.
// See https://apps.dtic.mil/dtic/tr/fulltext/u2/a439796.pdf
// Also view https://eugeneleeslover.com/USN-GUNS-AND-RANGE-TABLES/OP-770-1.html

using System;
using System.Collections.Generic;

namespace PointMassBallistics
{
    class RungeKutta4
	{
		private double BC;
		private readonly double g = 32.17405;
		private readonly double[] GarveN = {
			2, 3, 5, 3, 2, 1.7, 1.55 };
		private readonly double[] log10K = {
			5.66989 - 10, 2.77344 - 10, 6.80187 - 20,
			2.98090 - 10, 6.11926 - 10, 7.09620 - 10, 7.60905 - 10 };
		private readonly double[] K = new double[7];
		private int zone;

		static private double Density(double y)
		{
			return Math.Pow(10, -0.00001372 * y);
		}

		static private int ComputeIndex(double v)
		{
			int index;

			if (v > 3600)
				index = 6;

			else if (v > 2600 && v <= 3600)
				index = 5;

			else if (v > 1800 && v <= 2600)
				index = 4;

			else if (v > 1370 && v <= 1800)
				index = 3;

			else if (v > 1230 && v <= 1370)
				index = 2;

			else if (v > 790 && v <= 1230)
				index = 1;

			else
				index = 0;

			return index;
		}

		public double MayevskiRetardation(double v, int zone)
		{
			// See Exterior Ballistics 1935 by Ernest Edward Herrmann
			// Garve function

			return K[zone] * Math.Pow(v, GarveN[zone]);
		}

		private double Vx(double syn, double vxn, double vyn)
		{ 
			double v = Math.Sqrt(vxn * vxn + vyn * vyn);
			zone = ComputeIndex(v);
			double E = Density(syn) * MayevskiRetardation(v, zone) / BC;
			return -E * vxn / v;
        }

		private double Vy(double syn, double vxn, double vyn)
		{
			double v = Math.Sqrt(vxn * vxn + vyn * vyn);
			zone = ComputeIndex(v);
			double E = Density(syn) * MayevskiRetardation(v, zone) / BC;
			return -E * vyn / v - g;
		}

		static private double Sx(double vxn)
        {
			return vxn;
        }

		static private double Sy(double vyn)
		{
			return vyn;
		}

		public void Solve(
			double t0, double t1,
			double vx0, double vy0,
			double sx0, double sy0,
			double BC, int nSteps, ref List<double> lt,
			ref List<double> lvx, ref List<double> lvy,
			ref List<double> lsx, ref List<double> lsy)
		{ 
			double k1, k2, k3, k4;
			double l1, l2, l3, l4;
			double m1, m2, m3, m4;
			double n1, n2, n3, n4;
			double h = (t1 - t0) / nSteps, tn = t0;
			double vxn = vx0, vyn = vy0, sxn = sx0, syn = sy0;
			int n = 1;

			for (int i = 0; i < log10K.Length; i++)
				K[i] = Math.Pow(10, log10K[i]);

			this.BC = BC;

			lt.Add(tn);
			lvx.Add(vxn);
			lvy.Add(vyn);
			lsx.Add(sxn / 3);
			lsy.Add(syn);

			while (true)
			{
				tn = t0 + n * h;
				k1 = h * Vx(syn, vxn, vyn);
				l1 = h * Vy(syn, vxn, vyn);
				m1 = h * Sx(vxn);
				n1 = h * Sy(vyn);

				k2 = h * Vx(syn + 0.5 * n1, vxn + 0.5 * k1, vyn + 0.5 * l1);
				l2 = h * Vy(syn + 0.5 * n1, vxn + 0.5 * k1, vyn + 0.5 * l1);
				m2 = h * Sx(vxn + 0.5 * m1);
				n2 = h * Sy(vyn + 0.5 * n1);

				k3 = h * Vx(syn + 0.5 * n2, vxn + 0.5 * k2, vyn + 0.5 * l2);
				l3 = h * Vy(syn + 0.5 * n2, vxn + 0.5 * k2, vyn + 0.5 * l2);
				m3 = h * Sx(vxn + 0.5 * m2);
				n3 = h * Sy(vyn + 0.5 * n2);

				k4 = h * Vx(syn + n3, vxn + k3, vyn + l3);
				l4 = h * Vy(syn + n3, vxn + k3, vyn + l3);
				m4 = h * Sx(vxn + m3);
				n4 = h * Sy(vyn + n3);

				vxn = vx0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
				vyn = vy0 + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;
				sxn = sx0 + (m1 + 2 * m2 + 2 * m3 + m4) / 6.0;
				syn = sy0 + (n1 + 2 * n2 + 2 * n3 + n4) / 6.0;

				vx0 = vxn;
				vy0 = vyn;
				sx0 = sxn;
				sy0 = syn;
				
				n++;

				lt.Add(tn);
				lvx.Add(vxn);
				lvy.Add(vyn);
				lsx.Add(sxn / 3);
				lsy.Add(syn);
				
				if (syn <= 1.0e-2)
					break;
			}
		}
	}
}

Double Spring-Mass System

// Solves the following system of first order
// ordinary differential equations. Formulas
// are from "Elementary Numerical Analysis:
// An Algorithmic Approach 3rd Edition" by S.
// D. Conte & Carl de Boor (c) 1980 8.12 page 398.
// Extended from two to four equations.
// See "Ordinary Differential Equations
// from Calculus to Dynamical Systems"
// by Virginia W. Noonburg for the exact
// solution, view pages 167 - 170.

using System.Collections.Generic;

namespace DoubleSpring
{
	class RungeKutta4
	{
		private double omega2;

		static private double Dx1(double x3)
        {
			return x3;
        }

		static private double Dx2(double x4)
		{
			return x4;
		}

		private double Dx3(double x1, double x2)
		{
			return -2 * omega2 * x1 + omega2 * x2;
		}

		private double Dx4(double x1, double x2)
		{
			return omega2 * x1 - omega2 * x2;
		}

		public void Solve(
			double t0, double t1,
			double x10, double x20,
			double x30, double x40,
			double omega, int nSteps, ref List<double> lt,
			ref List<double> lx1, ref List<double> lx2,
			ref List<double> lx3, ref List<double> lx4)
		{
			double k1, k2, k3, k4;
			double l1, l2, l3, l4;
			double m1, m2, m3, m4;
			double n1, n2, n3, n4;
			double h = (t1 - t0) / nSteps, tn = t0;
			double x1n = x10, x2n = x20, x3n = x30, x4n = x40;
			int n = 1;

			omega2 = omega * omega;

			lt.Add(tn);
			lx1.Add(x1n);
			lx2.Add(x2n);
			lx3.Add(x3n);
			lx4.Add(x4n);

			while (tn <= t1)
			{
				tn = t0 + n * h;
				k1 = h * Dx1(x3n);
				l1 = h * Dx2(x4n);
				m1 = h * Dx3(x1n, x2n);
				n1 = h * Dx4(x1n, x2n);

				k2 = h * Dx1(x3n + 0.5 * k1);
				l2 = h * Dx2(x4n + 0.5 * l1);
				m2 = h * Dx3(x1n + 0.5 * m1, x2n + 0.5 * n1);
				n2 = h * Dx4(x1n + 0.5 * m1, x2n + 0.5 * n1);

				k3 = h * Dx1(x3n + 0.5 * k2);
				l3 = h * Dx2(x4n + 0.5 * l2);
				m3 = h * Dx3(x1n + 0.5 * m2, x2n + 0.5 * n2);
				n3 = h * Dx4(x1n + 0.5 * m2, x2n + 0.5 * n2);

				k4 = h * Dx1(x3n + k3);
				l4 = h * Dx2(x4n + l3);
				m4 = h * Dx3(x1n + m3, x2n + n3);
				n4 = h * Dx4(x1n + m3, x2n + n3);

				x1n = x10 + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
				x2n = x20 + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;
				x3n = x30 + (m1 + 2 * m2 + 2 * m3 + m4) / 6.0;
				x4n = x40 + (n1 + 2 * n2 + 2 * n3 + n4) / 6.0;

				x10 = x1n;
				x20 = x2n;
				x30 = x3n;
				x40 = x4n;

				n++;

				lt.Add(tn);
				lx1.Add(x1n);
				lx2.Add(x2n);
				lx3.Add(x3n);
				lx4.Add(x4n);
			}
		}
	}
}

Two-Tank Mixing Problem C++ Implementation Using Runge-Kutta 4 for a System of Two Ordinary Differential Equations Implementation by James Pate Williams, Jr.

// Two-tank Mixing Problem
// x'(t) = f(t, x, y)
// y'(t) = g(t, x, y)
// x(0) = 5, y(0) = 1
// See “Ordinary Differential Equations
// from Calculus to Dynamical Systems”
// by Virginia W. Noonburg for the exact
// solution, view pages 165 – 167. Also
// see “Elementary Numerical Analysis:
// An Algorithmic Approach Third
// Edition” by S. D. Conte and Carl de
// Boor pages 398 – 491 for the Runge-
// Kutta-4 equations.

#include "RK4System.h"

void RK4System::Solve(
	real t0, real t1,
	real x0, real  z0,
	real(*f)(real, real, real),
	real(*g)(real, real, real),
	int nSteps, vector<real>& tv,
	vector<real>& xv, vector<real>& yv)
{
	real k1, k2, k3, k4, l1, l2, l3, l4;
	real h, tn, xn, yn, xn1, yn1;
	
	h = (t1 - t0) / nSteps;
	tn = t0;
	xn = x0;
	yn = z0;
	tv.push_back(tn);
	xv.push_back(xn);
	yv.push_back(yn);

	for (unsigned int n = 1; n <= nSteps; n++)
	{
		tn = t0 + n * h;
		k1 = h * f(tn, xn, yn);
		l1 = h * g(tn, xn, yn);
		k2 = h * f(tn + 0.5 * h, xn + 0.5 * k1, yn + 0.5 * l1);
		l2 = h * g(tn + 0.5 * h, xn + 0.5 * k1, yn + 0.5 * l1);
		k3 = h * f(tn + 0.5 * h, xn + 0.5 * k2, yn + 0.5 * l2);
		l3 = h * g(tn + 0.5 * h, xn + 0.5 * k2, yn + 0.5 * l2);
		k4 = h * f(tn + h, xn + k3, yn + l3);
		l4 = h * g(tn + h, xn + k3, yn + l3);
		xn1 = xn + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
		yn1 = yn + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;
		xn = xn1;
		yn = yn1;
		tv.push_back(tn);
		xv.push_back(xn);
		yv.push_back(yn);
	}
}
#pragma once
// Solves the following system of first order
// ordinary differential equations. Formulas
// are from "Elementary Numerical Analysis:
// An Algorithmic Approach" by S. D. Conte &
// Carl de Boor (c) 1980 8.12 page 398.
// x' = f(t, x, y)
// y' = g(t, x, y)

#include <vector>
using namespace std;

typedef long double real;

class RK4System
{
public:
	static void Solve(
		real t0, real t1,
		real x0, real  y0,
		real(*f)(real, real, real),
		real(*g)(real, real, real),
		int nSteps, vector<real>& tv,
		vector<real>& xv, vector<real>& yv);
};

Brute Force Elliptic Curve Point Counting Algorithm Implementation by James Pate Williams, Jr.

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <chrono>
#include <vector>
using namespace std;

typedef long long ll;

typedef struct ecPoint
{
	ll x, y;
} ECPOINT, *PECPOINT;

typedef struct ecPointOrder
{
	ll order;
	ECPOINT pt;
} ECPOINTORDER, * PECPOINTORDER;

int JACOBI(ll a, ll n)
{
	int e = 0, s;
	ll a1, b = a, m, n1;

	if (a == 0) return 0;
	if (a == 1) return 1;
	while ((b & 1) == 0)
	{
		b >>= 1;
		e++;
	}
	a1 = b;
	m = n % 8;
	if (!(e & 1)) s = 1;
	else if (m == 1 || m == 7) s = +1;
	else if (m == 3 || m == 5) s = -1;
	if (n % 4 == 3 && a1 % 4 == 3) s = -s;
	if (a1 != 1) n1 = n % a1; else n1 = 1;
	return s * JACOBI(n1, a1);
}

ll ExpMod(ll x, ll b, ll n)
/* returns x ^ b mod n */
{
	ll a = 1LL, s = x;

	if (b == 0)
		return 1LL;
	while (b != 0) {
		if (b & 1l) a = (a * s) % n;
		b >>= 1;
		if (b != 0) s = (s * s) % n;
	}
	if (a < 0)
		a += n;
	return a;
}

ll ExtendedEuclidean(ll b, ll n)
{
	ll b0 = b, n0 = n, t = 1, t0 = 0, temp, q, r;

	q = n0 / b0;
	r = n0 - q * b0;

	while (r > 0) {
		temp = t0 - q * t;
		if (temp >= 0) temp = temp % n;
		else temp = n - (-temp % n);
		t0 = t;
		t = temp;
		n0 = b0;
		b0 = r;
		q = n0 / b0;
		r = n0 - q * b0;
	}
	if (b0 != 1) return 0;
	else return t % n;
}

ll Weierstrass(ll a, ll b, ll x, ll p)
{
	return ((((x * x) % p) * x) % p + (a * x) % p + b) % p;
}

ll EPoints(
	bool print,
	ll a, ll b, ll p,
	vector<ECPOINT>& e)
	/* returns the number of points on the elliptic
	curve y ^ 2 = x ^ 3 + ax + b mod p */
{
	ll count = 0, m = (p + 1) / 4, x, y;
	ECPOINT pt{};

	if (p % 4 == 3) {
		for (x = 0; x < p; x++) {
			y = Weierstrass(a, b, x, p);
			if (JACOBI(y, p) != -1) {
				y = ExpMod(y, m, p);
				if (print)
				{
					cout << "(" << setw(2) << x << ", ";
					cout << y << ") " << endl;
				}
				pt.x = x;
				pt.y = y;
				e.push_back(pt);
				count++;
				y = -y % p;
				if (y < 0) y += p;
				if (y != 0)
				{
					if (print)
					{
						cout << "(" << setw(2) << x << ", ";
						cout << y << ") " << endl;
					}
					pt.x = x;
					pt.y = y;
					e.push_back(pt);
					count++;
				}
				if (print && count % 5 == 0)
					cout << endl;
			}
		}
		if (print && count % 5 != 0)
			cout << endl;
	}
	return count;
}

void Add(ll a, ll p, ECPOINT P,
	ECPOINT Q, ECPOINT& R)
	/* elliptic curve point partial addition */
{
	ll i, lambda;

	if (P.x == Q.x && P.y == 0 && Q.y == 0) {
		R.x = 0;
		R.y = 1;
		return;
	}
	if (P.x == Q.x && P.y == p - Q.y) {
		R.x = 0;
		R.y = 1;
		return;
	}
	if (P.x == 0 && P.y == 1) {
		R = Q;
		return;
	}
	if (Q.x == 0 && Q.y == 1) {
		R = P;
		return;
	}
	if (P.x != Q.x) {
		i = Q.x - P.x;
		if (i < 0) i += p;
		i = ExtendedEuclidean(i, p);
		lambda = ((Q.y - P.y) * i) % p;
	}
	else {
		i = ExtendedEuclidean((2 * P.y) % p, p);
		lambda = ((3 * P.x * P.x + a) * i) % p;
	}
	if (lambda < 0) lambda += p;
	R.x = (lambda * lambda - P.x - Q.x) % p;
	R.y = (lambda * (P.x - R.x) - P.y) % p;
	if (R.x < 0) R.x += p;
	if (R.y < 0) R.y += p;
}

void Multiply(
	ll a, ll k, ll p,
	ECPOINT P,
	ECPOINT& R)
{
	ECPOINT S;

	R.x = 0;
	R.y = 1;
	S = P;
	while (k != 0) {
		if (k & 1) Add(a, p, R, S, R);
		k >>= 1;
		if (k != 0) Add(a, p, S, S, S);
	}
}

ll Order(ll a, ll p, ECPOINT P)
{
	ll order = 1;
	ECPOINT Q = P, R{};

	do {
		order++;
		Add(a, p, P, Q, R);
		Q = R;
	} while (R.x != 0 && R.y != 1);
	return order;
}

const int PrimeSize = 10000000;
bool sieve[PrimeSize];

void PopulateSieve() {

	// sieve of Eratosthenes

	int c, inc, i, n = PrimeSize - 1;

	for (i = 0; i < n; i++)
		sieve[i] = false;

	sieve[1] = false;
	sieve[2] = true;

	for (i = 3; i <= n; i++)
		sieve[i] = (i & 1) == 1 ? true : false;

	c = 3;

	do {
		i = c * c;
		inc = c + c;

		while (i <= n) {
			sieve[i] = false;
			i += inc;
		}

		c += 2;
		while (!sieve[c])
			c++;
	} while (c * c <= n);
}

ll Partition(
	vector<ECPOINTORDER>& a, ll n, ll lo, ll hi)
{
	ll pivotIndex = lo + (hi - lo) / 2;
	ECPOINTORDER po = a[(unsigned int)pivotIndex];
	ECPOINTORDER x = po;
	ECPOINTORDER t = x;

	a[(unsigned int)pivotIndex] = a[(unsigned int)hi];
	a[(unsigned int)hi] = t;

	ll storeIndex = lo;

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

	t = a[(unsigned int)storeIndex];
	a[(unsigned int)storeIndex] = a[(unsigned int)hi];
	a[(unsigned int)hi] = t;

	return storeIndex;
}

static void DoQuickSort(
	vector<ECPOINTORDER>& a, ll n, ll p, ll r)
{
	if (p < r)
	{
		ll q = Partition(a, n, p, r);

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

void QuickSort(vector<ECPOINTORDER>& a, ll n)
{
	DoQuickSort(a, n, 0, n - 1);
}

int main()
{
	PopulateSieve();

	while (true)
	{
		bool noncyclic = false;
		int option;
		ll a, b, p;

		cout << "Enter 1 to continue or 0 to quit = ";
		cin >> option;

		if (option == 0)
			break;

		cout << "Weierstrass parameter a = ";
		cin >> a;
		cout << "Weierstrass parameter b = ";
		cin >> b;
		cout << "Enter the prime parameter p = ";
		cin >> p;

		if (!sieve[p])
		{
			cout << "p is not a prime number. " << endl;
			continue;
		}

		if (p % 4 != 3)
		{
			cout << "p mod 4 must be 3 for this algorithm." << endl;
			continue;
		}

		bool print = false, enumerate = false;
		vector<ECPOINT> e;
		vector<ECPOINTORDER> eOrder;

		cout << "Enumerate points (0 for no or 1 for yes)? ";
		cin >> option;

		enumerate = option == 1;
		auto time0 = chrono::high_resolution_clock::now();
		ll count = EPoints(print, a, b, p, e);
		auto time1 = chrono::high_resolution_clock::now();
		auto elapsed = time1 - time0;
		ll h = count + 1;
		ll maxOrder = 0;
		ll minOrder = h;
		ll order = -1;
		ECPOINT P{}, Q{};

		if (enumerate)
			cout << "The points (x, y) on E and their orders are:" << endl;

		for (unsigned int i = 0; i < count; i++)
		{
			ECPOINTORDER ptOrder{};

			order = Order(a, p, e[i]);
			ptOrder.order = order;
			ptOrder.pt = e[i];
			eOrder.push_back(ptOrder);

			if (order < h)
				noncyclic = true;
			if (order < minOrder) {
				P = e[i];
				minOrder = order;
			}
			if (order > maxOrder) {
				Q = e[i];
				maxOrder = order;
			}
		}
		if (enumerate)
		{
			QuickSort(eOrder, count);

			for (unsigned int i = 0; i < eOrder.size(); i++)
			{
				cout << "(" << eOrder[i].pt.x << ", ";
				cout << eOrder[i].pt.y << ") ";
				cout << eOrder[i].order << "\t";
				if ((i + 1) % 5 == 0)
					cout << endl;
			}
			if (count % 5 != 0)
				cout << endl;
		}
		cout << "#E = " << h << endl;
		cout << "Noncyclic = " << noncyclic << endl;
		cout << "Minimum order = " << minOrder << endl;
		cout << "Maximum order = " << maxOrder << endl;
		cout << "(" << P.x << ", " << P.y << ")" << endl;
		cout << "(" << Q.x << ", " << Q.y << ")" << endl;

		long long runtime = chrono::duration_cast<chrono::milliseconds>(elapsed).count();
		if (runtime != 0)
			cout << "Point counting runtime in milliseconds = " << runtime << endl;
		else
		{
			runtime = chrono::duration_cast<chrono::microseconds>(elapsed).count();
			cout << "Point counting runtime in microseconds = " << runtime << endl;
		}
	}
	return 0;
}

Stinson Exercise 5.6 Galois Field

/*
Author:  Pate Williams (c) 1997

Exercise
"5.6 The field GF(2 ^ 5) can be constructed as
Z_2[x]/(x ^ 5 + x ^ 2 + 1). Perform the following
computations in this field.
(a) Compute (x ^ 4 + x ^ 2) * (x ^ 3 + x + 1).
(b) Using the Extended Euclidean algorithm
compute (x ^ 3 + x ^ 2) ^ - 1.
(c) Using the square and multiply algorithm,
compute x ^ 25." -Douglas R. Stinson-
See "Cryptography: Theory and Practice" by
Douglas R. Stinson page 201.
*/

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

#define SIZE 32l

void poly_mul(long m, long n, long *a, long *b, long *c, long *p)
{
	long ai, bj, i, j, k, sum;

	*p = m + n;
	for (k = 0; k <= *p; k++) {
		sum = 0;
		for (i = 0; i <= k; i++) {
			j = k - i;
			if (i > m) ai = 0; else ai = a[i];
			if (j > n) bj = 0; else bj = b[j];
			sum += ai * bj;
		}
		c[k] = sum;
	}
}

void poly_div(long m, long n, long *u, long *v,
	long *q, long *r, long *p, long *s)
{
	long j, jk, k, nk, vn = v[n];

	for (j = 0; j <= m; j++) r[j] = u[j];
	if (m < n) {
		*p = 0, *s = m;
		q[0] = 0;
	}
	else {
		*p = m - n, *s = n - 1;
		for (k = *p; k >= 0; k--) {
			nk = n + k;
			q[k] = r[nk] * pow(vn, k);
			for (j = nk - 1; j >= 0; j--) {
				jk = j - k;
				if (jk >= 0)
					r[j] = vn * r[j] - r[nk] * v[j - k];
				else
					r[j] = vn * r[j];
			}
		}
		while (*p > 0 && q[*p] == 0) *p = *p - 1;
		while (*s > 0 && r[*s] == 0) *s = *s - 1;
	}
}

void poly_exp_mod(long degreeA, long degreem, long n,
	long modulus, long *A, long *m, long *s,
	long *ds)
{
	int zero;
	long dp, dq, dx = degreeA, i;
	long p[SIZE], q[SIZE], x[SIZE];

	*ds = 0, s[0] = 1;
	for (i = 0; i <= dx; i++) x[i] = A[i];
	while (n > 0) {
		if ((n & 1) == 1) {
			/* s = (s * x) % m; */
			poly_mul(*ds, dx, s, x, p, &dp);
			poly_div(dp, degreem, p, m, q, s, &dq, ds);
			for (i = 0; i <= *ds; i++) s[i] %= modulus;
			zero = s[*ds] == 0, i = *ds;
			while (i > 0 && zero) {
				if (zero) *ds = *ds - 1;
				zero = s[--i] == 0;
			}
		}
		n >>= 1;
		/* x = (x * x) % m; */
		poly_mul(dx, dx, x, x, p, &dp);
		poly_div(dp, degreem, p, m, q, x, &dq, &dx);
		for (i = 0; i <= dx; i++) x[i] %= modulus;
		zero = x[dx] == 0, i = dx;
		while (i > 0 && zero) {
			if (zero) dx--;
			zero = x[--i] == 0;
		}
	}
}

void poly_copy(long db, long *a, long *b, long *da)
/* a = b */
{
	long i;

	*da = db;
	for (i = 0; i <= db; i++) a[i] = b[i];
}

int poly_Extended_Euclidean(long db, long dn,
	long *b, long *n,
	long *t, long *dt)
{
	int nonzero;
	long db0, dn0, dq, dr, dt0 = 0, dtemp, du, i;
	long b0[SIZE], n0[SIZE], q[SIZE];
	long r[SIZE], t0[SIZE], temp[SIZE], u[SIZE];

	*dt = 0;
	poly_copy(dn, n0, n, &dn0);
	poly_copy(db, b0, b, &db0);
	t0[0] = 0;
	t[0] = 1;
	poly_div(dn0, db0, n0, b0, q, r, &dq, &dr);
	nonzero = r[0] != 0;
	for (i = 1; !nonzero && i <= dr; i++)
		nonzero = r[i] != 0;
	while (nonzero) {
		poly_mul(dq, *dt, q, t, u, &du);
		if (dt0 < du)
			for (i = dt0 + 1; i <= du; i++) t0[i] = 0;
		for (i = 0; i <= du; i++)
			temp[i] = t0[i] - u[i];
		dtemp = du;
		poly_copy(*dt, t0, t, &dt0);
		poly_copy(dtemp, t, temp, dt);
		poly_copy(db0, n0, b0, &dn0);
		poly_copy(dr, b0, r, &db0);
		poly_div(dn0, db0, n0, b0, q, r, &dq, &dr);
		nonzero = r[0] != 0;
		for (i = 1; !nonzero && i <= dr; i++)
			nonzero = r[i] != 0;
	}
	if (db0 != 0 && b0[0] != 1) return 0;
	return 1;
}

void poly_mod(long da, long p, long *a, long *new_da)
{
	int zero;
	long i;

	for (i = 0; i <= da; i++) {
		a[i] %= p;
		if (a[i] < 0) a[i] += p;
	}
	zero = a[da] == 0;
	for (i = da - 1; zero && i >= 0; i--) {
		da--;
		zero = a[i] == 0;
	}
	*new_da = da;
}

void poly_write(char *label, long da, long *a)
{
	long i;

	printf("%s", label);
	for (i = da; i >= 0; i--)
		printf("%ld ", a[i]);
	printf("\n");
}

int main(void)
{
	long da = 4, db = 3, dc = 3, dd = 5;
	long de, dq, dr, ds, p = 2;
	long a[5] = { 0, 0, 1, 0, 1 };
	long b[4] = { 1, 1, 0, 1 };
	long c[4] = { 0, 0, 1, 1 };
	long d[6] = { 1, 0, 1, 0, 0, 1 };
	long e[SIZE], q[SIZE], r[SIZE], s[SIZE];

	poly_write("A = ", da, a);
	poly_write("B = ", db, b);
	poly_write("C = ", dc, c);
	poly_write("D = ", dd, d);
	poly_mul(da, db, a, b, e, &de);
	poly_mod(de, p, e, &de);
	poly_div(de, dd, e, d, q, r, &dq, &dr);
	poly_mod(dq, p, q, &dq);
	poly_mod(dr, p, r, &dr);
	poly_write("A * B mod 2 = ", de, e);
	poly_write("A * B / D mod 2 = ", dq, q);
	poly_write("A * B mod D, 2  = ", dr, r);
	if (!poly_Extended_Euclidean(dc, dd, c, d, e, &de))
		printf("*error*\in poly_Extended_Euclidean\n");
	poly_mod(de, p, e, &de);
	poly_write("C ^ - 1 mod D = ", de, e);
	poly_mul(dc, de, c, e, s, &ds);
	poly_mod(ds, p, s, &ds);
	poly_div(ds, dd, s, d, q, r, &dq, &dr);
	poly_mod(dr, p, r, &dr);
	poly_write("C * C ^ - 1 mod D, 2 = ", dr, r);
	dq = 1;
	q[0] = 0;
	q[1] = 1;
	poly_exp_mod(dq, dd, 7, p, q, d, s, &ds);
	poly_mod(ds, p, s, &ds);
	poly_write("x ^  7 mod D, 2 = ", ds, s);
	poly_exp_mod(dq, dd, 9, p, q, d, s, &ds);
	poly_mod(ds, p, s, &ds);
	poly_write("x ^  9 mod D, 2 = ", ds, s);
	poly_exp_mod(dq, dd, 10, p, q, d, s, &ds);
	poly_mod(ds, p, s, &ds);
	poly_write("x ^ 10 mod D, 2 = ", ds, s);
	poly_exp_mod(dq, dd, 25, p, q, d, s, &ds);
	poly_mod(ds, p, s, &ds);
	poly_write("x ^ 25 mod D, 2 = ", ds, s);
	return 0;
}

More Elliptic Curve Point Counting Results Using the Shanks-Mestre Algorithm by James Pate Williams, Jr.

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.

Elliptic Curve Point Counting Using the Shanks-Mestre Algorithm in C# and C++ by James Pate Williams, Jr.

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#.

Artificial Neural Network Learns the Exclusive OR (XOR) Function by James Pate Williams, Jr.

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";
            }
        }
    }
}

Comparison of Euler Method, Runge-Kutta 2, Runge-Kutta 4, and Taylor Series Solutions of a Single Ordinary Differential Equation

#include "RungeKutta4.h"
#include <iomanip>
#include <iostream>
using namespace std;

int main()
{
	int nSteps;
	vector<Point> solution1, solution2;
	vector<Point> solution3, solution4;

	while (true)
	{
		cout << "n-steps (zero to exit app) = ";
		cin >> nSteps;

		if (nSteps == 0)
			break;

		cout << "n\t\txn\tyn2p(xn)\tyn4p(xn)\tEuler\t\tTaylor\t\tDFDX2\t\tDFDX4\r\n";
		RungeKutta4::ComputeEuler(nSteps, 1, 2, -1, solution1);
		RungeKutta4::ComputeRK2(nSteps, 1, 2, -1, solution2);
		RungeKutta4::ComputeRK4(nSteps, 1, 2, -1, solution3);
		RungeKutta4::ComputeTaylor(nSteps, 1, 2, -1, solution4);

		for (int i = 0; i <= nSteps; i++)
		{
			POINT pt1 = solution1[i];
			POINT pt2 = solution2[i];
			POINT pt3 = solution3[i];
			POINT pt4 = solution4[i];
			cout << setw(2) << pt1.n << "\t";
			cout << setprecision(8) << fixed << pt1.x << "\t";
			cout << setprecision(8) << fixed << pt3.y << "\t";
			cout << setprecision(8) << fixed << pt4.y << "\t";
			cout << setprecision(8) << fixed << pt1.y << "\t";
			cout << setprecision(8) << fixed << pt2.y << "\t";
			cout << setprecision(8) << fixed << pt2.derivative << "\t";
			cout << setprecision(8) << fixed << pt3.derivative << "\r\n";
		}
	}
}
#pragma once
#include <vector>
using namespace std;

typedef struct Point
{
	int n;
	long double derivative, x, y;
} POINT, *PPOINT;

class RungeKutta4
{
private:
	static long double FX(
		long double x,
		long double y);
	static long double DFDX(
		long double x,
		long double y,
		long double yp);
public:
	static void ComputeEuler(
		int nSteps,
		long double x0,
		long double x1,
		long double y0,
		vector<Point>& pts);
	static void ComputeRK2(
		int nSteps,
		long double x0,
		long double x1,
		long double y0,
		vector<Point>& pts);
	static void ComputeRK4(
		int nSteps,
		long double x0,
		long double x1,
		long double y0,
		vector<Point>& pts);
	static void ComputeTaylor(
		int nSteps,
		long double x0,
		long double x1,
		long double y0,
		vector<Point>& pts);
};
// 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 20, 2023

#include "RungeKutta4.h"

long double RungeKutta4::FX(
	long double x, long double y)
{
	long double x2 = x * x;
	long double term1 = 1.0 / x2;
	long double term2 = -y / x;
	long double term3 = -y * y;
	return term1 + term2 + term3;
}

long double RungeKutta4::DFDX(
	long double x,
	long double y,
	long double yp)
{
	long double x3 = x * x * x;
	long double term1 = -2.0 / x3;
	long double term2 = -yp / x;
	long double term3 = y / x / x;
	long double term4 = -2.0 * y * yp;
	return term1 + term2 + term3 + term4;
}

void RungeKutta4::ComputeEuler(
	int nSteps,
	long double x0,
	long double x1,
	long double y0,
	vector<Point>& pts)
{
	Point pt0 = {};
	pt0.n = 0;
	pt0.x = x0;
	pt0.y = y0;
	double derivative = FX(x0, y0);
	pt0.derivative = derivative;
	pts.push_back(pt0);

	if (nSteps == 1)
		return;

	long double h = (x1 - x0) / nSteps;
	long double xn = x0;
	long double yn = y0;

	for (int n = 1; n <= nSteps; n++)
	{
		xn = x0 + n * h;
		long double yn = y0 + h * FX(xn, y0);
		derivative = FX(xn, y0);
		pt0.n = n;
		pt0.derivative = derivative;
		pt0.x = xn;
		pt0.y = yn;
		pts.push_back(pt0);
		y0 = yn;
	}
}

void RungeKutta4::ComputeRK2(
	int nSteps,
	long double x0,
	long double x1,
	long double y0,
	vector<Point>& pts)
{
	Point pt0 = {};
	pt0.n = 0;
	pt0.x = x0;
	pt0.y = y0;
	double derivative = FX(x0, y0);
	pt0.derivative = derivative;
	pts.push_back(pt0);

	if (nSteps == 1)
		return;

	long double h = (x1 - x0) / nSteps;
	long double xn = x0;
	long double yn = y0;

	for (int n = 1; n <= nSteps; n++)
	{
		long double k1 = h * FX(xn, yn);
		long double k2 = h * FX(xn + h, yn + k1);
		yn += 0.5 * (k1 + k2);
		xn = x0 + n * h;
		derivative = FX(xn, yn);
		pt0.n = n;
		pt0.derivative = derivative;
		pt0.x = xn;
		pt0.y = yn;
		pts.push_back(pt0);
	}
}

void RungeKutta4::ComputeRK4(
	int nSteps,
	long double x0,
	long double x1,
	long double y0,
	vector<Point>& pts)
{
	Point pt0 = {};
	pt0.n = 0;
	pt0.x = x0;
	pt0.y = y0;
	double derivative = FX(x0, y0);
	pt0.derivative = derivative;
	pts.push_back(pt0);

	if (nSteps == 1)
		return;

	long double h = (x1 - x0) / nSteps;
	long double xn = x0;
	long double yn = y0;

	for (int n = 1; n <= nSteps; n++)
	{
		long double k1 = h * FX(xn, yn);
		long double k2 = h * FX(xn + h / 2, yn + k1 / 2);
		long double k3 = h * FX(xn + h / 2, yn + k2 / 2);
		long double k4 = h * FX(xn + h, yn + k3);
		xn = x0 + n * h;
		yn = yn + (k1 + k2 + k3 + k4) / 6.0;
		derivative = FX(xn, yn);
		pt0.n = n;
		pt0.derivative = derivative;
		pt0.x = xn;
		pt0.y = yn;
		pts.push_back(pt0);
	}
}

void RungeKutta4::ComputeTaylor(
	int nSteps,
	long double x0,
	long double x1,
	long double y0,
	vector<Point>& pts)
{
	Point pt0 = {};
	pt0.n = 0;
	pt0.x = x0;
	pt0.y = y0;
	double derivative = FX(x0, y0);
	pt0.derivative = derivative;
	pts.push_back(pt0);

	if (nSteps == 1)
		return;

	long double h = (x1 - x0) / nSteps;
	long double xn = x0;
	long double yn = y0;
	long double yp = FX(x0, y0);

	for (int n = 1; n <= nSteps; n++)
	{
		xn = x0 + n * h;
		yp = FX(xn, y0);
		long double yn = y0 + h * (FX(xn, y0) + h * DFDX(xn, y0, yp) / 2);
		derivative = FX(xn, y0);
		pt0.n = n;
		pt0.derivative = derivative;
		pt0.x = xn;
		pt0.y = yn;
		pts.push_back(pt0);
		y0 = yn;
	}
}