PRIVATE: Polynomials

Files
file	fnft__poly_chirpz.h
	Fast evaluation of polynomials using the Chirp Z-transform.
file	fnft__poly_eval.h
	Evaluation of polynomials and their derivatives using Horner's method.
file	fnft__poly_fmult.h
	Fast multiplication of n polynomials.
file	fnft__poly_roots_fasteigen.h
	Fast root finding of polynomials.
file	fnft__poly_roots_fftgridsearch.h
	Approximation of polynomial roots on the unit circle using gridding.
file	fnft__poly_specfact.h
	Spectral factorization of polynomials.

Functions
FNFT_INT	fnft__poly_chirpz (const FNFT_UINT deg, FNFT_COMPLEX const *const p, const FNFT_COMPLEX A, const FNFT_COMPLEX W, const FNFT_UINT M, FNFT_COMPLEX *const result)
	Fast evaluation of a polynomial on a spiral in the complex plane. More...
FNFT_INT	fnft__poly_eval (const FNFT_UINT deg, FNFT_COMPLEX const *const p, const FNFT_UINT nz, FNFT_COMPLEX *const z)
	Evaluation of polynomials. More...
FNFT_INT	fnft__poly_evalderiv (const FNFT_UINT deg, FNFT_COMPLEX const *const p, const FNFT_UINT nz, FNFT_COMPLEX *const z, FNFT_COMPLEX *const deriv)
	Evaluation of polynomials and their derivatives. More...
FNFT_INT	fnft__poly_fmult_two_polys_len (const FNFT_UINT deg)
	Length of the FFTs used by fnft__poly_fmult_two_polys. More...
FNFT_INT	fnft__poly_fmult_two_polys (const FNFT_UINT deg, FNFT_COMPLEX const *const p1, FNFT_COMPLEX const *const p2, FNFT_COMPLEX *const result, fnft__fft_wrapper_plan_t plan_fwd, fnft__fft_wrapper_plan_t plan_inv, FNFT_COMPLEX *const buf0, FNFT_COMPLEX *const buf1, FNFT_COMPLEX *const buf2, const FNFT_UINT mode)
	Multiplies two polynomials. More...
FNFT_INT	fnft__poly_fmult_two_polys2x2 (const FNFT_UINT deg, FNFT_COMPLEX const *const p1_11, const FNFT_UINT p1_stride, FNFT_COMPLEX const *const p2_11, const FNFT_UINT p2_stride, FNFT_COMPLEX *const result_11, const FNFT_UINT result_stride, fnft__fft_wrapper_plan_t plan_fwd, fnft__fft_wrapper_plan_t plan_inv, FNFT_COMPLEX *const buf0, FNFT_COMPLEX *const buf1, FNFT_COMPLEX *const buf2, const FNFT_UINT mode_offset)
	Multiplies two 2x2 matrices of polynomials. More...
FNFT_UINT	fnft__poly_fmult_numel (const FNFT_UINT deg, const FNFT_UINT n)
	Number of elements that the input p to fnft__poly_fmult should have. More...
FNFT_INT	fnft__poly_fmult (FNFT_UINT *const d, FNFT_UINT n, FNFT_COMPLEX *const p, FNFT_INT *const W_ptr)
	Fast multiplication of multiple polynomials of same degree. More...
FNFT_UINT	fnft__poly_fmult2x2_numel (const FNFT_UINT deg, const FNFT_UINT n)
	Number of elements that the inputs p and result to fnft__poly_fmult2x2 should have. More...
FNFT_INT	fnft__poly_fmult2x2 (FNFT_UINT *d, FNFT_UINT n, FNFT_COMPLEX *const p, FNFT_COMPLEX *const result, FNFT_INT *const W_ptr)
	Fast multiplication of multiple 2x2 matrix-valued polynomials of the same degree. More...
FNFT_INT	fnft__poly_roots_fasteigen (const FNFT_UINT deg, FNFT_COMPLEX const *const p, FNFT_COMPLEX *const roots)
	Fast computation of polynomial roots. More...
