* Complex polynomial evaluation
@ 2007-12-11 23:04 Frank Reininghaus
2007-12-12 16:48 ` Brian Gough
0 siblings, 1 reply; 5+ messages in thread
From: Frank Reininghaus @ 2007-12-11 23:04 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 1379 bytes --]
Hi,
I've written a bit of code that evaluates polynomials with real or
complex coefficients for complex variables. There was some discussion
about this a couple of years ago, see
http://sourceware.org/ml/gsl-discuss/2004-q2/msg00040.html
and follow-ups.
I've followed the example of gsl_poly_eval () in gsl_poly.h and
poly/eval.c and wrote two functions:
gsl_poly_complex_eval () evaluates a polynomial with real coefficients
for a complex variable,
gsl_complex_poly_complex_eval () does the same with complex coefficients.
I'm not quite sure if the choice of names is appropriate, but I think
that it makes sense to have both functions because one often needs to
evaluate real polynomials for complex values (like Taylor series
expansions for some functions).
If you are interested in including this in GSL, you can use the patches
for gsl_poly.h and eval.c which I have attached (an alternative might be
to create a new header file gsl_complex_poly.h). If this is the case, I
could also write a bit of documentation and code for some test cases to
be included in 'make check'. So far, I've checked the results for many
random polynomials with Maxima. Moreover, I've checked with random data
that both functions yield the same results for real polynomials and that
their results coincide with gsl_poly_eval () if the variable is also real.
Regards,
Frank
[-- Attachment #2: patch_eval.c --]
[-- Type: text/x-csrc, Size: 1198 bytes --]
35a36,68
>
> gsl_complex
> gsl_poly_complex_eval(const double c[], const int len, const gsl_complex x)
> {
> int i;
> gsl_complex ans;
> GSL_SET_COMPLEX (&ans, c[len-1], 0.0);
> for(i=len-1; i>0; i--) {
> /* The following three lines are equivalent to
> ans = gsl_complex_add_real (gsl_complex_mul (x, ans), c[i-1]);
> but faster */
> double tmp = c[i-1] + GSL_REAL (x) * GSL_REAL (ans) - GSL_IMAG (x) * GSL_IMAG (ans);
> GSL_SET_IMAG (&ans, GSL_IMAG (x) * GSL_REAL (ans) + GSL_REAL (x) * GSL_IMAG (ans));
> GSL_SET_REAL (&ans, tmp);
> }
> return ans;
> }
>
> gsl_complex
> gsl_complex_poly_complex_eval(const gsl_complex c[], const int len, const gsl_complex x)
> {
> int i;
> gsl_complex ans = c[len-1];
> for(i=len-1; i>0; i--) {
> /* The following three lines are equivalent to
> ans = gsl_complex_add (c[i-1], gsl_complex_mul (x, ans));
> but faster */
> double tmp = GSL_REAL (c[i-1]) + GSL_REAL (x) * GSL_REAL (ans) - GSL_IMAG (x) * GSL_IMAG (ans);
> GSL_SET_IMAG (&ans, GSL_IMAG (c[i-1]) + GSL_IMAG (x) * GSL_REAL (ans) + GSL_REAL (x) * GSL_IMAG (ans));
> GSL_SET_REAL (&ans, tmp);
> }
> return ans;
> }
[-- Attachment #3: patch_gsl_poly.h --]
[-- Type: text/x-chdr, Size: 1562 bytes --]
44a45,46
>
> /* real polynomial, real x */
46a49,54
> /* real polynomial, complex x */
> gsl_complex gsl_poly_complex_eval (const double c [], const int len, const gsl_complex x);
>
> /* complex polynomial, complex x */
> gsl_complex gsl_complex_poly_complex_eval (const gsl_complex c [], const int len, const gsl_complex x);
>
57a66,100
>
> extern inline
> gsl_complex
> gsl_poly_complex_eval(const double c[], const int len, const gsl_complex x)
> {
> int i;
> gsl_complex ans;
> GSL_SET_COMPLEX (&ans, c[len-1], 0.0);
> for(i=len-1; i>0; i--) {
> /* The following three lines are equivalent to
> ans = gsl_complex_add_real (gsl_complex_mul (x, ans), c[i-1]);
> but faster */
> double tmp = c[i-1] + GSL_REAL (x) * GSL_REAL (ans) - GSL_IMAG (x) * GSL_IMAG (ans);
> GSL_SET_IMAG (&ans, GSL_IMAG (x) * GSL_REAL (ans) + GSL_REAL (x) * GSL_IMAG (ans));
> GSL_SET_REAL (&ans, tmp);
> }
> return ans;
> }
>
> extern inline
> gsl_complex
> gsl_complex_poly_complex_eval(const gsl_complex c[], const int len, const gsl_complex x)
> {
> int i;
> gsl_complex ans = c[len-1];
> for(i=len-1; i>0; i--) {
> /* The following three lines are equivalent to
> ans = gsl_complex_add (c[i-1], gsl_complex_mul (x, ans));
> but faster */
> double tmp = GSL_REAL (c[i-1]) + GSL_REAL (x) * GSL_REAL (ans) - GSL_IMAG (x) * GSL_IMAG (ans);
> GSL_SET_IMAG (&ans, GSL_IMAG (c[i-1]) + GSL_IMAG (x) * GSL_REAL (ans) + GSL_REAL (x) * GSL_IMAG (ans));
> GSL_SET_REAL (&ans, tmp);
> }
> return ans;
> }
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: Complex polynomial evaluation
2007-12-11 23:04 Complex polynomial evaluation Frank Reininghaus
@ 2007-12-12 16:48 ` Brian Gough
2007-12-13 20:05 ` Frank Reininghaus
0 siblings, 1 reply; 5+ messages in thread
From: Brian Gough @ 2007-12-12 16:48 UTC (permalink / raw)
To: Frank Reininghaus; +Cc: gsl-discuss
At Wed, 12 Dec 2007 00:04:05 +0100,
Frank Reininghaus wrote:
> If you are interested in including this in GSL, you can use the patches
> for gsl_poly.h and eval.c which I have attached (an alternative might be
> to create a new header file gsl_complex_poly.h). If this is the case, I
> could also write a bit of documentation and code for some test cases to
> be included in 'make check'.
Yes - that sounds good. Thanks.
--
Brian Gough
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: Complex polynomial evaluation
2007-12-12 16:48 ` Brian Gough
@ 2007-12-13 20:05 ` Frank Reininghaus
2007-12-14 13:40 ` Brian Gough
0 siblings, 1 reply; 5+ messages in thread
From: Frank Reininghaus @ 2007-12-13 20:05 UTC (permalink / raw)
To: gsl-discuss
[-- Attachment #1: Type: text/plain, Size: 694 bytes --]
Hi,
I've written some documentation for my complex polynomial evaluation
functions and some test cases to be checked by 'make check'. I've
created patches containing all my modifications to the current CVS
version with
diff gsl/poly/ [my directory with new version of 'poly'] > patch-poly
diff gsl/doc/ [my directory with new version of 'doc'] > patch-doc
Let me know if there is a way I could create patch files that are more
convenient to you.
patch-doc also contains a correction for a typo in the macro @inlinefns
in doc/gsl-ref.texi. I think it should be "Inline versions of these
*functions* are used...", not "Inline versions of these *function* are
used..."
Regards,
Frank
[-- Attachment #2: patch-doc --]
[-- Type: text/plain, Size: 1533 bytes --]
Common subdirectories: gsl/doc/CVS and /home/frank/tmp/doc-with-complex-polynomials/CVS
Common subdirectories: gsl/doc/examples and /home/frank/tmp/doc-with-complex-polynomials/examples
diff gsl/doc/gsl-ref.texi /home/frank/tmp/doc-with-complex-polynomials/gsl-ref.texi
127c127
< Inline versions of these function are used when @code{HAVE_INLINE} is defined.
---
> Inline versions of these functions are used when @code{HAVE_INLINE} is defined.
diff gsl/doc/poly.texi /home/frank/tmp/doc-with-complex-polynomials/poly.texi
25,26c25
< @deftypefun double gsl_poly_eval (const double @var{c}[], const int @var{len}, const double @var{x})
< This function evaluates the polynomial
---
> The functions described here evaluate the polynomial
29c28,39
< Horner's method for stability. @inlinefn{}
---
> Horner's method for stability. @inlinefns{}
>
> @deftypefun double gsl_poly_eval (const double @var{c}[], const int @var{len}, const double @var{x})
> This function evaluates a polynomial with real coefficients for the real variable @var{x}.
> @end deftypefun
>
> @deftypefun gsl_complex gsl_poly_complex_eval (const double @var{c}[], const int @var{len}, const gsl_complex @var{x})
> This function evaluates a polynomial with real coefficients for the complex variable @var{x}.
> @end deftypefun
>
> @deftypefun gsl_complex gsl_complex_poly_complex_eval (const gsl_complex @var{c}[], const int @var{len}, const gsl_complex @var{x})
> This function evaluates a polynomial with complex coefficients for the complex variable @var{x}.
[-- Attachment #3: patch-poly --]
[-- Type: text/plain, Size: 5231 bytes --]
Common subdirectories: gsl/poly/CVS and /home/frank/tmp/poly-with-complex-evaluation/CVS
diff gsl/poly/eval.c /home/frank/tmp/poly-with-complex-evaluation/eval.c
35a36,68
>
> gsl_complex
> gsl_poly_complex_eval(const double c[], const int len, const gsl_complex x)
> {
> int i;
> gsl_complex ans;
> GSL_SET_COMPLEX (&ans, c[len-1], 0.0);
> for(i=len-1; i>0; i--) {
> /* The following three lines are equivalent to
> ans = gsl_complex_add_real (gsl_complex_mul (x, ans), c[i-1]);
> but faster */
> double tmp = c[i-1] + GSL_REAL (x) * GSL_REAL (ans) - GSL_IMAG (x) * GSL_IMAG (ans);
> GSL_SET_IMAG (&ans, GSL_IMAG (x) * GSL_REAL (ans) + GSL_REAL (x) * GSL_IMAG (ans));
> GSL_SET_REAL (&ans, tmp);
> }
> return ans;
> }
>
> gsl_complex
> gsl_complex_poly_complex_eval(const gsl_complex c[], const int len, const gsl_complex x)
> {
> int i;
> gsl_complex ans = c[len-1];
> for(i=len-1; i>0; i--) {
> /* The following three lines are equivalent to
> ans = gsl_complex_add (c[i-1], gsl_complex_mul (x, ans));
> but faster */
> double tmp = GSL_REAL (c[i-1]) + GSL_REAL (x) * GSL_REAL (ans) - GSL_IMAG (x) * GSL_IMAG (ans);
> GSL_SET_IMAG (&ans, GSL_IMAG (c[i-1]) + GSL_IMAG (x) * GSL_REAL (ans) + GSL_REAL (x) * GSL_IMAG (ans));
> GSL_SET_REAL (&ans, tmp);
> }
> return ans;
> }
diff gsl/poly/gsl_poly.h /home/frank/tmp/poly-with-complex-evaluation/gsl_poly.h
44a45,46
>
> /* real polynomial, real x */
46a49,54
> /* real polynomial, complex x */
> gsl_complex gsl_poly_complex_eval (const double c [], const int len, const gsl_complex x);
>
> /* complex polynomial, complex x */
> gsl_complex gsl_complex_poly_complex_eval (const gsl_complex c [], const int len, const gsl_complex x);
>
57a66,100
>
> extern inline
> gsl_complex
> gsl_poly_complex_eval(const double c[], const int len, const gsl_complex x)
> {
> int i;
> gsl_complex ans;
> GSL_SET_COMPLEX (&ans, c[len-1], 0.0);
> for(i=len-1; i>0; i--) {
> /* The following three lines are equivalent to
> ans = gsl_complex_add_real (gsl_complex_mul (x, ans), c[i-1]);
> but faster */
> double tmp = c[i-1] + GSL_REAL (x) * GSL_REAL (ans) - GSL_IMAG (x) * GSL_IMAG (ans);
> GSL_SET_IMAG (&ans, GSL_IMAG (x) * GSL_REAL (ans) + GSL_REAL (x) * GSL_IMAG (ans));
> GSL_SET_REAL (&ans, tmp);
> }
> return ans;
> }
>
> extern inline
> gsl_complex
> gsl_complex_poly_complex_eval(const gsl_complex c[], const int len, const gsl_complex x)
> {
> int i;
> gsl_complex ans = c[len-1];
> for(i=len-1; i>0; i--) {
> /* The following three lines are equivalent to
> ans = gsl_complex_add (c[i-1], gsl_complex_mul (x, ans));
> but faster */
> double tmp = GSL_REAL (c[i-1]) + GSL_REAL (x) * GSL_REAL (ans) - GSL_IMAG (x) * GSL_IMAG (ans);
> GSL_SET_IMAG (&ans, GSL_IMAG (c[i-1]) + GSL_IMAG (x) * GSL_REAL (ans) + GSL_REAL (x) * GSL_IMAG (ans));
> GSL_SET_REAL (&ans, tmp);
> }
> return ans;
> }
diff gsl/poly/test.c /home/frank/tmp/poly-with-complex-evaluation/test.c
33a34,35
> /* Polynomial evaluation */
>
52a55,108
> {
> gsl_complex x, y;
> double c[1] = {0.3};
> GSL_SET_REAL (&x, 0.75);
> GSL_SET_IMAG (&x, 1.2);
> y = gsl_poly_complex_eval (c, 1, x);
>
> gsl_test_rel (GSL_REAL (y), 0.3, eps, "y.real, gsl_poly_complex_eval ({0.3}, 0.75 + 1.2i)");
> gsl_test_rel (GSL_IMAG (y), 0.0, eps, "y.imag, gsl_poly_complex_eval ({0.3}, 0.75 + 1.2i)");
> }
>
> {
> gsl_complex x, y;
> double c[4] = {2.1, -1.34, 0.76, 0.45};
> GSL_SET_REAL (&x, 0.49);
> GSL_SET_IMAG (&x, 0.95);
> y = gsl_poly_complex_eval (c, 4, x);
>
> gsl_test_rel (GSL_REAL (y), 0.3959143, eps, "y.real, gsl_poly_complex_eval ({2.1, -1.34, 0.76, 0.45}, 0.49 + 0.95i)");
> gsl_test_rel (GSL_IMAG (y), -0.6433305, eps, "y.imag, gsl_poly_complex_eval ({2.1, -1.34, 0.76, 0.45}, 0.49 + 0.95i)");
> }
>
> {
> gsl_complex x, y;
> gsl_complex c[1];
> GSL_SET_REAL (&c[0], 0.674);
> GSL_SET_IMAG (&c[0], -1.423);
> GSL_SET_REAL (&x, -1.44);
> GSL_SET_IMAG (&x, 9.55);
> y = gsl_complex_poly_complex_eval (c, 1, x);
>
> gsl_test_rel (GSL_REAL (y), 0.674, eps, "y.real, gsl_complex_poly_complex_eval ({0.674 - 1.423i}, -1.44 + 9.55i)");
> gsl_test_rel (GSL_IMAG (y), -1.423, eps, "y.imag, gsl_complex_poly_complex_eval ({0.674 - 1.423i}, -1.44 + 9.55i)");
> }
>
> {
> gsl_complex x, y;
> gsl_complex c[4];
> GSL_SET_REAL (&c[0], -2.31);
> GSL_SET_IMAG (&c[0], 0.44);
> GSL_SET_REAL (&c[1], 4.21);
> GSL_SET_IMAG (&c[1], -3.19);
> GSL_SET_REAL (&c[2], 0.93);
> GSL_SET_IMAG (&c[2], 1.04);
> GSL_SET_REAL (&c[3], -0.42);
> GSL_SET_IMAG (&c[3], 0.68);
> GSL_SET_REAL (&x, 0.49);
> GSL_SET_IMAG (&x, 0.95);
> y = gsl_complex_poly_complex_eval (c, 4, x);
>
> gsl_test_rel (GSL_REAL (y), 1.82462012, eps, "y.real, gsl_complex_poly_complex_eval ({-2.31 + 0.44i, 4.21 - 3.19i, 0.93 + 1.04i, -0.42 + 0.68i}, 0.49 + 0.95i)");
> gsl_test_rel (GSL_IMAG (y), 2.30389412, eps, "y.imag, gsl_complex_poly_complex_eval ({-2.31 + 0.44i, 4.21 - 3.19i, 0.93 + 1.04i, -0.42 + 0.68i}, 0.49 + 0.95i)");
> }
>
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: Complex polynomial evaluation
2007-12-13 20:05 ` Frank Reininghaus
@ 2007-12-14 13:40 ` Brian Gough
2007-12-17 22:10 ` Brian Gough
0 siblings, 1 reply; 5+ messages in thread
From: Brian Gough @ 2007-12-14 13:40 UTC (permalink / raw)
To: Frank Reininghaus; +Cc: gsl-discuss
At Thu, 13 Dec 2007 21:05:13 +0100,
Frank Reininghaus wrote:
> Let me know if there is a way I could create patch files that are more
> convenient to you.
Can you send them with diff -u, thanks.
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: Complex polynomial evaluation
2007-12-14 13:40 ` Brian Gough
@ 2007-12-17 22:10 ` Brian Gough
0 siblings, 0 replies; 5+ messages in thread
From: Brian Gough @ 2007-12-17 22:10 UTC (permalink / raw)
To: Frank Reininghaus; +Cc: gsl-discuss
At Fri, 14 Dec 2007 13:39:57 +0000,
Brian Gough wrote:
>
> At Thu, 13 Dec 2007 21:05:13 +0100,
> Frank Reininghaus wrote:
> > Let me know if there is a way I could create patch files that are more
> > convenient to you.
>
> Can you send them with diff -u, thanks.
>
Thanks, I've added these into the CVS repository now.
--
Brian Gough
^ permalink raw reply [flat|nested] 5+ messages in thread
end of thread, other threads:[~2007-12-17 22:10 UTC | newest]
Thread overview: 5+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2007-12-11 23:04 Complex polynomial evaluation Frank Reininghaus
2007-12-12 16:48 ` Brian Gough
2007-12-13 20:05 ` Frank Reininghaus
2007-12-14 13:40 ` Brian Gough
2007-12-17 22:10 ` Brian Gough
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).