Category: Computer Science

Multiple Integration Using Simpson’s Rule – Calculation of the Ground State of the Non-Relativistic Helium Atom by James Pate Williams, Jr.

The six-dimensional Cartesian coordinate wavefunction is calculated by the C# method:

public double Psi(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double exp1 = Math.Exp(-alpha[0] * r1);
            double exp2 = Math.Exp(-alpha[1] * r2);
            double exp3 = Math.Exp(-alpha[2] * r12);

            return exp1 * exp2 * exp3;
        }

        public double Psi2(double[] x, double[] alpha)
        {
            double psi = Psi(x, alpha);

            return psi * psi;
        }

The wavefunction normalization method is:

public double Normalize(double[] alpha, int nSteps)
        {
            double[] lower = new double[6];
            double[] upper = new double[6];
            int[] steps = new int[6];

            lower[0] = 0.001;
            lower[1] = 0.001;
            lower[2] = 0.001;
            lower[3] = 0.001;
            lower[4] = 0.001;
            lower[5] = 0.001;

            upper[0] = 10.0;
            upper[1] = 10.0;
            upper[2] = 10.0;
            upper[3] = 10.0;
            upper[4] = 10.0;
            upper[5] = 10.0;

            for (int i = 0; i < 6; i++)
                steps[i] = nSteps;

            double norm = Math.Sqrt(integ.Integrate(
                lower, upper, alpha, Psi2, 6, steps));

            return norm;
        }

The two kinetic energy integrands are encapsulated in the following two methods:

private double Integrand1(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double term = -2.0 * alpha[0] / r1 + alpha[0] * alpha[0];
            double mul1 = Math.Exp(-2.0 * alpha[0] * r1);
            double mul2 = Math.Exp(-2.0 * alpha[1] * r2);
            double mul3 = Math.Exp(-2.0 * alpha[2] * r12);

            return N * N * term * mul1 * mul2 * mul3;
        }

        private double Integrand2(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double term = -2.0 * alpha[1] / r1 + alpha[1] * alpha[1];
            double mul1 = Math.Exp(-2.0 * alpha[0] * r1);
            double mul2 = Math.Exp(-2.0 * alpha[1] * r2);
            double mul3 = Math.Exp(-2.0 * alpha[2] * r12);

            return N * N * term * mul1 * mul2 * mul3;
        }

The two potential energy terms use the following two integrands:

private double Integrand3(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double term = 1.0 / r1;
            double mul =
                Math.Exp(-2.0 * alpha[0] * r1) *
                Math.Exp(-2.0 * alpha[1] * r2) *
                Math.Exp(-2.0 * alpha[2] * r12);

            return -N * N * Z * term * mul;
        }

        private double Integrand4(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));
            double term = 1.0 / r2;
            double mul =
                Math.Exp(-2.0 * alpha[0] * r1) *
                Math.Exp(-2.0 * alpha[1] * r2) *
                Math.Exp(-2.0 * alpha[2] * r12);

            return -N * N * Z * term * mul;
        }

The electron-electron repulsion integrand is given by:

private double Integrand5(double[] x, double[] alpha)
        {
            double r1 = Math.Sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
            double r2 = Math.Sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5]);
            double r12 = Math.Sqrt(Math.Pow(x[0] - x[3], 2.0) +
                Math.Pow(x[1] - x[4], 2.0) + Math.Pow(x[2] - x[5], 2.0));

            if (r12 == 0)
                r12 = 0.01;

            double term = 1.0 / r12;
            double mul =
                Math.Exp(-2.0 * alpha[0] * r1) *
                Math.Exp(-2.0 * alpha[1] * r2) *
                Math.Exp(-2.0 * alpha[2] * r12);

            return N * N * term * mul;
        }

The ground state non-relativistic energy is computed by the method:

