Date: Wed, 8 Nov 2017 01:49:32 +0000 From: "Montgomery-Smith, Stephen" <stephen@missouri.edu> To: "freebsd-numerics@freebsd.org" <freebsd-numerics@freebsd.org> Subject: Re: cpow and clog Message-ID: <4330a817-ae31-6336-7aa6-3e9981b02ab7@missouri.edu> In-Reply-To: <20171106204121.GB37361@troutmask.apl.washington.edu> References: <20171106194937.GA87725@freebird> <20171106204121.GB37361@troutmask.apl.washington.edu>
index | next in thread | previous in thread | raw e-mail
On 11/06/2017 02:41 PM, Steve Kargl wrote: > On Mon, Nov 06, 2017 at 08:49:43PM +0100, Michael Danilov wrote: >> Hello, >> >> I would like to have some feedback on my attempt to import OpenBSD >> code for cpow and clog: >> >> https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=221341 >> https://bugs.freebsd.org/bugzilla/attachment.cgi?id=187693 The problem with this implementation of clog is that they will produce very inaccurate results when |z| is close to 1. A proper implementation is given in the paper: T. E. Hull, Thomas F. Fairgrieve, Ping Tak Peter Tang, “Implementing Complex Elementary Functions Using Exception Handling”, ACM Transactions on Mathematical Software, Volume 20 (1994), pp 215–244, Corrigenda, p553. They use a precision which is twice that of the required precision. For example, if you want log(sqrt(x^2+y^2)) accurate to about 17 digits, you will need to compute it to 34 digits. If all you care about is getting a a good relative error for the total answer (rather than insisting on a good relative error for the real and imaginary parts separately), then for |x+I*y| close to one, compute the real part using 0.5*logp1((x-1)*(x+1)+y*y) >> >> What happened to the alternative implementation mentioned in the thread below? > > bde has an implementation of clog[fl]. He may someday > commit it. I don't know if anyone ever worked on cpow[fl]. > I stopped working on powl and tgammal when I returned my > commit bit due to differences with "higher-ranking" committers. > bde's implementation was fantastic. He managed to use two double variables to represent a single quadruple precision real (using a technique which I believe is standard, but it was new to me). Why he never committed it, I don't know. People on this list take way too long to commit excellent software. For example, I had code for cacos, catan, etc, ready for a long time before someone else stepped in and committed it. And I tested it to death. I even found a bug in T. E. Hull, Thomas F. Fairgrieve, Ping Tak Peter Tang, “Implementing the complex arcsin and arccosine functions using exception handling”, ACM Transactions on Mathematical Software, Volume 23 (1997) pp 299–335.help
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?4330a817-ae31-6336-7aa6-3e9981b02ab7>