Blog Entry © Sunday, July 27, 2025, A Bit of Programming Nostalgia Prime Number Related Programs by James Williams, Jr.

/*
  Author:  Pate Williams c 1995

  The following program is a solution to problem 18.15
  in Pascalgorithms by Edwin D. Reilly and Francis D.
  Federighi page 627. The program uses Simpson's rule
  to calculate the number of primes less than or equal
  a given number.
*/

#include <math.h>
#include <stdio.h>

typedef double real;

static real f(real x)
{
	return(1.0 / log(x));
}

static real simpson(int n, real a, real b)
{
	int i;
	real evensum, h, oddsum, twoh, x;

	if (n % 2 == 1) n = n - 1;
	h = (b - a) / n;
	twoh = h + h;
	x = h + a;
	oddsum = 0.0;
	for (i = 1; i <= n / 2; i++)
	{
		oddsum += f(x);
		x = twoh + x;
	}
	x = twoh + a;
	evensum = 0.0;
	for (i = 1; i <= n / 2 - 1; i++)
	{
		evensum += f(x);
		x = twoh + x;
	}
	return(h / 3.0 * (f(a) + f(b) + 4.0 * oddsum + 2.0 * evensum));
}

int main(void)
{
	int i, n, Nmaximum = 0, Nminimum = 0, Nstep;

	printf("n = ");			scanf_s("%d", &n);
	printf("N minimum = "); scanf_s("%d", &Nminimum);
	printf("N maximum = "); scanf_s("%d", &Nmaximum);
	printf("N step = ");	scanf_s("%d", &Nstep);
	printf("\n");
	printf("----------------------------------------\n");
	printf("Min\t\tMax\t\tprimes\n");
	printf("----------------------------------------\n");
	for (i = Nminimum; i <= Nmaximum; i += Nstep)
	{
		printf("%8d\t%8d\t%8.0lf\n", Nminimum, i + Nstep,
			simpson(n, Nminimum, i + Nstep));
	}
	printf("----------------------------------------\n");
	return(0);
}
n = 1024
N minimum = 0
N maximum = 10000000
N step = 1000000

----------------------------------------
Min             Max             primes
----------------------------------------
       0         1000000           78551
       0         2000000          148923
       0         3000000          216788
       0         4000000          283122
       0         5000000          348361
       0         6000000          412754
       0         7000000          476461
       0         8000000          539590
       0         9000000          602224
       0        10000000          664424
       0        11000000          726239
----------------------------------------

D:\PrimeCounter\x64\Release\PrimeCounter.exe (process 51884) exited with code 0 (0x0).
Press any key to close this window . . .
/*
  Author:  Pate Williams c 1995

  The following is a translation of the Pascal program
  sieve found in Pascalgorithms by Edwin D. Reilly and
  Francis D. Federighi page 652. This program uses sets
  to represent the sieve (see C Programming Language An
  Applied Perspective by Lawrence Miller and Alec Qui-
  lici pages 160 - 162).
*/

#include <math.h>
#include <stdio.h>

#define _WORD_SIZE 32
#define _VECT_SIZE 524288
#define SET_MIN    0
#define SET_MAX    16777215

typedef unsigned long SET[_VECT_SIZE];
typedef long ELEMENT;
typedef unsigned long LONG;

SET set;

static int get_bit_pos(int* long_ptr, int* bit_ptr,
	ELEMENT element)
{
	*long_ptr = element / _WORD_SIZE;
	*bit_ptr = element % _WORD_SIZE;
	return(element >= SET_MIN && element <= SET_MAX);
}

static void set_bit(ELEMENT element, int inset)
{
	int bit, word;

	if (get_bit_pos(&word, &bit, element))
	{
		if (inset > 0)
			set[word] |= (01 << bit);
		else
			set[word] &= ~(01 << bit);
	}
}

static int get_bit(ELEMENT element)
{
	int bit, word;

	return(get_bit_pos(&word, &bit, element) ?
		(set[word] >> bit) & 01 : 0);
}

static void set_Add(ELEMENT element)
{
	set_bit(element, 1);
}

static void set_Del(ELEMENT element)
{
	set_bit(element, 0);
}

static int set_Mem(ELEMENT element)
{
	return get_bit(element);
}

