Status of expl logl

[Stephen Montgomery-Smith]

On 08/14/2012 05:15 PM, Bruce Evans wrote:
> On Tue, 14 Aug 2012, Steve Kargl wrote:
>
>> On Tue, Aug 14, 2012 at 02:39:29PM -0500, Stephen Montgomery-Smith wrote:
>>> On 08/14/2012 01:35 PM, Steve Kargl wrote:
>>>> On Tue, Aug 14, 2012 at 10:52:57AM -0700, Steve Kargl wrote:
>>>>> On Tue, Aug 14, 2012 at 12:37:16PM -0500, Stephen Montgomery-Smith
>>>>> wrote:
>>>>>> Are people working on expl, logl and log1pl?
>
> -rw-r--r--  1 bde  wheel  37470 Aug 10 00:08 ld128/s_logl.c
> -rw-r--r--  1 bde  wheel  32946 Aug 10 00:07 ld80/s_logl.c
> -rw-r--r--  1 bde  wheel  32056 Aug 10 00:06 s_log.c
> -rw-r--r--  1 bde  wheel  30053 Aug 10 00:05 s_logf.c
>
> These are mostly large comments and larger tables.  Each file implements
> log, log10, log2, log1p for the given precision.  Internal accuracy is
> about 7 extra bits of precision.  Speed on core2 with optimal CFLAGS is
> 30-70 cycles depending on the precision.
>
>>> ...
>>> So I am looking through src/e_pow.c.
>>>
>>> It seems to me that the constants L1, L2, L3, etc, are 3/5, 3/7, 3/9,
>>> etc, but not exactly these constants.  So they must have used some
>>> process where they jiggled the constants around, perhaps using trial and
>>> error, to get a few extra ulp.  Is that right?
>
> Hopefully not by trial and error.  There is the Remes algorithm.

Just looked it up.  It looks very nice.

>>> Also, I am trying to see what P1, P2, P3, etc are.  They seem to be
>>> related to the factorial (maybe a power series related to exp(x)), but I
>>> must admit that I am not getting it.
>
> These must be for exp().  Indeed, they are identical with the Pn's in
> e_exp.c.  They are essentially Bernoulli numbers.  One generating
> function for Bernoulli numbers with a certain normalization is
>
>      z/(exp(z) - 1) := sum(n = 0, Inf, B[n]/n!*z^n)
>
> and e_exp.c transforms exp(x) to essentially x/(exp(x) - 1).  Applying
> Remes to them makes Pn not quite a nice fraction.

Yes.  They are the coefficients of 2z/(exp(z)-1).  Thank you.