The curve in this best fit example is the function y(x) = exp(x).
Curve Fitting Using Orthogonal Polynomials in C# by James Pate Williams, Jr.
Derivation of de Broglie’s Wave-Particle Duality Principle of Quantum Mechanics and Macroscopic Exercise by James Pate Williams, Jr.
Euler’s Method and Runge-Kutta 4 Algorithm for Numerically Solving an Ordinary Differential Equation by James Pate Williams, Jr.
I added the Runge-Kutta 4 algorithm found in Numerical Analysis: An Algorithmic Approach Third Edition by S. D. Conte and Carl de Boor. I also added a multistep method, the Adams-Bashforth Method.
Estimated Babe Ruth 1921 Homerun Parameters by James Pate Williams, Jr.
“Babe Ruth is generally considered the owner of the record for the longest home run in MLB history with a 575-foot bomb launched at Navin Field in Detroit in 1921.” – https://www.msn.com/en-us/sports/mlb/what-is-the-longest-home-run-in-mlb-history/ar-AA1dGwlZ
As you can see, I estimated the pitch velocity at 90 miles per hour and Babe Ruth’s (Sultan of Swing) at 90 miles per hour also. My analytic calculations yield a range of the baseball’s trajectory as about 576 feet.
Added the Pseudo-Random Number Generators: Mersenne Twister and Lagged Fibonacci to My C# Application by James Pate Wiliams, Jr. with the Help of the Internet
// "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;
}
}
}
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;
}
Numerical Explorations 