public double Energy(double[] alpha, int nSteps, int Z)
        {
            double[] lower = new double[6];
            double[] upper = new double[6];
            int[] steps = new int[6];

            lower[0] = lower[1] = lower[2] = lower[3] = lower[4] = lower[5] = 0.001;
            upper[0] = upper[1] = upper[2] = upper[3] = upper[4] = upper[5] = 10.0;
            steps[0] = steps[1] = steps[2] = steps[3] = steps[4] = steps[5] = nSteps;

            N = Normalize(alpha, nSteps);

            this.Z = Z;

            double integ1 = integ.Integrate(lower, upper, alpha, Integrand1, 6, steps);
            double integ2 = integ.Integrate(lower, upper, alpha, Integrand2, 6, steps);
            double integ3 = integ.Integrate(lower, upper, alpha, Integrand3, 6, steps);
            double integ4 = integ.Integrate(lower, upper, alpha, Integrand4, 6, steps);
            double integ5 = integ.Integrate(lower, upper, alpha, Integrand5, 6, steps);
            
            return integ1 + integ2 - integ3 - integ4 + integ5;
        }

Using trial and error we calculate Alpha as: 0.535139999999

public double Integrate(double[] lower, double[] upper, double[] alpha,
            Func<double[], double[], double> f, int n, int[] steps)
        {
            double p = 1;
            double[] h = new double[n];
            double[] h2 = new double[n];
            double[] s = new double[n];
            double[] t = new double[n];
            double[] x = new double[n];
            double[] w = new double[n];

            for (int i = 0; i < n; i++)
            {
                h[i] = (upper[i] - lower[i]) / steps[i];
                h2[i] = 2.0 * h[i];
            }

            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < n; j++)
                    x[j] = lower[j];
                
                x[i] = lower[i] + h[i];

                for (int j = 1; j < steps[i]; j += 2)
                {
                    s[i] += f(x, alpha);
                    x[i] += h2[i];
                }

                for (int j = 0; j < n; j++)
                    x[j] = lower[j];

                x[i] = lower[i] + h2[i];

                for (int j = 2; j < steps[i]; j += 2)
                {
                    t[i] += f(x, alpha);
                    x[i] += h2[i];
                }

                w[i] = h[i] * (f(lower, alpha) + 4 * s[i] + 2 * t[i] + f(upper, alpha)) / 3.0;
            }

            for (int i = 0; i < n; i++)
                p *= w[i];

            return p;
        }

C++ Searching: Binary and Linear Implemented by James Pate Williams, Jr.

#pragma once

// Binary searches from the following website
// https://www.geeksforgeeks.org/binary-search/

typedef long long ll;

class Search
{
public:
	int IterativeBinarySearch(
		ll A[], ll x, int low, int high);
	int RecursiveBinarySearch(
		ll A[], ll x, int low, int high);
	int IterativeLinearSearch(
		ll A[], ll x, int low, int high);
};

#include "Search.h"

int Search::IterativeBinarySearch(
	ll A[], ll x, int low, int high)
{
	do
	{
		int middle = low + (high - low) / 2;

		if (x == A[middle])
			return middle;

		else if (x > A[middle])
			low = middle + 1;

		else if (x < A[middle])
			high = middle - 1;

	} while (high >= low);

	return -1;
}

int Search::RecursiveBinarySearch(
	ll A[], ll x, int low, int high)
{
	if (high >= low)
	{
		int middle = low + (high - low) / 2;

		if (x == A[middle])
			return middle;

		else if (x > A[middle])
			return RecursiveBinarySearch(
				A, x, middle + 1, high);

		else if (x < A[middle])
			return RecursiveBinarySearch(
				A, x, low, middle - 1);
	}
	
	return -1;
}

int Search::IterativeLinearSearch(
	ll A[], ll x, int low, int high)
{
	for (int i = low; i <= high; i++)
		if (x == A[i])
			return i;

	return -1;
}

// Searching.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include "Search.h"
#include <chrono>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;
using chrono::duration_cast;
using chrono::nanoseconds;

ll A[1000001];

ll RandomLongLong()
{
	ll lo = rand();
	ll hi = rand();

	return (lo << 31) | hi;
}

void GenerateArray(ll A[], int n, int seed)
{
	for (int i = 1; i <= n; i++)
		A[i] = RandomLongLong();
}

