// GauntCoefficients.cpp (c) Monday, November 4, 2024
// by James Pate Williams, Jr., BA, BS, MSWE, PhD
// Computes the Gaunt angular momentum coefficients
// Also the Wigner-3j symbols are calculated
// https://en.wikipedia.org/wiki/3-j_symbol
// https://doc.sagemath.org/html/en/reference/functions/sage/functions/wigner.html#
// https://www.geeksforgeeks.org/factorial-large-number/
#include <iostream>
using namespace std;
typedef long double real;
real pi;
// iterative n-factorial function
real Factorial(int n)
{
real factorial = 1;
for (int i = 2; i <= n; i++)
factorial *= i;
if (n < 0)
factorial = 0;
return factorial;
}
real Delta(int lt, int rt)
{
return lt == rt ? 1.0 : 0.0;
}
real Wigner3j(
int j1, int j2, int j3,
int m1, int m2, int m3)
{
real delta = Delta(m1 + m2 + m3, 0) *
powl(-1.0, j1 - j2 - m3);
real fact1 = Factorial(j1 + j2 - j3);
real fact2 = Factorial(j1 - j2 + j3);
real fact3 = Factorial(-j1 + j2 + j3);
real denom = Factorial(j1 + j2 + j3 + 1);
real numer = delta * sqrt(
fact1 * fact2 * fact3 / denom);
real fact4 = Factorial(j1 - m1);
real fact5 = Factorial(j1 + m1);
real fact6 = Factorial(j2 - m2);
real fact7 = Factorial(j2 + m2);
real fact8 = Factorial(j3 - m3);
real fact9 = Factorial(j3 + m3);
real sqrt1 = sqrtl(
fact4 * fact5 * fact6 * fact7 * fact8 * fact9);
real sumK = 0;
int K = (int)fmaxl(0, fmaxl((real)j2 - j3 - m1,
(real)j1 - j3 + m2));
int N = (int)fminl((real)j1 + j2 - j3,
fminl((real)j1 - m1, (real)j2 + m2));
for (int k = K; k <= N; k++)
{
real f0 = Factorial(k);
real f1 = Factorial(j1 + j2 - j3 - k);
real f2 = Factorial(j1 - m1 - k);
real f3 = Factorial(j2 + m2 - k);
real f4 = Factorial(j3 - j2 + m1 + k);
real f5 = Factorial(j3 - j1 - m2 + k);
sumK += powl(-1.0, k) / (f0 * f1 * f2 * f3 * f4 * f5);
}
return numer * sqrt1 * sumK;
}
real GauntCoefficient(
int l1, int l2, int l3, int m1, int m2, int m3)
{
real factor = sqrtl(
(2.0 * l1 + 1.0) *
(2.0 * l2 + 1.0) *
(2.0 * l3 + 1.0) /
(4.0 * pi));
real wigner1 = Wigner3j(l1, l2, l3, 0, 0, 0);
real wigner2 = Wigner3j(l1, l2, l3, m1, m2, m3);
return factor * wigner1 * wigner2;
}
int main()
{
pi = 4.0 * atanl(1.0);
cout << "Gaunt(1, 0, 1, 1, 0, 0) = ";
cout << GauntCoefficient(1, 0, 1, 1, 0, 0);
cout << endl;
cout << "Gaunt(1, 0, 1, 1, 0, -1) = ";
cout << GauntCoefficient(1, 0, 1, 1, 0, -1);
cout << endl;
real number = -1.0 / 2.0 / sqrtl(pi);
cout << "-1.0 / 2.0 / sqrt(pi) = ";
cout << number << endl;
return 0;
}
Tag: C/C++
The programming language I utilized.
Blog Entry (c) Tuesday, July 23, 2024, by James Pate Williams, Jr. Mueller’s Method for Finding the Complex and/or Real Roots of a Complex and/or Real Polynomial
I originally implemented this algorithm in FORTRAN IV in the Summer Quarter of 1982 at the Georgia Institute of Technology. I was taking a course named “Scientific Computing I” taught by Professor Gunter Meyer. I made a B in the class. Later in 2015 I re-implemented the recipe in C# using Visual Studio 2008 Professional. VS 2015 did not have support for complex numbers nor large integers. In December of 2015 I upgraded to Visual Studio 2015 Professional which has support for big integers and complex numbers. I used Visual Studio 2019 Community version for this project. Root below should be function.
Degree (0 to quit) = 2
coefficient[2].real = 1
coefficient[2].imag = 0
coefficient[1].real = 1
coefficient[1].imag = 0
coefficient[0].real = 1
coefficient[0].imag = 0
zero[0].real = -5.0000000000e-01 zero[0].imag = 8.6602540378e-01
zero[1].real = -5.0000000000e-01 zero[1].imag = -8.6602540378e-01
root[0].real = 0.0000000000e+00 root[0].imag = -2.2204460493e-16
root[1].real = 3.3306690739e-16 root[1].imag = -7.7715611724e-16
Degree (0 to quit) = 3
coefficient[3].real = 1
coefficient[3].imag = 0
coefficient[2].real = 0
coefficient[2].imag = 0
coefficient[1].real = -18.1
coefficient[1].imag = 0
coefficient[0].real = -34.8
coefficient[0].imag = 0
zero[0].real = -2.5026325486e+00 zero[0].imag = -8.3036679880e-01
zero[1].real = -2.5026325486e+00 zero[1].imag = 8.3036679880e-01
zero[2].real = 5.0052650973e+00 zero[2].imag = 2.7417672687e-15
root[0].real = 0.0000000000e+00 root[0].imag = 1.7763568394e-15
root[1].real = 3.5527136788e-14 root[1].imag = -1.7763568394e-14
root[2].real = 2.8421709430e-14 root[2].imag = 1.5643985575e-13
Degree (0 to quit) = 5
coefficient[5].real = 1
coefficient[5].imag = 0
coefficient[4].real = 2
coefficient[4].imag = 0
coefficient[3].real = 3
coefficient[3].imag = 0
coefficient[2].real = 4
coefficient[2].imag = 0
coefficient[1].real = 5
coefficient[1].imag = 0
coefficient[0].real = 6
coefficient[0].imag = 0
zero[0].real = -8.0578646939e-01 zero[0].imag = 1.2229047134e+00
zero[1].real = -8.0578646939e-01 zero[1].imag = -1.2229047134e+00
zero[2].real = 5.5168546346e-01 zero[2].imag = 1.2533488603e+00
zero[3].real = 5.5168546346e-01 zero[3].imag = -1.2533488603e+00
zero[4].real = -1.4917979881e+00 zero[4].imag = 1.8329656063e-15
root[0].real = 8.8817841970e-16 root[0].imag = 4.4408920985e-16
root[1].real = -2.6645352591e-15 root[1].imag = -4.4408920985e-16
root[2].real = 8.8817841970e-16 root[2].imag = 1.7763568394e-15
root[3].real = 3.4638958368e-14 root[3].imag = -1.4210854715e-14
root[4].real = 8.8817841970e-16 root[4].imag = 2.0710031449e-14
Blog Entry Friday, July 19, 2024, Easy Internet Math “Puzzle” (c) James Pate Williams, Jr.
#include <math.h>
#include <iostream>
using namespace std;
long double f(long double x)
{
return powl(8.0, x) - powl(2.0, x) -
2.0 * (powl(6.0, x) - powl(3.0, x));
}
long double g(long double x)
{
return powl(8.0, x) * logl(8.0) - powl(2.0, x) * logl(2.0) -
2.0 * (powl(6.0, x) * logl(6.0) - powl(3.0, x) * logl(3.0));
}
long double Newton(long double x, int maxIts, int& iterations)
{
long double x0 = x;
long double x1 = 0.0;
iterations = 0;
while (true) {
long double dx = 0.0;
long double fx = f(x0);
long double gx = g(x0);
x1 = x0 - fx / gx;
dx = fabsl(x1 - x0);
iterations++;
if (dx < 1.0e-15)
break;
if (fabsl(fx) < 1.0e-15)
break;
if (iterations == maxIts)
break;
x0 = x1;
}
return x1;
}
int main() {
int iterations = 0, maxIts;
long double x0 = 0.0, x1 = 0.0;
while (true) {
cout << "x0 = ";
cin >> x0;
if (x0 == 0)
break;
cout << "maximum iterations = ";
cin >> maxIts;
x1 = Newton(x0, maxIts, iterations);
cout << "x1 = " << x1 << endl;
cout << "iterations = ";
cout << iterations << endl;
}
return 0;
}
Blog Entry Monday, June 24, 2024 (c) James Pate Williams, Jr. Computing Binomial Coefficients and Pascal’s Triangle in the C Language
Enter n (<= 18) below:
5
Enter k (<= 18) below:
0
1 1
Enter n (<= 18) below:
5
Enter k (<= 18) below:
1
5 5
Enter n (<= 18) below:
5
Enter k (<= 18) below:
2
10 10
Enter n (<= 18) below:
0
Enter n (<= 18) below:
0
Pascal's Triangle:
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
1 7 21 35 35 21 7 1
1 8 28 56 70 56 28 8 1
1 9 36 84 126 126 84 36 9 1
1 10 45 120 210 252 210 120 45 10 1
1 11 55 165 330 462 462 330 165 55 11 1
1 12 66 220 495 792 924 792 495 220 66 12 1
1 13 78 286 715 1287 1716 1716 1287 715 286 78 13 1
1 14 91 364 1001 2002 3003 3432 3003 2002 1001 364 91 14 1
1 15 105 455 1365 3003 5005 6435 6435 5005 3003 1365 455 105 15 1
1 16 120 560 1820 4368 8008 11440 12870 11440 8008 4368 1820 560 120 16 1
1 17 136 680 2380 6188 12376 19448 24310 24310 19448 12376 6188 2380 680 136 17 1
1 18 153 816 3060 8568 18564 31824 43758 48620 43758 31824 18564 8568 3060 816 153 18 1
C:\Users\james\source\repos\BinomialCoefficeint\Debug\BinomialCoefficeint.exe (process 40028) exited with code 0.
Press any key to close this window . . .
// BinomialCoefficient.c (c) Monday, June 24, 2024
// by James Pate Williams, Jr. BA, BS, MSwE, PhD
#include <stdio.h>
#include <stdlib.h>
typedef long long ll;
ll** Binomial(ll n)
{
ll** C = (ll**)calloc(n + 1, sizeof(ll*));
if (C == NULL)
exit(-1);
for (int i = 0; i < n + 1; i++)
{
C[i] = (ll*)calloc(n + 1, sizeof(ll));
if (C[i] == NULL)
exit(-1);
}
if (n >= 0)
{
C[0][0] = 1;
}
if (n >= 1)
{
C[1][0] = 1;
C[1][1] = 1;
}
if (n >= 2)
{
for (int i = 2; i <= n; i++)
{
for (int j = 2; j <= n; j++)
{
C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
}
}
}
return C;
}
ll Factorial(ll n)
{
ll fact = 1;
if (n > 1)
{
for (int i = 2; i <= n; i++)
fact = i * fact;
}
return fact;
}
ll BC(ll n, ll k)
{
return Factorial(n) / (Factorial(n - k) * Factorial(k));
}
int main()
{
int i = 0, j = 0;
ll** C = Binomial(20);
while (1)
{
char buffer[256] = { '\0' };
printf_s("Enter n (<= 18) below:\n");
scanf_s("%s", buffer, sizeof(buffer));
printf_s("\n");
ll n = atoll(buffer);
if (n == 0)
break;
printf_s("Enter k (<= 18) below:\n");
scanf_s("%s", buffer, sizeof(buffer));
printf_s("\n");
ll k = atoll(buffer);
printf_s("%lld\t%lld\n\n", C[n + 2][k + 2], BC(n, k));
}
printf_s("Pascal's Triangle:\n\n");
for (i = 2; i <= 20; i++)
{
for (j = 2; j <= 20; j++)
if (C[i][j] != 0)
printf_s("%5lld ", C[i][j]);
printf_s("\n");
}
for (i = 0; i <= 20; i++)
free(C[i]);
free(C);
}

