cpow and clog

Tue Nov 7 19:24:57 UTC 2017

On Tue, Nov 07, 2017 at 05:53:56PM +0000, dimpase at gmail.com wrote:
> On Tue, Nov 07, 2017 at 07:31:49AM -0800, Steve Kargl wrote:
> > On Tue, Nov 07, 2017 at 10:38:59AM +0000, dimpase at gmail.com wrote:
> > > On Mon, Nov 06, 2017 at 12:41:21PM -0800, Steve Kargl wrote:
> > > > 
> > > > How have you tested the NetBSD and/or OpenBSD code?  What is the
> > > > quality?  Have the long double clogl and cpowl been tested on both
> > > > ld80 and ld128 hardware?  See FreeBSD's lib/msun/src/math_private.h
> > > > for a discussion of possible issues of using I from complex.h in this
> > > > code.
> > > 
> > > I would like to point out that various FreeBSD ports already contain
> > > implementations of the functions in question.
> > 
> > What is the quality of those implementations?  Have you tested the
> > code?
> Do you have a testsuite for these functions?

Not for these functions.  That's why I'm asking if someone 
has tested them.   I'm still working on non-complex functions,
I can actually use in my *real* job.

> Perhaps it can be found somewhere, do you know?

You could ask the Netbsd and Openbsd projects if there
are tests available or you could attempt to write a
test.

Google 'n1256.pdf'.  Look at 7.3.7.2 and Annex G.  clog
has a branch cut along the negative real axis.  If one is
close to the branch cut, does clog return the correct
answer or does it (occasionally) jump across the branch cut.
clog needs to satisfy a number conditions for exceptional
values (See Annex G).  Does the NetBSd/OpenBSD code meet
those conditions.

> We are able to make bits of code needed for conformance to pass
> our tests in Sagemath
> on FreeBSD 11 (x86_64 only, although I probably can get FreeBSD running
> on aarch64 somehow) with various clang versions, and with gcc6. 
> Basically I am talking about these functions from the Cephes library:
> http://files.sagemath.org/spkg/upstream/cephes/index.html
> 
> > Have you read the comment in math_private.h concerning the
> > use of I in computations?
> I presume you refer to
> 
> "...x+I*y is
>  * currently unusable in general since gcc introduces many overflow,
>  * underflow, sign and efficiency bugs..."
> 
> We are building FreeBSD with clang, and not gcc, I don't know how
> relevant it is, for clang is a fast-moving target.

Last time I checked, clang had the same issue as gcc.
The issue was reported to clang developers 7 years ago.

https://bugs.llvm.org/show_bug.cgi?id=8532

> By the way, https://wiki.freebsd.org/Numerics still mentiones FreeBSD 10
> as "current"...

David Das use to keep the wiki up-to-date.  I haven't seen das@
in a long time.  Feel free to fix the wiki.

> > > Sorry for being blunt, but IMHO the attitude on this list appears to be
> > > to let the numerics stack on FreeBSD die a slow death.
> > 
> > You're talking to the one person who has spent decades trying to 
> > improve libm through bug fixes and implementing missing functionality.
> > 
> > > Indeed, most people hate to reinvent the wheel. It's
> > > really no fun at all to scramble to get these missing implementations
> > > somehow, there are certainly much better ways to use one's time and
> > > brainpower.  On this list people prefer to point at some private code in
> > > uncertain shape, and hope that somehow by some magic FreeBSD will have
> > > the best humanely possible implementation of the complex transcendental
> > > functions... 
> > 
> > On this list, which is mostly Bruce and me it seems, people care about
> > the quality of the code.  People, which is mostly Bruce and me, spend
> > quite a bit of time benchmarking proposed patches go ensure that users
> > aren't fed wrong results.
> > 
> > > Why don't you first of all try to provide *some* reasonably
> > > working implementation (thus allowing porters not to have to reinvent
> > > this wheel, badly, for $n$-th time over, and then having *fun* making
> > > sure the tools know where to get these functions), and only then try to
> > > improve it?
> > 
> > Who is 'you' here?  I recently provided implementations of sinpi[fl],
> > cospi[fl], and tanpi[fl] for ISO/IEC TS 18661-4.  Patches were submitted
> > to this list (check the archive).  
> 
> Well, this involves browsing entries of
> https://lists.freebsd.org/pipermail/freebsd-numerics/
> Would you mind at least mention year and month?

Google "FreeBSD sinpi" the top two results are

218514 – [LIBM] implementations of sinpi[fl], cospi ... - FreeBSD Bugzilla
https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=218514