void GenerateArrayMod(ll A[], int n, int mod, int seed)
{
	Search search;

	for (int i = 1; i <= n; i++)
	{
		while (true)
		{
			ll Ai = RandomLongLong() % mod;

			if (search.IterativeLinearSearch(
				A, Ai, 1, i) == -1)
			{
				A[i] = Ai;
				break;
			}
		}
	}
}

int Partition(ll* A, int p, int r)
{
	int q = p;
	ll t;

	for (int u = p; u <= r - 1; u++)
	{
		if (A[u] <= A[r])
		{
			t = A[q];
			A[q] = A[u];
			A[u] = t;
			q++;
		}
	}

	t = A[q];
	A[q] = A[r];
	A[r] = t;

	return q;
}

void RunQuickSort(ll* A, int p, int r)
{
	if (p < r)
	{
		int q = Partition(A, p, r);

		RunQuickSort(A, p, q - 1);
		RunQuickSort(A, q + 1, r);
	}
}

void RunSearches(
	Search search,
	ll A[],
	int n, int seed)
{
	int index[4];
	ll x = RandomLongLong() % n;

	GenerateArrayMod(A, n, 2LL * n, seed);
	RunQuickSort(A, 1, n);
	
	cout << "key = " << x << endl;

	if (n <= 25)
	{		for (int i = 1; i <= n; i++)
			cout << A[i] << " ";

		cout << endl;
	}

	cout << "Runtimes in nanoseconds" << endl;
	cout << "n" << "\t" << "Iter" << "\t";
	cout << "Recur" << "\t" << "Linear" << endl;
	cout << setw(5) << n << "\t";

	auto start1 = chrono::steady_clock::now();
	index[0] = search.IterativeBinarySearch(A, x, 1, n);
	auto final1 = chrono::steady_clock::now();

	cout << setw(5) << duration_cast<nanoseconds>(final1 - start1).count();
	cout << "\t";

	auto start2 = chrono::steady_clock::now();
	index[1] = search.RecursiveBinarySearch(A, x, 1, n);
	auto final2 = chrono::steady_clock::now();

	cout << setw(5) << duration_cast<nanoseconds>(final2 - start2).count();
	cout << "\t";

	auto start3 = chrono::steady_clock::now();
	index[2] = search.IterativeLinearSearch(A, x, 1, n);
	auto final3 = chrono::steady_clock::now();

	cout << setw(5) << duration_cast<nanoseconds>(final3 - start3).count();
	cout << endl;
	cout << "Search indicies" << endl;

	for (int i = 0; i < 3; i++)
		cout << "index[" << i << "] = " << index[i] << endl;
}

int main()
{
	cout << "Comparison of Iterative Binary Search," << endl;
	cout << "Recursive Binary Search, and Linear Search" << endl;

	while (true)
	{
		int n, seed;
		Search search;

		cout << "n = ";
		cin >> n;
		
		if (n == 0)
			break;

		cout << "seed = ";
		cin >> seed;
		srand(seed);

		RunSearches(search, A, n, seed);
	}
}

C++ std::sort Comparison by James Pate Williams, Jr.

We compare selection sort, merge sort, quick sort, and std::sort. The runtime orders were quick sort < merge sort < std::sort < selection sort.

#pragma once
class Sort
{
public:
	void RunSelectionSort(long long* A, int p, int r);
	void RunMergeSort(
		long long* A,
		long long* B,
		long long* C,
		int p, int r);
	void RunQuickSort(long long* A, int p, int r);
private:
	void Merge(
		long long* A,
		long long* B,
		long long* C,
		int p, int q, int r);
	int Partition(long long* A, int lo, int hi);
};

#include <limits.h>
#include "Sort.h"