FNFT_INT	fnft__poly_roots_fftgridsearch (const FNFT_UINT deg, FNFT_COMPLEX const *const p, FNFT_UINT *const M_ptr, FNFT_REAL const *const PHI, FNFT_COMPLEX *const roots)
	Unit circle roots of a polynomial via grid search. More...
FNFT_INT	fnft__poly_roots_fftgridsearch_paraherm (const FNFT_UINT deg, FNFT_COMPLEX const *const p, FNFT_UINT *const M_ptr, FNFT_REAL const *const PHI, FNFT_COMPLEX *const roots)
	Unit circle roots of a parahermitian Laurent polynomial via grid search. More...
FNFT_INT	fnft__poly_specfact (const FNFT_UINT deg, FNFT_COMPLEX const *const poly, FNFT_COMPLEX *const result, const FNFT_UINT oversampling_factor, const FNFT_INT kappa)
	Spectral factorization of polynomial. More...

Detailed Description

Function Documentation

fnft__poly_chirpz()

FNFT_INT fnft__poly_chirpz	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	p,
		const FNFT_COMPLEX	A,
		const FNFT_COMPLEX	W,
		const FNFT_UINT	M,
		FNFT_COMPLEX *const	result
	)

Fast evaluation of a polynomial on a spiral in the complex plane.

This routine implements the Chirp Z-transform. Given a polynomial

\[ p(z)=p_0+p_1 z^1+p_2 z^2+...+p_{deg} z^{deg} \]

and complex numbers A and W, this routine evaluates the polynomial \( p(z) \) at the M points \( z=1/w_k \), where

\[ w_k=AW^{-m}, \quad m=0,1,\dots,M-1, \]

using only \( O\{(N+M)\log(N+M)\} \) floating point operations.

See also

https://doi.org/10.1109/TAU.1969.1162034

Parameters

[in]	deg	Degree of the polynomial
[in]	p	Array containing the deg+1 coefficients of the polynomial in descending order (i.e., \( p_{deg}, p_{deg-1}, \dots, p_1, p_0 \)).
[in]	A	First constant defining the spiral at which the polynomial will be evaluated.
[in]	W	Second constant defining the spiral at which the polynomial will be evaluated
[in]	M	Number of points at which the polynomial will be evaluated.
[out]	result	Array of M points. Will be filled with \( p(1/w_1),...,p(1/w_M) \).

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_eval()

FNFT_INT fnft__poly_eval	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	p,
		const FNFT_UINT	nz,
		FNFT_COMPLEX *const	z
	)

Evaluation of polynomials.

Evaluates a polynomial

\[ p(z)=p_0+p_1 z^1+p_2 z^2+...+p_{deg} z^{deg} \]

at nz points

\[ z=z_1,...,z_{nz} \]

in the complex plane using Horner's method.

See also

http://cnx.org/contents/hnzXCQRy@7/Horners-Method-for-Evaluating-

Parameters

[in]	deg	Degree of the polynomial
[in]	p	Array containing the deg+1 coefficients of the polynomial in descending order (i.e., \( p_{deg}, p_{deg-1}, \dots, p_1, p_0 \)).
[in]	nz	Number of points z at which the polynomial should be evaluated
[in,out]	z	Array of nz points \(z_1,\dots,z_{nz}\) at which the polynomial should be evaluated. These values will be overwritten with the corresponding values \(p(z_1),\dots,p(z_{nz})\) of the polynomial.

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_evalderiv()

FNFT_INT fnft__poly_evalderiv	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	p,
		const FNFT_UINT	nz,
		FNFT_COMPLEX *const	z,
		FNFT_COMPLEX *const	deriv
	)

Evaluation of polynomials and their derivatives.

Evaluates a polynomial

\[ p(z)=p_0+p_1 z^1+p_2 z^2+...+p_{deg} z^{deg} \]

and its derivative

