cpow and clog

Montgomery-Smith, Stephen stephen at missouri.edu
Wed Nov 8 01:50:44 UTC 2017


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.


More information about the freebsd-numerics mailing list