From mboxrd@z Thu Jan 1 00:00:00 1970 X-Msuck: nntp://news.gmane.org/gmane.linux.lib.musl.general/13878 Path: news.gmane.org!.POSTED.blaine.gmane.org!not-for-mail From: Szabolcs Nagy Newsgroups: gmane.linux.lib.musl.general Subject: Re: ATANH Date: Wed, 27 Feb 2019 11:38:24 +0100 Message-ID: <20190227103824.GD21289@port70.net> References: Reply-To: musl@lists.openwall.com Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Injection-Info: blaine.gmane.org; posting-host="blaine.gmane.org:195.159.176.226"; logging-data="39901"; mail-complaints-to="usenet@blaine.gmane.org" User-Agent: Mutt/1.10.1 (2018-07-13) To: musl@lists.openwall.com Original-X-From: musl-return-13893-gllmg-musl=m.gmane.org@lists.openwall.com Wed Feb 27 12:45:52 2019 Return-path: Envelope-to: gllmg-musl@m.gmane.org Original-Received: from mother.openwall.net ([195.42.179.200]) by blaine.gmane.org with smtp (Exim 4.89) (envelope-from ) id 1gywbn-00074P-AW for gllmg-musl@m.gmane.org; Wed, 27 Feb 2019 11:38:39 +0100 Original-Received: (qmail 10178 invoked by uid 550); 27 Feb 2019 10:38:36 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Original-Received: (qmail 10157 invoked from network); 27 Feb 2019 10:38:36 -0000 Mail-Followup-To: musl@lists.openwall.com Content-Disposition: inline In-Reply-To: Xref: news.gmane.org gmane.linux.lib.musl.general:13878 Archived-At: * Damian McGuckin [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