/* NUMroots.c
 *
 * Copyright (C) 1992-2002 Paul Boersma
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 * pb 1999/10/23
 * pb 2002/03/07 GPL
 */

#include "NUM.h"

#define MR 8
#define MT 10
#define MAXIT (MT*MR)
#define MACRO_NUMfindRoot_laguerre(complex,real,epsilon) \
	int iter, j; \
	real abx, abp, abm, err; \
	complex dx, x1, b, d, f, g, h, sq, gp, gm, g2; \
	static real frac [MR + 1] = { 0.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0 }; \
	for (iter = 1; iter <= MAXIT; iter ++) { \
		b = a [degree]; \
		err = complex##_abs (b); \
		d = f = complex##_create (0.0, 0.0); \
		abx = complex##_abs (*x); \
		for (j = degree - 1; j >= 0; j --) { \
			f = complex##_add (complex##_mul (*x, f), d); \
			d = complex##_add (complex##_mul (*x, d), b); \
			b = complex##_add (complex##_mul (*x, b), a [j]); \
			err = complex##_abs (b) + abx * err; \
		} \
		err *= epsilon; \
		if (complex##_abs (b) <= err) return 1; \
		g = complex##_div (d, b); \
		g2 = complex##_mul (g, g); \
		h = complex##_sub (g2, complex##_rmul (2.0, complex##_div (f, b))); \
		sq = complex##_sqrt (complex##_rmul ((real) (degree - 1), complex##_sub (complex##_rmul ((real) degree, h), g2))); \
		gp = complex##_add (g, sq); \
		gm = complex##_sub (g, sq); \
		abp = complex##_abs (gp); \
		abm = complex##_abs (gm); \
		if (abp < abm) gp = gm; \
		dx = (abp > abm ? abp : abm) > 0.0 ? \
			  complex##_div (complex##_create ((real) degree, 0.0), gp) : \
			  complex##_rmul (exp (log (1 + abx)), complex##_create (cos ((real) iter), sin ((real) iter))); \
		x1 = complex##_sub (*x, dx); \
		if (x -> re == x1. re && x -> im == x1. im) return 1; \
		if (iter % MT) *x = x1; \
		else *x = complex##_sub (*x, complex##_rmul (frac [iter / MT], dx)); \
	} \
	return 0;   /* Too many iterations. */

int NUMfindRoot_laguerre_f (fcomplex a [], long degree, fcomplex *x) {
	MACRO_NUMfindRoot_laguerre (fcomplex, float, 1.0e-7)
}
int NUMfindRoot_laguerre_d (dcomplex a [], long degree, dcomplex *x) {
	MACRO_NUMfindRoot_laguerre (dcomplex, double, 1.0e-15)
}
#undef MR
#undef MT
#undef MAXIT

#define MACRO_NUMfindRoots(complex,real,letter,epsilon) \
	int i, j, jj; \
	static long maximumDegree = 0; \
	static complex *ad; \
	if (degree < 1) return 1; \
	if (ad == NULL || degree > maximumDegree) { \
		NUM##letter##cvector_free (ad, 0); \
		ad = NUM##letter##cvector (0, degree); \
		if (ad == NULL) return 0;   /* Out of memory. */ \
		maximumDegree = degree; \
	} \
	NUM##letter##cvector_copyElements (a, ad, 0, degree); \
	for (j = degree; j >= 1; j --) { \
		complex root = complex##_create (0.0, 0.0), b; \
		NUMfindRoot_laguerre_##letter (ad, j, & root); \
		/* \
		 * Detect any real-valued root. \
		 */ \
		if (fabs (root. im) <= 2.0 * epsilon * fabs (root. re)) { \
			root. im = 0.0; \
		} \
		roots [j] = root; \
		b = ad [j]; \
		for (jj = j - 1; jj >= 0; jj --) { \
			complex temp = ad [jj]; \
			ad [jj] = b; \
			b = complex##_add (complex##_mul (root, b), temp); \
		} \
	} \
	/* Polish. */ \
	for (j = 1; j <= degree; j ++) { \
		NUMfindRoot_laguerre_##letter (a, degree, & roots [j]); \
	} \
	/* Sort. */ \
	for (j = 2; j <= degree; j ++) { \
		complex root = roots [j]; \
		for (i = j - 1; i >= 1; i --) { \
			if (roots [i]. re <= root. re) break; \
			roots [i + 1] = roots [i]; \
		} \
		roots [i + 1] = root; \
	} \
	return 1;

int NUMfindRoots_f (fcomplex a [], long degree, fcomplex roots []) {
	MACRO_NUMfindRoots (fcomplex, float, f, 2.0e-6)
}
int NUMfindRoots_d (dcomplex a [], long degree, dcomplex roots []) {
	MACRO_NUMfindRoots (dcomplex, double, d, 2.0e-14)
}

/* End of file NUMroots.c */