The x -axis is the number to be tested, the y-axis is prime number bound for factoring, and the z-axis is the runtime in seconds.
I corrected my powering modulo a prime routine. I added Pollard’s p – 1 factoring method and Shanks-Mestre elliptic curve point counting algorithm.
number to be tested or 0 to quit:
10000019
number of primes in factor base:
10000
Prime sieving time = 3.220000
N[0] = 10000019
a = 7838973
b = 2449531
m = 9995356
q = 356977
P = (9786147, 3226544)
P1 = (0, 1)
P2 = (5887862, 8051455)
N[1] = 356977
a = 45561
b = 178451
m = 357946
q = 178973
P = (80627, 163299)
P1 = (0, 1)
P2 = (52101, 282559)
N[2] = 178973
a = 135281
b = 76426
m = 178996
q = 73
P = (10238, 98035)
P1 = (0, 1)
P2 = (46702, 94326)
number is proven prime
runtime in seconds = 35.471000
number to be tested or 0 to quit:
10015969
number of primes in factor base:
10000
Prime sieving time = 3.424000
N[0] = 10015969
a = 6613193
b = 3951715
m = 10013908
q = 2503477
P = (998314, 8329764)
P1 = (0, 1)
P2 = (6944357, 1053776)
N[1] = 2503477
a = 1175442
b = 379813
m = 2505736
q = 293
P = (646462, 1631861)
P1 = (0, 1)
P2 = (1477980, 88719)
number is proven prime
runtime in seconds = 5.612000
number to be tested or 0 to quit:
99997981
number of primes in factor base:
10000
Prime sieving time = 4.152000
N[0] = 99997981
a = 34129462
b = 80482974
m = 100001414
q = 181
P = (19305995, 40493835)
P1 = (0, 1)
P2 = (33828245, 72969559)
number is proven prime
runtime in seconds = 11.500000
number to be tested or 0 to quit:
100001819
number of primes in factor base:
100000
Prime sieving time = 3.218000
N[0] = 100001819
a = 2694060
b = 17329746
m = 100008102
q = 5569
P = (124594, 14596756)
P1 = (0, 1)
P2 = (32514144, 56926555)
number is proven prime
runtime in seconds = 76.301000
number to be tested or 0 to quit:
100005317
number of primes in factor base:
100000
Prime sieving time = 3.269000
N[0] = 100005317
a = 45478318
b = 328034
m = 99988256
q = 3124633
P = (62548529, 30179124)
P1 = (0, 1)
P2 = (70379514, 76899689)
N[1] = 3124633
a = 2605576
b = 1809212
m = 3127654
q = 503
P = (1236288, 2081401)
P1 = (0, 1)
P2 = (2264479, 2583693)
number is proven prime
runtime in seconds = 459.979000
number to be tested or 0 to quit:
100000007
number of primes in factor base:
100000
Prime sieving time = 3.209000
N[0] = 100000007
a = 50593669
b = 72502607
m = 100005736
q = 2053
P = (72365335, 69885097)
P1 = (0, 1)
P2 = (55023241, 20078454)
number is proven prime
runtime in seconds = 163.705000
number to be tested or 0 to quit:
100014437
number of primes in factor base:
100000
Prime sieving time = 3.919000
N[0] = 100014437
a = 49955472
b = 45482796
m = 100024160
q = 263
P = (41650735, 8652103)
P1 = (0, 1)
P2 = (53790105, 37282431)
number is proven prime
runtime in seconds = 12.915000
prime number sieve creation
time in seconds = 3.483000
number to be factored or 0 to quit:
2111222333
1 11 1 p
2 17 1 p
3 11289959 1 p
factoring time in seconds = 0.063000
number to be factored or 0 to quit:
1234567890
1 2 1 p
2 3 2 p
3 5 1 p
4 3607 1 p
5 3803 1 p
factoring time in seconds = 0.133000
number to be factored or 0 to quit:
2^30+0
prime powers are not allowed
number to be factored or 0 to quit:
0
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.
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;
}
The slow computations used 999,999,999 terms. I seem to recall from my first numerical analysis (Scientific Computing I) course in the Mathematics Department at the Georgia Institute of Technology with Professor Gunter Meyer in the Summer of 1982 that computing a truncated infinite series is more accurate to start with the smallest terms.
// RiemannZetaFunctionWin32Console.c (c) Saturday, July 13-15, 2024
// James Pate Williams, Jr. some computations use 999,999,998 terms
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
typedef signed long long sll;
// https://en.wikipedia.org/wiki/Particular_values_of_the_Riemann_zeta_function#cite_note-7
long double EvenZeta(int n)
{
sll A[12] = { 0ll,
6ll, 90ll, 945ll, 9450ll, 93555ll, 638512875ll,
18243225ll, 325641566250ll, 38979295480125ll,
1531329465290625ll, 13447856940643125ll };
sll B[12] = { 0ll, 1ll, 1ll, 1ll, 1ll, 1ll, 691ll, 2ll,
3617ll, 43867ll, 174611ll, 155366ll };
long double pi = 4.0 * atan(1.0);
return (long double)B[n / 2] * powl(pi, n) / (long double)A[n / 2];
}
long double Zeta(long s, long *terms)
{
int n = 0;
long double sum = 0.0;
*terms = 0;
for (n = 1000000000; n >= 2; n--)
{
sum += 1.0 / powl(n, s);
(*terms)++;
}
return sum + 1.0;
}
void PrintEvenZetaValue(int n)
{
clock_t time0 = clock();
long double ez = EvenZeta(n);
clock_t time1 = clock();
double runtime = ((double)time1 - time0) / CLOCKS_PER_SEC;
printf_s("Runtime in seconds to compute Zeta(%ld) = %16.15Lf: %Lf\n",
n, ez, runtime);
}
void PrintZetaValue(int n)
{
clock_t time0 = clock();
long terms = 0;
long double zeta = Zeta(n, &terms);
clock_t time1 = clock();
double runtime = ((double)time1 - time0) / CLOCKS_PER_SEC;
printf_s("Runtime in seconds to compute Zeta(%ld) = %16.15Lf terms = %ld: time = %Lf\n",
n, zeta, terms, runtime);
}
int main()
{
PrintEvenZetaValue(2);
PrintEvenZetaValue(4);
PrintEvenZetaValue(6);
PrintEvenZetaValue(8);
PrintZetaValue(2);
PrintZetaValue(4);
PrintZetaValue(6);
PrintZetaValue(8);
PrintZetaValue(3);
PrintZetaValue(5);
PrintZetaValue(7);
PrintZetaValue(9);
return 0;
}
I was performing an Internet search for prime counting functions, and I ran across the following PDF:
It was made available online in April of 2021. That made me remember the opening line of Geoffery Chaucer’s “Canterbury Tales”:
Whan that Aprill with his shoures soote
I seem to recall that in my English literature textbook ‘with’ was replaced with the German ‘mit’. Middle English was related to the old French and German languages I seem to recall.
Numerical Explorations 