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):
#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;
}
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.
I became fascinated with secret key cryptography as a child. Later, as an adult, in around 1979, I started creating crude symmetric cryptographic algorithms. I became further enthralled with cryptography and number theory in 1996 upon reading Applied Cryptography, Second Edition: Protocols, Algorithms, and Source Code in C by Bruce Schneier and later the Handbook of Applied Cryptography by Alfred J. Menezes, Paul C. van Oorschot, and Scott A. Vanstone. After implementing many of the algorithms in both tomes, I communicated my results to two of the authors namely Bruce Schneier and Professor Alfred J. Menezes. In 1997 I developed a website devoted to constraint satisfaction problems and their solutions, cryptography, and number theory. I posted legal C and C++ source code. Professor Menezes advertised my website along with his treatise. See the following blurb:
In the spirit of my twin scientific infatuations, I offer yet another C integer factoring implementation utilizing the Free Large Integer Package (known more widely as lip) which was created by Arjen K. Lenstra (now a Professor Emeritus). This implementation includes Henri Cohen’s Trial Division algorithm, the Brent-Cohen-Pollard rho method, the Cohen-Pollard p – 1 stage 1 method, and the Lenstra lip Elliptic Curve Method. If I can get the proper authorization, I will later post the source code.
total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
1
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.014000 seconds:
2
5 ^ 2
41
397
2113
enter number below:
0
total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
2
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 1.531000 seconds:
2
5 ^ 2
41
397
2113
415878438361
3630105520141
enter number below:
0
total time required for initialization: 0.055000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
3
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.066000 seconds:
2
5 ^ 2
41
838861
415878438361
3630105520141
enter number below:
0
total time required for initialization: 0.056000 seconds
enter number below:
2^111+2
== Menu ==
1 Trial Division
2 Pollard-Brent-Cohen rho
3 p - 1 Pollard-Cohen
4 Lenstra's Elliptic Curve Method
5 Pollard-Lenstra rho
4
2596148429267413814265248164610050
number is composite
factors:
total time required factoring: 0.013000 seconds:
2
5
205
838861
415878438361
3630105520141
enter number below:
0
The sieve of Eratosthenes handles primes up to an upper bound of 100,000,000. The number of primes is 5,761,455. Below are a few examples runs of the app. I have also created C source code for trial division and other factoring algorithms that use Professor Emeritus Arjen K. Lenstra’s Free Large Integer Package also known as lip.
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);
}
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;
}
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;
}
Numerical Explorations 