\[ p'(z)=\frac{dp}{dz}(z)=p_1+2p_2 z^1+...+deg~p_{deg} z^{deg-1} \]

at nz points

\[ z=z_1,...,z_{nz} \]

in the complex plane using Horner's method.

See also

http://cnx.org/contents/hnzXCQRy@7/Horners-Method-for-Evaluating-

Parameters

[in]	deg	Degree of the polynomial
[in]	p	Array containing the deg+1 coefficients of the polynomial in descending order (i.e., \( p_{deg}, p_{deg-1}, \dots, p_1, p_0 \)).
[in]	nz	Number of points z at which the polynomial should be evaluated
[in,out]	z	Array of nz points \(z_1,\dots,z_{nz}\) at which the polynomial and its derivative should be evaluated. These values will be overwritten with the corresponding values \(p(z_1),\dots,p(z_{nz})\) of the polynomial.
[out]	deriv	Array of length nz. It will be filled with the values \(p'(z_1),\dots,p'(z_{nz})\) of the derivative.

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_fmult()

FNFT_INT fnft__poly_fmult	(	FNFT_UINT *const	d,
		FNFT_UINT	n,
		FNFT_COMPLEX *const	p,
		FNFT_INT *const	W_ptr
	)

Fast multiplication of multiple polynomials of same degree.

Fast multiplication of n polynomials of degree d. Their coefficients are stored in the array p and will be overwritten. If W_ptr != NULL, the result has been normalized by a factor 2^W. Upon exit, W has been stored in *W_ptr.

Parameters

[in,out]	d	Upon entry, degree of the input polynomials. Upon exit, degree of their product.
[in]	n	Number of polynomials.
[in,out]	p	Complex valued array with m entries, where m is determined using fnft__poly_fmult_numel. Upon entry, the first (*d+1)*n elements of this array contain the coefficients of the polynomials. Upon exit, the first *d+1 elements contain the result.
[in]	W_ptr	Pointer to normalization flag.

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_fmult2x2()

FNFT_INT fnft__poly_fmult2x2	(	FNFT_UINT *	d,
		FNFT_UINT	n,
		FNFT_COMPLEX *const	p,
		FNFT_COMPLEX *const	result,
		FNFT_INT *const	W_ptr
	)

Fast multiplication of multiple 2x2 matrix-valued polynomials of the same degree.

Fast multiplication of n 2x2 matrix-valued polynomials of degree d. Their coefficients are stored in the array p and will be overwritten. If W_ptr != NULL, the result has been normalized by a factor 2^W. Upon exit, W has been stored in *W_ptr.

Parameters

[in]	d	Pointer to a FNFT_UINT containing the degree of the polynomials.
[in]	n	Number of 2x2 matrix-valued polynomials.
[in,out]	p	Complex valued array which holds the coefficients of the polynomials being multiplied. Should be of length m*(*d+1), where m is obtained using fnft__poly_fmult2x2_numel. WARNING: p is overwritten.
[out]	result	Complex valued array that holds the result of the multiplication. Should be of the same size as p.
[in]	W_ptr	Pointer to normalization flag.

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_fmult2x2_numel()

FNFT_UINT fnft__poly_fmult2x2_numel	(	const FNFT_UINT	deg,
		const FNFT_UINT	n
	)

Number of elements that the inputs p and result to fnft__poly_fmult2x2 should have.

Specifies how much memory (in number of elements) the user needs to allocate for the inputs p and result of the routine fnft__poly_fmult2x2.

Parameters

[in]	deg	Degree of the polynomials
[in]	n	Number of polynomials

Returns

A number m. The inputs p and result to fnft__poly_fmult2x2 should be a array with m entries.

fnft__poly_fmult_numel()

FNFT_UINT fnft__poly_fmult_numel	(	const FNFT_UINT	deg,
		const FNFT_UINT	n
	)

Number of elements that the input p to fnft__poly_fmult should have.

Specifies how much memory (in number of elements) the user needs to allocate for the input p of the routine fnft__poly_fmult.

Parameters

[in]	deg	Degree of the polynomials
[in]	n	Number of polynomials

Returns

A number m. The input p to fnft__poly_fmult should be a array with m entries.

fnft__poly_fmult_two_polys()

FNFT_INT fnft__poly_fmult_two_polys	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	p1,
		FNFT_COMPLEX const *const	p2,
		FNFT_COMPLEX *const	result,
		fnft__fft_wrapper_plan_t	plan_fwd,
		fnft__fft_wrapper_plan_t	plan_inv,
		FNFT_COMPLEX *const	buf0,
		FNFT_COMPLEX *const	buf1,
		FNFT_COMPLEX *const	buf2,
		const FNFT_UINT	mode
	)

