Blog Entry Monday, July 15, 2024, (c) James Pate Williams, Jr. Some Values of the Riemann Zeta Function Computed with Fast and Very Slow Algorithms

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.

Runtime in seconds to compute Zeta FunctionsDownload
// 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;
}

Blog Entry Friday, July 12, 2024, © James Pate Williams, Jr.

I was performing an Internet search for prime counting functions, and I ran across the following PDF:

1905.09818 (arxiv.org)

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

  1. General Prologue | Harvard’s Geoffrey Chaucer Website

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.

Blog Entry Wednesday, July 10, 2024, © James Pate Williams, Jr. My Dual Interests in Cryptography and Number Theory

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

Blog Entry Monday, July 8, 2024, (c) James Pate Williams, Jr. Relatively Fast 64-bit Trial Division Using Henri Cohen’s Algorithm

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.

Long Long Trial Division C Source CodeDownload

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;
}

Blog Entry Tuesday, June 18, 2024 (c) James Pate Williams, Jr. FreeLIP Computation of Euler Numbers and Tangent Numbers

FreeLIP is a free large integer package solely created by Professor Emeritus Arjen K. Lenstra of the Number Field Sieve fame. He developed FreeLIP while he was an employee of AT&T – Lucent in the late 1980s. His copyright notice in the header file, lip.h, states copyright from 1989 to 1997. I have been using this excellent number theoretical package since the late 1990s. See the paper by Donald E. Knuth and Thomas J. Buckholtz for the formula for Tangent Numbers. I can’t remember where I got the Euler Numbers recurrence relation. I wrote a C# application in 2015 for computing Euler Numbers. The code below is in the vanilla C computer language. Excellent resources for the Euler and tangent numbers also known as zag numbers are:

https://oeis.org/A122045

https://oeis.org/A000182

Click to access S0025-5718-1967-0221735-9.pdf

euler-and-tangent-numbers-1Download
fleulernumbers-c-source-codeDownload
fltangentnumbersconsole-c-source-code-1Download