/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Elementary Numerical Analysis: An
* Algorithmic Approach" by S. D. Conte and Carl
* de Boor Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1
*/
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
//#define DEBUG
static _Lcomplex HornersMethod(_Lcomplex coeff[], _Lcomplex z, int degree) {
int i = 0;
_Lcomplex c = coeff[degree];
for (i = degree; i >= 1; i--) {
_Lcomplex product = _LCmulcc(c, z);
c._Val[0] = product._Val[0] + coeff[i - 1]._Val[0];
c._Val[1] = product._Val[1] + coeff[i - 1]._Val[1];
}
#ifdef DEBUG
_Lcomplex sum = { 0 };
for (i = 0; i <= degree; i++) {
_Lcomplex term =
_LCmulcc(cpowl(z, _LCbuild(i, 0.0)), coeff[i]);
sum._Val[0] += term._Val[0];
sum._Val[1] += term._Val[1];
}
long double delta = fabsl(cabsl(c) - cabsl(sum));
if (delta > 1.0e-12)
exit(-5);
#endif
return c;
}
static void comdiv(
long double xr, long double xi,
long double yr, long double yi,
long double* zr, long double* zi)
{
long double h, d;
if (fabs(yi) < fabs(yr)) {
if (yi == 0.0) {
*zr = xr / yr;
*zi = xi / yr;
}
else {
h = yi / yr;
d = h * yi + yr;
*zr = (xr + h * xi) / d;
*zi = (xi - h * xr) / d;
}
}
else {
h = yr / yi;
d = h * yr + yi;
*zr = (xr * h + xi) / d;
*zi = (xi * h - xr) / d;
}
}
#ifdef DEBUG
static _Lcomplex MyComplexDivide(_Lcomplex numer, _Lcomplex denom) {
long double norm2 =
denom._Val[0] * denom._Val[0] +
denom._Val[1] * denom._Val[1];
_Lcomplex result = { 0 };
result._Val[0] = (
numer._Val[0] * denom._Val[0] +
numer._Val[1] * denom._Val[1]) / norm2;
result._Val[1] = (
numer._Val[1] * denom._Val[0] -
numer._Val[0] * denom._Val[1]) / norm2;
return result;
}
#endif
static _Lcomplex ComplexDivide(_Lcomplex numer, _Lcomplex denom) {
_Lcomplex result = { 0 };
comdiv(
numer._Val[0], numer._Val[1],
denom._Val[0], denom._Val[1],
&result._Val[0], &result._Val[1]);
#ifdef DEBUG
_Lcomplex myResult = MyComplexDivide(numer, denom);
long double delta = fabsl(cabsl(result) - cabsl(myResult));
if (delta > 1.0e-12)
exit(-6);
#endif
return result;
}
static int Deflate(
_Lcomplex coeff[], _Lcomplex zero,
_Lcomplex* fzero, _Lcomplex* fzrdfl,
_Lcomplex zeros[], int i, int* count,
int degree) {
_Lcomplex denom = { 0 };
(*count)++;
*fzero = HornersMethod(coeff, zero, degree);
*fzrdfl = *fzero;
if (i < 1)
return 0;
for (int j = 1; j <= i; j++) {
denom._Val[0] = zero._Val[0] - zeros[j - 1]._Val[0];
denom._Val[1] = zero._Val[1] - zeros[j - 1]._Val[1];
if (cabsl(denom) == 0.0) {
zeros[i] = _LCmulcr(zero, 1.001);
return 1;
}
else
*fzrdfl = ComplexDivide(*fzrdfl, denom);
}
return 0;
}
static void Mueller(
_Lcomplex coeff[], _Lcomplex zeros[],
double epsilon1, double epsilon2,
int degree, int fnreal, int maxIts, int n, int nPrev) {
double eps1 = max(epsilon1, 1.0e-12);
double eps2 = max(epsilon2, 1.0e-20);
int count = 0, i = 0;
_Lcomplex c = { 0 };
_Lcomplex denom = { 0 };
_Lcomplex divdf1 = { 0 };
_Lcomplex divdf2 = { 0 };
_Lcomplex divdf1p = { 0 };
_Lcomplex fzero = { 0 };
_Lcomplex fzr = { 0 };
_Lcomplex fzdfl = { 0 };
_Lcomplex fzrdfl = { 0 };
_Lcomplex fzrprv = { 0 };
_Lcomplex four = _LCbuild(4.0, 0.0);
_Lcomplex h = { 0 };
_Lcomplex hprev = { 0 };
_Lcomplex sqr = { 0 };
_Lcomplex zero = { 0 };
_Lcomplex p5 = _LCbuild(0.5, 0.0);
_Lcomplex zeropp5 = { 0 };
_Lcomplex zeromp5 = { 0 };
_Lcomplex diff = { 0 };
_Lcomplex tadd = { 0 };
_Lcomplex tmul = { 0 };
_Lcomplex umul = { 0 };
_Lcomplex vmul = { 0 };
for (i = nPrev; i < n; i++) {
count = 0;
Label1:
zero = zeros[i];
h = p5;
zeropp5._Val[0] = zero._Val[0] + p5._Val[0];
zeropp5._Val[1] = zero._Val[1] + p5._Val[1];
if (Deflate(
coeff, zeropp5, &fzr, &divdf1p,
zeros, i, &count, degree))
goto Label1;
zeromp5._Val[0] = zero._Val[0] - p5._Val[0];
zeromp5._Val[1] = zero._Val[1] - p5._Val[1];
if (Deflate(
coeff, zeromp5, &fzr, &fzrprv,
zeros, i, &count, degree))
goto Label1;
hprev._Val[0] = -1.0;
hprev._Val[1] = 0.0;
diff._Val[0] = fzrprv._Val[0] - divdf1p._Val[0];
diff._Val[1] = fzrprv._Val[1] - divdf1p._Val[1];
if (cabsl(hprev) == 0)
exit(-2);
divdf1p = ComplexDivide(diff, hprev);
if (Deflate(
coeff, zero, &fzr, &fzrdfl,
zeros, i, &count, degree))
goto Label1;
Label2:
diff._Val[0] = fzrdfl._Val[0] - fzrprv._Val[0];
diff._Val[1] = fzrdfl._Val[1] - fzrprv._Val[1];
if (cabsl(h) == 0)
exit(-3);
divdf1 = ComplexDivide(diff, h);
diff._Val[0] = divdf1._Val[0] - divdf1p._Val[0];
diff._Val[1] = divdf1._Val[1] - divdf1p._Val[1];
tadd._Val[0] = h._Val[0] + hprev._Val[0];
tadd._Val[1] = h._Val[1] + hprev._Val[1];
if (cabsl(tadd) == 0)
exit(-3);
divdf2 = ComplexDivide(diff, tadd);
hprev = h;
divdf1p = divdf1;
tmul = _LCmulcc(h, divdf2);
c._Val[0] = divdf1._Val[0] + tmul._Val[0];
c._Val[1] = divdf1._Val[1] + tmul._Val[1];
tmul = _LCmulcc(c, c);
umul = _LCmulcc(four, fzrdfl);
vmul = _LCmulcc(umul, divdf2);
sqr._Val[0] = tmul._Val[0] - vmul._Val[0];
sqr._Val[1] = tmul._Val[1] - vmul._Val[1];
if (fnreal && sqr._Val[0] < 0.0)
{
sqr._Val[0] = 0.0;
sqr._Val[1] = 0.0;
}
sqr = csqrtl(sqr);
if ((c._Val[0] * sqr._Val[0] + c._Val[1] * sqr._Val[1]) < 0.0) {
denom._Val[0] = c._Val[0] - sqr._Val[0];
denom._Val[1] = c._Val[1] - sqr._Val[1];
}
else {
denom._Val[0] = c._Val[0] + sqr._Val[0];
denom._Val[1] = c._Val[1] + sqr._Val[1];
}
if (cabsl(denom) <= 0.0)
{
denom._Val[0] = 1.0;
denom._Val[1] = 0.0;
}
if (cabsl(denom) == 0)
exit(-4);
tmul = _LCmulcr(fzrdfl, -2.0);
h = ComplexDivide(tmul, denom);
fzrprv = fzrdfl;
zero._Val[0] = zero._Val[0] + h._Val[0];
zero._Val[1] = zero._Val[1] + h._Val[1];
if (count > maxIts)
goto Label4;
Label3:
if (Deflate(
coeff, zero, &fzr, &fzrdfl,
zeros, i, &count, degree))
goto Label1;
if (cabsl(h) < eps1 * cabsl(zero))
goto Label4;
if (max(cabsl(fzr), cabsl(fzdfl)) < eps2)
goto Label4;
if (cabsl(fzrdfl) >= 10.0 * cabsl(fzrprv)) {
h = _LCmulcr(h, 0.5);
zero._Val[0] = zero._Val[0] - h._Val[0];
zero._Val[1] = zero._Val[1] - h._Val[1];
goto Label3;
}
else
goto Label2;
Label4:
zeros[i] = zero;
}
}
int main(void)
{
double epsilon1 = 1.0e-15;
double epsilon2 = 1.0e-15;
int degree = 0, fnreal = 0, i = 0, maxIts = 1000;
int n = 0, nPrev = 0;
while (1) {
_Lcomplex* coeff = NULL;
_Lcomplex* zeros = NULL;
printf_s("Degree (0 to quit) = ");
scanf_s("%d", °ree);
if (degree == 0)
break;
n = degree;
coeff = calloc(degree + 1, sizeof(_Lcomplex));
if (coeff == NULL)
exit(-1);
zeros = calloc(n, sizeof(_Lcomplex));
if (zeros == NULL)
exit(-1);
for (i = degree; i >= 0; i--) {
printf_s("coefficient[%d].real = ", i);
scanf_s("%Lf", &coeff[i]._Val[0]);
printf_s("coefficient[%d].imag = ", i);
scanf_s("%Lf", &coeff[i]._Val[1]);
}
Mueller(
coeff, zeros, epsilon1,
epsilon2, degree, fnreal,
maxIts, n, nPrev);
printf_s("\n");
for (i = 0; i < degree; i++) {
printf_s("zero[%d].real = %17.10e\t", i, zeros[i]._Val[0]);
printf_s("zero[%d].imag = %17.10e\n", i, zeros[i]._Val[1]);
}
printf_s("\n");
for (i = 0; i < degree; i++) {
_Lcomplex func = HornersMethod(coeff, zeros[i], degree);
printf_s("func[%d].real = %17.10e\t", i, func._Val[0]);
printf_s("func[%d].imag = %17.10e\n", i, func._Val[1]);
}
printf_s("\n");
free(coeff);
free(zeros);
}
return 0;
}
Category: Complex Functions
Blog Entry © Sunday, November 9, 2025, by James Pate Williams, Jr. Hydrogenic Wavefunctions, Radial Probability Functions, Distribution Functions, and First Moment Integrals
Blog Entry © Wednesday, October 8, 2025, Quadratic and Cubic Equation Solver by James Pate Williams, Jr.
Blog Entry (c) Wednesday, August 13, 2025, by James Pate Williams, Jr. Exercises from an Online Textbook
#include <complex>
#include <vector>
class CmpLinearAlgebra
{
public:
static void CmpPrintMatrix(
int m, int n,
std::vector<std::vector<std::complex<double>>>& Ac);
static void CmpAddition(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpSubtraction(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpMultiply(
size_t m, size_t n, size_t p,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpAnticommutator(
size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpCommutator(
size_t n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::vector<std::complex<double>>>& B,
std::vector<std::vector<std::complex<double>>>& C);
static void CmpAdjoint(
size_t m, size_t n,
std::vector<std::vector<std::complex<double>>>& Ac,
std::vector<std::vector<std::complex<double>>>& Ad);
static std::complex<double> CmpDeterminant(
bool& failure, int n,
std::vector<std::vector<std::complex<double>>>& A);
static bool CmpGaussianElimination(
int m, int n,
std::vector<std::vector<std::complex<double>>>& A,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpGaussianFactor(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<size_t>& pivot);
static bool CmpGaussianSolution(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpSubstitution(
int m, int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::complex<double>>& b,
std::vector<std::complex<double>>& x,
std::vector<size_t>& pivot);
static bool CmpInverse(
int n, std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& Mi);
static void CmpCharPolyAndAdjoint(
int n,
std::vector<std::vector<std::complex<double>>>& C,
std::vector<std::vector<std::complex<double>>>& I,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& adjoint,
std::vector<std::complex<double>>& a);
static void CmpMatrixKernel(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& X,
size_t& r);
static void CmpMatrixImage(
int m, int n,
std::vector<std::vector<std::complex<double>>>& M,
std::vector<std::vector<std::complex<double>>>& N,
std::vector<std::vector<std::complex<double>>>& X,
int rank);
};
#include <vector>
class DblLinearAlgebra
{
public:
static void DblPrintMatrix(
int m, int n, std::vector<std::vector<double>>& A);
static void DblAddition(
size_t m, size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblSubtraction(
size_t m, size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblMultiply(
size_t m, size_t n, size_t p,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblAnticommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static void DblCommutator(
size_t n,
std::vector<std::vector<double>>& A,
std::vector<std::vector<double>>& B,
std::vector<std::vector<double>>& C);
static double DblDeterminant(
bool& failure, int n,
std::vector<std::vector<double>>& A);
static bool DblGaussianElimination(
int m, int n, std::vector<std::vector<double>>& A,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblGaussianFactor(
int n, std::vector<std::vector<double>>& M,
std::vector<size_t>& pivot);
static bool DblGaussianSolution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblSubstitution(
int n, std::vector<std::vector<double>>& M,
std::vector<double>& b, std::vector<double>& x,
std::vector<size_t>& pivot);
static bool DblInverse(
int n, std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& A);
static void DblCharPolyAndAdjoint(
int n,
std::vector<std::vector<double>>& C,
std::vector<std::vector<double>>& I,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& adjoint,
std::vector<double>& a);
static void DblMatrixKernel(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& X,
size_t& r);
static void DblMatrixImage(
int m, int n,
std::vector<std::vector<double>>& M,
std::vector<std::vector<double>>& N,
std::vector<std::vector<double>>& X,
int rank);
};
// Exercises from "Modern Quantum Chemistry An Introduction to Advanced
// Electronic Structure Theory" by Attila Szabo and Neil S. Ostlund
// https://chemistlibrary.wordpress.com/wp-content/uploads/2015/02/modern-quantum-chemistry.pdf
// Program (c) Tuesday, August 12, 2025 by James Pate Williams, Jr.
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include "DblLinearAlgebra.h"
#include "CmpLinearAlgebra.h"
int main()
{
double AArcb[3][3] = { { 2, 3, -1 }, { 4, 4, -3 }, { -2, 3, -1 } };
double AArso[3][3] = { { 1, 1, 0 }, { 1, 2, 2 }, { 0, 2, -1 } };
double BBrso[3][3] = { { 1, -1, 1 }, { -1, 0, 0 }, { 1, 0, 1} };
double BBr[3][3] = { { 1, -1, 1 }, { -1 , 0, 0 }, { 1, 0, 1 } };
double AAcr[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
double AAci[3][3] = { { 1, 1, 2 }, { 3, 0, 1 }, { 0, 2, 4 } };
double BBcr[3][3] = { { 1, 0, 1 }, { 1 , 1, 0 }, { 0, 1, 1 } };
double BBci[3][3] = { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
int m = 3, n = 3, p = 3;
std::vector<double> br(3);
std::vector<size_t> pivot(3);
std::vector<std::vector<double>> Arcb(3, std::vector<double>(3));
std::vector<std::vector<double>> Arso(3, std::vector<double>(3));
std::vector<std::vector<double>> Brso(3, std::vector<double>(3));
std::vector<std::vector<double>> Br(3, std::vector<double>(3));
std::vector<std::vector<double>> Cr(3, std::vector<double>(3));
std::vector<std::vector<double>> Ai(3, std::vector<double>(3));
std::vector<std::vector<double>> Ari(3, std::vector<double>(3));
std::vector<std::vector<std::complex<double>>> Ac(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Bc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Cc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Dc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Ec(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Fc(3,
std::vector<std::complex<double>>(3));
std::vector<std::vector<std::complex<double>>> Gc(3,
std::vector<std::complex<double>>(3));
for (int i = 0; i < m; i++)
{
for (int j = 0; j < p; j++)
{
Arcb[i][j] = AArcb[i][j];
Arso[i][j] = AArso[i][j];
Brso[i][j] = BBrso[i][j];
Ac[i][j]._Val[0] = AAcr[i][j];
Ac[i][j]._Val[1] = AAci[i][j];
}
}
for (int i = 0; i < p; i++)
{
for (int j = 0; j < n; j++)
{
Br[i][j] = BBr[i][j];
Bc[i][j]._Val[0] = BBcr[i][j];
Bc[i][j]._Val[1] = BBci[i][j];
}
}
DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Br, Cr);
std::cout << "Ar * Br = Cr Conte & de Boor" << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
std::cout << "Ac * Bc = Cc" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
std::cout << std::endl;
// Exercise 1.2
std::cout << "Exercise 1.2 page 5 Commutator" << std::endl;
DblLinearAlgebra::DblCommutator(3, Arso, Brso, Cr);
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
std::cout << "Exercise 1.2 page 5 Anticommutator" << std::endl;
DblLinearAlgebra::DblAnticommutator(3, Arso, Brso, Cr);
DblLinearAlgebra::DblPrintMatrix(3, 3, Cr);
std::cout << std::endl;
CmpLinearAlgebra::CmpAdjoint(3, 3, Cc, Dc);
std::cout << "Exercise 1.3 page 6 Cc adjoint" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Dc);
std::cout << std::endl;
CmpLinearAlgebra::CmpAdjoint(3, 3, Ac, Ec);
CmpLinearAlgebra::CmpAdjoint(3, 3, Bc, Fc);
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Fc, Ec, Gc);
std::cout << "Exercise 1.3 page 6 Bc adjoint * Ac adjoint" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Gc);
std::cout << std::endl;
std::cout << "Ar matrix" << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Arcb);
bool inv = DblLinearAlgebra::DblInverse(n, Arcb, Ai);
std::cout << std::endl;
std::cout << "Ar Conte & de Boor inverse = " << inv << std::endl;
DblLinearAlgebra::DblPrintMatrix(3, 3, Ai);
std::cout << std::endl;
std::cout << "Ar * Ar inverse" << std::endl;
DblLinearAlgebra::DblMultiply(3, 3, 3, Arcb, Ai, Ari);
DblLinearAlgebra::DblPrintMatrix(3, 3, Ari);
std::cout << std::endl;
std::cout << "Ac" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Ac);
std::cout << std::endl;
inv = CmpLinearAlgebra::CmpInverse(3, Ac, Bc);
std::cout << "Ac inverse = " << inv << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Bc);
CmpLinearAlgebra::CmpMultiply(3, 3, 3, Ac, Bc, Cc);
std::cout << std::endl;
std::cout << "Ac * AC inverse" << std::endl;
CmpLinearAlgebra::CmpPrintMatrix(3, 3, Cc);
}
Blog Entry © Sunday, August 10, 2025, First-Order Perturbation Treatment of the Helium Atom by James Pate Williams, Jr.
Approximation of the Ground-State Total Energy of a Beryllium Atom © Sunday, March 30 to Tuesday April 1, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Computer Science
Blog Entry © Thursday, March 27, 2025, by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Lithium (Li, Z = 3) Total Ground-State Energy Numerical Experiments
Blog Entry © Tuesday, March 25, 2025, by James Pate Williams, Jr. Hydrogen Radial Wavefunctions and Related Functions
Live Person-to-Person Tutoring
Blog Entry © Wednesday, November 6, 2024, by James Pate Williams, Jr. Particle in a Finite Spherical Three-Dimensional Potential Well
The Bessel functions are from A Numerical Library in C for Scientists and Engineers (c) 1995 by H. T. Lau, PhD.
Numerical Explorations 