void Sort::Merge(
    long long* A,
    long long* B,
    long long* C,
    int p, int q, int r)
{
    int n1 = q - p + 1, n2 = r - q;

    for (int i = p, j = 1; i <= q; i++, j++)
        B[j] = A[i];

    for (int i = q + 1, j = 1; i <= r; i++, j++)
        C[j] = A[i];

    B[n1 + 1] = C[n2 + 1] = LLONG_MAX;

    int ii = 1, jj = 1;

    for (int k = p; k <= r; k++)
    {
        if (B[ii] <= C[jj])
        {
            A[k] = B[ii];
            ii++;
        }

        else
        {
            if (B[ii] > C[jj])
            {
                A[k] = C[jj];
                jj++;
            }
        }
    }
}

int Sort::Partition(long long* A, int p, int r)
{
    int q = p;
    long long t;

    for (int u = p; u <= r - 1; u++)
    {
        if (A[u] <= A[r])
        {
            t = A[q];
            A[q]= A[u];
            A[u] = t;
            q++;
        }
    }

    t = A[q];
    A[q] = A[r];
    A[r] = t;

    return q;
}

void Sort::RunMergeSort(
    long long* A,
    long long* B,
    long long* C,
    int p, int r)
{
    int q = (p + r) / 2;

    if (p < r)
    {
        RunMergeSort(A, B, C, p, q);
        RunMergeSort(A, B, C, q + 1, r);
        Merge(A, B, C, p, q, r);
    }
}

void Sort::RunQuickSort(long long* A, int p, int r)
{
    if (p < r)
    {
        int q = Partition(A, p, r);

        RunQuickSort(A, p, q - 1);
        RunQuickSort(A, q + 1, r);
    }
}

void Sort::RunSelectionSort(long long* A, int p, int r)
{
    for (int i = p; i <= r - 1; i++)
    {
        for (int j = i + 1; j <= r; j++)
        {
            if (A[i] > A[j])
            {
                long long t = A[i];

                A[i] = A[j];
                A[j] = t;
            }
        }
    }
}

#include "Sort.h"
#include <stdlib.h>
#include <algorithm>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;
using chrono::duration_cast;
using chrono::milliseconds;

long long A[16], B[18], C[18];
long long AA[100002], BB[100002], CC[100002];

long long RandomLongLong()
{
	long long lo = rand();
	long long hi = rand();

	return (lo << 31) | hi;
}

void GenerateArray(long long A[], int n, int seed)
{
	srand(seed);

	for (int i = 1; i <= n; i++)
		A[i] = RandomLongLong();
}

void GenerateArrayMod(long long A[], int n, int mod, int seed)
{
	srand(seed);

	for (int i = 1; i <= n; i++)
		A[i] = RandomLongLong() % mod;
}

void RunSorts(
	Sort sort,
	long long A[],
	long long B[],
	long long C[],
	int n, int seed)
{
	cout << setw(5) << n << "\t";

	auto start1 = chrono::steady_clock::now();
	sort.RunSelectionSort(A, 1, n);
	auto final1 = chrono::steady_clock::now();

	cout << setw(5) << duration_cast<milliseconds>(final1 - start1).count();
	cout << "\t";

	GenerateArray(A, n, seed);

	auto start2 = chrono::steady_clock::now();
	sort.RunMergeSort(A, B, C, 1, n);
	auto final2 = chrono::steady_clock::now();

	cout << setw(5) << duration_cast<milliseconds>(final2 - start2).count();
	cout << "\t";

	GenerateArray(A, n, seed);

	auto start3 = chrono::steady_clock::now();
	sort.RunQuickSort(A, 1, n);
	auto final3 = chrono::steady_clock::now();

	cout << setw(5) << duration_cast<milliseconds>(final3 - start3).count();
	cout << "\t";

	GenerateArray(A, n, seed);
	vector<long long> V;

	for (int i = 0; i < n; i++)
		V.push_back(A[i + 1]);

	auto start4 = chrono::steady_clock::now();
	std::sort(V.begin(), V.end());
	auto final4 = chrono::steady_clock::now();

	cout << setw(5) << duration_cast<milliseconds>(final4 - start4).count();
	cout << endl;
}