static void primes(long n)
{
	long c, i, inc, k;
	double x;

	set_Add(2);
	for (i = 3; i <= n; i++)
		if ((i + 1) % 2 == 0)
			set_Add(i);
		else
			set_Del(i);
	c = 3;
	do
	{
		i = c * c;
		inc = c + c;
		while (i <= n)
		{
			set_Del(i);
			i = i + inc;
		}
		c += 2;
		while (set_Mem(c) == 0) c += 1;
	} while (c * c <= n);
	k = 0;
	for (i = 2; i <= n; i++)
		if (set_Mem(i) == 1) k++;
	x = n / log(n) - 5.0;
	x = x + exp(1.0 + 0.15 * log(n) * sqrt(log(n)));
	printf("%8ld\t%8ld\t%8.0lf\n", n, k, x);
}

int main(void)
{
	long n = 100L;

	printf("----------------------------------------\n");
	printf("n\t\tprimes\t\ttheory\n");
	printf("----------------------------------------\n");
	do
	{
		primes((int)n);
		n = 10L * n;
	} while (n < (long)SET_MAX);
	printf("----------------------------------------\n");
	return(0);
}
----------------------------------------
n               primes          theory
----------------------------------------
     100              25              29
    1000             168             181
   10000            1229            1261
  100000            9592            9634
 1000000           78498           78396
10000000          664579          665060
----------------------------------------

D:\Sieve\x64\Release\Sieve.exe (process 60092) exited with code 0 (0x0).
Press any key to close this window . . .

Blog Entry (c) Friday, April 11, 2025, by James Pate Williams, Jr. Multiplication and Division of Finite Power Series

x = 0.25
N = 5
Series Cosine(x)  = 0.968912421711
C++ cos(x)        = 0.968912421711
Series sine(x)    = 0.247403959255
C++ sin(x)        = 0.247403959255
Series Tangent(x) = 0.255341921221
C++ tan(x)        = 0.255341921221
Series sin(2x)    = 0.479425538604
C++ sin(2x)       = 0.479425538604
C++ 2sin(x)cos(x) = 0.479425538604
End app ? y = yes = n
x = 0.5
N = 5
Series Cosine(x)  = 0.877582561890
C++ cos(x)        = 0.877582561890
Series sine(x)    = 0.479425538604
C++ sin(x)        = 0.479425538604
Series Tangent(x) = 0.546302489844
C++ tan(x)        = 0.546302489844
Series sin(2x)    = 0.841470984807
C++ sin(2x)       = 0.841470984808
C++ 2sin(x)cos(x) = 0.841470984808
End app ? y = yes = y
x = 0.75
N = 5
Series Cosine(x)  = 0.731688868808
C++ cos(x)        = 0.731688868874
Series sine(x)    = 0.681638760020
C++ sin(x)        = 0.681638760023
Series Tangent(x) = 0.931596460023
C++ tan(x)        = 0.931596459944
Series sin(2x)    = 0.997494986509
C++ sin(2x)       = 0.997494986604
C++ 2sin(x)cos(x) = 0.997494986604
End app ? y = yes = n
 
x = 1.00
N = 5
Series Cosine(x)  = 0.540302303792
C++ cos(x)        = 0.540302305868
Series sine(x)    = 0.841470984648
C++ sin(x)        = 0.841470984808
Series Tangent(x) = 1.557407730344
C++ tan(x)        = 1.557407724655
Series sin(2x)    = 0.909297423159
C++ sin(2x)       = 0.909297426826
C++ 2sin(x)cos(x) = 0.909297426826
End app ? y = yes = n
x = 1.25
N = 5
Series Cosine(x)  = 0.315322332275
C++ cos(x)        = 0.315322362395
Series sine(x)    = 0.948984616456
C++ sin(x)        = 0.948984619356
Series Tangent(x) = 3.009569952151
C++ tan(x)        = 3.009569673863
Series sin(2x)    = 0.598472085108
C++ sin(2x)       = 0.598472144104
C++ 2sin(x)cos(x) = 0.598472144104
End app ? y = yes = n
x = 1.50
N = 5
Series Cosine(x)  = 0.070736934117
C++ cos(x)        = 0.070737201668
Series sine(x)    = 0.997494955682
C++ sin(x)        = 0.997494986604
Series Tangent(x) = 14.101472846329
C++ tan(x)        = 14.101419947172
Series sin(2x)    = 0.141119469924
C++ sin(2x)       = 0.141120008060
C++ 2sin(x)cos(x) = 0.141120008060
End app ? y = yes = y

// DivMulPowerSeries.cpp (c) Tuesday, April 8, 2025
// by James Pate Williams, Jr.
// https://math.libretexts.org/Bookshelves/Calculus/Calculus_(OpenStax)/10%3A_Power_Series/10.02%3A_Properties_of_Power_Series
// https://en.wikipedia.org/wiki/Formal_power_series

