Skip site navigation (1)Skip section navigation (2)
Date:      Sat, 20 Oct 2012 18:32:14 +1100 (EST)
From:      Bruce Evans <brde@optusnet.com.au>
To:        Bruce Evans <brde@optusnet.com.au>
Cc:        svn-src-head@freebsd.org, svn-src-all@freebsd.org, src-committers@freebsd.org, Warner Losh <imp@freebsd.org>
Subject:   Re: svn commit: r241755 - head/lib/msun/src
Message-ID:  <20121020181443.X1561@besplex.bde.org>
In-Reply-To: <20121020150917.I1095@besplex.bde.org>
References:  <201210192246.q9JMkm4R092929@svn.freebsd.org> <20121020150917.I1095@besplex.bde.org>

next in thread | previous in thread | raw e-mail | index | archive | help
On Sat, 20 Oct 2012, Bruce Evans wrote:

> On Fri, 19 Oct 2012, Warner Losh wrote:
>
>> Log:
>>  Document the methods used to compute logf. Taken and edited from the
>>  double version, with adaptations for the differences between it and
>>  the float version.
>
> Please back this out.
>
> This was intentionally left out.  It is painful to maintain large
> comments that should be identical for all precisions.  Double precision
> is primary, and only comments that depend on the precision should be
> given in other precisions, except for short comments to to right of
> coulde whose non-duplication would just give larger diffs.  If you
> want to duplicate them, then they would often have to be quadruplicated
> in 4 files for the 4 supported precisions:
> - double precision
> - float precision
> - long double with MANT_DIG = 64 and 16 other bits
> - long double with MANT_DIG = 113 and 15 other bits
> Only 4 now, but there could easily be variations for long doubles.

PS: in long double trig functions, I handle this carefully by giving
references to the file containing the main description of the algorithm.
Old Cygnus float functions are no so careful, and I didn't want to churn
them for this.  Example from k_cosl.c:

% #include <sys/cdefs.h>
% __FBSDID("$FreeBSD: src/lib/msun/ld80/k_cosl.c,v 1.1 2008/02/17 07:32:14 das Exp $");
% 
% /*
%  * ld80 version of k_cos.c.  See ../src/k_cos.c for most comments.
%  */

Most details are here.  They are described in the same style as the
e_log.c comments that you duplicated.  Here they are in the kernel
instead of in the main file.  Another problem is that they really
do belong in the kernel, but are harder to read there.  We already
have kernels for log and logf (k_log.h and k_logf.c).  These
unfortunately have to duplicate the code, but they are careful to
not duplicate comments and carefully point to the main file for the
comments.  The kernels are only used for log10*() and log2*(), so
strictly whatever the kernels do doesn't affect log() and logf().
You didn't duplicate the comments into these files.  If you did,
there would already be 4 or 6 copies of them without even long double
precision:
- e_log.c, e_logf.c
- k_log.h, k_logf.c
- e_log10.c, e_log10f.c, e_log2.c, e_log2f.c (if we don't use kernels
   and/or don't put the comments in the kernels).

% 
% #include "math_private.h"
% 
% /*
%  * Domain [-0.7854, 0.7854], range ~[-2.43e-23, 2.425e-23]:
%  * |cos(x) - c(x)| < 2**-75.1
%  *
%  * The coefficients of c(x) were generated by a pari-gp script using
%  * a Remez algorithm that searches for the best higher coefficients
%  * after rounding leading coefficients to a specified precision.

This documents some precision-dependent details too verbosely.  It goes
without saying that we generate minimax polynomials by some method.
If we want to document the method(s), then it should be once per
method and not placed in the source files whose polynomials happened
to be generated using the methods.

%  *
%  * Simpler methods like Chebyshev or basic Remez barely suffice for
%  * cos() in 64-bit precision, because we want the coefficient of x^2
%  * to be precisely -0.5 so that multiplying by it is exact, and plain
%  * rounding of the coefficients of a good polynomial approximation only
%  * gives this up to about 64-bit precision.  Plain rounding also gives
%  * a mediocre approximation for the coefficient of x^4, but a rounding
%  * error of 0.5 ulps for this coefficient would only contribute ~0.01
%  * ulps to the final error, so this is unimportant.  Rounding errors in
%  * higher coefficients are even less important.
%  *
%  * In fact, coefficients above the x^4 one only need to have 53-bit
%  * precision, and this is more efficient.  We get this optimization
%  * almost for free from the complications needed to search for the best
%  * higher coefficients.
%  */

This documents in excessive detail the limits of some methods.  When it
was written, I didn't know the limits very well.  This is essentially
a way of saying "don't document your methods here, else it would expand
to be as verbose as this comment, if not to a book.  I don't want to
read this book in every source file".

Bruce



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20121020181443.X1561>