Multiplies two polynomials.

Fast multiplication of two polynomials using the FFT.

Parameters

	deg	Degree of the polynomials.
[in]	p1	Array of coefficients for the first polynomial. If p1 is NULL, then buf1 is expected to contain the FFT of the zero-padded coefficient vector of the first polynomial instead. The purpose of this feature is that the buf2 returned by another run of this function can be reused, thus saving an FFT.
[in]	p2	Array of coefficients for the second polynomial. If p2 is NULL, then buf2 is expected to contain the FFT of the zero-padded coefficient vector of the second polynomial instead.
[out]	result	Array in which the coefficients of the result are stored.
	plan_fwd	Plan generated by fnft__fft_wrapper_create_plan for a forward FFT of the length returned by fnft__poly_fmult_two_polys_len.
	plan_inv	Plan generated by fnft__fft_wrapper_create_plan for an inverse FFT of the length returned by fnft__poly_fmult_two_polys_len.
[in,out]	buf0	Buffer of the same length as the FFTs. Must be allocated allocated and free by the user using fnft__fft_wrapper_malloc and fnft__fft_wrapper_free, respectively.
[in,out]	buf1	See buf0.
[in,out]	buf2	See buf0. Upon exit, buf2 contains the FFT of the zero-padded coefficient vector of the second polynomial.
[in]	mode	If mode==0, then the product of the two polynomials is stored in result. This is the standard mode. The following modes are useful to speed-up the multiplication of polynomial matrices. If mode==1, then the product is added to the values currently stored in result. If mode==3, then the FFT of the product of the (zero-padded) polynomials is stored in result. If mode==4, then the FFT of the two (zero-padded) polynomials is added to the values currently stored in result, an inverse FFT is applied to these values, and the outcome is finally stored in result. IMPORTANT: The modes 3 and 4 require that the array result is large enough to store an FFT of the size given above.

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_fmult_two_polys2x2()

FNFT_INT fnft__poly_fmult_two_polys2x2	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	p1_11,
		const FNFT_UINT	p1_stride,
		FNFT_COMPLEX const *const	p2_11,
		const FNFT_UINT	p2_stride,
		FNFT_COMPLEX *const	result_11,
		const FNFT_UINT	result_stride,
		fnft__fft_wrapper_plan_t	plan_fwd,
		fnft__fft_wrapper_plan_t	plan_inv,
		FNFT_COMPLEX *const	buf0,
		FNFT_COMPLEX *const	buf1,
		FNFT_COMPLEX *const	buf2,
		const FNFT_UINT	mode_offset
	)

Multiplies two 2x2 matrices of polynomials.

Fast multiplication of two 2x2 matrices of polynomials using the FFT.

Parameters

	deg	Degree of the polynomials.