[PATCH] lang/python27 -- Fix namespace collision - FreeBSD lists
https://lists.freebsd.org/pipermail/freebsd-ports/2017-July/109675.html

It is also not too difficult to click on the various *recent*
months in the mail archive, but since you ask.

https://lists.freebsd.org/pipermail/freebsd-numerics/2017-April/thread.html
https://lists.freebsd.org/pipermail/freebsd-numerics/2017-May/thread.html

> > Bruce is the only person to comment on the patch.  I took
> > the time to fix the lang/python* ports, which break with this
> > patch, and sent the patches to freebsd-ports list where the
> >  patches are lost in time.
> 
> Nowadays one would rather expect a git branch with the changes, not
> patches, as you might know.

I don't use git.  

cd /usr/ports
svn diff lang/python27 > python27.diff
svn diff lang/python36 > python36.diff

Mail patches to freebsd-ports and freebsd-python.  Watch patches wither.

> > Check out the history of msun/src/imprecise.c.  This is a 'workaround' 
> > committed by David Chisnall so that FreeBSD could claim C++ conformance.
> > David has never lifted a finger to fix that kludge.  I am the person
> > who has clean up that mess (until I stopped working on powl and tgammal).  
> > This is exactly the the situation that people on this list are trying to 
> > avoid when someone wants to simply "provide *some* reasonably working
> > implementation".
>
> See, it was great work, it has allowed him to publish stuff, e.g.:
> http://www.cl.cam.ac.uk/research/security/ctsrd/pdfs/201403-asiabsdcon2014-llvmbsd.pdf
> :-)
> 
> His homepage however says that he's "Former Core Team member" of FreeBSD.
> If https://wiki.freebsd.org/Numerics is the right place to look at, this
> issue ought to be mentioned there, don't you think so?

I objected when it was first suggested.  However, I wasn't
a member of Core, so you see were those objections went.

Feel free to update the wiki.

> > Instead of complaining about the people on this list, which is mostly
> > Bruce and me, why not help.  So, I'll repeat myself here.  What is the
> > quality of those implementations?  Have you tested the code?  Have you
> > read the comment in math_private.h concerning the use of I in computations?
> 
> While I do have some time to spend on this in the coming year or so, I
> would like to understand how the developent process works in this
> particular case,

For libm stuff, the development process in the recent past has been
Steve (or Stephen Montegomery or Peter Jeremy) writes code for a
function.  Steve does a significant amount of testing.  Steve submits
patch for freebsd-numerics.  Bruce notices patch and grabs it.  Bruce
tests patch, fixes patch, tests fixed patch, and then sends Steve a
long email with suggested improvement and/or issues with the original
code.  Steve understands about 10% of Bruce's email and tries to fix
the code.  Steve then does significant testing of the new code.  Steve
creates a new patch and submits it to FreeBSD-numerics.  Rinse-and-repeat.

It should also be noted that  I no longer have a commit bit and Bruce
rarely commits, so patch to libm lnger in bugzilla.

> who have rights to commit changes, how tests are run,
> whether there is any CI system to be used, etc.
> It seems that one has to start by setting up a set of tests for these
> functions.

At a minimum grab n1256.pdf.  Annex G has

G.6.3.2 The clog functions

* clog(conj(z)) = conj(clog(z)).
* clog(-0 + i*0) returns -inf + i*pi and raises the divide-by-zero
    floating-point exception.
* clog(+0 + i*0) returns -inf + i*0 and raises the divide-by-zero 
    floating-point exception.
* clog(x + i*inf) returns +inf + i*pi/2, for finite x.
* clog(x + i*NaN) returns NaN + i*NaN and optionally raises the 
    invalid floating-point exception, for finite x.
* clog(-inf + i*y) returns +inf + i*pi, for finite positive-signed y.
* clog(+inf + i*y) returns +inf + i*0, for finite positive-signed y.
* clog(-inf + i*inf) returns +inf + i*3*pi/4.
* clog(+inf + i*inf) returns +inf + i*pi/4.
* clog(+-inf + i*NaN) returns +inf + i*NaN.
* clog(NaN + iy) returns NaN + i*NaN and optionally raises the 
    invalid floating-point exception, for finite y.
* clog(NaN + i*inf) returns +inf + i*NaN.
* clog(NaN + i*NaN) returns NaN + i*NaN.

Do the NetBSD/OpenBSD meet the above conditions?  For finite
x+I*y, what is the accuracy of the return results? 

-- 
Steve