From: Frank Reininghaus <frank78ac@googlemail.com>
To: gsl-discuss@sourceware.org
Subject: Complex polynomial evaluation
Date: Tue, 11 Dec 2007 23:04:00 -0000 [thread overview]
Message-ID: <475F1765.5020300@googlemail.com> (raw)
[-- 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;
> }
next reply other threads:[~2007-12-11 23:04 UTC|newest]
Thread overview: 5+ messages / expand[flat|nested] mbox.gz Atom feed top
2007-12-11 23:04 Frank Reininghaus [this message]
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
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=475F1765.5020300@googlemail.com \
--to=frank78ac@googlemail.com \
--cc=gsl-discuss@sourceware.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
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).