From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 8695 invoked by alias); 11 Dec 2007 23:04:25 -0000 Received: (qmail 8684 invoked by uid 22791); 11 Dec 2007 23:04:22 -0000 X-Spam-Check-By: sourceware.org Received: from nf-out-0910.google.com (HELO nf-out-0910.google.com) (64.233.182.185) by sourceware.org (qpsmtpd/0.31) with ESMTP; Tue, 11 Dec 2007 23:04:12 +0000 Received: by nf-out-0910.google.com with SMTP id b11so1553698nfh for ; Tue, 11 Dec 2007 15:04:09 -0800 (PST) Received: by 10.82.121.15 with SMTP id t15mr5498856buc.1197414248162; Tue, 11 Dec 2007 15:04:08 -0800 (PST) Received: from ?192.168.178.21? ( [87.189.21.92]) by mx.google.com with ESMTPS id j9sm12944849mue.2007.12.11.15.04.06 (version=TLSv1/SSLv3 cipher=RC4-MD5); Tue, 11 Dec 2007 15:04:06 -0800 (PST) Message-ID: <475F1765.5020300@googlemail.com> Date: Tue, 11 Dec 2007 23:04:00 -0000 From: Frank Reininghaus User-Agent: Thunderbird 2.0.0.6 (X11/20071022) MIME-Version: 1.0 To: gsl-discuss@sourceware.org Subject: Complex polynomial evaluation Content-Type: multipart/mixed; boundary="------------060301020202040703020808" Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2007-q4/txt/msg00039.txt.bz2 This is a multi-part message in MIME format. --------------060301020202040703020808 Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Content-length: 1379 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 --------------060301020202040703020808 Content-Type: text/x-csrc; name="patch_eval.c" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="patch_eval.c" Content-length: 1198 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; > } --------------060301020202040703020808 Content-Type: text/x-chdr; name="patch_gsl_poly.h" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="patch_gsl_poly.h" Content-length: 1562 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; > } --------------060301020202040703020808--