#include <iostream>
#include <iomanip>
#include <string>
#include <vector>

std::vector<double> Multiplication(
	std::vector<double> c,
	std::vector<double> d,
	int N) {
	std::vector<double> e(N + 1);

	for (int n = 0; n <= N; n++) {
		double sum = 0.0;

		for (int k = 0; k <= n; k++) {
			sum += c[k] * d[n - k];
		}

		e[n] = sum;
	}

	return e;
}

std::vector<double> Division(
	std::vector<double> a,
	std::vector<double> b,
	int N) {
	double a0 = a[0];
	std::vector<double> c(N + 1);

	for (int n = 0; n <= N; n++) {
		double sum = 0.0;

		for (int k = 1; k <= n; k++) {
			sum += a[k] * c[n - k];
		}

		c[n] = (b[n] - sum) / a0;
	}

	return c;
}

double Factorial(int n) {
	double nf = 1.0;

	for (int i = 2; i <= n; i++) {
		nf *= i;
	}

	return nf;
}

std::vector<double> Cosine(double x, int N, double& fx) {
	std::vector<double> series(N + 1);

	fx = 0.0;

	for (int n = 0; n <= N; n++) {
		int argument = 2 * n;
		double coeff = pow(-1, n) / Factorial(argument);
		series[n] = coeff;
		fx += coeff * pow(x, argument);
	}

	return series;
}

std::vector<double> Sine(double x, int N, double& fx) {
	std::vector<double> series(N + 1);
	
	fx = 0.0;

	for (int n = 0; n <= N; n++) {
		int argument = 2 * n + 1;
		double coeff = pow(-1, n) / Factorial(argument);
		series[n] = coeff;
		fx += coeff * pow(x, argument);
	}

	return series;
}

std::vector<double> Tangent(double x, int N, double& fx) {
	double fc = 0.0, fs = 0.0;
	std::vector<double> seriesC = Cosine(x, N, fc);
	std::vector<double> seriesS = Sine(x, N, fs);
	std::vector<double> seriesT = Division(seriesS, seriesC, N);
	fx = fs / fc;

	return seriesT;
}

std::vector<double> Sine2x(double x, int N, double& fx) {
	double fc = 0.0, fs = 0.0;
	std::vector<double> seriesC = Cosine(x, N, fc);
	std::vector<double> seriesS = Sine(x, N, fs);
	std::vector<double> series2 = Multiplication(seriesS, seriesC, N);
	fx = 2.0 * fs * fc;

	return series2;
}

int main()
{
	while (true) {
		char line[128] = { };
		std::cout << "x = ";
		std::cin.getline(line, 127);
		std::string str1(line);
		double x = std::stod(str1);
		std::cout << "N = ";
		std::cin.getline(line, 127);
		std::string str2(line);
		int N = std::stoi(str2);
		double cx = 0.0, sx = 0.0, tx = 0.0, xx = 0.0;
		std::vector<double> cSeries = Cosine(x, N, cx);
		std::vector<double> sSeries = Sine(x, N, sx);
		std::vector<double> tSeries = Tangent(x, N, tx);
		std::vector<double> xSeries = Sine2x(x, N, xx);
		std::cout << std::fixed << std::setprecision(12);
		std::cout << "Series Cosine(x)  = " << cx << std::endl;
		std::cout << "C++ cos(x)        = " << cos(x) << std::endl;
		std::cout << "Series sine(x)    = " << sx << std::endl;
		std::cout << "C++ sin(x)        = " << sin(x) << std::endl;
		std::cout << "Series Tangent(x) = " << tx << std::endl;
		std::cout << "C++ tan(x)        = " << tan(x) << std::endl;
		std::cout << "Series sin(2x)    = " << xx << std::endl;
		std::cout << "C++ sin(2x)       = " << sin(x + x) << std::endl;
		std::cout << "C++ 2sin(x)cos(x) = " << 2.0 * sin(x) * cos(x);
		std::cout << std::endl;
		std::cout << "End app ? y = yes = ";
		std::cin.getline(line, 127);

		if (line[0] == 'Y' || line[0] == 'y') {
			break;
		}
	}

	return 0;
}

Blog Entry © Sunday, March 29, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Slater Determinant Coefficients for Z = 2 to 4