[in]	p1_11	Array of deg+1 coefficients for the upper left polynomial p1_11(z) of the first matrix p1(z).
[in]	p1_stride	The deg+1 coefficients of the other polynomials p1_12(z), p1_21(z) and p1_22(z) in the first matrix are expected to be stored at p1_11+p1_stride, p1_11+2*p1_stride and p1+3*p1_stride, respectively.
[in]	p2_11	Array of deg+1 coefficients for the upper left polynomial p2_11(z) of the second matrix p2(z).
[in]	p2_stride	The deg+1 coefficients of the other polynomials p2_12(z), p2_21(z) and p2_22(z) in the first matrix are expected to be stored at p2_11+p2_stride, p2_11+2*p2_stride and p2+3*p2_stride, respectively.
[out]	result_11	Array in which the 2*deg+1 coefficients of the upper left polynomial result_11(z) of the product result(z)=p1(z)*p2(z) will be stored.
[in]	result_stride	The 2*deg+1 coefficients of the other polynomials result_12(z), result_21(z) and result_22(z) in the resulting matrix are will be stored at result_11+result_stride, result_11+2*result_stride and result+3*result_stride, respectively.
	plan_fwd	Plan generated by fnft__fft_wrapper_create_plan for a forward FFT of the length returned by fnft__poly_fmult_two_polys_len.
	plan_inv	Plan generated by fnft__fft_wrapper_create_plan for an inverse FFT of the length returned by fnft__poly_fmult_two_polys_len.
[in,out]	buf0	Buffer of the same length as the FFTs. Must be allocated allocated and free by the user using fnft__fft_wrapper_malloc and fnft__fft_wrapper_free, respectively.
[in,out]	buf1	See buf0.
[in,out]	buf2	See buf0.
[in]	mode_offset	Set to 2 if the arrays result_11+n*result_stride, where n=0,1,2,3, are long enough to store the FFT's. This allows improve performance by storing some intermediate values there. Otherwise, set to 0.

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_fmult_two_polys_len()

FNFT_INT fnft__poly_fmult_two_polys_len

(

const FNFT_UINT

deg

)

Length of the FFTs used by fnft__poly_fmult_two_polys.

Parameters

deg	The degree of the polynomials to be multiplied.

Returns

The length of the (inverse) FFT (in number of complex elements) required by fnft__poly_fmult_two_polys.

fnft__poly_roots_fasteigen()

FNFT_INT fnft__poly_roots_fasteigen	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	p,
		FNFT_COMPLEX *const	roots
	)

Fast computation of polynomial roots.

This routine compute the roots of a polynomial

\[ p(z)=p_0+p_1 z^1+p_2 z^2+...+p_{deg} z^{deg} \]

using only \( O\{ deg^2 \}\) floating point operations.

See also

https://arxiv.org/abs/1611.02435v2

Parameters

[in]	deg	Degree of the polynomial
[in]	p	Array containing the deg+1 coefficients of the polynomial in descending order (i.e., \( p_{deg}, p_{deg-1}, \dots, p_1, p_0 \)).
[out]	roots	Array of deg+1 points. Will be filled with the roots of \( p(z) \).

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_roots_fftgridsearch()

FNFT_INT fnft__poly_roots_fftgridsearch	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	p,
		FNFT_UINT *const	M_ptr,
		FNFT_REAL const *const	PHI,
		FNFT_COMPLEX *const	roots
	)

Unit circle roots of a polynomial via grid search.

This routine approximates the unit circle roots of a polynomial

\[ p(z)=p_0 + p_{1}z^1 + \dots p_{deg} + p_N z^{deg}, \]

where the star denote the complex conjugate. (That is, it looks for the solutions of \( p(z)=0 \) that satisfy \( |z|=1 \).) The polynomial \( p(z) \) is evaluated on the grid \( z=r_k e^{j\theta_m} \), where

\[ \theta_m=\Phi_0 + m \epsilon, \ \ \ \ m=0,1,\dots,M-1, \ \epsilon=\frac{\Phi_1 - \Phi_0}{M-1}\]

and

\[ r_k = 1 - k\epsilon, \ \ \ \ k=-1,0,1. \]

The evaluations are carried out using the Chirp transform. The grid points on the unit circle where the absolute value of the polynomial is lower than or equal to the absolute value at the eight neighboring grid points are candidates for a root. (This is the first stage of the Lindsey-Fox algorithm.) The polynomial is then approximated locally around the center grid point and the root of this linerization is computed. The root of the linearization is kept if it is not too far from the grid point. The routine gurantees only linear convergence.

