* ATANH
@ 2019-02-27 2:59 Damian McGuckin
2019-02-27 10:38 ` ATANH Szabolcs Nagy
0 siblings, 1 reply; 2+ messages in thread
From: Damian McGuckin @ 2019-02-27 2:59 UTC (permalink / raw)
To: musl
The comments for this routine say:
* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5)
There is a point where atanh(x) can be approximated by just x to machine
precision. This is where the exponent of x is less than some number, or
where x itself is less than some number.
In MUSL, for the 80-bit version, directly from the code, this is
2^(-LDBL_MANT_DIG/2) = 0x1.0p-32;
Interestingly, this same number, is used double atanh and float atanhf.
Note that in fdlibm, 0x1.0p-28 is used for all types
In fact, isn't
x * (1 + x^2/3) == x (after roundoff)
if
x < 2^(p/2),
i.e. 0x1.0p-26 for double, 0x1.0p-12 for float.
Note that for 80-bit wide, it is 0x1.0p-32 agreeing with MUSL!!
Regards - Damian
Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: ATANH
2019-02-27 2:59 ATANH Damian McGuckin
@ 2019-02-27 10:38 ` Szabolcs Nagy
0 siblings, 0 replies; 2+ messages in thread
From: Szabolcs Nagy @ 2019-02-27 10:38 UTC (permalink / raw)
To: musl
* Damian McGuckin <damianm@esi.com.au> [2019-02-27 13:59:41 +1100]:
> The comments for this routine say:
>
> * atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5)
>
> There is a point where atanh(x) can be approximated by just x to machine
> precision. This is where the exponent of x is less than some number, or
> where x itself is less than some number.
>
> In MUSL, for the 80-bit version, directly from the code, this is
>
> 2^(-LDBL_MANT_DIG/2) = 0x1.0p-32;
>
> Interestingly, this same number, is used double atanh and float atanhf.
>
> Note that in fdlibm, 0x1.0p-28 is used for all types
that's because musl assumes i386 uses FLT_EVAL_METHOD==2
so essentially evaluates everything to long double
precision.
fdlibm in the various BSDs assumes FLT_EVAL_METHOD==0
on all targets i think.
(in c99 the return statement is not required to drop
excess precision, so the different result is observable)
but i think either way is fine, it's not going to cause
huge errors.
note that the largest possible threshold may not be the
best choice: you want that branch to be correctly
predicted, so e.g. if 0x1p-40 input happens a few times,
but inputs below 0x1p-50 are very rare then a x < 0x1p-50
check may be better in practice: a 0x1p-40 input won't
break the branch predictor trained for the common case.
>
> In fact, isn't
>
> x * (1 + x^2/3) == x (after roundoff)
>
> if
> x < 2^(p/2),
>
> i.e. 0x1.0p-26 for double, 0x1.0p-12 for float.
>
> Note that for 80-bit wide, it is 0x1.0p-32 agreeing with MUSL!!
>
> Regards - Damian
>
> Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
> Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
> Views & opinions here are mine and not those of any past or present employer
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2019-02-27 10:38 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2019-02-27 2:59 ATANH Damian McGuckin
2019-02-27 10:38 ` ATANH Szabolcs Nagy
Code repositories for project(s) associated with this public inbox
https://git.vuxu.org/mirror/musl/
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).