Blog Entry Wednesday, August 14, 2024 (c) James Pate Williams, Jr. Goldwasser-Kilian Primality Test
The Goldwasser-Kilian Primality proving algorithm was the first method to utilize elliptic curves to generate primality proving certificates. What follows is a file of two certificates and the single precision C source code.
Blog Entry (c) Saturday, August 10, 2024, by James Pate Williams, Jr. The LLL Lattice Reduction Algorithm
Blog Entry August 9, 2024, (c) James Pate Williams, Jr. Another Online Math Problem
The problem is to find the real root of the equation: f(x)=x^(x^8)-8=0. I use the Newton-Raphson method, a root finding algorithm. A first guess is x = 2. The solution is: x = 1.2968395547, f(x) = -2.6645353e-15. I compute the necessary derivative using central-finite differences with a step size of h = 2/10000.
Blog Entry August 1, 2024, (c) James Pate Williams, Jr. Online Math Problem
Blog Entry Thursday, July 25, 2024, (c) James Pate Williams, Jr. Newton’s Method for Finding the Real Roots of a Real Polynomial
degree (0 to quit) = 5
a[5] = 1
a[4] = -15.5
a[3] = 77.5
a[2] = -155
a[1] = 124
a[0] = -32
x[0] = 0.45
x[1] = 0.9
x[2] = 1.8
x[3] = 3.6
x[4] = 7.2
iterations = 31
root[0] = 5.0000000000e-01
root[1] = 1.0000000000e+00
root[2] = 2.0000000000e+00
root[3] = 4.0000000000e+00
root[4] = 8.0000000000e+00
func[0] = 0.0000000000e+00
func[1] = 0.0000000000e+00
func[2] = 0.0000000000e+00
func[3] = 0.0000000000e+00
func[4] = 0.0000000000e+00
degree (0 to quit) =
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 Sunday, July 21, 2024 (c) James Pate Williams, Jr. Another Easy Internet Mathematics Problem
x = 1.2679491924e+00 y = 4.7320508076e+00
f = 0.0000000000e+00 g = 0.0000000000e+00
iterations = 100
legend: f = x + y – 6, g = x * y – 6
The solution was found via my Win32 C application whose source code is presented below (the method is the Newton iteration for systems of linear and/or non-linear equations):
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 Tuesday, July 16, 2024, Prime Number Zeta Function Values for s = 2 and Primes < 1,000,000 (c) James Pate Williams, Jr.
Runtime in seconds to build sieve: 4.462000
N = <= 1000000 17
s > 1 = 2
prime zeta = 0.449282078106864
N = <= 1000000 257
s > 1 = 2
prime zeta = 0.452175428603043
N = <= 1000000 65537
s > 1 = 2
prime zeta = 0.452247336475390
N = <= 1000000 100003
s > 1 = 2
prime zeta = 0.452247368830138
N = <= 1000000 900001
s > 1 = 2
prime zeta = 0.452247415884861
N = <= 1000000 0
C:\Users\james\source\repos\PrimeZetaFunction\Debug\PrimeZetaFunction.exe (process 20584) exited with code 0.
Press any key to close this window . . .
/*
* PrimeZetaFunction.c (c) Monday, July 15, 2024
* by James Pate Williams, Jr.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define BITS_PER_LONG 32
#define BITS_PER_LONG_1 31
#define MAX_SIEVE 100000000
#define MAX_PRIME_INDEX 5761454
#define SIEVE_SIZE (MAX_SIEVE / BITS_PER_LONG + 1)
long prime[MAX_PRIME_INDEX], sieve[SIEVE_SIZE];
long 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)
{
long 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);
}
long double primeZetaFunction(long N, long s)
{
long n, p;
long double sum = 0.0L;
for (n = N; n >= 0; n--)
{
p = prime[n];
sum += 1.0 / powl((long double)p, (long double)s);
}
return sum;
}
int main()
{
long N = 0, i = 0, p = 2, s = 0;
double runtime = 0.0;
clock_t time0 = clock(), time1 = 0;
Sieve(MAX_SIEVE, sieve);
for (i = 0; i <= MAX_PRIME_INDEX; i++) {
while (!get_bit(p, sieve)) p++;
prime[i] = p++;
}
time1 = clock();
runtime = ((double)time1 - time0) / CLOCKS_PER_SEC;
printf_s("Runtime in seconds to build sieve: %Lf\n", runtime);
for (;;) {
printf_s("N = <= %ld ", 1000000);
scanf_s("%ld", &N);
if (N == 0)
break;
printf_s("s > 1 = ");
scanf_s("%ld", &s);
printf_s("prime zeta = %16.15Lf\n", primeZetaFunction(N, s));
}
return 0;
}
Numerical Explorations 