void TestSorts(
	Sort sort,
	long long A[],
	long long B[],
	long long C[],
	int n, int seed)
{
	cout << "Testing Sorting Algorthms" << endl;
	cout << "Selection Sort" << endl;

	GenerateArrayMod(A, n, 100000, 1);

	for (int i = 1; i <= n; i++)
		cout << setw(5) << A[i] << " ";

	cout << endl;

	sort.RunSelectionSort(A, 1, n);

	for (int i = 1; i <= n; i++)
		cout << setw(5) << A[i] << " ";

	cout << endl << endl;
	cout << "Merge Sort" << endl;

	GenerateArrayMod(A, n, 100000, 1);

	for (int i = 1; i <= n; i++)
		cout << setw(5) << A[i] << " ";

	cout << endl;

	sort.RunMergeSort(A, B, C, 1, n);

	for (int i = 1; i <= n; i++)
		cout << setw(5) << A[i] << " ";

	cout << endl << endl;
	cout << "Quick Sort" << endl;

	GenerateArrayMod(A, n, 100000, 1);

	for (int i = 1; i <= n; i++)
		cout << setw(5) << A[i] << " ";

	cout << endl;

	sort.RunQuickSort(A, 1, n);

	for (int i = 1; i <= n; i++)
		cout << setw(5) << A[i] << " ";

	cout << endl << endl;
	cout << "std::sort" << endl;

	GenerateArrayMod(A, n, 100000, 1);
	vector<long long> V;

	for (int i = 1; i <= n; i++)
	{
		cout << setw(5) << A[i] << " ";
		V.push_back(A[i]);
	}

	cout << endl;

	std::sort(V.begin(), V.end());

	for (size_t i = 0; i < V.size(); i++)
		cout << setw(5) << V[i] << " ";

	cout << endl << endl;
}

int main(int argc, char** argv)
{
	Sort sort;
	
	TestSorts(sort, A, B, C, 15, 1);

	cout << "Runtimes in Milliseconds" << endl;

	cout << "    n" << "\tSelect" << "\t Merge" << "\t Quick\tstd::sort" << endl;

	for (int n = 10000; n <= 60000; n += 10000)
	{
		RunSorts(sort, AA, BB, CC, n, 1);
	}
	
	return 0;
}

C++ Low Density Subset Sum Problem Solver Application Using the L3 Lattice Basis Reduction Algorithm Implemented by James Pate Williams, Jr. Utilizing the “Handbook of Applied Cryptography” Edited by Alfred J. Menezes Among Others

See the website https://cacr.uwaterloo.ca/hac/about/chap3.pdf

Atkin’s Primality Test by James Pate Williams, Jr.

Back in 2022 I reimplemented Henri Cohen’s Atkin’s Primality Test algorithm. This test makes use of an elliptic curve analog of Pocklington’s theorem. I restate the theorem utilized from Henri Cohen’s “A course in Computational Algebraic Number Theory” on pages 467 to 468: “Proposition 9.2.1. Let N be an integer coprime to 6 and different from 1, and E be an elliptic curve modulo N. Assume that we know an integer m a point P contained on the elliptic curve satisfying the following conditions. (1) There exists a prime divisor q of m such that q  >  (N^1/4 + 1) ^ 2 (2) m * P = O_E = (0 : 1 : ). (3) (m / q) * P = (x : y : t) with t contained in (Z/NZ)*. Then N is prime.” I used C# and Microsoft’s BigInteger class. I have not been able to prove numbers greater than 14 decimal digits to be prime. I am recoding the algorithm in C++ which limits me to 19 decimal digits since 2 ^ 63 – 1 = 9,223,372,036,854,775,807 (Int64).

Henri Cohen’s Brent’s Modification of the Pollard rho Factorization Method Implemented by James Pate Williams, Jr.

I implemented the algorithm in C using unsigned 64-bit integers. This method is good for integers of around 18 decimal digits in length. For comparison I tested against my blazingly fast Shor’s classical factoring algorithm which works on arbitrarily large integers of around 50 to 60 or more decimal digits.

Henri Cohen’s Trial Division Algorithm Implemented by James Pate Williams, Jr.

I implemented this algorithm for C long (32-bit) number in 1997, C# big integers in approximately 2015 or 2016, and C unsigned long long (64-bit) integers in March, 2023.

