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) =
Category: Numerical Analysis
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;
}
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.
// 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 Sunday, June 16, 2024 (c) James Pate Williams, Jr. Chapter 4 Matrices and Systems of Linear Equations from a Textbook by S. D. Conte and Carl de Boor
Blog Entry Friday, June 14, 2024 (c) James Pate Williams, Jr.
For the last week or so I have been working my way through Chapter 3 The Solution of Nonlinear Equations found in the textbook “Numerical Analysis: An Algorithmic Approach” by S. D. Conte and Carl de Boor. I also used some C source code from “A Numerical Library in C for Scientists and Engineers” by H. T. Lau, PhD. I implemented twenty examples and exercises from the previously mentioned chapter.
Numerical Explorations 