From daf5550c7808e77884f4eb405fe14566cd8e1bdc Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 16 Oct 2012 10:27:05 +0200 Subject: [PATCH] more helper functions, and some documenting --- doc/source/fmpcb.rst | 5 + doc/source/fmpcb_poly.rst | 166 ++++++++++++++++++++++++ doc/source/index.rst | 2 + fmpcb_poly.h | 4 + fmpcb_poly/refine_roots_durand_kerner.c | 8 -- fmpcb_poly/root_inclusion.c | 5 - fmpcb_poly/set2_fmpq_poly.c | 10 ++ fmpcb_poly/set2_fmprb_poly.c | 18 +++ fmpcb_poly/validate_roots.c | 14 -- 9 files changed, 205 insertions(+), 27 deletions(-) create mode 100644 doc/source/fmpcb.rst create mode 100644 doc/source/fmpcb_poly.rst diff --git a/doc/source/fmpcb.rst b/doc/source/fmpcb.rst new file mode 100644 index 00000000..8ddebee6 --- /dev/null +++ b/doc/source/fmpcb.rst @@ -0,0 +1,5 @@ +fmpcb.h -- complex numbers represented as pairs of real balls +=============================================================================== + +documentation tbd + diff --git a/doc/source/fmpcb_poly.rst b/doc/source/fmpcb_poly.rst new file mode 100644 index 00000000..bc89e2b9 --- /dev/null +++ b/doc/source/fmpcb_poly.rst @@ -0,0 +1,166 @@ +fmpcb_poly.h -- polynomials of complex numbers +=============================================================================== + +Types, macros and constants +------------------------------------------------------------------------------- + +.. type:: fmpcb_poly_struct + +.. type:: fmpcb_poly_t + + Contains a pointer to an array of coefficients (coeffs), the used + length (length), and the allocated size of the array (alloc). + + An *fmpcb_poly_t* is defined as an array of length one of type + *fmpcb_poly_struct*, permitting an *fmpcb_poly_t* to + be passed by reference. + +Memory management +------------------------------------------------------------------------------- + +.. function:: void fmpcb_poly_init(fmpcb_poly_t poly) + + Initializes the polynomial for use, setting it to the zero polynomial. + +.. function:: void fmpcb_poly_clear(fmpcb_poly_t poly) + + Clears the polynomial, deallocating all coefficients and the + coefficient array. + +.. function:: void fmpcb_poly_fit_length(fmpcb_poly_t poly, long len) + + Makes sures that the coefficient array of the polynomial contains at + least *len* initialized coefficients. + +.. function:: void _fmpcb_poly_set_length(fmpcb_poly_t poly, long len) + + Directly changes the length of the polynomial, without allocating or + deallocating coefficients. The value shold not exceed the allocation length. + +.. function:: void _fmpcb_poly_normalise(fmpcb_poly_t poly) + + Strips any trailing coefficients which are identical to zero. + +.. function:: void fmpcb_poly_swap(fmpcb_poly_t poly1, fmpcb_poly_t poly2) + + Swaps *poly1* and *poly2* efficiently. + + +Basic properties and manipulation +------------------------------------------------------------------------------- + +.. function:: long fmpcb_poly_length(const fmpcb_poly_t poly) + + Returns the length of the polynomial. + +.. function:: void fmpcb_poly_zero(fmpcb_poly_t poly) + + Sets *poly* to the zero polynomial. + +.. function:: void fmpcb_poly_one(fmpcb_poly_t poly) + + Sets *poly* to the constant polynomial 1. + +Input and output +------------------------------------------------------------------------------- + +.. function:: void fmpcb_poly_printd(const fmpcb_poly_t poly, long digits) + + Prints the polynomial as an array of coefficients, printing each + coefficient using *fmprb_printd*. + +Conversions +------------------------------------------------------------------------------- + +.. function:: void fmpcb_poly_set_fmprb_poly(fmpcb_poly_t poly, const fmprb_poly_t re) + +.. function:: void fmpcb_poly_set2_fmprb_poly(fmpcb_poly_t poly, const fmprb_poly_t re, const fmprb_poly_t im) + +.. function:: void fmpcb_poly_set_fmpq_poly(fmpcb_poly_t poly, const fmpq_poly_t re, long prec) + +.. function:: void fmpcb_poly_set2_fmpq_poly(fmpcb_poly_t poly, const fmpq_poly_t re, const fmpq_poly_t im, long prec) + + Sets *poly* to the given real polynomial *re* plus the polynomial *im* + multiplied by the imaginary unit. + + +Evaluation +------------------------------------------------------------------------------- + +.. function:: void _fmpcb_poly_evaluate(fmpcb_t res, const fmpcb_struct * f, long len, const fmpcb_t a, long prec) + +.. function:: void fmpcb_poly_evaluate(fmpcb_t res, const fmpcb_poly_t f, const fmpcb_t a, long prec) + + Evaluates the polynomial using Horner's rule. + + +Derivatives +------------------------------------------------------------------------------- + +.. function:: void _fmpcb_poly_derivative(fmpcb_struct * res, const fmpcb_struct * poly, long len, long prec) + +.. function:: void fmpcb_poly_derivative(fmpcb_poly_t res, const fmpcb_poly_t poly, long prec) + + Sets *res* to the derivative of *poly*. + + +Root-finding +------------------------------------------------------------------------------- + +.. function:: void _fmpcb_poly_root_inclusion(fmpcb_t r, const fmpcb_t m, const fmpcb_struct * poly, const fmpcb_struct * polyder, long len, long prec); + + Given any complex number `m`, and a nonconstant polynomial `f` and its + derivative `f'`, sets *r* to an interval centered on `m` that is + guaranteed to contain at least one root of `f`. + Such an interval is obtained by taking a ball of radius `|f(m)/f'(m)| n` + where `n` is the degree of `f`. + +.. function:: long _fmpcb_poly_validate_roots(fmpcb_struct * roots, const fmpcb_struct * poly, long len, long prec) + + Given a list of approximate roots of the input polynomial, this + function sets a rigorous bounding interval for each root, and determines + which roots are isolated from all the other roots. + It then rearranges the list of roots so that the isolated roots + are at the front of the list, and returns the count of isolated roots. + + If the return value equals the degree of the polynomial, then all + roots have been found. If the return value is smaller, all the + remaining output intervals are guaranteed to contain roots, but + it is possible that not all of the polynomial's roots are contained + among them. + +.. function:: void _fmpcb_poly_refine_roots_durand_kerner(fmpcb_struct * roots, const fmpcb_struct * poly, long len, long prec) + + Refines the given roots simultaneously using a single iteration + of the Durand-Kerner method. The radius of each root is set to an + approximation of the correction, giving a rough estimate of its error (not + a rigorous bound). + +.. function:: long _fmpcb_poly_find_roots(fmpcb_struct * roots, const fmpcb_struct * poly, const fmpcb_struct * initial, long len, long maxiter, long prec) + +.. function:: long fmpcb_poly_find_roots(fmpcb_struct * roots, const fmpcb_poly_t poly, const fmpcb_struct * initial, long maxiter, long prec) + + Attempts to compute all the roots of the given nonzero polynomial *poly* + using a working precision of *prec* bits. If *n* denotes the degree of *poly*, + the function writes *n* approximate roots with rigorous error bounds to + the preallocated array *roots*, and returns the number of + roots that are isolated. + + If the return value equals the degree of the polynomial, then all + roots have been found. If the return value is smaller, all the output + intervals are guaranteed to contain roots, but it is possible that + not all of the polynomial's roots are contained among them. + + The roots are computed numerically by performing several steps with + the Durand-Kerner method and terminating if the estimated accuracy of + the roots approaches the working precision or if the number + of steps exceeds *maxiter*, which can be set to zero in order to use + a default value. Finally, the approximate roots are validated rigorously. + + Initial values for the iteration can be provided as the array *initial*. + If *initial* is set to *NULL*, default values `(0.4+0.9i)^k` are used. + + The polynomial is assumed to be squarefree. If there are repeated + roots, the iteration is likely to find them (with low numerical accuracy), + but the error bounds will not converge as the precision increases. + diff --git a/doc/source/index.rst b/doc/source/index.rst index e9c45a9a..50c53c32 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -30,6 +30,8 @@ Module documentation fmprb.rst fmprb_poly.rst fmprb_mat.rst + fmpcb.rst + fmpcb_poly.rst fmpz_holonomic.rst Credits and references diff --git a/fmpcb_poly.h b/fmpcb_poly.h index 566dbb60..ee66a1a3 100644 --- a/fmpcb_poly.h +++ b/fmpcb_poly.h @@ -90,8 +90,12 @@ void _fmpcb_poly_derivative(fmpcb_struct * res, const fmpcb_struct * poly, long void fmpcb_poly_derivative(fmpcb_poly_t res, const fmpcb_poly_t poly, long prec); +void fmpcb_poly_set_fmprb_poly(fmpcb_poly_t poly, const fmprb_poly_t re); + void fmpcb_poly_set2_fmprb_poly(fmpcb_poly_t poly, const fmprb_poly_t re, const fmprb_poly_t im); +void fmpcb_poly_set_fmpq_poly(fmpcb_poly_t poly, const fmpq_poly_t re, long prec); + void fmpcb_poly_set2_fmpq_poly(fmpcb_poly_t poly, const fmpq_poly_t re, const fmpq_poly_t im, long prec); void _fmpcb_poly_root_inclusion(fmpcb_t r, const fmpcb_t m, diff --git a/fmpcb_poly/refine_roots_durand_kerner.c b/fmpcb_poly/refine_roots_durand_kerner.c index acd39010..767f143a 100644 --- a/fmpcb_poly/refine_roots_durand_kerner.c +++ b/fmpcb_poly/refine_roots_durand_kerner.c @@ -135,14 +135,6 @@ _fmpcb_poly_evaluate_mid(fmpcb_t res, const fmpcb_struct * f, long len, fmpcb_clear(t); } - -/* -Refines the given roots simultaneously using a single iteration -of the Durand-Kerner method. The radius of each root is set to an -approximation of the correction, giving a rough estimate of its error (not -a rigorous bound). -*/ - void _fmpcb_poly_refine_roots_durand_kerner(fmpcb_struct * roots, const fmpcb_struct * poly, long len, long prec) diff --git a/fmpcb_poly/root_inclusion.c b/fmpcb_poly/root_inclusion.c index 1d5c4b19..bec8513c 100644 --- a/fmpcb_poly/root_inclusion.c +++ b/fmpcb_poly/root_inclusion.c @@ -25,11 +25,6 @@ #include "fmpcb_poly.h" -/* -Given any complex number m, and a polynomial f and its derivative f', -sets r to an interval centered on m that is guaranteed to contain -at least one root of f. Assumes len > 1. -*/ void _fmpcb_poly_root_inclusion(fmpcb_t r, const fmpcb_t m, const fmpcb_struct * poly, diff --git a/fmpcb_poly/set2_fmpq_poly.c b/fmpcb_poly/set2_fmpq_poly.c index 3d627214..1f0553bd 100644 --- a/fmpcb_poly/set2_fmpq_poly.c +++ b/fmpcb_poly/set2_fmpq_poly.c @@ -25,6 +25,16 @@ #include "fmpcb_poly.h" +void +fmpcb_poly_set_fmpq_poly(fmpcb_poly_t poly, const fmpq_poly_t re, long prec) +{ + fmprb_poly_t t; + fmprb_poly_init(t); + fmprb_poly_set_fmpq_poly(t, re, prec); + fmpcb_poly_set_fmprb_poly(poly, t); + fmprb_poly_clear(t); +} + void fmpcb_poly_set2_fmpq_poly(fmpcb_poly_t poly, const fmpq_poly_t re, const fmpq_poly_t im, long prec) { diff --git a/fmpcb_poly/set2_fmprb_poly.c b/fmpcb_poly/set2_fmprb_poly.c index 8be119d5..d8fc7c21 100644 --- a/fmpcb_poly/set2_fmprb_poly.c +++ b/fmpcb_poly/set2_fmprb_poly.c @@ -25,6 +25,24 @@ #include "fmpcb_poly.h" +void +fmpcb_poly_set_fmprb_poly(fmpcb_poly_t poly, const fmprb_poly_t re) +{ + long i, len; + + len = fmprb_poly_length(re); + + fmpcb_poly_fit_length(poly, len); + + for (i = 0; i < len; i++) + { + fmprb_set(fmpcb_realref(poly->coeffs + i), re->coeffs + i); + fmprb_zero(fmpcb_imagref(poly->coeffs + i)); + } + + _fmpcb_poly_set_length(poly, len); +} + void fmpcb_poly_set2_fmprb_poly(fmpcb_poly_t poly, const fmprb_poly_t re, const fmprb_poly_t im) { diff --git a/fmpcb_poly/validate_roots.c b/fmpcb_poly/validate_roots.c index 4a5d41b2..3eeb6969 100644 --- a/fmpcb_poly/validate_roots.c +++ b/fmpcb_poly/validate_roots.c @@ -25,20 +25,6 @@ #include "fmpcb_poly.h" -/* -Given a list of approximate roots of the input polynomial, this -function sets a rigorous bounding interval for each root, and determines -which roots are isolated from all the other roots. -It then rearranges the list of roots so that the isolated roots -are at the front of the list, and returns the count of isolated roots. - -If the return value equals the degree of the polynomial, then all -roots have been separated. If the return value is smaller, all the -remaining intervals are guaranteed to contain roots, but -it is possible that not all of the polynomial's roots are contained -among them. -*/ - long _fmpcb_poly_validate_roots(fmpcb_struct * roots, const fmpcb_struct * poly, long len, long prec)