/*
Author:  Pate Williams (c) 1997 - 2023

"Algorithm 8.1.1 (Trial Division). We assume given
a table of prime numbers p[1] = 2, p[2] = 3,..., p[k],
with k > 3, an array t <- [6, 4, 2, 4, 4, 6, 2], and
an index j such that if p[k] mod 30 is equal to
1, 7, 11, 13, 17, 19, 23 or 29 then j is equal to
0, 1, 2, 3, 4, 5, 6 or 7 respectively. Finally, we
give ourselves an upper bound B such that B >= p[k],
essentially to avoid spending too much time. Then
given a positive integer N, this algorithm tries to
factor (or split N) and if it fails, N will be free
of prime factors less than or equal to B."
-Henri Cohen-
See "A Course in Computational Algebraic Number
Theory" by Henri Cohen page 420.
*/

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

#define BITS_PER_LONG   32
#define BITS_PER_LONG_1 31
#define MAX_SIEVE 100000000L
#define LARGEST_PRIME 99999989L
#define SIEVE_SIZE (MAX_SIEVE / BITS_PER_LONG)

typedef unsigned long long ull;

struct factor { int expon;  ull prime; };
long *prime, sieve[SIEVE_SIZE];

int get_bit(long i, long *sieve) {
	long b = i % BITS_PER_LONG;
	long c = i / BITS_PER_LONG;

	return (sieve[c] >> (BITS_PER_LONG_1 - b)) & 1;
}

void set_bit(long i, long v, long *sieve) {
	long b = i % BITS_PER_LONG;
	long c = i / BITS_PER_LONG;
	long mask = 1 << (BITS_PER_LONG_1 - b);

	if (v == 1)
		sieve[c] |= mask;
	else
		sieve[c] &= ~mask;
}

void Sieve(long n, long *sieve) {
	int c, i, inc;

	set_bit(0, 0, sieve);
	set_bit(1, 0, sieve);
	set_bit(2, 1, sieve);
	for (i = 3; i <= n; i++)
		set_bit(i, i & 1, sieve);
	c = 3;
	do {
		i = c * c, inc = c + c;
		while (i <= n) {
			set_bit(i, 0, sieve);
			i += inc;
		}
		c += 2;
		while (!get_bit(c, sieve)) c++;
	} while (c * c <= n);
}

void get_bits(ull b, int *number, int *bits) {
	int i = 0;
	do {
		bits[i++] = b % 2;
		b = b / 2;
	} while (b > 0);
	*number = i;
}

 ull pow_mod(
	 ull x,
	 ull b, 
	 ull n) {
	// Figure 4.4 in "Cryptography Theory and Practice"
	// by Douglas R. Stinson page 127
	int bits[64], i, l;
	ull z = 1, s = 0;

	get_bits(b, &l, bits);
	s = bits[l - 1];
	for (i = l - 2; i >= 0; i--)
		s = 2 * s + bits[i];
	if (b != s) {
		printf("Error in pow_mod function\n");
		exit(-1);
	}
	for (i = l - 1; i >= 0; i--) {
		z = (z * z) % n;
		if (bits[i] == 1)
			z = (z * x) % n;
	}
	return z;
}

ull rand_ull(void) {
	ull r = 0;
	for (int i = 0; i < 64; i += 15 /*30*/) {
		r = r * ((ull)RAND_MAX + 1) + rand();
	}
	return r;
}

int Miller_Rabin(ull n, int t) {
	// returns 1 for prime
	// returns 0 for composite
	// "Handbook of Applied Cryptography"
	// Alfred J. Menezes among others
	// 4.24 Algorithm page 139
	int i, j, s = 0;
	ull a, n1 = n - 1, r, y;

	r = n1;
	while ((r % 2) == 0) {
		s++;
		r /= 2;
	}
	for (i = 0; i < t; i++)
	{
		while (true)
		{
			a = rand_ull() % n;
			
			if (a >= 2 && a <= n - 2)
				break;
		}
		y = pow_mod(a, r, n);
		if (y != 1 && y != n1) {
			j = 1;
			while (j <= s - 1 && y != n1) {
				y = (y * y) % n;
				if (y == 1)
					return 0;
				j++;
			}
			if (y == n1)
				return 0;
		}
	}
	return 1;
}