See also

https://doi.org/10.1109/MSP.2003.1253552

poly_roots_fftgridsearch_paraherm

poly_chirpz

Parameters

[in]	deg	The degree of the polynomial.
[in]	p	Array containing the deg+1 coefficients of the Laurent polynomial in descending order (i.e., \( p_{deg}, p_{deg-1}, \dots, p_{1}, p_{0} \)).
[in,out]	M_ptr	Upon entry, *M_ptr contains the desired value for the number of points M. Upon return, *M_ptr has been overwritten with the number of detected roots.
[in]	PHI	Array with two entries, \( \Phi_0 \) and \( \Phi_1 \). The first value should be lower than the second one.
[out]	roots	Array of M points. Will be filled with the detected roots. (The number of detected roots is stored in *M_ptr.)

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_roots_fftgridsearch_paraherm()

FNFT_INT fnft__poly_roots_fftgridsearch_paraherm	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	p,
		FNFT_UINT *const	M_ptr,
		FNFT_REAL const *const	PHI,
		FNFT_COMPLEX *const	roots
	)

Unit circle roots of a parahermitian Laurent polynomial via grid search.

This routine approximates the unit circle roots of a para-hermitian Laurent polynomial

\[ p(z)=p_{-N} + p_{-N+1} + \dots p_{N-1} + p_N, \ \ \ \ p(z)=p^*(1/z^*), \]

where the star denote the complex conjugate. (That is, it looks for the solutions of \( p(z)=0 \) that satisfy \( |z|=1 \).) The Laurent polynomial \( p(z) \) is evaluated on the grid \( z=e^{j\theta_m} \), where

\[ \theta_m=\Phi_0 + m \epsilon, \ \ \ \ m=0,1,\dots,M-1, \ \epsilon=\frac{\Phi_1 - \Phi_0}{M-1}.\]

The evaluations are carried out using the Chirp transform. Since the Laurent polynomial is para-hermitian, it will be real on the unit circle. If the sign of the Laurent polynomial changes between two consequtive points on the grid, a root is detected. The position of the roots is determined using linear interpolation.

See also

poly_roots_fftgridsearch

poly_chirpz

Parameters

[in]	deg	This should be the value of \( 2N+1 \).
[in]	p	Array containing the deg+1 coefficients of the Laurent polynomial in descending order (i.e., \( p_{N}, p_{N-1}, \dots, p_{-N+1}, p_{-N} \)).
[in,out]	M_ptr	Upon entry, *M_ptr contains the desired value for the number of points M. Upon return, *M_ptr has been overwritten with the number of detected roots.
[in]	PHI	Array with two entries, \( \Phi_0 \) and \( \Phi_1 \). The first value should be lower than the second one.
[out]	roots	Array of M points. Will be filled with the detected roots. (The number of detected roots is stored in *M_ptr.)

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.

fnft__poly_specfact()

FNFT_INT fnft__poly_specfact	(	const FNFT_UINT	deg,
		FNFT_COMPLEX const *const	poly,
		FNFT_COMPLEX *const	result,
		const FNFT_UINT	oversampling_factor,
		const FNFT_INT	kappa
	)

Spectral factorization of polynomial.

Computes the spectral factor of given polynomial with no zeros inside the unit circle.

Parameters

[in]	deg	Degree of the polynomial.
[in]	poly	Array of coefficients for the polynomial.
[out]	result	Array in which the coefficients of the result are stored.
[in]	oversampling_factor	Oversampling factor is used to improve conditioning of the spectral factorization problem.
[in]	kappa	kappa =+1 for the focusing nonlinear Schroedinger equation, =-1 for the defocusing one.

Returns

FNFT_SUCCESS or one of the FNFT_EC_... error codes defined in fnft_errwarn.h.