From octave-sources-request at bevo dot che dot wisc dot edu Mon Oct 4 16:53:42 1999 Subject: function cerf(z) - calculate erf(z) for complex z From: John Smith To: octave-sources at bevo dot che dot wisc dot edu Date: Mon, 4 Oct 1999 22:52:37 +0100 (BST) This is a bit of C++ to compile into a dynamically loadable octave function that calculates erf(z) for complex z. The erf(z) present in Octave-2.0.14 is restricted to real z only. The resulting function I have called cerf (Complex erf). It would be best to put the rest of this message into a file cerf.cc, and compile it with mkoctfile cerf.cc At least under some variant of linux, and gcc 2.8.1 it compiles. This function was previously offered to octave-sources in November 1998. Someone tried to use it, and told me of a bug which added a spurious 2 to the answer for some negative z. So this is a slightly revised copy. Good luck.... John ------------------------------------------------- /* Copyright (C) 1998, 1999 John Smith This file is part of Octave. Or it might be one day.... Octave is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. Octave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Octave; see the file COPYING. If not, write to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ // Put together by John Smith john at arrows dot demon dot co dot uk, // using ideas by others. // // Calculate erf(z) for complex z. // Three methods are implemented; which one is used depends on z. // // The code includes some hard coded constants that are intended to // give about 14 decimal places of accuracy. This is appropriate for // 64-bit floating point numbers. // // Oct 1999: Fixed a typo that in // const Complex cerf_continued_fraction( const Complex z ) // that caused erroneous answers for erf(z) where real(z) negative // #include #include #include "f77-fcn.h" // // Abramowitz and Stegun: (eqn: 7.1.14) gives this continued // fraction for erfc(z) // // erfc(z) = sqrt(pi).exp(-z^2). 1 1/2 1 3/2 2 5/2 // --- --- --- --- --- --- ... // z + z + z + z + z + z + // // This is evaluated using Lentz's method, as described in the narative // of Numerical Recipes in C. // // The continued fraction is true providing real(z)>0. In practice we // like real(z) to be significantly greater than 0, say greater than 0.5. // const Complex cerfc_continued_fraction( const Complex z ) { double tiny = 1e-20 ; // a small number, large enough to calculate 1/tiny double eps = 1e-15 ; // large enough so that 1.0+eps > 1.0, when using // the floating point arithmetic // // first calculate z+ 1/2 1 // --- --- ... // z + z + Complex f(z) ; Complex C(f) ; Complex D(0.0) ; Complex delta ; double a ; a = 0.0 ; do { a = a + 0.5 ; D = z + a*D ; C = z + a/C ; if (D.real() == 0.0 && D.imag() == 0.0) D = tiny ; D = 1.0 / D ; delta = (C * D) ; f = f * delta ; } while (abs(1.0-delta) > eps ) ; // // Do the first term of the continued fraction // f = 1.0 / f ; // // and do the final scaling // f = f * exp(-z*z)/ sqrt(M_PI) ; return f ; } const Complex cerf_continued_fraction( const Complex z ) { // warning("cerf_continued_fraction:"); if (z.real() > 0) return 1.0 - cerfc_continued_fraction( z ) ; else return -1.0 + cerfc_continued_fraction( -z ) ; } // // Abramawitz and Stegun, Eqn. 7.1.5 gives a series for erf(z) // good for all z, but converges faster for smallish abs(z), say abs(z)<2. // const Complex cerf_series( const Complex z ) { double tiny = 1e-20 ; // a small number compared with 1. // warning("cerf_series:"); Complex sum(0.0) ; Complex term(z) ; Complex z2(z*z) ; for (int n=0; n<3 || abs(term) > abs(sum)*tiny; n++) { sum = sum + term / (2*n+1) ; term = -term * z2 / (n+1) ; } return sum * 2.0 / sqrt(M_PI) ; } // // Numerical Recipes quotes a formula due to Rybicki for evaluating // Dawson's Integral: // // exp(-x^2) integral exp(t^2).dt = 1/sqrt(pi) lim sum exp(-(z-n.h)^2) / n // 0 to x h->0 n odd // // This can be adapted to erf(z). // const Complex cerf_rybicki( const Complex z ) { // warning("cerf_rybicki:"); double h = 0.2 ; // numerical experiment suggests this is small enough // // choose an even n0, and then shift z->z-n0.h and n->n-h. // n0 is chosen so that real((z-n0.h)^2) is as small as possible. // int n0 = 2*(int) (floor( z.imag()/(2*h) + 0.5 )) ; Complex z0( 0.0, n0*h ) ; Complex zp(z-z0) ; Complex sum(0.0,0.0) ; // // limits of sum chosen so that the end sums of the sum are // fairly small. In this case exp(-(35.h)^2)=5e-22 // // for (int np=-35; np<=35; np+=2) { Complex t( zp.real(), zp.imag()-np*h) ; Complex b( exp(t*t) / (np+n0) ) ; sum += b ; } sum = sum * 2 * exp(-z*z) / M_PI ; return Complex(-sum.imag(), sum.real()) ; } const Complex cerf( const Complex z ) { // // Use the method appropriate to size of z - // there probably ought to be an extra option for NaN z, or infinite z // // if (abs(z) < 2.0) return cerf_series( z ) ; else if (abs(z.real()) < 0.5) return cerf_rybicki( z ) ; else return cerf_continued_fraction( z ) ; } DEFUN_DLD (cerf, args, , "erf(z)\n\ \n\ Compute the error function erf(z) for any complex values of z.\n\ \n\ ") { octave_value arg = args(0) ; ComplexMatrix cm = arg.complex_matrix_value() ; int nr = cm.rows() ; int nc = cm.columns() ; ComplexMatrix ans( nr, nc ) ; for (int i=0; i