int trial_division(ull *N, long *prime,
	struct factor *f, int *count) {
	int found;
	ull B, d;
	int e, i, j = 0, k, m, n;
	int t[8] = { 6, 4, 2, 4, 2, 4, 6, 2 };
	int table[8] = { 1, 7, 11, 13, 17, 19, 23, 29 };
	ull l, r;

	if (*N <= 5) {
		*count = 1;
		f[0].expon = 1;
		switch (*N) {
		case 1: f[0].prime = 1; break;
		case 2: f[0].prime = 2; break;
		case 3: f[0].prime = 3; break;
		case 4: f[0].expon = 2; f[0].prime = 2; break;
		case 5: f[0].prime = 5; break;
		}
		return 1;
	}
	*count = 0;
	B = LARGEST_PRIME;
	i = -1, l = (ull)sqrt((double)*N), m = -1;
L2:
	m++;
	if (i == MAX_SIEVE) {
		i = j - 1;
		goto L5;
	}
	d = prime[m];
	k = d % 30;
	found = 0;
	for (n = 0; !found && n < 8; n++) {
		found = k == table[n];
		if (found) j = n;
	}
L3:
	r = *N % d;
	if (r == 0) {
		e = 0;
		do {
			e++;
			*N /= d;
		} while (*N % d == 0);
		f[*count].expon = e;
		f[*count].prime = d;
		*count = *count + 1;
		if (*N == 1) return 1;
		goto L2;
	}
	if (d >= l) {
		if (*N > 1) {
			f[*count].expon = 1;
			f[*count].prime = *N;
			*count = *count + 1;
			*N = 1;
			return 1;
		}
	}
	else if (i < 0) goto L2;
L5:
	i = (i + 1) % 8;
	d = d + t[i];
	if (d > B) return 0;
	goto L3;
}

int main(void) {
	char e_buffer[256], n_buffer[256] = { 0 };
	long k, e_length, n_length = 0, valid;
	double time_spent;
	int count;
	long i, p = 2;
	ull largest_prime;
	long prime_count = 0;
	ull N, sqrtN;
	struct factor f[32];
	clock_t begin, end;

	begin = clock();
	Sieve(MAX_SIEVE, sieve);
	for (p = 2; p < MAX_SIEVE; p++) {
		while (!get_bit(p, sieve))
			p++;
		prime_count++;
	}
	prime_count--;
	prime = (long*)malloc(prime_count * sizeof(long));
	i = 0;
	p = 2;
	while (i < prime_count)
	{
		while (!get_bit(p, sieve))
			p++;
		prime[i++] = p++;
	}
	largest_prime = LARGEST_PRIME;
	printf("Largest prime = %I64u\n", largest_prime);
	end = clock();
	time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
	printf("Time to create prime number sieve in seconds %lf\n", time_spent);
	for (;;) {
		printf("N = ");
		scanf_s("%s", e_buffer, 256);
		e_length = strlen(e_buffer);
		n_length = 0;
		valid = 1;
		for (k = 0; valid && k < e_length; k++) {
			if (e_buffer[k] >= '0' &&
				e_buffer[k] <= '9' ||
				e_buffer[k] == ',') {
				if (e_buffer[k] != ',') {
					n_buffer[n_length++] = e_buffer[k];
				}
			}
			else {
				valid = 0;
			}
		}
		if (!valid) {
			printf("Invalid character must in set {0, 1, ..., 9, ',')\n");
			continue;
		}
		else {
			n_buffer[n_length] = '\0';
			N = atoll(n_buffer);
			sqrtN = (ull)sqrt((double)N);
			if (sqrtN > largest_prime) {
				printf("Number is too large!\n");
				printf("Square root %I64u must be < %I64u\n", sqrtN, largest_prime);
				continue;
			}
		}
		if (N <= 0) break;
		trial_division(&N, prime, f, &count);
		printf("Non-trivial factors:\n");
		for (i = 0; i < count; i++) {
			printf("%I64u", f[i].prime);
			if (f[i].expon > 1) printf(" ^ %ld\n", f[i].expon);
			else printf("\n");
		}
		if (N != 1)
		{
			if (Miller_Rabin(N, 20) == 1)
				printf("Residue is a prime number\n");
			else
				printf("Last factor is composite\n");
		}
		else
			printf("Number was completely factored\n");
	}
	return 0;
}