Blog Entry Sunday, June 23, 2024 (c) James Pate Williams, Jr.
The object of this C Win32 application is to find a multiple of 9 with its digits summing to a multiple of 9 also. The first column below is a multiple of 9 whose digits sum to 9 also. The second column is the sum of digits found in the column one number. The last column is the first column divided by 9.
Enter PRNG seed:
1
Enter number of bits (4 to 16):
4
9 9 1
Enter number of bits (4 to 16):
5
27 9 3
Enter number of bits (4 to 16):
6
45 9 5
Enter number of bits (4 to 16):
7
117 9 13
Enter number of bits (4 to 16):
8
252 9 28
Enter number of bits (4 to 16):
0
C:\Users\james\source\repos\CProductOf9Console\Debug\CProductOf9Console.exe (process 23280) exited with code 0.
Press any key to close this window . . .
Enter PRNG seed:
1
Enter number of bits (4 to 16):
9
369 18 41
Enter number of bits (4 to 16):
10
846 18 94
Enter number of bits (4 to 16):
11
1080 9 120
Enter number of bits (4 to 16):
12
3015 9 335
Enter number of bits (4 to 16):
13
5040 9 560
Enter number of bits (4 to 16):
14
10350 9 1150
Enter number of bits (4 to 16):
15
30870 18 3430
Enter number of bits (4 to 16):
16
57798 36 6422
Enter number of bits (4 to 16):
0
// CProductOf9Console.c (c) Sunday, June 23, 2024
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
char nextStr[256], numbStr[256];
void ConvertToString(int number, int radix)
{
int i = 0;
while (number > 0)
{
nextStr[i++] = (char)(number % radix + '0');
number /= radix;
}
nextStr[i++] = '\0';
_strrev(nextStr);
}
int Sum(int next)
{
long sum = 0;
ConvertToString(next, 10);
for (int i = 0; i < (int)strlen(nextStr); i++)
sum += (long)nextStr[i] - '0';
if (sum % 9 == 0 && sum != 0)
return sum;
return -1;
}
long GetNext(int numBits, int* next)
{
long hi = 0, lo = 0, nine = 0;
nextStr[0] = '\0';
numbStr[0] = '\0';
if (numBits == 4)
{
while (1)
{
*next = 9 * (long)(rand() % 16);
if (*next != 0 && *next >= 8 && *next < 16)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 5)
{
while (1)
{
*next = 9 * (long)(rand() % 32);
if (*next >= 16 && *next < 32)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 6)
{
while (1)
{
*next = 9 * (long)(rand() % 64);
if (*next >= 32 && *next < 64)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 7)
{
while (1)
{
*next = 9 * (long)(rand() % 128);
if (*next >= 64 && *next < 128)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 8)
{
while (1)
{
*next = 9 * (long)(rand() % 256);
if (*next >= 128 && *next < 256)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 9)
{
while (1)
{
*next = 9 * (long)(rand() % 512);
if (*next >= 256 && *next < 512)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 10)
{
while (1)
{
*next = 9 * (long)(rand() % 1024);
if (*next >= 512 && *next < 1024)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 11)
{
while (1)
{
*next = 9 * (long)(rand() % 2048);
if (*next >= 1024 && *next < 2048)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 12)
{
while (1)
{
*next = 9 * (long)(rand() % 4096);
if (*next >= 2048 && *next < 4096)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 13)
{
while (1)
{
*next = 9 * (long)(rand() % 8192);
if (*next >= 4096 && *next < 8192)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 14)
{
while (1)
{
*next = 9 * (long)(rand() % 16384);
if (*next >= 8192 && *next < 16384)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 15)
{
while (1)
{
*next = 9 * (long)(rand() % 32768);
if (*next >= 16384 && *next < 32768)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
else if (numBits == 16)
{
while (1)
{
*next = 9 * (long)(rand() % 65536);
if (*next >= 32768 && *next < 65536)
{
nine = Sum(*next);
if (nine % 9 == 0)
return nine;
}
}
}
return -1;
}
int main()
{
char buffer[256] = { '\0' };
long seed = 0;
printf_s("Enter PRNG seed:\n");
scanf_s("%s", buffer, sizeof(buffer));
seed = atol(buffer);
srand((unsigned int)seed);
while (1)
{
int next = 0, nine = 0, numberBits = 0;
printf_s("Enter number of bits (4 to 16):\n");
scanf_s("%s", buffer, sizeof(buffer));
numberBits = atol(buffer);
if (numberBits == 0)
break;
if (numberBits < 4 || numberBits > 16)
{
printf_s("illegal number of bits must >= 4 and <= 16\n");
continue;
}
nine = GetNext(numberBits, &next);
if (nine == -1)
{
printf_s("illegal result, try again\n");
continue;
}
printf_s("%5ld\t%5ld\t%5ld\n", next, nine, next / 9);
}
return 0;
}
Blog Entry (c) Friday, June 21, 2024, by James Pate Williams, Jr. Comparison of Two Prime Number Sieves
First the C++ results:
Limit = 1000000
Number of primes <= 1000000 78498
Milliseconds taken by Sieve of Atkin: 12
Number of primes <= 1000000 78498
Milliseconds taken by Sieve of Eratosthenes: 14
Limit = 10000000
Number of primes <= 10000000 664579
Milliseconds taken by Sieve of Atkin: 159
Number of primes <= 10000000 664579
Milliseconds taken by Sieve of Eratosthenes: 204
Limit = 100000000
Number of primes <= 100000000 5761455
Milliseconds taken by Sieve of Atkin: 1949
Number of primes <= 100000000 5761455
Milliseconds taken by Sieve of Eratosthenes: 2343
Limit = 0
Next, we have the Java results:
C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 1000000 0
number of primes less than equal 1000000 = 78498
total computation time in seconds = 0.008
C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 10000000 0
number of primes less than equal 10000000 = 664579
total computation time in seconds = 0.098
C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 1000000 0
number of primes less than equal 1000000 = 78498
total computation time in seconds = 0.011
C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 10000000 0
number of primes less than equal 10000000 = 664579
total computation time in seconds = 0.151
C:\WINDOWS\system32>java -jar k:\SieveOfAtkin\build\Debug\SieveOfAtkin.jar 100000000 0
number of primes less than equal 100000000 = 5761455
total computation time in seconds = 1.511
C:\WINDOWS\system32>java -jar k:\SieveOfEratosthenes\build\Debug\SieveOfEratosthenes.jar 100000000 0
number of primes less than equal 100000000 = 5761455
total computation time in seconds = 1.995
Notice that the Java application outperforms the C++ application.
// PrimeSieveComparison.cpp (c) Friday, June 21, 2024
// by James Pate Williams, Jr.
//
// SieveOfAtkin.java
// SieveOfAtkin
//
// Created by James Pate Williams, Jr. on 9/29/07.
// Copyright (c) 2007 James Pate Williams, Jr. All rights reserved.
//
// SieveOfEratosthenes.java
// SieveOfEratosthenes
//
// Created by James Pate Williams, Jr. on 9/29/07.
// Copyright (c) 2007 James Pate Williams, Jr. All rights reserved.
//
#include <math.h>
#include <iostream>
#include <chrono>
using namespace std::chrono;
using namespace std;
const int Maximum = 100000000;
bool sieve[Maximum + 1];
void SieveOfAtkin(int limit)
{
auto start = high_resolution_clock::now();
int e, k, n, p, x, xx3, xx4, y, yy;
int primeCount = 2, sqrtLimit = (int)sqrt(limit);
for (n = 5; n <= limit; n++)
sieve[n] = false;
for (x = 1; x <= sqrtLimit; x++) {
xx3 = 3 * x * x;
xx4 = 4 * x * x;
for (y = 1; y <= sqrtLimit; y++) {
yy = y * y;
n = xx4 + yy;
if (n <= limit && (n % 12 == 1 || n % 12 == 5))
sieve[n] = !sieve[n];
n = xx3 + yy;
if (n <= limit && n % 12 == 7)
sieve[n] = !sieve[n];
n = xx3 - yy;
if (x > y && n <= limit && n % 12 == 11)
sieve[n] = !sieve[n];
}
}
for (n = 5; n <= sqrtLimit; n++) {
if (sieve[n]) {
e = 1;
p = n * n;
while (true) {
k = e * p;
if (k > limit)
break;
sieve[k] = false;
e++;
}
}
}
for (n = 5; n <= limit; n++)
if (sieve[n])
primeCount++;
auto stop = high_resolution_clock::now();
auto duration = duration_cast<milliseconds>(stop - start);
std::cout << "Number of primes <= " << limit << ' ';
std::cout << primeCount << endl;
std::cout << "Milliseconds taken by Sieve of Atkin: "
<< duration.count() << endl;
}
void SieveOfEratosthenes(int limit)
{
auto start = high_resolution_clock::now();
int i = 0, k = 0, n = 0, nn = 0;
int primeCount = 0, sqrtLimit = (int)sqrt(limit);
// initialize the prime number sieve
for (n = 2; n <= limit; n++)
sieve[n] = true;
// eliminate the multiples of n
for (n = 2; n <= sqrtLimit; n++)
for (i = 2; i <= n - 1; i++)
sieve[i * n] = false;
// eliminate squares
for (n = 2; n <= sqrtLimit; n++) {
if (sieve[n]) {
k = 0;
nn = n * n;
i = nn + k * n;
while (i <= limit) {
sieve[i] = false;
i = nn + k * n;
k++;
}
}
}
primeCount = 0;
for (n = 2; n <= limit; n++)
if (sieve[n])
primeCount++;
auto stop = high_resolution_clock::now();
auto duration = duration_cast<milliseconds>(stop - start);
std::cout << "Number of primes <= " << limit << ' ';
std::cout << primeCount << endl;
std::cout << "Milliseconds taken by Sieve of Eratosthenes: "
<< duration.count() << endl;
}
int main()
{
while (true)
{
int limit = 0;
std::cout << "Limit = ";
cin >> limit;
if (limit == 0)
break;
SieveOfAtkin(limit);
SieveOfEratosthenes(limit);
}
return 0;
}
Classical Shor’s Algorithm Versus J. M. Pollard’s Factoring with Cubic Integers
We tried to factor the following numbers with each algorithm: 11^3+2, 2^33+2, 5^15+2, 2^66+2, 2^72+2, 2^81+2, 2^101+2, 2^129+2, and 2^183+2. Shor’s algorithm fully factored all of the numbers. Factoring with cubic integers fully factored all numbers except 2^66+2, 2^71+2, 2^129+2, and 2^183+2.


















Typical full output from factoring with cubic integers:
A-Solutions = 973 B-Solutions = 234 Known Eqs = 614 Solutions = 1821 Rows = 1821 Columns = 1701 Kernel rank = 423 Sieved = 326434 Successes0 = 200863 Successes1 = 47073 Successes2 = 2708 Successes3 = 973 Successes4 = 1735 2417851639229258349412354 - 25 DDs 2 p 65537 p 414721 p 44479210368001 p Sets = 189 #Factor Base 1 = 501 #Factor Base 2 = 868 FactB1 time = 00:00:00.000 FactB2 time = 00:00:05.296 Sieve time = 00:00:17.261 Kernel time = 00:00:06.799 Factor time = 00:00:02.327 Total time = 00:00:31.742
A-solutions have no large prime. B-solutions have a large prime between B0 and B1 exclusively which is this case is between 3272 and 50000 exclusively. The known equations are between the rational primes and the cubic primes and their associates of the form p = 6k + 1 that have -2 as a cubic residue. There are 81 rational primes of the form and 243 cubic primes but we keep many other associates of the cubic primes so more a and b pairs are successfully algebraically factored. In out case the algebraic factor base has 868 members. The rational prime factor base also includes the negative unit -1. The kernel rank is the number of independent columns in the matrix. The number of dependent sets is equal to columns – rank which is this case 1701 – 423 = 1278. The number of (a, b) pairs sieved is 326434. Successes0 is the pairs that have gcd(a, b) = 1. Successes1 is the number of (a, b) pairs such that a+b*r is B0-smooth or can be factored by the first 500 primes and the negative unit. r is equal to 2^27. Successes2 is the number of (a, b) pairs whose N[a, b] = a^2-2*b^3 can be factored using the norms of the algebraic primes. Successes3 is the number of A-solutions that are algebraically and rationally smooth. Successes4 is the number of B-solutions without combining to make the count modulo 2 = 0. Successes3 + Successes4 should equal Successes2 provided all proper algebraic primes and their associates are utilized.
Note factoring with cubic integers is very fickle with respect to parameter choice.
Root Finding Algorithms by James Pate Williams, BA, BS, MSwE, PhD
We designed and implemented a C# application that uses the following root finding algorithms:
- Bisection Method
- Brent’s Method
- Newton’s Method
- Regula Falsi
https://en.wikipedia.org/wiki/Bisection_method
https://en.wikipedia.org/wiki/Brent%27s_method
https://en.wikipedia.org/wiki/Newton%27s_method
https://en.wikipedia.org/wiki/False_position_method

















The source code files are displayed below as Word files:
Roots of Small Degree Polynomials with Real Coefficients by James Pate Williams, BA, BS, MSwE, PhD
We designed and implemented quadratic formula, cubic formula, and quartic formula solvers using the formulas in the Wikipedia articles:
https://en.wikipedia.org/wiki/Quadratic_formula
https://en.wikipedia.org/wiki/Cubic_function
https://en.wikipedia.org/wiki/Quartic_function
We tested our C# implementation against:
https://keisan.casio.com/exec/system/1181809415
http://www.wolframalpha.com/widgets/view.jsp?id=3f4366aeb9c157cf9a30c90693eafc55
https://keisan.casio.com/exec/system/1181809416
Here are screenshots of the C# application:








C# source code files for the application:
The Bailey-Borwein-Plouffe Formula for Calculating the First n Digits of Pi
The Bailey-Borwein-Plouffe formula for determining the digits of pi was discovered in 1995. This formula has been utilized to find the exact digits of pi to many decimal places.
I recently re-implemented my legacy C and FreeLIP program that utilized the BBP formula. The new C# application uses a homegrown big unsigned decimal number package that includes methods for +, -, *, / operators and an exponentiation (power) function. I used short integers (16-bit signed integers) to represent the individual digits of the number in any base whose square can be expressed as a positive short integer. That includes the decimal base 10 and hexadecimal base 16. For this application the base was chosen to be 10. Also, included was a n-digits of pi function that used the C# language’s built-in BigInteger data type. Below are some screen shots of the program in action.


As you can easily see the BigInteger implementation is an order of magnitude faster that the BigDecimal version (actually around 27+ times faster).
Last, we include a link to a PDF containing data comparing calculations performed on a Intel based desktop versus an AMD based laptop.