Category: Root Finding Algorithms

/*
* 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", &degree);

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