Romberg Integration of the Logarithmic Integral by James Pate Williams, Jr.

/*
Author:  Pate Williams (c) January 22, 1995
The following program is a translation of the FORTRAN
subprogram found in Elementary Numerical Analysis by
S. D. Conte and Carl de Boor pages 343-344. The program
uses Romberg extrapolation to find the integral of a
function.
*/

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

typedef long double real;

real t[10][10];

real f(real x)
{
	return(expl(-x * x));
}

real Li(real x)
{
	return(1.0 / logl(x));
}

real romberg(
	real a, real b, int start, int row,
	real(*f)(real))
{
	int i, k, m;
	real h, ratio, sum;

	m = start;
	h = (b - a) / m;
	sum = 0.5 * (f(a) + f(b));
	if (m > 1)
		for (i = 1; i < m; i++)
			sum += f(a + i * h);
	t[1][1] = h * sum;
	if (row < 2) return(t[1][1]);
	for (k = 2; k <= row; k++)
	{
		h = 0.5 * h;
		m *= 2;
		sum = 0.0;
		for (i = 1; i <= m; i += 2)
			sum += f(a + h * i);
		t[k][1] = 0.5 * t[k - 1][1] + sum * h;
		for (i = 1; i < k; i++)
		{
			t[k - 1][i] = t[k][i] - t[k - 1][i];
			t[k][i + 1] = t[k][i] - t[k - 1][i] /
				(powl(4.0, (real)i) - 1.0);
		}
	}
	if (row < 3) return(t[2][2]);
	for (k = 1; k <= row - 2; k++)
		for (i = 1; i <= k; i++)
		{
			if (t[k + 1][i] == 0.0)
				ratio = 0.0;
			else
				ratio = t[k][i] / t[k + 1][i];
			t[k][i] = ratio;
		}
	return(t[row][row - 1]);
}

int main(void)
{
	long double experimental = 0.7468241328124271;
	int row = 8;

	printf("first test\n");
	printf("Romberg integration of f(x) = exp(- x * x)\n");
	printf("Internet value = 0.7468241328124271\n");
	printf("from website = https://www.integral-calculator.com/\n");
	printf("Approximation\t\tPercent Difference\tSteps\tRuntime (s)\n");
	for (long steps = 100000; steps <= 900000; steps += 100000)
	{
		clock_t start = clock();
		long double trial = romberg(0.0, 1.0, steps, row, f);
		clock_t endTm = clock();
		long double runtime = ((long double)endTm - start) /
			CLOCKS_PER_SEC;
		long double pd = 100.0 * fabsl(experimental - trial) /
			(0.5 *(experimental + trial));
		printf("%16.15lf\t%16.15le\t%6d\t%4.3lf\n", 
			trial, pd, steps, runtime);
	}
	printf("second test\n");
	printf("Romberg integration of Li(x) = 1 / ln x between 2 and some n\n");
	printf("approximates the number of prime numbers between 2 and n\n");
	printf("n\t\tLi(n)\t\tSteps\tTime (s)\n");
	for (long n = 10; n <= 1000000000; n *= 10)
	{
		long steps = 1000000;

		clock_t start = clock();
		long double trial = romberg(2.0, n, steps, row, Li);
		clock_t endTm = clock();
		long double runtime = ((long double)endTm - start) /
			CLOCKS_PER_SEC;
		printf("%12ld\t%12.0lf\t%7ld\t%4.3lf\n",
			n, trial, steps, runtime);
	}
	return(0);
}