Kody źródłowe/Metoda Laguerre'a

C

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

double horner(double x, int n, double *P)
{
	double result = P[n];
	for (int i = n - 1; i >= 0; i--)
	{
		result = result*x + P[i];
	}
	return result;
}

double dhorner(double x, int n, double *P)
{
	double result = P[n] * n;
	for (int i = n - 1; i >= 1; i--)
	{
		result = result*x + P[i]*i;
	}
	return result;
}

double d2horner(double x, int n, double *P)
{
	double result = P[n] * n * (n-1);
	for (int i = n - 1; i >= 2; i--)
	{
		result = result*x + P[i] * i *(i-1);
	}
	return result;
}

/* complex
mul: (a + bi)(c + di) = (ac - bd) + (bc + ad)i
div (a + bi)/(c + di) = ((ac + bd) + (bc - ad)i)/(c^2 + d^2)
*/
void hornerZ(double &resRe, double &resIm, double xRe, double xIm, int n, double *P)
{
	resRe = P[n];
	resIm = 0;
	for (int i = n - 1; i >= 0; i--)
	{
		double origResRe = resRe;
		resRe = (resRe*xRe - resIm*xIm);
		resIm = resIm*xRe + origResRe*xIm;
		resRe += P[i];
	}
}

void dhornerZ(double &resRe, double &resIm, double a, double b, int n, double *P)
{
	resRe = P[n] * n;
	resIm = 0;
	for (int i = n - 1; i >= 1; i--)
	{
		double origResRe = resRe;
		resRe = (resRe*a - resIm*b);
		resIm = resIm*a + origResRe*b;
		resRe += P[i] * i;
	}
}

void d2hornerZ(double &resRe, double &resIm, double a, double b, int n, double *P)
{
	resRe = P[n] * n * (n - 1);
	resIm = 0;
	for (int i = n - 1; i >= 2; i--)
	{
		double origResRe = resRe;
		resRe = (resRe*a - resIm*b);
		resIm = resIm*a + origResRe*b;
		resRe += P[i] * i *(i - 1);
	}
}

double Laguerre(double x, int n, double *P)
{
	double maxA = 0,absA;
	for (int i = 0; i <= n; i++)
	{
		absA = fabs(P[i]);
		if (absA > maxA) maxA = absA;
	}
	while (true)
	{
		double h = horner(x, n, P);
		if (fabs(h) < fabs(DBL_EPSILON*maxA)) break;
		double dh = dhorner(x, n, P);
		double G = dh / h;
		double H = G*G - d2horner(x, n, P) / h;
		double v = (n - 1)*(n*H - G*G);
		if (v < 0) return NAN;
		double s = sqrt(v);
		double denom1 = G + s;
		double denom2 = G - s;
		double denom;
		if (fabs(denom1) > fabs(denom2))
			denom = denom1;
		else
			denom = denom2;
		double a = n / denom;
		if (fabs(a) < 2 * DBL_EPSILON*fabs(x))
		{
			x = x - 0.5*a;
			break;
		}
		x = x - a;
	}
	return x;
}

int LaguerreZ(double &resRe, double &resIm, double xRe, double xIm, int n, double *P)
{
	int loopCnt = 0;
	double maxA = 0, absA;
	for (int i = 0; i <= n; i++)
	{
		absA = fabs(P[i]);
		if (absA > maxA) maxA = absA;
	}
	while (true)
	{
		double hRe, hIm;
		hornerZ(hRe, hIm, xRe, xIm, n, P);
		if (sqrt(hRe*hRe + hIm*hIm) < DBL_EPSILON*maxA) break;
		double dhRe, dhIm;
		dhornerZ(dhRe, dhIm, xRe, xIm, n, P);
		double GRe, GIm;
		double denomh = hRe*hRe + hIm*hIm;
		GRe = (dhRe*hRe + dhIm*hIm) / denomh;
		GIm = (dhIm*hRe - dhRe*hIm) / denomh;
		double G2Re, G2Im;
		G2Re = GRe*GRe - GIm*GIm;
		G2Im = 2 * GRe * GIm;
		double d2hRe, d2hIm;
		d2hornerZ(d2hRe, d2hIm, xRe, xIm, n, P);
		double GRe2, GIm2;
		GRe2 = (d2hRe*hRe + d2hIm*hIm) / denomh;
		GIm2 = (d2hIm*hRe - d2hRe*hIm) / denomh;
		double HRe, HIm;
		HRe = G2Re - GRe2;
		HIm = G2Im - GIm2;
		double vRe, vIm;
		vRe = (n - 1)*(n*HRe - G2Re);
		vIm = (n - 1)*(n*HIm - G2Im);
		double modv = sqrt(vRe*vRe + vIm*vIm);
		double sRe, sIm;
		sRe = sqrt((modv + vRe) / 2);
		sIm = sqrt((modv - vRe) / 2);
		if (vIm < 0) sIm = -sIm;
		double denom1Re, denom1Im;
		denom1Re = GRe + sRe;
		denom1Im = GIm + sIm;
		double denom2Re, denom2Im;
		denom2Re = GRe - sRe;
		denom2Im = GIm - sIm;
		double denomRe, denomIm;
		if (denom1Re*denom1Re + denom1Im*denom1Im > denom2Re*denom2Re + denom2Im*denom2Im)
		{
			denomRe = denom1Re;
			denomIm = denom1Im;
		}
		else
		{
			denomRe = denom2Re;
			denomIm = denom2Im;
		}
		double aRe, aIm;
		double denoma = denomRe*denomRe + denomIm*denomIm;
		aRe = n*denomRe / denoma;
		aIm = -n*denomIm / denoma;
		if (sqrt(aRe*aRe + aIm*aIm) < 2 * DBL_EPSILON*sqrt(xRe*xRe + xIm*xIm))
		{
			xRe -= 0.5*aRe;
			xIm -= 0.5*aIm;
			break;
		}
		xRe -= aRe;
		xIm -= aIm;
		loopCnt++;
	}
	resRe = xRe;
	resIm = xIm;
	return loopCnt;
}

const int degP = 6;
double P[degP + 1] = { 36,-48,-3,17,2,-5,1 }; //(x^2 + 3·x + 3)(x-1)(x-2)(x-2)(x-3)
int main()
{	
	double a,b;	
	LaguerreZ(a, b, 0, 0, degP, P);	
}