Enter the atomic number Z (2 to 6 or 0 to quit): 2
2       1       1       +       a(1)b(2)
1       0       0       -       a(2)b(1)
# Even Permutations = 1
Enter the atomic number Z (2 to 6 or 0 to quit): 3
6       3       1       +       a(1)b(2)c(3)
5       2       0       -       a(1)b(3)c(2)
4       2       0       -       a(2)b(1)c(3)
3       1       1       +       a(2)b(3)c(1)
2       1       1       +       a(3)b(1)c(2)
1       0       0       -       a(3)b(2)c(1)
# Even Permutations = 3
Enter the atomic number Z (2 to 6 or 0 to quit): 4
24      12      0       +       a(1)b(2)c(3)d(4)
23      11      1       -       a(1)b(2)c(4)d(3)
22      11      1       -       a(1)b(3)c(2)d(4)
21      10      0       +       a(1)b(3)c(4)d(2)
20      10      0       +       a(1)b(4)c(2)d(3)
19      9       1       -       a(1)b(4)c(3)d(2)
18      9       1       -       a(2)b(1)c(3)d(4)
17      8       0       +       a(2)b(1)c(4)d(3)
16      8       0       +       a(2)b(3)c(1)d(4)
15      7       1       -       a(2)b(3)c(4)d(1)
14      7       1       -       a(2)b(4)c(1)d(3)
13      6       0       +       a(2)b(4)c(3)d(1)
12      6       0       +       a(3)b(1)c(2)d(4)
11      5       1       -       a(3)b(1)c(4)d(2)
10      5       1       -       a(3)b(2)c(1)d(4)
9       4       0       +       a(3)b(2)c(4)d(1)
8       4       0       +       a(3)b(4)c(1)d(2)
7       3       1       -       a(3)b(4)c(2)d(1)
6       3       1       -       a(4)b(1)c(2)d(3)
5       2       0       +       a(4)b(1)c(3)d(2)
4       2       0       +       a(4)b(2)c(1)d(3)
3       1       1       -       a(4)b(2)c(3)d(1)
2       1       1       -       a(4)b(3)c(1)d(2)
1       0       0       +       a(4)b(3)c(2)d(1)
# Even Permutations = 12
Enter the atomic number Z (2 to 6 or 0 to quit):
// AOPermutations.cpp : This file contains the 'main' function.
// Program execution begins and ends there.
// Copyright (c) Saturday, March 29, 2025
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Signs of the atomic orbitals in a Slater Determinant

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>

int main()
{
    char alpha[] = { 'a', 'b', 'c', 'd', 'e', 'f' }, line[128] = {};
    int factorial[7] = { 1, 1, 2, 6, 24, 120, 720 };

    while (true)
    {
        int col = 0, counter = 0, row = 0, sign = 1, t = 0, Z = 0, zfact = 0;
        int numberEven = 0;
        std::cout << "Enter the atomic number Z (2 to 6 or 0 to quit): ";
        std::cin.getline(line, 127);
        std::string str(line);
        Z = std::stoi(str);

        if (Z == 0)
        {
            break;
        }

        if (Z < 2 || Z > 6)
        {
            std::cout << "Illegal Z, please try again" << std::endl;
            continue;
        }

        zfact = factorial[Z];

        std::vector<char> orb(Z);
        std::vector<int> tmp(Z), vec(Z);

        for (int i = 0; i < Z; i++)
        {
            orb[i] = alpha[i];
            vec[i] = i + 1;
        }

        do
        {
            for (int i = 0; i < (int)vec.size(); i++)
            {
                tmp[i] = vec[i];
            }

            t = 0;

            do
            {
                t++;
            } while (std::next_permutation(tmp.begin(), tmp.end()));

            std::cout << t << '\t' << t / 2 << '\t';
            std::cout << (t / 2 & 1) << '\t';

            if (Z == 2 || Z == 3)
            {
                if ((t / 2 & 1) == 0)
                {
                    std::cout << "-\t";
                }

                else
                {
                    std::cout << "+\t";
                    numberEven++;
                }
            }

            else
            {
                if ((t / 2 & 1) == 1)
                {
                    std::cout << "-\t";
                }

                else
                {
                    std::cout << "+\t";
                    numberEven++;
                }
            }

            for (int i = 0; i < Z; i++)
            {
                std::cout << orb[i] << '(' << vec[i] << ')';
            }

            row++;
            std::cout << std::endl;

            if (zfact != 2 && row == zfact)
            {
                std::cout << std::endl;
                break;
            }

            row %= Z;
        } while (std::next_permutation(vec.begin(), vec.end()));

        std::cout << "# Even Permutations = ";
        std::cout << numberEven << std::endl;
    }

    return 0;
}